RNAseqからSNV解析 - kuroの覚え書き

さて、ちょっと変則的にRNAseqデータからsnpを検出すべく、あれこれやっているわけだが、
普通にゲノムデータでSNV解析するならBWAでマッピングしてsamtools mpileupと行くのだろうが、ここをHisat2でマッピングしてmpileupと持っていこうと思う。
で、sambambaは結局うまく動いてくれなかったので、普通にsamtoolsで順番にやっていくことにしてみた。

$ samtools mpileup -uf "$ref" "$output_dir"/"$samp"_"$exp"_sort.bam | bcftools call -c -v -o "$output_dir"/"$samp"_"$exp".vcf.txt
$ vcfutils.pl varFilter "$output_dir"/"$samp"_"$exp".vcf.txt >"$output_dir"/"$samp"_"$exp"_filter.vcf.txt

ディレクトリ構造はさておき、パイプライン的にはこんな感じで処理をしてみた。
この処理によって出来上がるvcfファイルを理解してみたい。
vcfのカラムは

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT

このような並びになっている。
vcfファイルのヘッダを見てみると

##ALT=<ID=*,Description="Represents allele(s) other than observed.">
##INFO=<ID=INDEL,Number=0,Type=Flag,Description="Indicates that the variant is an INDEL.">
##INFO=<ID=IDV,Number=1,Type=Integer,Description="Maximum number of reads supporting an indel">
##INFO=<ID=IMF,Number=1,Type=Float,Description="Maximum fraction of reads supporting an indel">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Raw read depth">
##INFO=<ID=VDB,Number=1,Type=Float,Description="Variant Distance Bias for filtering splice-site artefacts in RNA-seq data (bigger is better)",Version="3">
##INFO=<ID=RPB,Number=1,Type=Float,Description="Mann-Whitney U test of Read Position Bias (bigger is better)">
##INFO=<ID=MQB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality Bias (bigger is better)">
##INFO=<ID=BQB,Number=1,Type=Float,Description="Mann-Whitney U test of Base Quality Bias (bigger is better)">
##INFO=<ID=MQSB,Number=1,Type=Float,Description="Mann-Whitney U test of Mapping Quality vs Strand Bias (bigger is better)">
##INFO=<ID=SGB,Number=1,Type=Float,Description="Segregation based metric.">
##INFO=<ID=MQ0F,Number=1,Type=Float,Description="Fraction of MQ0 reads (smaller is better)">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##INFO=<ID=AF1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele frequency (assuming HWE)">
##INFO=<ID=AF2,Number=1,Type=Float,Description="Max-likelihood estimate of the first and second group ALT allele frequency (assuming HWE)">
##INFO=<ID=AC1,Number=1,Type=Float,Description="Max-likelihood estimate of the first ALT allele count (no HWE assumption)">
##INFO=<ID=MQ,Number=1,Type=Integer,Description="Root-mean-square mapping quality of covering reads">
##INFO=<ID=FQ,Number=1,Type=Float,Description="Phred probability of all samples being the same">
##INFO=<ID=PV4,Number=4,Type=Float,Description="P-values for strand bias, baseQ bias, mapQ bias and tail distance bias">
##INFO=<ID=G3,Number=3,Type=Float,Description="ML estimate of genotype frequencies">
##INFO=<ID=HWE,Number=1,Type=Float,Description="Chi^2 based HWE test P-value based on G3">
##INFO=<ID=DP4,Number=4,Type=Integer,Description="Number of high-quality ref-forward , ref-reverse, alt-forward and alt-reverse bases">

このように説明されている。
確からしさのスコアも多数記載されているが、この後の抽出で使いそうなところを抜き出してみる。

７つ目のカラムまでがvcfの基本情報

8つ目のカラムINFOにいろいろ記載されているらしい。
ID=INDEL　　INDELであるかどうか、違う場合はこのタグは出てこない。
ID=DP　　read depth
ID=AF1　　最も多いALTの割合（refに対してaltがホモだと１）
ID=AC1　　ALTがホモだと２、ヘテロだと１のように見えるが、、、
ID=DP4　　クオリティーの高いreadの本数、左からrefのfwd、refのrev、altのfwd、altのrevの順

9と10のカラムには

GT:PL	1/1:255,33,0

こんな感じに入っているが
９番目のGTに相当するのが10番目の1/1の部分
PLが255,33,0の部分となっている。
GTはgenotypeなのでここが1/1だと変異ホモ、0/1だとヘテロということだ。今回は変異のない部分はcallしないようにしていたので0/0の正常ホモはない。
例えばref=G alt=Aだった場合GG, GA, AAの3つのgenotypeが存在するはずだが、PLはそれぞれのもっともらしさを表しているらしく、小さいほど（0に近いほど）もっともらしい、ということらしい。なのでこの例だと3つ目が０となっているのでAAの変異ホモが最もこのポジションを表している、ということになるらしい。