k-mer
現在のアセンブラで主流となってる「ド・ブラン（de Brujin1）グラフ」を使ったアセンブリで使用される、配列断片長の事。
で、ド・ブラングラフに関しては、
を読め。で話が終了。
が、インフォマティシャン以外の人に向けて、もう、ほんのちょっとだけ実例が欲しい場合の資料として。
本職の人は読む必要が無い、つか読まないで。
適当な例。
以下は、 E.coli K-12 の適当な部分を適当に持ってきた配列。
A: GCTTTTTTTTTCGACCAAAGGTAACGAGGTAAC
B: AGGTAACAACCATGCGAGTGTTGAAGTTC
で、コンセンサスが判らないものとして、 A, B のリードが得られてると仮定。
今回は、 5mer で 1bp スライドのスライディングウィンドウ的な断片を作成。それを並べてグラフを作成してアセンブリする。
目次と、配列断片を得る
$ function get_frag(){
perl -le 'while(1){ $a = substr $ARGV[0], $b ++ , 5 ; last if length $a < 5 ; print join "\t", $b, $a }' "$1"
}
$ get_frag GCTTTTTTTTTCGACCAAAGGTAACGAGGTAAC
1 GCTTT
2 CTTTT
3 TTTTT
4 TTTTT
5 TTTTT
6 TTTTT
7 TTTTT
8 TTTTC
9 TTTCG
10 TTCGA
11 TCGAC
12 CGACC
13 GACCA
14 ACCAA
15 CCAAA
16 CAAAG
17 AAAGG
18 AAGGT
19 AGGTA
20 GGTAA
21 GTAAC
22 TAACG
23 AACGA
24 ACGAG
25 CGAGG
26 GAGGT
27 AGGTA
28 GGTAA
29 GTAAC
$ get_frag AGGTAACAACCATGCGAGTGTTGAAGTTC
1 AGGTA
2 GGTAA
(ommitted)
24 AAGTT
25 AGTTC
A, B それぞれの目次があるので1から順にリード内の並びは再構成出来る。
断片のテーブルを作成する
上記、 A, B のコマンド結果を適当な1ファイルに保存。
$ perl -Mvars=%h -lane '$h{$F[1]} ++ }{ print join "\t", $h{$_}, $_ for sort { $h{$b} <=> $h{$a} } keys %h ' FILE
5 TTTTT
3 AGGTA
3 GGTAA
3 GTAAC
1 AAAGG
1 AACAA
1 AACCA
(omitted)
指定長の断片の出現頻度が得られる2。
これを、目次と照らし合わせて3、オイラー経路として解析する、、、と。
グラフを描く
例なので目次を便りに簡略化。
1 GCTTT
2 CTTTT
3 TTTTT # polyT
4 TTTTT #
5 TTTTT #
6 TTTTT #
7 TTTTT #
8 TTTTC
↓
18 AAGGT
19 AGGTA # 3
20 GGTAA # 3
21 GTAAC # 3
22 TAACG
↓
26 GAGGT
27 AGGTA # 3
28 GGTAA # 3
29 GTAAC # 3
1 AGGTA # 3
2 GGTAA # 3
3 GTAAC # 3
4 TAACA
↓
25 AGTTC
まあ、見たまんまな訳だけど、一応整理して並べると
- A の 1~18 までは（polyT はあるものの）目次通り => A_frag1 ( GCTTTTTTTTTCGACCAA )
- A の 19,20,21 => segment ( AGG )
- A の 22~26 => A_frag2 ( TAACG )
- B の 4~25 => B_frag ( TAACAACCATGCGAGTGTTGAA )
A_frag1
↓
+----> segment ----> B_frag
| ↓
+----- A_frag2
で、並び順は、
A_frag1 => segment => A_frag2 => segment => B_frag
になる。
よって答えは、経路の最後の 4 塩基を足して
GCTTTTTTTTTCGACCAA . AGG . TAACG . AGG . TAACAACCATGCGAGTGTTGAA + GTTC
アライメント結果は、
GCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTC
A: GCTTTTTTTTTCGACCAAAGGTAACGAGGTAAC
B: AGGTAACAACCATGCGAGTGTTGAAGTTC
で、ぐだぐだと述べて来た、この 5mer を k-mer = 5 と表記する訳。
k-mer = 20 なら 20塩基だし、 50 なら 50塩基。
ちょっと昔の次世代シーケンサーは、リード長が 100bp にも満たなかったから、k-mer 20 辺りでやるのが当たり前とか言われてたけど、最近 300bp 読めてるんだから、当然、もっとデカい値でのアセンブルが最適なものも出てくるよね。って話が続くんだけど、、、