Pindel output についてのメモ
The third line and further show the sequence of the read ... the position of the mapped half of the paired-end read, the mapping quality of the mapped read
というのはVariant Call する際に用いたリードのペアになっているリード（Pindel のアルゴリズム上、アンカーになっているはずのリード）がレファレンス配列上にマップされた位置のことである。これが誤っている場合があるのでメモしておく。
Pindel では、ペアエンドで塩基を読むとき中央付近の塩基が両サイドから二度読まれることを考慮していないようである。例えばペアエンドで100塩基ずつ読んでいる場合、フラグメントサイズが200を下回ると不具合が起きる。config ファイルでフラグメント長の平均値を入力するが、このとき200以下を入力すると走らない。
これと同様のことがoutput ファイルにおける上記エラーにも起きている。ペアリードのマップ位置が200以上離れていない時には、なぜだかはわからないが、confファイルに入力したフラグメント平均長の数値を足してマップ位置をずらして出力するようである。
東大医科研のスパコンはデフォルトのC コンパイラがICC
最適化の関係でしょうか？
GCCでないとコンパイルできないものもあるんですが、ICCだからダメだと言われるんじゃなくて「GCCのバージョンが古い」と怒られることが多く、gcc --version とやっても最新バージョンなのになー、と首を傾げることになります。例えばこんな感じ
$ ./configure : : configure: error: Your version of gcc does not support the 'std::tr1' standard. Recommended gcc version is 4.1.2 or later. Please use a newer gcc version, or try to download the pre-compiled binaries from the fastx-toolkit website. $ gcc --version gcc (GCC) 4.3.2 Copyright (C) 2008 Free Software Foundation, Inc. This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
環境変数CC を書き換えて使いましょう。
CASAVAで用いるSampleSheetには空白を入れてはいけない
題の通り。このことはよく知られていると思うんだけど、どのエラーが出た時にサンプルシートの空白を疑ったらいいのかなかなかわからない。僕の環境では以下のエラーが出たから、もし検索でここが引っかかった人はサンプルシートを確認するといいと思う。
CASAVA 1.8.2 configureBclToFastqにて
*** [Temp/L003_R2_demux_summary.xml] エラー 1 make: *** [Temp/L007_R1_demux_summary.xml] エラー 1 make: *** [Temp/L007_R2_demux_summary.xml] エラー 1 make: *** [Temp/L002_R1_demux_summary.xml] エラー 1 make: *** [Temp/L002_R2_demux_summary.xml] エラー 1 make: *** [Temp/L005_R1_demux_summary.xml] エラー 1 make: *** [Temp/L005_R2_demux_summary.xml] エラー 1 make: *** [Temp/L004_R1_demux_summary.xml] エラー 1 make: *** [Temp/L004_R2_demux_summary.xml] エラー 1 make: *** [Temp/L001_R1_demux_summary.xml] エラー 1 make: *** [Temp/L001_R2_demux_summary.xml] エラー 1 Command exited with non-zero status 2
しかし、空白が入っただけで動かなくなるなんて、大昔のロケットの制御プログラムみたいだね。
bam_header_read EOF marker is absent. The input is probably truncated.
これも海外から見る人が多そうなので英語も書いておこう。
$ samtools index foo.sorted.bam [bam_header_read] EOF marker is absent. The input is probably truncated. [bam_header_read] invalid BAM binary header (this is not a BAM file). Segmentation fault (core dump)
.bamファイルが壊れているらしい。
.samはまともなはずなのにいつからこうなったか。
What broke foo.sorted.bam?
I compared the bam files before and after sorting.
$samtools flagstat foo.bam 8000000 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates 7910957 + 0 mapped (98.89%:-nan%) 8000000 + 0 paired in sequencing 4000000 + 0 read1 4000000 + 0 read2 7859822 + 0 properly paired (98.25%:-nan%) 7899729 + 0 with itself and mate mapped 11228 + 0 singletons (0.14%:-nan%) 15109 + 0 with mate mapped to a different chr 14075 + 0 with mate mapped to a different chr (mapQ>=5) $samtools flagstat foo.sorted.bam [bam_header_read] EOF marker is absent. The input is probably truncated. [bam_header_read] invalid BAM binary header (this is not a BAM file). [bam_flagstat_core] Truncated file? Continue anyway. 0 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates 0 + 0 mapped (-nan%:-nan%) 0 + 0 paired in sequencing 0 + 0 read1 0 + 0 read2 0 + 0 properly paired (-nan%:-nan%) 0 + 0 with itself and mate mapped 0 + 0 singletons (-nan%:-nan%) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)
Sorting process would have had trouble.
Although the cause of this phenomenon is still unclear, resorting bam fixed this trouble.
どうもsortのときに問題が起きているらしいな。sortのログを見ても特におかしなことは起きていないのだけれども。
結局sortをかけ直したらなおった。
$samtools sort foo.bam foo.sorted $ samtools flagstat foo.sorted.bam 8000000 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 duplicates 7910957 + 0 mapped (98.89%:-nan%) 8000000 + 0 paired in sequencing 4000000 + 0 read1 4000000 + 0 read2 7859822 + 0 properly paired (98.25%:-nan%) 7899729 + 0 with itself and mate mapped 11228 + 0 singletons (0.14%:-nan%) 15109 + 0 with mate mapped to a different chr 14075 + 0 with mate mapped to a different chr (mapQ>=5)
fastx インストール中に困ったことなど
pipe_fitter.c: In function ‘pipe_close’: pipe_fitter.c:30:6: error: variable "i" set but not used [-Werror=unused-but-set-variable]
Makefile を見ると-Wunused(使われない変数があったら警告を出す)と-Werror(警告があったらエラーを出す)が両方立っているので、片方を消せば通る。ただ、Makefileはconfigureで生成されているのでconfigureから直すこと。
fastx_uncollapser.cpp:31:39: Fatal error: gtextutils/stream_wrapper.h: No such file or directory
fastx toolkit はlibgtextutilsを必要としているから、fastxのサイトにはlibgtextutilsもおいてあるのだけれど、それだけでは足りないみたい。上記stream_wrapper.hというのはlibgtextutils-develの中にあるらしく、それを落としてくれば良い。ただ、パッケージ管理ソフト使える人ならそっちの方が早いよ。
# yum install libgtextutils-devel
あとはインストール先を確認して、PKG_CONFIG_PATHに指定してやれば良い。
$ rpm -ql libgtextutils-devel /usr/include/gtextutils /usr/include/gtextutils/gtextutils /usr/include/gtextutils/gtextutils/container_join.h /usr/include/gtextutils/gtextutils/exit_manip.h /usr/include/gtextutils/gtextutils/inbuf1.hpp /usr/include/gtextutils/gtextutils/natsort.h /usr/include/gtextutils/gtextutils/outbuf3.hpp /usr/include/gtextutils/gtextutils/pipe_fitter.h /usr/include/gtextutils/gtextutils/stream_wrapper.h /usr/include/gtextutils/gtextutils/string_tokenize.h /usr/include/gtextutils/gtextutils/strnatcmp.h /usr/include/gtextutils/gtextutils/text_line_reader.h /usr/include/gtextutils/gtextutils/tuple_parser.h /usr/lib64/libgtextutils.so /usr/lib64/pkgconfig/gtextutils.pc $export PKG_CONFIG_PATH=/usr/include/gtextutils/gtextutils/:/usr/include/gtextutils/
参考にさせていただいた記事：cygwinでfastx toolkitをインストール - バイオ系研究室PC管理担当のメモ
http://d.hatena.ne.jp/hornistyf/20120401/1333273058
.sam fileのbitwise flagについて "0x2 each segment properly aligned according to the aligner" 途中まで
0x02each segment properly aligned がどういう意味なのか分からないので調べている。
sam formatについての公式文書があまりにも頼りないので、samtoolsでいろいろとフィルターをかけながら調べてみた。
リード自身がunmappedなものとmateがunmappedなものについてはこのbitは意味のない値になるので除外しておく。
properly aligned を抜き出してくるフィルター
$samtools view -bh -f 0x02 foo.bam | samtools view -bh -F 0x04 - | samtools view -h -F 0x08 -
unproperly aligned を抜き出してくるフィルター
$ samtools view -bh -F 0x02 foo.bam | samtools view -bh -F 0x04 - | samtools view -h -F 0x08 -
これをベースにして両者の違いを調べ、bitflag0x02の意味を明らかにしていきたい。が、最初に言っておくと明らかにしきれていない。ごめんなさい。
まずはmateがどの染色体にmapされているかどうか調べる。異なる染色体にmapされていたら、properとは言えないだろう。
上記フィルターにさらに
awk '/^[^@]/{print $7}'|sort | uniq
をかましてmateのいる染色体を調べる。
properly:"="
unproperly:"1,10,11,12,13,14,15,16,17,18,19,2,3,4,5,6,7,8,9,=,X,Y"
each properly で異所的アラインメントがフィルターされていることは分かる。しかしフィルターされた同所的アラインメントはいったい何がunproperlyだったのだろうか。
まず考えられるのがアラインメントの向きである。
awk '/^[^@]/&&$7=="="{print $2}'|sort | uniq
proper
163:128,32,2,1
147:128,16,2,1
99:64,32,2,1
83:64,16,2,1
ここで、
128 the last segment in the template
64 the first segment in the template
32 SEQ of the next segment in the template being reversed
16 SEQ being reverse complemented
2 each segment properly aligned according to the aligner
1 template having multiple segments in sequencing
であるため、ペアの両方がリバース、ペアの両方がフォワードというリード対は取り除かれていることが分かる。
さて、unproperly paired の方はどうか。
177:128,32,16,1
161:128,32,1
145:128,16,1
129:128,1
113:64,32,16,1
97:64,32,1
81:64,16,1
65:64,1
177,129,113,65はアラインメントの向きが揃ってしまっているため不適と分かる。
他にペアエンド間の距離も調べたところ、properの方は2~500程度に収まっているのに対して、unproperでは数万に及ぶものもあった。しかしそれらを除外してもなお0.1%程度unproperなリードが残ってしまう。
ひとまずこの辺であきらめ、備忘録的にこいつを公開しておきます。誰かの役に立つかもしれないし（主に僕の役にしかたたないだろうけど）
ちなみにこのsamはbwaの割と最近のバージョン（2010年以降）で生成されました。人のデータをもらって解析してるから詳しいことは知りませんが、大してかわらんでしょう。
計画を立てよう
今週末はbwaのアルゴリズムの理解に費やそうと決めました。土曜日の午後から仕事の手を止めます。
