??最近一直忙期末匯報(bào)和腳本編寫，沒來的及接著往下寫文章惜姐，年前把這一塊寫了剪撬，然后再往下的分析流程就比較特異性了，做一步寫一步那種漓穿，相對于這種已經(jīng)常規(guī)化的流程（除了一些細(xì)節(jié)上的差異嗤军，別的都是大同小異的了），再往下的可能就不是很常規(guī)的分析手段了晃危，不同的實(shí)驗(yàn)室有不同的分析方法叙赚，希望能夠有大佬提點(diǎn)意見，多多交流僚饭。
“突變”區(qū)別
??對于WGS的數(shù)據(jù)震叮，在處理完bam文件以后，就是call variations了鳍鸵，之后所有的工作其實(shí)都是針對variations進(jìn)行分析苇瓣，那么說到變異呢，其實(shí)中文的“變異”在英文里對應(yīng)多個(gè)單詞偿乖，經(jīng)郴髯铮看到傻傻分不清，這里稍微區(qū)別一下贪薪。
??mutation：核苷酸序列的永久性改變（來源于ACMG）媳禁，在人群中小于1%或5%（來源于北大生物信息學(xué)公開課）
??polymorphism：在人群中頻率超過1%（來源于ACMG），或超過5%（來源于北大生物信息學(xué)公開課）
??variation/variant：以上兩個(gè)的總和（來源于北大生物信息學(xué)公開課）
變異類型
??我們一般所說的variations主要有四大類：SNV画切，INDEL竣稽，SV，CNV。
????SNV即單核苷酸位點(diǎn)變異（single nucleotide variants）
????INDEL即小片段插入缺失（insertion and deletion）
????SV即結(jié)構(gòu)變異（structural variation）
????CNV即拷貝數(shù)變異（copy number variation）
SNV與SNP
??說到SNV丧枪，可能大家也經(jīng)常看到SNP（single nucleotide polymorphism）庞萍，同樣是混淆使用拧烦，這倆其實(shí)還是有一些區(qū)別的，看到上面對變異的區(qū)分钝计，其實(shí)也能看出這倆的區(qū)別來：
??一般SNP是二態(tài)的恋博，SNV沒有這樣的限制，如果在一個(gè)物種中該單堿基變異的頻率達(dá)到一定水平（1%）就叫SNP私恬，而頻率未知（比如僅僅在一個(gè)個(gè)體中發(fā)現(xiàn)）就叫SNV债沮，SNV包含SNP。
變異簡述
??人基因組通常有4.1~5.0M的變異本鸣，但是99.9%都是由SNV和short indel造成疫衩。
??通常，一個(gè)人全基因組內(nèi)會有約 3.6~4.4 M 個(gè) SNVs荣德，絕大數(shù)（大于 95%）的高頻（群體中等位基因頻率大于 5%）的 SNP 在 dbSNP中有記錄闷煤，高頻的SNP一般都不是致病的主要突變位點(diǎn)。
??通常涮瞻，一個(gè)人全基因組內(nèi)會有約 600K 的 Indel（<50bp的插入缺失為small indel）鲤拿。
??編碼區(qū)或剪接位點(diǎn)處發(fā)生的插入缺失都可能會改變蛋白的翻譯。
SV
??結(jié)構(gòu)變異指的是在基因組上一些大的結(jié)構(gòu)性的變異署咽，比如大片段丟失（deletion）近顷、大片段插入（insertion）、大片段重復(fù)（duplication）宁否、拷貝數(shù)變異（copy number variants）窒升、倒位（inversion）、易位（translocation）家淤。一般來說异剥，結(jié)構(gòu)變異涉及的序列長度在1kb到3Mb之間。結(jié)構(gòu)變異普遍存在于人類基因組中絮重，是個(gè)人差異和一些疾病易感性的來源冤寿。結(jié)構(gòu)變異還可能導(dǎo)致融合基因的發(fā)生，一些癌癥已經(jīng)證實(shí)和結(jié)構(gòu)變異導(dǎo)致的基因融合事件有關(guān)青伤。
CNV
??拷貝數(shù)變異指的是基因組上大片段序列拷貝數(shù)的增加或者減少督怜，可分為缺失（deletion）和重復(fù)(duplication)兩種類型，是一種重要的分子機(jī)制狠角。CNV能夠?qū)е旅系聽栠z傳病與罕見疾病, 同時(shí)與包括癌癥在內(nèi)的復(fù)雜疾病相關(guān)号杠，因此對于染色體水平的缺失、擴(kuò)增的研究已經(jīng)成為疾病研究熱點(diǎn)。
以上數(shù)據(jù)在不同的數(shù)據(jù)庫或文獻(xiàn)上可能有所差異姨蟋，但相去不遠(yuǎn)屉凯，具體的變異分類以及分布就不詳述，可參考一些文獻(xiàn)眼溶，下面說說怎么進(jìn)行其中的SNV和indel的檢測悠砚。
??以下我將提供兩種分析方式的腳本（DeepVariant以及bcftools）和流程，為啥沒有GATK堂飞？因?yàn)槲矣X得黃老師的這篇GATK分析流程寫得已經(jīng)很好了灌旧，大家可以參考一下。
Bcftools
??與舊版的samtools+bcftools不同绰筛，作者為了避免bcftools和samtools的版本不同導(dǎo)致的不兼容枢泰，新版的bcftools可以自己完成call snv/indel的工作。
??使用bcftools進(jìn)行變異檢測铝噩，一般分為三部曲衡蚂，分別為三個(gè)模塊mpileup、call薄榛、filter讳窟，當(dāng)然，我們一般也不會分三步進(jìn)行操作敞恋，而是使用管道（pipeline）進(jìn)行編寫腳本丽啡，這樣能減少產(chǎn)生一些不必要的過程文件，同時(shí)提高自動化和效率硬猫，下面是我的實(shí)際使用腳本补箍。
$ bcftools mpileup --threads 12 -q 20 -Q 20 -Ou -f /your/path/of/reference /your/bamfile/after/sorted.merged.markdup | bcftools call --threads 12 -vm -Ov | bcftools filter --threads 12 -s FILTER -g 10 -G 10 -i "%QUAL>20 && DP>6 && MQ>50 && (DP4[2]+DP4[3])>4" > raw.tmp.vcf
## bcftools mpileup檢測變異;
# --threads線程數(shù)
# -q表示reads比對質(zhì)量選擇，MAPQ啸蜜，默認(rèn)0;
# -Q表示reads堿基對質(zhì)量選擇坑雅，默認(rèn)13
# -O表示輸出格式,u表示未壓縮bcf格式;
# -f參考序列位置
## bcftools call參數(shù)
# --threads線程;
# -v只輸出變異位點(diǎn);
# m為克服-c調(diào)用模型中已知的局限性(與-c沖突)而設(shè)計(jì)的多等位和罕見變異調(diào)用的替代模型;
# -O輸出文件格式v未壓縮vcf;
## bcftools filter篩選變異;-i只保留后面條件的;-s對不符合的變異打上標(biāo)簽
$ awk -F "\t" '{if($1~/#/){print}else if($7~/PASS/){print}}' raw.tmp.vcf > var.flt.vcf
# 將標(biāo)為低質(zhì)量的變異去掉
??有幾個(gè)點(diǎn)值得討論一下，首先是mpileup的-q參數(shù)衬横，這個(gè)和之前提到的samtools view的-q是一樣的裹粤，前面的文章有大篇幅說過這個(gè)MAPQ的數(shù)值，20翻譯過來的意思其實(shí)就是比對正確率99%蜂林。
??filter中的-s是軟過濾的意思遥诉，就是把不符合后面條件的variations打上標(biāo)簽，但不過濾掉噪叙；-g矮锈，-G這對參數(shù)是說，indel附近的indel或snp是不準(zhǔn)確的睁蕾，大多是假陽性苞笨，過濾掉债朵，這里我設(shè)的10bp，這個(gè)值還是比較合理的瀑凝；-i是保留后面符合條件的變異序芦，剛好和-e相反，兩者選其一粤咪，我這里用的-i編寫過濾表達(dá)式
????QUAL：基于Phred格式的表示ALT的質(zhì)量芝加，也可以理解為可靠性；可以理解為所call出來的變異位點(diǎn)的質(zhì)量值射窒。Q=-10lgP，Q表示質(zhì)量值将塑；P表示這個(gè)位點(diǎn)發(fā)生錯(cuò)誤的概率脉顿。因此，如果想把錯(cuò)誤率從控制在90%以上点寥，P的閾值就是1/10艾疟，那lg（1/10）=-1，Q=（-10）*（-1）=10敢辩。同理蔽莱，當(dāng)Q=20時(shí)，錯(cuò)誤率就控制在了0.01戚长。這個(gè)參數(shù)其實(shí)和mpileup的-Q是重復(fù)的盗冷。
????DP是堿基的覆蓋深度，一般很多公司和課題組會選擇10同廉，但我這里選擇的是6仪糖，本著寬進(jìn)嚴(yán)出的原則，保留更多的陽性variations迫肖，10也是沒有關(guān)系的锅劝。
????參考上一篇里的MAPQ分布可以看到剥纷，使用bwa mem后，其實(shí)大多數(shù)的reads的MAPQ都在60呢铆，這樣的話晦鞋，MQ最好也就是60，但是在40有一個(gè)小突起，對于MQ的篩選悠垛，大家的選擇可以酌情而定线定，50是一個(gè)大家使用較多的閾值點(diǎn)。
????DP4分別是正反鏈上REF和ALT的深度确买，我用“DP4[2]+DP4[3]>4”篩選ALT的深度至少是5的variations斤讥。
??那么，以上的腳本其實(shí)是符合單個(gè)樣本進(jìn)行分析的湾趾，但如果是家系樣本進(jìn)行分析芭商，其實(shí)是有問題的，因?yàn)樵赾all這一步的時(shí)候用-v參數(shù)只輸出變異位點(diǎn)搀缠，在獲得了三個(gè)人的vcf文件進(jìn)行merge的時(shí)候铛楣，對于一些變異（這些變異只存在于三個(gè)人的某一個(gè)或某兩個(gè)），你就不知道不存在的那個(gè)人身上艺普，是因?yàn)闊o變異還是沒有覆蓋到reads簸州。這不利于做trio分析。那么要克服這個(gè)問題只需要去掉-v參數(shù)即可歧譬。
??然后進(jìn)行merge岸浑，是在完成了一個(gè)家系的call snp/indel以后。
$ bgzip -c -f -@ 12 var.flt.vcf > var.vcf.gz
# 進(jìn)行文件壓縮;-c不改變內(nèi)容; -f強(qiáng)制輸出瑰步，存在就覆蓋;-@線程數(shù)
$ bcftools index -t var.vcf.gz
# 建立索引矢洲，merge需要; -t建立tbi格式索引
$ bcftools merge -Ov --force-samples -l file.list -o merge.var.vcf
# 進(jìn)行合成，-O輸出文件格式,v表示vcf格式文件;
# --force-samples缩焦，對于重名樣本強(qiáng)制合成;-o輸出文件
# -l包含文件名的文件兵钮，一行一個(gè)文件名，將先證者放在第一位
??之后可以用bcftools對結(jié)果做一下統(tǒng)計(jì)處理
## 統(tǒng)計(jì)結(jié)果plot-vcfstats
$ bcftools stats -F/your/path/of/reference -s - merge.var.vcf > merge.var.vcf.stats && \
$ plot-vcfstats merge.var.vcf.stats -p vars_output
# plot-vcfstats程序在bcftools下的misc目錄中
DeepVariant
??谷歌提供了分別適用于WGS和WES的腳本舌界，可供大家參考掘譬。我親測了一下這個(gè)軟件，速度有點(diǎn)慢……完全沒有谷歌自己描述的那么快呻拌，做為嘗鮮吧葱轩，把當(dāng)時(shí)的腳本放上來，這里要特別感謝師姐的指導(dǎo)藐握！雖然最后我也沒打算用這個(gè)軟件完成我的課題吧靴拱。
????https://github.com/google/deepvariant/tree/r0.7/scripts
??Deepvariant無需安裝，直接拉docker下來就行猾普，不然要是安裝袜炕，這個(gè)環(huán)境配置怕是要折磨死人的。
# docker安裝初家，僅針對ubuntu用戶
$ sudo apt-get install apt-transport-https ca-certificates curl software-properties-common
$ curl -fsSL https://download.docker.com/linux/ubuntu/gpg | sudo apt-key add -
$ sudo apt-key fingerprint 0EBFCD88
$ sudo add-apt-repository "deb [arch=amd64] https://download.docker.com/linux/ubuntu $(lsb_release -cs) stable"
$ sudo apt-get update
$ sudo apt-get install docker-ce
$ sudo docker run hello-world
$ sudo usermod -aG docker $USER
??拉取一下docker就可以使用deepvariant了
$ sudo docker pull dajunluo/deepvariant
$ sudo docker images
$ sudo docker run -it dajunluo/deepvariant:latest
$ sudo docker run -it --name deepvariant -v /you/path/of/refrence/dir/:/home/ref_hg19 -v /home/biowork/:/home/biowork dajunluo/deepvariant
# -v是把我本地的數(shù)據(jù)對接到docker環(huán)境里偎窘，這個(gè)請自行根據(jù)需要對接
$ nohup python /home/bin/make_examples.zip --mode calling --ref /you/path/of/refrence --reads //you/path/of/bam --examples example.gz > 1.log &
$ nohup python /home/bin/call_variants.zip --outfile call_variants_output.gz --examples example.gz --checkpoint /home/models/model.ckpt > 2.log &
$ nohup python /home/bin/postprocess_variants.zip --ref /you/path/of/refrence --infile call_variants_output.gz --outfile output.vcf.gz > 3.log &
??之后便可以用bcftools或者gatk進(jìn)行merge乌助，對于這一步可以參考上面的gatk鏈接或者bcftools步驟。
??到了這一步就可以完了么陌知？顯然不是他托，你用上述兩個(gè)方法或者參考黃老師的腳本用gatk得到的結(jié)果，如果你仔細(xì)看仆葡，就會發(fā)現(xiàn)赏参。
??這是什么鬼，這又是什么鬼沿盅，為什么有那么多的多等位位點(diǎn)把篓，你要是去統(tǒng)計(jì)，會發(fā)現(xiàn)這樣的多等位位點(diǎn)還挺多的腰涧，所以你還不能直接過濾纸俭。對于這些變異我們一般是不會考慮嵌合的，因?yàn)槠骄?0X的WGS是沒法檢測出嵌合體的南窗。那么對于這樣的位點(diǎn)怎么去考慮分析呢？
??比較便捷的方法就是用bcftools里面的norm工具了郎楼。
$ bgzip -c -f -@ 12 merge.var.vcf > merge.var.vcf.gz
# 進(jìn)行文件壓縮;-c不改變內(nèi)容; -f強(qiáng)制輸出万伤，存在就覆蓋;-@線程數(shù)
$ bcftools index -t merge.var.vcf.gz
# 建立索引, -t建立tbi格式索引
$ bcftools norm -Ov -m-any -f /you/path/of/refrence merge.var.vcf.gz > norm.vcf
??這樣，多等位位點(diǎn)就變成了二等位位點(diǎn)了呜袁，便于后續(xù)的分析敌买。今天的內(nèi)容就到這里了，下一篇就是注釋以及各種過濾了阶界。
??水平有限虹钮，要是存在什么錯(cuò)誤請?jiān)u論指出！請大家多多批評指正膘融，相互交流芙粱，共同成長，謝謝Ｑ跤场４号稀！
