Conseils: [Python] Restaurer et extraire au hasard un tableau à partir d'un fichier fasta

Il y a des moments où vous souhaitez échantillonner au hasard un grand nombre de séquences. Cette fois, je présenterai le code à utiliser dans un tel cas.

La méthode consiste simplement à ouvrir le fichier, à sélectionner un tableau avec des nombres aléatoires et à y accéder.

Accès en sélectionnant un tableau avec des nombres aléatoires

#Obtenir un tableau
seqnum=random.randint(1,seqnumMAX)
name = linecache.getline(infn, seqnum*2-1) #Nom de la séquence
seq = linecache.getline(infn, seqnum*2)　 #Tableau
#l'écriture
outfdl    =   open(outfn, 'w')
outfdl.writelines(name+seq)
outfdl.close()
linecache.clearcache() #Vider le cache

À ce stade, utilisez linecache.getline (nom du fichier d'entrée, nombre de lignes) pour lire. Cela optimisera en interne la lecture de la ligne spécifiée. La bioinformatique a tendance à gérer de gros fichiers de séquence, ce qui peut être utile.

Une extraction de restauration aléatoire peut être effectuée en bouclant autant de fois que vous le souhaitez pour acquérir ce processus.

Ce qui précède peut être résumé comme un processus approprié comme suit.

randsampleFasta.py

#!/usr/local/bin/python3
# -*- coding: utf-8 -*-

"""
Restaurez et extrayez au hasard un tableau à partir d'un fichier fasta.
"""
__author__  = "Kazuki Nakamae <[email protected]>"
__version__ = "0.00"
__date__    = "29 May 2017"

import sys
import mmap
import linecache
import random

def randsampleFasta(infn,outfn,n):
    """
Restaurez et extrayez au hasard un tableau à partir d'un fichier fasta.
    @param infn nom de fichier fasta à lire{string} 
    @param sortie sortie fasta nom de fichier{string} 
    @param n Nombre de tableaux à restaurer et extraire{int} 
    """

    infdl     =   open(infn, 'r')
    #Vérifiez le nombre de séquences dans le fichier d'entrée
    print('Checking how many sequences are in a FASTA format sequence file....(1/2)')
    buf = mmap.mmap(infdl.fileno(), 0, prot=mmap.PROT_READ)
    seqnumMAX= 0
    readline = buf.readline
    while readline():
        seqnumMAX += 1
    seqnumMAX = int(seqnumMAX / 2)
    infdl.close()
    print(str(seqnumMAX)+' sequences')
    #Restauration et extraction aléatoires de séquences
    print('random sampling from a sequence with replacement....(2/2)')
    outfdl    =   open(outfn, 'w')
    random.seed(a='hoge', version=2)    #Réglage des semences
    seqi=1
    while seqi<=n:
        seqnum=random.randint(1,seqnumMAX)
        name = linecache.getline(infn, seqnum*2-1)
        seq = linecache.getline(infn, seqnum*2)
        outfdl.writelines(name+seq)
        seqi += 1
    outfdl.close()
    linecache.clearcache() 
    print('done.')

if __name__ == '__main__':
    argvs = sys.argv  #Arguments de ligne de commande
    argc = len(argvs) #Nombre d'arguments

    if (argc != 4):   #Vérification des arguments
        print("USAGE : python3 randsampleFasta.py <INPUT.fa> <OUTPUT.fa> <NUMBER_OF_SEQUENCES>")
        quit()

    randsampleFasta(argvs[1],argvs[2],int(argvs[3]))
quit()

Exemple

Fichier à saisir

test.fa

>1
CCGTATTGGAAAGCTC
>2
AGGATTATCGGATACT
>3
ATCCGGACGGGGGGTT
>4
GACCTCGTTATCATCC
>5
AGTCAGGTTACCCGCA

Entrée sur Bash

contribution

python3 randsampleFasta.py test.fa out.fa 4

Sortie sur Bash

Sortie standard

Checking how many sequences are in a FASTA format sequence file....(1/2)
5 sequences
random sampling from a sequence with replacement....(2/2)
done.

Fichier de sortie

out.fa

>3
ATCCGGACGGGGGGTT
>2
AGGATTATCGGATACT
>4
GACCTCGTTATCATCC
>3
ATCCGGACGGGGGGTT

Utilisation du code

Tu es libre.