multiFASTA文件處理

Question

我很想知道是否有任何生物信息學工具可以處理multiFASTA文件，這些文件可以爲我提供諸如序列號，長度，核苷酸/氨基酸含量等信息，也可以自動繪製描述性圖表。 R生物傳感器解決方案或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，這是從Kent Source Tree的工具，雖然這需要一些更多的努力（你必須DLOAD和編譯）不僅僅是用grep ...這裏是一些示例輸出：

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