procesamiento de archivos multiFASTA

Question

Tenía curiosidad por saber si hay alguna herramienta bioinformática capaz de procesar un archivo multiFASTA que me dé información como cantidad de secuencias, longitud, contenido de nucleótidos/aminoácidos, etc. y tal vez dibujar automáticamente gráficos descriptivos. También funcionaría una solución R BIoconductor o un módulo BioPerl, pero no logré encontrar nada.

¿Me puede ayudar? Muchas gracias :-)

Answer 1

Algunas de las herramientas de relieve son una colección de pequeñas herramientas que pueden ayudarte.

• seqstats vuelve secuencia de longitud
• pepstats deben darle contenido de aminoácidos, etc. Algunas de las herramientas también ofrecen funciones de trazado. 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 son duplicado, puedo comprobar que tengo el mismo número al hacer esto: grep '^>' mySequences.fasta | sort | uniq | wc -l

Answer 2

También puede estar interesado en faSize, que es una herramienta de la Kent Source Tree, aunque esto requiere un poco más de esfuerzo (debe DLOAD y compilar) que simplemente usar grep ... aquí está una cierta 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 robusta biblioteca de python específicamente diseñada para manejar estas tareas llamada ed Biopython. En algunas líneas de código, puede acceder rápidamente a las respuestas de todas las preguntas anteriores. Aquí hay algunos ejemplos muy básicos, en su mayoría adaptados del enlace. También hay gráficas GC% de la placa de caldera y gráficos de longitud de secuencia en el 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(), '*'))