VarScan somatic の出力ファイル形式
なぜか研修中はこれをpileupファイルと呼びならわしていたのだけれど、いわゆるsamtoolsが吐くようなpileupファイルとは形式が違う。
可読なファイルだし、各コラムの意味はヘッダ見ればわかるけど、念のため調べてみた。
VarScanのプロジェクトページにあったのでメモがわりに貼っておきます。
Output
http://varscan.sourceforge.net/somatic-calling.html#somatic-output
VarScan creates two output files by default, one for SNVs (.snp) and one for indels (.indel). If the --validation flag is turned on, a third file (.validation) will also be generated containing all positions that were called. Output files have headers, and all share the same format:
Field Description
chrom Chromosome or reference name
position Position from pileup (1-based)
ref Reference base at this position
var Variant base seen in tumor
normal_reads1 Reads supporting reference in normal
normal_reads2 Reads supporting variant in normal
normal_var_freq Variant allele frequency in normal
normal_gt Consensus genotype call in normal
normal_reads1 Reads supporting reference in tumor
normal_reads2 Reads supporting variant in tumor
normal_var_freq Variant allele frequency in tumor
normal_gt Consensus genotype call in tumor
somatic_status Somatic status call (Germline, Somatic, LOH, or Unknown)
variant_p_value Variant p-value for Germline events
somatic_p_value Somatic p-value for Somatic/LOH events
tumor_reads1_plus Tumor reference-supporting reads on + strand
tumor_reads1_minus Tumor reference-supporting reads on - strand
tumor_reads2_plus Tumor variant-supporting reads on + strand
tumor_reads2_minus Tumor variant-supporting reads on - strand
p-valueが何を指して(なにを帰無仮説として)いるのかいまいち判然としない。
某シーケンサーが吐き出す設定ファイルを読んでいる。
なかなか面白く、使用法をレクチャーしてもらったときにはいまいち腑に落ちなかったことがわかってくる。例えば、画像の1ピクセルは実物の何mmになるように調整されている。とか、各光学フィルターによる位置のずれはどれくらい補正されているか。とか。
しかしここで詳細を書くのはためらわれる。レクチャーの資料ですら「取り扱いには注意してください」とあり、これは部外秘ということだと受け取っていいだろう。いわんや設定ファイルをや。だ。
もちろん企業だってお金をかけて開発しているのだから、ツールの中身についておおっぴらにされたくないだろうが、これのおかげでネット上にほとんど情報がないことになってしまい、トラブルシューティングはいったんサポートセンターに問い合わせをしなければならなくなる。これが結構手間で、企業のツールからフリーウェアに乗り換えてしまう原因になっていると思う。
GATK ERROR MESSAGE: Invalid GZIP header
何だか海外から検索してやってきてらっしゃる方もいるようなので、英語でも書いておこう。
java -jar GenomeAnalysisTK.jar -T DepthOfCoverage -R foo.fa -I bar.bam -L bar.bed -omitBaseOutput -ct 10 -ct 20 & ##### ERROR ------------------------------------------------------------------------------------------ ##### ERROR A BAM ERROR has occurred (version 1.5-0-g04cafff): ##### ERROR The invalid inputs must be corrected before the GATK can proceed ##### ERROR Please do not post this error to the GATK forum until you have followed the instructions below ##### ERROR ##### ERROR Please make sure that your BAM file is well-formed by running Picard's validator on it ##### ERROR (see http://picard.sourceforge.net/command-line-overview.shtml#ValidateSamFile for details) ##### ERROR Also, please ensure that your BAM index is not corrupted: delete the current one and regenerate it with 'samtools index' ##### ERROR Visit our wiki for extensive documentation http://www.broadinstitute.org/gsa/wiki ##### ERROR Visit our forum to view answers to commonly asked questions http://getsatisfaction.com/gsa ##### ERROR ##### ERROR MESSAGE: Invalid GZIP header ##### ERROR ------------------------------------------------------------------------------------------
Check if your index file is older than the bam file itself.
In my case, this error was completely fixed after reindexing the bam file.
扱ってるファイルでgzipなのってbamだけなんだけど、samtoolsが吐いたbamの形式についてあれこれ言われてもどうしようもない。
追記：どうやらindexが古かったようなので
samtools index bar.bam
したら通った。
picard/MarkDuplicatesはオプションの書式が少々違ってもエラーを返さない。
いろいろ不審な点がある。
| raw | after MarkDuplicates | After samtools -F 12 | |
| total | 199260064 + 0 | 199260064 + 0 | 168640351 + 0 |
| duplicates | 0 + 0 | 8801525 + 0 | 4705460 + 0 |
| mapped | 177408978 + 0 (89.03%:-nan%) | 177408978 + 0 (89.03%:-nan%) | 168640351 + 0 (100.00%:-nan%) |
| paired in sequencing | 199260064 + 0 | 199260064 + 0 | 168640351 + 0 |
| read1 | 99630032 + 0 | 99630032 + 0 | 84320175 + 0 |
| read2 | 99630032 + 0 | 99630032 + 0 | 84320176 + 0 |
| properly paired | 167347120 + 0 (83.98%:-nan%) | 167347120 + 0 (83.98%:-nan%) | 167347119 + 0 (99.23%:-nan%) |
| with itself and mate mapped | 168640351 + 0 | 168640351 + 0 | 168640351 + 0 |
| singletons | 8768627 + 0 (4.40%:-nan%) | 8768627 + 0 (4.40%:-nan%) | 0 + 0 (0.00%:-nan%) |
| with mate mapped to a different chr | 999173 + 0 | 999173 + 0 | 999173 + 0 |
| with mate mapped to a different chr (mapQ>=5) | 925569 + 0 | 925569 + 0 | 925569 + 0 |
追記：MarkDuplicates の REMOVE_DUPLICATESオプションを
"=true"
ではなく
"=ture"
と書いてしまっていたらしい。そういう時こそエラー出してくれよ！
samtools でPCR Duplicates flag(1024)を除去したところ、
samtools view -hbF 1024 foo.sorted.redup.bam
samtools flagstat
| total | 190458539 + 0 | |
| duplicates | 0 + 0 | |
| mapped (88.53%:-nan%) | 168607453 + 0 | |
| paired in sequencing | 190458539 + 0 | |
| read1 | 95082294 + 0 | |
| read2 | 95376245 + 0 | |
| properly paired (85.40%:-nan%) | 162655156 + 0 | |
| with itself and mate mapped | 163934891 + 0 | |
| singletons (2.45%:-nan%) | 4672562 + 0 | |
| with mate mapped to a different chr | 990781 + 0 | |
| with mate mapped to a different chr (mapQ>=5) | 917396 + 0 |
きちんとduplicatesが消えている。
.sam fileのoption タグ
XTはUniquely mappedかRepeatedly mappedかMate rescued(repeatedだけどpairがUnique)か。
NMはEdit distanceで、SNVやIn/Delの入った数。
XMはNMのうちSNVのみ数えたものだと思うが、細かいところはまだよくわからない。
MDは変異の入った位置を示す。
MD:Z:86T0A0
だったら86塩基目まではreference通りで87,88塩基目でSNV(referenceではTA。リードでどうなっているかはSEQフィールドを見れば良い。)
In/Delの位置はCIGAR stringsフィールドに書かれている。
2M2I86M
だったら3,4塩基目にInsertionあり。
AM SM のtemplate independent mapping quality というものの正体がわからん。
CIGAR strings
- S
Smith-Waterman アルゴリズムでは、リードを最初から最後までアラインメントするわけじゃない。
両端が切れたものがアラインされることもある。そんなときにはCIGAR string "S" でタグづけするのさ。（のさ？）
In Smith-Waterman alignment, a sequence may not be aligned from the first residue to the last one. Subsequences at the ends may be clipped off. We introduce operation 'S' to describe clipped alignment.
http://davetang.org/wiki/tiki-index.php?page=SAM
