我很想知道是否有任何生物信息学工具能够处理 multiFASTA 文件,为我提供序列数量、长度、核苷酸/氨基酸含量等信息,并可能自动绘制描述图。也可以使用 R BIOconductor 解决方案或 BioPerl 模块,但我没有找到任何东西。
你能帮助我吗?非常感谢 :-)
我很想知道是否有任何生物信息学工具能够处理 multiFASTA 文件,为我提供序列数量、长度、核苷酸/氨基酸含量等信息,并可能自动绘制描述图。也可以使用 R BIOconductor 解决方案或 BioPerl 模块,但我没有找到任何东西。
你能帮助我吗?非常感谢 :-)
一些浮雕工具是可以帮助您的小工具的集合。
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
您可能还对faSize感兴趣,它是Kent Source Tree中的一个工具,尽管这比仅使用 grep 需要更多的努力(您必须 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
应该注意(对于任何偶然发现这一点的人,就像我刚才所做的那样),有一个强大的 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(), '*'))
python 中的熨平板很棒。
import screed
for record in screed.open(fastafilename):
print(len(record.sequence))