processamento de arquivos multiFASTA
-
22-09-2019 - |
Pergunta
Fiquei curioso para saber se existe alguma ferramenta de bioinformática capaz de processar um arquivo multiFASTA, fornecendo informações como número de sequências, comprimento, conteúdo de nucleotídeos/aminoácidos, etc.e talvez desenhar gráficos descritivos automaticamente.Também uma solução R BIoconductor ou um módulo BioPerl serviriam, mas não consegui encontrar nada.
Pode me ajudar?Muito obrigado :-)
Solução
Algumas das ferramentas de relevo são uma coleção de pequenas ferramentas que podem ajudá -lo.
seqstats
Retorna o comprimento da sequênciapepstats
deve fornecer conteúdo aminoácido etc. Algumas das ferramentas também oferecem funções de plotagem. Muito conveniente.http://emboss.sourceforge.net/apps/release/5.0/emboss/apps/groups.html
Para contar o número de entradas da FASTA, eu uso: grep -c '^>' mySequences.fasta
.
Para garantir que nenhuma das entradas seja duplicada, verifico que recebo o mesmo número ao fazer isso: grep '^>' mySequences.fasta | sort | uniq | wc -l
Outras dicas
Você também pode estar interessado em Fasize, que é uma ferramenta do Árvore fonte de Kent, embora isso exija um pouco mais de esforço (você deve dar e compilar) do que apenas usar o Grep ... aqui está algum exemplo de saída:
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
Deve-se notar (para qualquer um que se depare com isso, como eu acabei de fazer) que existe uma biblioteca python robusta projetada especificamente para lidar com essas tarefas chamada Biopíton.Em algumas linhas de código, você pode acessar rapidamente as respostas para todas as perguntas acima.Aqui estão alguns exemplos muito básicos, em sua maioria adaptados do link.Existem gráficos padrão de GC% e gráficos de comprimento de sequência no tutorial também.
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(), '*'))