l'elaborazione di file MultiFASTA

Question

ero curioso di sapere se c'è qualche strumento di bioinformatica là fuori in grado di elaborare un file MultiFASTA avermi dato informazioni come il numero di sequenze, la lunghezza, nucleotide / contenuto di aminoacidi, ecc e magari disegnare automaticamente diagrammi descrittivi. Anche una soluzione R Bioconductor o un modulo Bioperl farebbe, ma l'ho fatto non riesce a trovare nulla.

Mi potete aiutare? Grazie molto: -)

Answer 1

Alcuni degli strumenti sono imprime una collezione di piccoli strumenti che possono aiutare a uscire.

• seqstats ritorna lunghezza della sequenza
• pepstats dovrebbe darvi contenuto di aminoacidi, ecc Alcuni degli strumenti offrono anche funzioni di stampa. Molto maneggevole. http://emboss.sourceforge.net/apps/release/5.0 /emboss/apps/groups.html

Per contare il numero di voci FASTA, io uso:  grep -c '^>' mySequences.fasta.

Per assicurarsi che nessuna delle voci sono duplicati, verifico che ottengo lo stesso numero quando si fa questo: grep '^>' mySequences.fasta | sort | uniq | wc -l

Answer 2

Si può anche essere interessati a faSize , che è uno strumento dal Kent Fonte Albero , anche se questo richiede un po 'più di sforzo (è necessario DLOAD e compilare) che usare grep ... ecco qualche esempio di output:

me@my-lab ~/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

Answer 3

Si segnala (per chiunque inciampo su questo, come ho appena fatto) che v'è una libreria Python robusta specificamente progettato per gestire questi compiti chiamati Biopython . In poche righe di codice, è possibile accedere rapidamente risposte per tutte le domande di cui sopra. Ecco alcuni esempi molto semplici, per lo più adattati da link. Esistono caldaia-piastra GC% grafici e grafici lunghezza sequenza nel tutorial anche.

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(), '*'))