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だけを取りたい場合は
```sh: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が低い物を除けば良い、

```sh: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 で簡単にできる

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

#bam をsortしてindexing

これもいつも通り

```sh: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毎に処理するなどして工夫する

```sh: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に変換する

```sh:conv_vcf.sh
bcftools  convert foo.bcf > foo.vcf
```

##デフォルト条件でcallingする

```sh: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情報をいれられる
以下例

```txt: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であれば以下の様な感じ

```txt:mono_ploidy.txt
* * * * 1
```

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

```sh:call_Vonly_bcftools.sh
bcftools  call \
-m \
-v \
--ploidy-file ./ploidy.txt \
foo.vcf > foo.called.v.vcf
```

##フィルタリングを行う
vcffilterを使ってdepthが浅いリードやマップクオリティが低い部分を除く

```sh: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で圧縮する

```sh:vcffilter.sh
bgzip foo.vcf 
```

tabixでindexをつける

```sh: tabix.sh
tabix foo.vcf.gz
```


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

```sh:conct.sh
vcf-concat *vcf.gz > allout.gz
mv allout.gz allout.vcf.gz
```

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

```sh: tabix2.sh
tabix allout.vcf.gz
```


##ほしい情報だけ取り出す
vcf-query で好きなフォーマットに変換できるので、R用に変換する

```sh: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では計算法が違うので
数値が異なる
＜作成中＞