multiFASTA ファイルの処理

StackOverflow https://stackoverflow.com/questions/1789307

質問

私は、配列数、長さ、ヌクレオチド/アミノ酸含有量などの情報を提供する multiFASTA ファイルを処理できるバイオインフォマティクス ツールがあるかどうか知りたいと思っていました。説明的なプロットを自動的に描画することもあります。R BIoconductor ソリューションや BioPerl モジュールも使えますが、何も見つかりませんでした。

手伝ってもらえますか？どうもありがとう ：-）

役に立ちましたか？

解決

一部のエンボスのツールを集めた小さなツールできます。

カウント数fasta作品に使っている: grep -c '^>' mySequences.fasta.

めていないかどうかの複製、まることを確認を取得しますことになっ grep '^>' mySequences.fasta | sort | uniq | wc -l

他のヒント

あなたも興味があるかもしれません ファサイズ, のツールです。 ケントソースツリー, ただし、これには 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

(私と同じように、これにつまずいた人のために) と呼ばれる、これらのタスクを処理するために特別に設計された堅牢な Python ライブラリがあることに注意してください。 バイオパイソン. 。数行のコードで、上記のすべての質問に対する答えにすぐにアクセスできます。以下に非常に基本的な例をいくつか示します。ほとんどがリンクから転用されたものです。このチュートリアルには、定型的な 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(), '*'))
ライセンス： CC-BY-SA と 帰属
所属していません StackOverflow
scroll top