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

あなたはまた、Kent Source TreeからツールですfaSize、に興味がある可能性があり：

は、エントリの確かどれも私はこれをやって同じ番号を取得することを確認し、重複しないようにするには

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

（私はちょうど行ったように、この時につまずく人のために）それは注意する必要があります。これは、ここではいくつかの出力例です...ただのgrepを使用するよりも（あなたがDLOADとコンパイルする必要があります）もう少し努力が必要これらのタスクを処理するために特別に設計された堅牢なPythonライブラリがありますed 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(), '*'))