k-mer カウントして、配列も出力するツール jellyfish、BFCounter (original) (raw)

2017.11追記

2018.03 -sフラグ修正

2019 5/14 リンク追加

2019 9/9 BFcounterインストール追記

2019 12/19 condaインストール追記、BFCounterインストール追記

k-merカウントを行うjellyfishと、k-merの全配列を書き出すBFCounterを紹介する。

Jellyfish

公式サイト

JELLYFISH - Fast, Parallel k-mer Counting for DNA

マニュアルPDF

http://www.cbcb.umd.edu/software/jellyfish/jellyfish-manual-1.1.pdf

Github

git clone https://github.com/gmarcais/Jellyfish.git
cd Jellyfish/
autoreconf -i
 ./configure --prefix=$HOME
make -j 8
sudo make install

#bioconda(link)
#version1系
conda install -c bioconda -y kmer-jellyfish==1.1.12-0

#version2 latest
conda install -c bioconda -y kmer-jellyfish

またはbrewでも導入できる。

ラン

k-merカウントする。

jellyfish count -m 31 -s 10000000 -t 20 input.fasta/fastq -o mer_counts.jf

• -m Length of mer
• -s Initial hash size
• -t Number of threads (1)

インプットファイルはfasta ,fastqとも対応している。

mer_counts.jfができる。中身はbinaryである。このファイルを以後の解析に使う。

簡単なstatistics

jellyfish stats mer_counts.jf

user$ jellyfish stats mer_counts.jf

Unique: 2413100

Distinct: 10097611

Total: 69409538

Max_count: 388

DIstinct k-merはたとえ複数回出現しても1回しか数えられないk-merで、Unique k-merは1回しか出現しないk-mer。Canonicalかどうかでカウント値は影響を受ける。

k-merの配列を全て出力する。

jellyfish dump mer_counts.jf > k31.fasta

中身を見てみる。

user-no-MacBook-Pro: user$ less -S k33.fa

>1

CCCAGGGTCTGGAACTATCGCCGGTGGAAAATATCGAAACGGAAGCGGGGCAAGGGGTTCAGGGTTGGTACCAAGGC

>15

AATGCCCTGGGGAGTAACTGAAAGGAAACGAAAAAATAATCTTCTTCCTATTCTATTCCTTCTGATTGAGTAACCAG

>1

CCTTGTTTGACAGTGCAGTCCGGGTTCCTAAAAGGTTGCAAGGGGATTTATTCATTCGCAGTACCCAAACGGACACC

>1

GGTGGAATAAATACCTGTGGTGAAGGTATTAACGCCAAAATATTGCACTGTGCCAAAATCGTTCAGAGTCCCCATCA

>1

AGAACGTCGTGGACGGCGTCCCCGTCGTGAGGGAACGGTGGAGGTCAAAGGGCGATCGGGAAGGGCGGCGGAACTAT

>16

TGGTGGACGGAAATCTTAAACATTGCCTAGAGAAAAGTTGGGAATTGACGGTGGTGGGGGCAGGGCAAGGAGTGGAA

>15

GGTTTCCTGGGCAAACTCGATGGAACCTGGACAATCGAGGAAATTTAAGCGCAAATTTTCATAATCAATCCCTGCCA

指定のk-merサイズの全塩基配列が出力されている。数値は登場回数である。重複させながら元の配列を分解しているので、元のfastqよりサイズははるかに大きくなる。

ヒストグラムデータを出力。

jellyfish histo mer_counts.jf > histogram.txt

出力の中身を見てみる。

user-no-MacBook-Pro: user$ less -S histogram.txt

1 3101726

2 8283

3 5424

4 12735

5 29386

6 61451

7 103516

8 164603

9 230629

10 291824

11 347466

Rに読み込んでグラフ化する。

純化されたバクテリアのシーケンスデータを分析する。

R #Rに入る
input <- read.table("histogram.txt") #読み込み
plot(input[2:100,],type="l") #2列目が頻度情報。

２つ目の山がゲノムのシングルコピーの領域に由来する33merのピークと思われる。ピーク座標を調べる。

input[0:20,]

13がピークと思われる（赤色）。

> input[0:20,]

V1 V2

1 1 136238639

2 2 3362948

3 3 446016

4 4 142167

5 5 122950

6 6 159056

7 7 216472

8 8 287219

9 9 362505

10 10 430060

11 11 487906

12 12 520707

13 13 531598

14 14 531210

15 15 504837

16 16 472747

17 17 435696

18 18 390053

19 19 342325

20 20 298075

もっとちゃんと調べる。

k15

k35

k55

k65

k77

k99

k127

kmergenuineでもbestなk-merを探す（参考にしたページ）。（kmergenuineはbrewで導入できる。）

kmergenie input.fq

fitting model to histograms to estimate best k

table of predicted num. of genomic k-mers: histograms.dat

recommended coverage cut-off for best k: 2

best k: 65

kmergenuineでは65-merがベストと判定された。

65-merの入力テーブルを使い、ゲノムサイズを推定する。

sum(as.numeric(input[12:1130,1]*input[12:1130,2]))/65

リピートがあると、その分だけ少ない値が出るはずである。

k-mer coverageの貧弱なものを除きたいならkhmerを使う。khmerはk-merグラフの分布から低頻度k-mer配列を除去したり様々な機能をもつde novo assemblyをアシストするツールである。（pipでインストール可能）

load-into-counting.py -T 20 -k 20 -x 3e6 output1 input.fastq

• T スレッド数
• k k-merサイズ
• x 大まかなゲノムサイズ。

filter-abund.py -C 2 -T 20 output1 input.fastq

注 khmrは別途記事にします。

BFCounter

BFCounterは k-merを数えて、全部の配列を出力するツール。

Github

https://github.com/pmelsted/BFCounter

ソースコードをダウンロードしてmakeする。

git clone https://github.com/pmelsted/BFCounter.git
cd BFCounter/

#ここではlarge k-merを扱えるようにしておく
make MAX_KMER_SIZE=300

パスを通しておく。

# ./BFCounter

Error: too few arguments

BFCounter BFCounter 1.0

A memory efficient k-mer counting program.

Usage: BFCounter [options] ...

Where can be one of:

count Counts the occurrences of k-mers in sequence files

dump Writes occurrences of k-mers into a tab-separated text file

cite Prints information for citing the paper

version Displays version number

ラン

BFCounter count -n 1000 -k 127 -o single_out -t 12 input.fq

• -n Estimated number of k-mers
• -k Size of k-mers, at most 301
• -t Number of threads to use (default 1)

出力はbinaryファイルである。ここからk-mer長の配列を抽出する。

BFCounter dump -i single_out -k 127 -o k-127_output

• -i Filename for input file from count
• -k Size of k-mers, at most 301
• -o Filename for output

出力の中身を見てみる。

このように全k-mer配列が出力されている。あとはスクリプトを使い自由にいじることができる。

このツールの機能はこれだけであるが、シンプルでなため使いやすい。

追記

１、シロイヌナズナゲノム (k-mer size 28)

$ jellyfish stats mer_counts.jf

Unique: 109610472

Distinct: 112718914

Total: 119479935

Max_count: 4778

28merという大きなサイズでゲノムのk-merを調べると、偶然２回以上出現する配列は限りなくゼロになる。そのため、ユニークなk-merの出現頻度は大半が１になり（下の図）、頻度1のユニークなk-merの数とゲノムのサイズはほぼ等しくなる。しかしリピートがあるゲノム（28merで調べているなら28以上完全に同一の配列がある）ならば、同じ28-merの配列が複数回出現するため、出現頻度1のユニークなk-merの合計と、全k-merのトータル（重複があるかどうかを問わず）に大きな差が生じる。リピート領域が多ければ、ヒストグラムの裾野の分布（肩）が右向きに広く（コピー数が多い）、高く（コピーの割合が多い）なる。

２、シーケンスデータのk-mer分布。

シロイヌナズナのillumina rawシーケンスデータ (k-mer size 28)

$ jellyfish stats mer_counts.jf

Unique: 241278948

Distinct: 475904703

Total: 3772597024

Max_count: 843576

シーケンスデータのk-merカウントを行うと、高いカバレッジでシーケンスしているため、ゲノムがリピートを持つかどうかに関わらず、ユニークなk-merはほぼゼロという結果が得られる（下の図）。それにもかかわらずユニークなk-merの数が非常に多いのは、シーケンスエラーが発生しているためである。平均クオリティQ30の100bpリードは、平均して10リードに１つしかシーケンスエラーを含まないが、それでも１つのシーケンスエラーが発生すると、28merでカウントしているなら28のユニークk-merが出現する（変化して他に配列と合致する可能性がゼロとして）。リアルデータだと、シーケンスエラーの発生はそこまで単純ではないため、おそらくはユニークk-mer数はもっと多くなる。参考リンク

図、左端はシーケンスエラー由来のピークで、ゲノムアセンブリに必要な真の情報は頻度5-15付近のピーク部分。リピートが多いゲノムのシーケンスデータでは、右に大きな第二ピークができる。

引用

Guillaume Marcais and Carl Kingsford, A fast, lock-free approach for efficient parallel counting of occurrences of k-mers.

Bioinformatics (2011) 27(6): 764-770 (first published online January 7, 2011) doi:10.1093/bioinformatics/btr011)

Efficient counting of k-mers in DNA sequences using a bloom filter.

BMC Bioinformatics. 2011 Aug 10;12:333. doi: 10.1186/1471-2105-12-333.

Melsted P1, Pritchard JK.

K-mer analysis and genome size estimate

関連