SNPデータからvcfを作って、手動でalt read のdepthなどを計算する方法

はじめに

解析に用いたサンプルについて、核間や、strain内で多型がある可能性がでたので
そのような1/1に分離しないSNPをざっくり検証してみる。

使用ツールのバージョン

uniq-hitだけを取り出す

bwaの場合

bwaを使った場合はmultiple-hit のリードについては、mapqは計算されない（＝０になる, ２つアライメントされている場合は3らしい）。ちなみにマップが全くない場合は255 らしい。

uniq hit の場合は以下のタグがつけられる
'XT:A:U'
なのでこれを使ってuniq-hitだけを取りたい場合は

uniq_bam.sh

samtools view reads.bam | grep 'XT:A:U' | samtools view -bS -T referenceSequence.fa - > reads.uniqueMap.bam

これ以外にNHを使った方法がある
http://yuifu.github.io/remove-multi-reads/

あと　　'XT:A:U'と　mapq = 0が共存する場合もあるらしいので（それぞれの算定方法が違うため）注意
https://www.biostars.org/p/18128/

bwa-mem/bowtie/bowtie2の場合

bwa-memではこのタグはもう使われない。
またbowtie系はもともとこのタグをつかわない

なので
そのときは以下の様にしてmapqが低い物を除けば良い、

uniq_bam.bowtie.sh

samtools view -b -q 4 foo.bam > foo.filtered.bam

-q int でint以下のmapqのリードを除ける

https://www.biostars.org/p/76893/
https://gist.github.com/davfre/8596159

PCR duplicate を除く

samtools で簡単にできる

rmPCR.sh

samtools rmdup -s --reference ../ref.fa \
bwa.sh.s.filteredQ10.bam \
bwa.sh.s.filteredQ10.PCR_rem.bam

bam をsortしてindexing

これもいつも通り

sort_index.bam.sh

samtools sort -@ 20 bwa.sh.s.filteredQ10.PCR_rem.bam > bwa.sh.s.filteredQ10.PCR_rem.s.bam
samtools index -@ 20  bwa.sh.s.filteredQ10.PCR_rem.s.bam

samtoolsでコーリングを行う

mpilup

時間がかかるのでcontig毎に処理するなどして工夫する

mpileup_multi.sh

# !/bin/sh
  
# $ -S /bin/bash
# $ -cwd
# $ -v PATH
# $ -t 1-500

samtools mpileup -r unitig_${SGE_TASK_ID} \ #-r contig名を指定することで、その部分だけ処理する
-f ../Ref.fa \
-u bwa.sh.s.bam > bwa.sh.s.bam.unitig_${SGE_TASK_ID}.bcf #bcfファイルが出来る

vcfに変換する

conv_vcf.sh

bcftools  convert foo.bcf > foo.vcf

デフォルト条件でcallingする

call_bcftools.sh

bcftools  call \
-m \
--ploidy-file ./ploidy.txt \
foo.vcf > foo.called.vcf

-mで マルチアレルやrare-variantに対応したコーリングを行う。
alternative model for multiallelic and rare-variant calling (conflicts with -c)

--ploidy-fileには、contigｔごとのploydy情報をいれられる
以下例

ploidy.txt

X 1 60000 M 1
X 2699521 154931043 M 1
Y 1 59373566 M 1
Y 1 59373566 F 0
MT 1 16569 M 1
MT 1 16569 F 1
*  * *     M 2
*  * *     F 2

タブかスペース区切りで
contig名　start end sex 0/1/2
といれていく
sexを指定すれば、bam毎に（サンプル毎に）挙動を変えられる。
＊がワイルドカードとして使える

monoploidであれば以下の様な感じ

mono_ploidy.txt

* * * * 1

-vオプションを使うとvariantが検出されたサイトだけを出力する

call_Vonly_bcftools.sh

bcftools  call \
-m \
-v \
--ploidy-file ./ploidy.txt \
foo.vcf > foo.called.v.vcf

フィルタリングを行う

vcffilterを使ってdepthが浅いリードやマップクオリティが低い部分を除く

vcffilter.sh

vcffilter -f "DP > 50 & DP < 300 &  MQ > 30" foo.vcf > foo.filtered.vcf

以上の例ではcalling後のdepthが 50以上300以下でmap quality が30以上のものを選択した

結果をマージ

contig毎に処理した結果をvcftoolsを使ってマージする

まずbgzipで圧縮する

vcffilter.sh

bgzip foo.vcf 

tabixでindexをつける

tabix.sh

tabix foo.vcf.gz

vcf-concat(vcftoolsの一部）でマージする

conct.sh

vcf-concat *vcf.gz > allout.gz
mv allout.gz allout.vcf.gz

再度インデックスをつける

tabix2.sh

tabix allout.vcf.gz

ほしい情報だけ取り出す

vcf-query で好きなフォーマットに変換できるので、R用に変換する

conct.sh

vcf-query allout.vcf.gz   -f '%CHROM,%POS,%REF,%ALT,%QUAL,%%INFO/DP,%INFO/DP4\n' >  allout.vcf.gz.txt

DPはdepthでDP4はその詳細を以下のコンマ区切りで与える

1. forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles

qualityをパスしたものだけがカウントされているが、DPとDP4では計算法が違うので
数値が異なる
＜作成中＞