bioinformatics - 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

Kent Source TreeのツールであるfaSizeにも興味があるかもしれませんが、これには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

Answer 3

Biopythonと呼ばれるこれらのタスクを処理するために特別に設計された堅牢な 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(), '*'))