multiFASTA文件处理

Question

我是好奇的知道如果有任何生物信息学工具有能够处理一multiFASTA文件给我的相关信息等数量的序列长度，核苷酸/氨基酸的内容，等等。也许动画描述性的地块。还R使用bioconductor解决方案或一个BioPerl模块会怎么做但是我没找到任何东西。

你能帮帮我吗？非常感谢:-)

Answer 1

一些浮雕工具是可以帮助您的小工具的集合。

• seqstats 返回序列长度
• pepstats 应该给你氨基酸含量等。一些工具还提供绘图功能。非常便利。http://emboss.sourceforge.net/apps/release/5.0/emboss/apps/groups.html

为了计算 fasta 条目的数量，我使用： grep -c '^>' mySequences.fasta.

为了确保所有条目都不重复，我检查执行此操作时是否获得相同的数字： grep '^>' mySequences.fasta | sort | uniq | wc -l

Answer 2

你也可能会感兴趣 faSize, ，这是一个工具，从 肯特源树, 虽然这需要一点更多的努力(你必须dload和编纂)不仅仅是使用查询...这里是一些例子产出：

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

应该指出的（任何人绊倒在此，就像我只是做了），有一个强大的Python库专门用来处理这些任务叫的 Biopython 。在几行代码，你可以为所有的上述问题迅速获得答案。这里有一些非常基本的例子，大多是从链接，其适于。有样板GC％图和序列长度的曲线图在教程还

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