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 :-)

Foi útil?

Solução

Algumas das ferramentas de relevo são uma coleção de pequenas ferramentas que podem ajudá -lo.

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(), '*'))
Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top