2009-11-24 12 views
7

Ich war neugierig zu wissen, ob es ein Bioinformatik-Tool da draußen gibt, das eine multiFASTA-Datei verarbeiten kann, die mir Informationen wie Anzahl der Sequenzen, Länge, Nucleotid-/Aminosäuregehalt usw. gibt und vielleicht automatisch beschreibende Plots zeichnet. Auch eine R BIoconductor-Lösung oder ein BioPerl-Modul würde reichen, aber ich habe nichts gefunden.multiFASTA Dateiverarbeitung

Können Sie mir helfen? Vielen Dank :-)

Antwort

7

Einige der Prägewerkzeuge sind eine Sammlung kleiner Werkzeuge, die Ihnen helfen können.

Anzahl von fasta Einträge zu zählen, die ich benutze: grep -c '^>' mySequences.fasta.

Um sicherzustellen, dass keiner der Einträge ist doppelt zu machen, überprüfe ich, dass ich die gleiche Nummer erhalten, wenn dies zu tun: grep '^>' mySequences.fasta | sort | uniq | wc -l

2

Sie können auch daran interessiert sein faSize, die ein Werkzeug aus den Kent Source Tree ist, obwohl dies ein wenig mehr Aufwand erfordert (Sie DLOAD und kompilieren müssen) als nur grep ... hier ist ein Beispiel einer Ausgabe:

[email protected] ~/data $ time faSize myfile.fna 
215400419 bases (104761 N's 215295658 real 215295658 upper 0 lower) in 731620 sequences in 1 files 
Total size: mean 294.4 sd 138.5 min 30 (F5854LK02GG895) max 1623 (F5854LK01AHBEH) median 307 
N count: mean 0.1 sd 0.4 
U count: mean 294.3 sd 138.5 
L count: mean 0.0 sd 0.0 
%0.00 masked total, %0.00 masked real 

real 0m3.710s 
user 0m3.541s 
sys  0m0.164s 
0

Es sollte beachtet werden (für jedermann auf diese stolpern, wie ich gerade), dass Es gibt eine robuste Python-Bibliothek, die speziell für diese Aufgaben entwickelt wurde ed Biopython. In einigen Codezeilen können Sie schnell auf Antworten für alle oben genannten Fragen zugreifen. Hier sind einige sehr einfache Beispiele, die größtenteils aus dem Link stammen. Es gibt GC-Graphen im GC-Diagramm und Sequenzlängengraphen im Tutorial.

In [1]: from Bio import SeqIO 

In [2]: allSeqs = [seq_record for seq_record in SeqIO.parse('/home/kevin/stack/ls_orchid.fasta', """fasta""")] 

In [3]: allSeqs[0] 
Out[3]: SeqRecord(seq=Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', SingleLetterAlphabet()), id='gi|2765658|emb|Z78533.1|CIZ78533', name='gi|2765658|emb|Z78533.1|CIZ78533', description='gi|2765658|emb|Z78533.1|CIZ78533 C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA', dbxrefs=[]) 

In [4]: len(allSeqs) #number of unique sequences in the file 
Out[4]: 94 

In [5]: len(allSeqs[0].seq) # call len() on each SeqRecord.seq object 
Out[5]: 740 

In [6]: A_count = allSeqs[0].seq.count('A') 
     C_count = allSeqs[0].seq.count('C') 
     G_count = allSeqs[0].seq.count('G') 
     T_count = allSeqs[0].seq.count('T') 

     ​print A_count # number of A's 

     144 

In [7]: allSeqs[0].seq.count("AUG") # or count how many start codons 
Out[7]: 0 

In [8]: allSeqs[0].seq.translate() # translate DNA -> Amino Acid 
Out[8]: Seq('RNKVSVGEPAEGSLMRPWNKRSSESGGPVYSAHRGHCSRGDPDLLLGRLGSVHG...*VY', HasStopCodon(ExtendedIUPACProtein(), '*'))