procesamiento de archivos multiFASTA

Question

Yo tenía curiosidad por saber si hay alguna herramienta bioinformática ahí fuera capaz de procesar un archivo multiFASTA darme informaciones como el número de secuencias, la longitud, el nucleótido / contenido de aminoácidos, etc., y tal vez llamar automáticamente parcelas descriptivos. También una solución Bioconductor R o un módulo BioPerl lo haría, pero no lograron encontrar nada.

¿Me puede ayudar? Muchas gracias: -)

Answer 1

Algunas de las herramientas de estampado son un conjunto de pequeñas herramientas que pueden ayudarle a cabo.

• longitud de la secuencia vuelve seqstats
• pepstats debe darle aminoácido contenido etc. Algunas de las herramientas también ofrecen tablas de funciones. Muy útil. http://emboss.sourceforge.net/apps/release/5.0 /emboss/apps/groups.html

Para contar el número de entradas FASTA, que utilizo:  grep -c '^>' mySequences.fasta.

Para asegurarse de que ninguna de las entradas están duplicados, puedo comprobar que tengo el mismo número al hacer esto: grep '^>' mySequences.fasta | sort | uniq | wc -l

Answer 2

Es posible que también esté interesado en faSize , que es una herramienta de la Kent árbol de código fuente , aunque esto requiere un poco más de esfuerzo (debe DLOAD y compilación) que simplemente usar grep ... aquí hay alguna salida de ejemplo:

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

Debe tenerse en cuenta (para cualquier persona tropezar con esto, como acabo de hacer) que hay una biblioteca de Python robusta diseñada específicamente para manejar estas tareas llamados Biopython . En unas pocas líneas de código, puede acceder rápidamente respuestas para todas las preguntas anteriores. Estos son algunos ejemplos muy básicos, sobre todo adaptados desde el enlace. Hay caldera de placa GC% gráficos y gráficos de longitud de secuencia en el tutorial también.

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