2.1. bwa建立索引文件
bwa index-a bwtsw-p zm437/home/chaim/disk/zm437/Zea_mays.AGPv4.dna.toplevel.fa
2.3 BWA的mem的使用，好用快速一步到位兵拢。
bwa mem-t24-M-P-R'@RG\tID:2447-20\tSM:2447-20\tLB:WES\tPL:Illumina'zm437../origin/fastp_out/2447-20_L1.R1.clean.fastq.gz../origin/fastp_out/2447-20_L1.R2.clean.fastq.gz>2447-20_L1.sam&
sam轉bam
samtools view -Sb -@40 new_ys_am43.sam -o new_ys_am43.bam
加頭部
java -Xmx256g -jar /public/home/zzp/biosoft/picard.jar AddOrReplaceReadGroups I=am43green_gy14.bam O=nam43green_gy14.bam LB=WGS PL=illumina ID=A00151 SM=am43green PU=am43green
java -Xmx256g -jar /public/home/zzp/biosoft/picard.jar AddOrReplaceReadGroups I=am43_yellow_gy14.bam O=nam43_yellow_gy14.bam LB=WGS PL=illumina ID=A00129 SM=am43yellow PU=am43yellow
ID = Read group identifier
每一個read group 獨有的ID，每一對reads 均有一個獨特的ID辙喂，可以自定義命名局齿；
PL = Platform
測序平臺罢浇；ILLUMINA, SOLID, LS454, HELICOS and PACBIO讼积，不區(qū)分大小寫；
SM = sample
reads屬于的樣品名烦衣；SM要設定正確篮赢，因為GATK產(chǎn)生的VCF文件也使用這個名字;
LB = DNA preparation library identifier
對一個read group的reads進行重復序列標記時，需要使用LB來區(qū)分reads來自那條lane;有時候琉挖，同一個庫可能在不同的lane上完成測序;為了加以區(qū)分启泣，
同一個或不同庫只要是在不同的lane產(chǎn)生的reads都要單獨給一個ID. 一般無特殊說明，成對兒read屬于同一庫示辈，可自定義寥茫，比如：library1
若是忘記添加read group信息還以通過 AddOrReplaceReadGroups 添加
gatk AddOrReplaceReadGroups -I .bam -O .add.bam -LB library1 -PL illumina -PU pl1 -SM name##參數(shù)-I Input file(BAM or SAM or a GA4GH url)；-O Output file(BAM or SAM)矾麻；-LB Read-Group library纱耻；-PL Read-Group platform(e.g. ILLUMINA,SOLID)；-PU Read-Group platform unit(eg. run barcode)险耀；-SM Read-Group sample name
ID str：輸入reads集ID號弄喘；LB：read集文庫名；PL：測序平臺（illunima或solid）甩牺；PU：測序平臺下級單位名稱（run的名稱）蘑志；SM：樣本名稱。
bam排序
samtools sort -@ 40 -m 100G new_ys_am43.bam ./new_ys_am43.sorted.bam
去重
java -jar /public/home/zzp/newam43/picard.jar MarkDuplicates? INPUT=new_ys_am43.sorted.bam.bam OUTPUT=new_ys_am43.sorted.bam.bam.repeatmark.bam METRICS_FILE=new_ys_am43.sorted.bam.bam.metrics
建索引
samtools index??new_ys_am43.sorted.bam.bam.repeatmark.bam
會生成bwa比對的參考基因組文件的INDEX贬派。 以及GATK所需要的dict.這個時候要把GATK的dict 該一個名稱急但，比如：mv Athaliana_447_TAIR10.fa.dict Athaliana_447_TAIR10.dict。 不然下邊GATKcall SNP 會報錯
1.1生成Gvcf文件
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar HaplotypeCaller -R Gy14_genome_v2.fa -I am43ys_gy14.reordered.sort.repeatmark.bam --native-pair-hmm-threads 56 --emit-ref-confidence GVCF -O am43ys_gy14.reordered.sort.repeatmark.g.vcf
本地跑GATK4的時候發(fā)現(xiàn)搞乏，它的參數(shù)設置中竟然沒有多線程的參數(shù)波桩，這和GATK3有點大不同，因為在GATK3中我們可以用-nt或者-nct來設定多線程運行请敦，但是GATK4卻沒有類似的參數(shù)镐躲，這是為啥呢？
早有研究者問過GATK的團隊侍筛，官方的回答簡單粗暴萤皂，沒有！如果要使用多線程來跑GATK4流程勾笆，那么就在本地節(jié)點配置好Spark敌蚜，那樣子就可以用GATK4中的Spark功能模塊（比如： HaplotypeCallerSpark桥滨、BaseRecalibratorSpark窝爪、ApplyBQSRSpark弛车、SortSamSpark、MarkDuplicatesSpark）運行蒲每，讓Spark來完成多線程》柞耍現(xiàn)在GATK4中有蠻多的應用采用了Spark的技術，大家可以使用下面的代碼查看：
/path/of/your/local/gatk4/gatk--list
跑過GATK3的同學們應該了解過邀杏，以前GATK3中的多線程功能的效果其實并不好贫奠，而且還容易出問題⊥可能也是由于這方面的原因唤崭，GATK團隊這一次在GATK4中就干脆放棄了通過自己團隊實現(xiàn)多線程的想法，直接使用現(xiàn)成的Spark來完成多線程使用的調(diào)度脖律。另外谢肾，值得一提的是在GATK4中跑并行任務的最好做法是采用WDL和Cromwell相結合的方式。
話雖如此小泉，但GATK團隊實際上還是留下了唯一的一個例外芦疏！那就是HaplotypeCaller中最消耗計算資源的模塊——pariHMM，這個是可以本地單獨多線程的微姊！通過“–native-pair-hmm-threads”這個參數(shù)來設置酸茴，它默認是4，功能有些隱蔽兢交！
HaplotypeCaller的應用有兩種做法薪捍，差別只在于要不要在中間生成一個gVCF： （1）直接進行HaplotypeCaller，這適合于單樣本配喳，或者那種固定樣本數(shù)量的情況飘诗，也就是只執(zhí)行一次HaplotypeCaller。如果增加一個樣本就得重新運行這個HaplotypeCaller界逛，而這個時候算法需要重新去讀取所有人的BAM文件昆稿，浪費大量時間； （2）每個樣本先各自生成gVCF息拜，然后再進行群體joint-genotype溉潭。gVCF全稱是genome VCF，是每個樣本用于變異檢測的中間文件少欺，格式類似于VCF喳瓣，它把joint-genotype過程中所需的所有信息都記錄在這里面，文件無論是大小還是數(shù)據(jù)量都遠遠小于原來的BAM文件赞别。這樣一旦新增加樣本也不需要再重新去讀取所有人的BAM文件了畏陕，只需為新樣本生成一份gVCF，然后重新執(zhí)行這個joint-genotype就行了仿滔。
推薦使用第二種惠毁，變異檢測不是一個樣本的事情犹芹，有越多的同類樣本放在一起joint calling結果將會越準確，而如果樣本足夠多的話鞠绰，在低測序深度的情況下也同樣可以獲得完整并且準確的結果腰埂。
1.2合并之前生成的GVCF文件到一個文件。
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar CombineGVCFs -R Gy14_genome_v2.fa --variant am43ys_gy14.reordered.sort.repeatmark.g.vcf --variant /public/home/zzp/am43_bsa/gy14/am43_yellow_gy14.g.vcf.gz --variant /public/home/zzp/am43_bsa/gy14/am43_green_gy14.g.vcf.gz -O am43_all.g.vcf.gz
1.3 GenotypeGVCFs
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar GenotypeGVCFs -R Gy14_genome_v2.fa -V am43_all.g.vcf.gz -O am43_all.vcf.gz
2.select SNP
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar SelectVariants -R Gy14_genome_v2.fa -O am43_all_snp.vcf --variant am43_all.vcf.gz --select-type-to-include SNP 2>select_snp.err
3.select indel
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar SelectVariants -R Gy14_genome_v2.fa -O am43_all_indel.vcf --variant am43_all.vcf.gz --select-type-to-include INDEL 2>select_indel.err
4.1 filter SNP（變異過濾蜈膨，硬過濾屿笼。）
質控的含義和目的是指通過一定的標準，最大可能地剔除假陽性的結果翁巍，并盡可能地保留最多的正確數(shù)據(jù)驴一。 第一種方法 GATK VQSR，它通過機器學習的方法利用多個不同的數(shù)據(jù)特征訓練一個模型（高斯混合模型）對變異數(shù)據(jù)進行質控灶壶，使用VQSR需要具備以下兩個條件： 第一蛔趴，需要一個精心準備的已知變異集，它將作為訓練質控模型的真集例朱。比如孝情，Hapmap、OMNI洒嗤，1000G和dbsnp等這些國際性項目的數(shù)據(jù)箫荡，這些可以作為高質量的已知變異集。 第二渔隶，要求新檢測的結果中有足夠多的變異羔挡，不然VQSR在進行模型訓練的時候會因為可用的變異位點數(shù)目不足而無法進行。
此方法要求新檢測的結果中有足夠多的變異间唉，不然VQSR在進行模型訓練的時候會因為可用的變異位點數(shù)目不足而無法進行绞灼。可能很多非人的物種在完成變異檢測之后沒法使用GATK VQSR的方法進行質控呈野，一些小panel低矮、外顯子測序，由于最后的變異位點不夠被冒，也無法使用VQSR军掂。全基因組分析或多個樣本的全外顯子組分析適合用此方法。
第二種方法通過過濾指標過濾昨悼。?QualByDepth（QD）：QD是變異質量值（Quality）除以覆蓋深度（Depth）得到的比值蝗锥。 FisherStrand (FS)：FS是一個通過Fisher檢驗的p-value轉換而來的值，它要描述的是測序或者比對時對于只含有變異的read以及只含有參考序列堿基的read是否存在著明顯的正負鏈特異性（Strand bias率触，或者說是差異性） StrandOddsRatio (SOR)：對鏈特異（Strand bias）的一種描述. RMSMappingQuality (MQ)：MQ這個值是所有比對至該位點上的read的比對質量值的均方根. MappingQualityRankSumTest (MQRankSum) ReadPosRankSumTest (ReadPosRankSum)?通過過濾指標過濾
如何理解硬過濾的指標和閾值的計算
考慮到SNP和Indel在判斷指標和閾值方面的思路是一致的终议，因此沒必要重復說。所以，下面我只以SNP為例子穴张，告訴大家設定閾值的思路细燎。強調(diào)一下，為了更具有通用價值陆馁，這些閾值是借用NA12878（來自GIAB）的高深度數(shù)據(jù)進行計算獲得的，所以如果你的數(shù)據(jù)（或者物種）相對比較特殊（不是哺乳動物）合愈，那么就不建議直接套用了叮贩，但可以依照類似的思路去尋找新閾值。
QualByDepth（QD）
QD是變異質量值（Quality）除以覆蓋深度（Depth）得到的比值佛析。這里的變異質量值就是VCF中QUAL的值——用來衡量變異的可靠程度益老，這里的覆蓋深度是這個位點上所有含有變異堿基的樣本的覆蓋深度之和，通俗一點說寸莫，就是這個值可以通過累加每個含變異的樣本（GT為非0/0的樣本）的覆蓋深度（VCF中每個樣本里面的DP）而得到捺萌。舉個例子：
11429249?.CT1044.77?. ?.GT:AD:DP:GQ:PL?0/1:48,15:63:99:311,0,1644?0/0:47,0:47:99:392,0,0?1/1:0,76:76:99:3010,228,0
這個位點是1:1429249，VCF格式膘茎，但我把FILTER和INFO的信息省略了桃纯，它的變異質量值QUAL=1044.77。我們可以從中看到一共有三個樣本披坏，其中一個是雜合變異（GT=0/1）态坦，一個純合的非變異（GT=0/0），最后一個是純合的變異（GT=1/1）棒拂。每個樣本的覆蓋深度都在其各自的DP域上伞梯，分別是：63，47和76帚屉。按照定義谜诫，這個位點的QD值就應該等于質量值除以另外兩個含有變異的樣本的深度之和（排除中間GT=0/0這個不含變異的樣本），也就是:
QD?=?1044.77?/ (63+76) =?7.516
QD這個值描述的實際上就是單位深度的變異質量值攻旦，也可以理解為是對變異質量值的一個歸一化喻旷，QD越高一般來說變異的可信度也越高。在質控的時候牢屋，相比于QUAL或者DP（深度）來說掰邢，QD是一個更加合理的值。因為我們知道伟阔，原始的變異質量值實際上與覆蓋的read數(shù)目是密切相關的辣之，深度越高的位點QUAL一般都是越高的，而任何一個測序數(shù)據(jù)皱炉，都不可避免地會存在局部深度不均的情況怀估，如果直接使用QUAL或者DP都會很容易因為覆蓋深度的差異而帶來有偏的質控結果。
在上面『執(zhí)行硬過濾』的例子里面，我們看到認為好的SNP變異多搀，QD的值不能低于2歧蕉，但問題是為什么是2，而不是3或者其它數(shù)值呢康铭？
要回答這個問題惯退，我們可以通過利用NA12878 VQSR質控之后的變異數(shù)據(jù)和原始的變異數(shù)據(jù)來進行比較，并把它說明白从藤。
首先催跪，我們可以先把所有變異位點的QD值都提取出來，然后畫一個密度分布圖（Y軸代表的是對應QD值所占總數(shù)的比例夷野，而不是個數(shù)）懊蒸，看看QD值的總體分布情況（如下圖，來自NA12878的數(shù)據(jù)）悯搔。
從這個圖里骑丸，我們可以看到QD的范圍主要集中在0~40之間。同時妒貌，我們可以明顯地看到有兩個峰值（QD=12和QD=32）通危。這兩個峰所反映的恰恰是雜合變異和純合變異的QD值所集中的地方。這里大家可以思考一下灌曙，哪一個是代表雜合變異的峰黄鳍，哪一個是代表純合變異的峰呢？
回答是平匈，第一個峰（QD=12）代表雜合框沟，而第二峰（QD=32）代表純合，為什么呢增炭？因為對于純合變異來說忍燥，貢獻于質量值的read是雜合變異的兩倍，同樣深度的情況下隙姿，QD會更大梅垄。對于大多數(shù)的高深度測序數(shù)據(jù)來說，QD的分布和上面的情況差不多输玷，因此這個分布具有一定的代表性队丝。
接著，我們同時畫出VQSR之后所有可信變異（FILTER=Pass）和不可信變異的QD分布圖欲鹏，如下机久，淺綠色代表可信變異的QD分布圖，淺紅色代表不可信變異的QD分布圖赔嚎。
你可以看到膘盖，大多數(shù)Fail的變異胧弛，都集中在左邊的低QD區(qū)域，而且紅波峰恰好是QD=2的地方侠畔，這就是為什么硬過濾時設置QD>2的原因了结缚。
可是在上面的圖里，我想你也看到了软棺，有很多Fail的變異它的QD還是很高的红竭，有些甚至高于30，通過這樣的硬過濾參數(shù)所得到的結果中就會包含這部分本該要過濾掉的壞變異喘落；而同樣的茵宪，在低QD（<2）區(qū)域其實也有一些是好的變異，但是卻被過濾掉了揖盘。這其實也是硬過濾的一大弊端眉厨，它不會像VQSR那樣锌奴，通過多個不同維度的數(shù)據(jù)構造合適的高維分類模型兽狭，從而能夠更準確地區(qū)分好和壞的變異，而僅僅是一刀切鹿蜀。
當你理解了上面有關QD的計算和閾值選擇的過程之后箕慧，要弄懂后面的指標和閾值也就容易了，因為用的也都是同樣的思路茴恰。
FisherStrand（FS）
FS是一個通過Fisher檢驗的p-value轉換而來的值颠焦，它要描述的是測序或者比對時對于只含有變異的read以及只含有參考序列堿基的read是否存在著明顯的正負鏈特異性（Strand bias，或者說是差異性）往枣。這個差異反應了測序過程不夠隨機伐庭，或者是比對算法在基因組的某些區(qū)域存在一定的選擇偏向。如果測序過程是隨機的分冈，比對是沒問題的圾另，那么不管read是否含有變異，以及是否來自基因組的正鏈或者負鏈雕沉，只要是真實的它們就都應該是比較均勻的集乔，也就是說，不會出現(xiàn)鏈特異的比對結果坡椒，F(xiàn)S應該接近于零扰路。
這里多說一句，在VCF的INFO中有時除了FS之外倔叼，有時你還會看到SB或者SOR汗唱。它們實際上是從不同的層面對鏈特異的現(xiàn)象進行描述。只不過SB給出的是原始各鏈比對數(shù)目丈攒，而FS則是對這些數(shù)據(jù)做了精確Fisher檢驗渡嚣；SOR原理和FS類似，但是用了不同的統(tǒng)計檢驗方法計算差異程度，它更適合于高覆蓋度的數(shù)據(jù)识椰。
與QD一樣绝葡，我們先來看一下質控前所有變異的FS總體密度分布圖（如下）。很明顯與QD相比腹鹉，F(xiàn)S的范圍更加的大藏畅，從0到好幾百的都有。不過從圖中也可以看出功咒，絕大部分的變異還是在100以下的愉阎。
下面這一個圖則是經(jīng)過VQSR之后，畫出來的FS分布圖力奋。跟上面的QD一樣榜旦，淺綠色代表好變異，淺紅色代表壞變異景殷。我們可以看到溅呢，大部分好變異的FS都集中在0~10之間，而且壞變異的峰值在60左右的位置上猿挚，因此過濾的時候咐旧，我們把FS設置為大于60。其實設置這么高的一個閾值是比較激進的（留下很多假變異）绩蜻，但是從圖中你也可以看到铣墨，不過設置得多低，我們總會保留下很多假的變異办绝，既然如此我們就干脆選擇盡可能保留更多好的變異伊约，然后祈禱可以通過『執(zhí)行硬過濾』里其他的閾值來過濾掉那些無法通過FS過濾的假變異。
StrandOddsRatio（SOR）
關于SOR在上面講到FS的時候孕蝉，我就在注釋里提及過了屡律。它同樣是對鏈特異（Strand bias）的一種描述，但是從上面我們也可以看到FS在硬過濾的時候并不是非常給力昔驱，而且由于很多時候read在外顯子區(qū)域末端的覆蓋存在著一定的鏈特異（這個區(qū)域的現(xiàn)象其實是正常的）疹尾，往往只有一個方向的read，這個時候該區(qū)域中如果有變異位點的話骤肛，那么FS通常會給出很差的分值纳本，這時SOR就能夠起到比較好的校正作用了。計算SOR所用的統(tǒng)計檢驗方法也與FS不同腋颠，它用的是symmetric odds ratio test繁成，數(shù)據(jù)是一個2×2的列聯(lián)表（如下），公式也十分簡單淑玫，我把公式進行了簡單的展開巾腕，從中可以清楚地看出面睛，它考慮的其實就是ALT和REF這兩個堿基的read覆蓋方向的比例是否有偏，如果完全無偏尊搬，那么應該等于1叁鉴。
sor?= (ref_fwd/ref_rev) / (alt_fwd/alt_rev) = (ref_fwd * alt_rev) / (ref_rev * alt_fwd)
OK，那么同樣的佛寿，我們先看一下這個值總體的密度分布情況（如下）幌墓。總的來說冀泻，這個分布明顯集中在0~3之間常侣，這也和我們的預期比較一致。不過也有比較明顯的長尾現(xiàn)象弹渔，這個時候我們也沒必要定下太過明確的閾值預期胳施，先看VQSR的分布結果。
下面這個圖就是在VQSR之后肢专，區(qū)分了好和壞變異之后舞肆，SOR的密度分布。很明顯鸟召，好的變異基本就在1附近胆绊。結合這個分布圖氨鹏，我們在上面的例子里把它的閾值定為3基本上也不會過損失好的變異了欧募，雖然單靠這個閾值還是會保留下不少假的變異，但是至少不合理的長尾部分可以被砍掉仆抵。
RMSMappingQuality（MQ）
MQ這個值是所有比對至該位點上的read的比對質量值的均方根（先平方跟继、再平均、然后開方镣丑，如下公式）舔糖。
它和平均值相比更能夠準確地描述比對質量值的離散程度。而且有意思的是莺匠，如果我們的比對工具用的是bwa mem金吗，那么按照它的算法，對于一個好的變異位點趣竣，我們可以預期摇庙，它的MQ值將等于60。
下面是所有未過濾的變異位點上MQ的密度分布圖遥缕∥捞唬基本上就只在60的地方存在一個很瘦很高的峰〉ハ唬可以說這是目前為止這幾個指標中圖形最為規(guī)則的了夕凝，在這個圖上宝穗，我們甚至就可以直接定出MQ的閾值了，比如所有小于50的就可以過濾掉了码秉。
但是逮矛，理性告訴我們還是要看一下VQSR的對比結果（下圖）。
你會發(fā)現(xiàn)似乎所有好的變異都緊緊集中在60旁邊了转砖，其它地方就都是假的變異了橱鹏，所以MQ的閾值設置為50也是合理的。但是同樣要注意到的地方是堪藐，60這個范圍實際上依然有假的變異位點在那里莉兰，我們把這個區(qū)域放大來看，如下圖礁竞，這里你就會發(fā)現(xiàn)其實假變異的密度分布圖也覆蓋到60這個范圍了糖荒。
考慮到篇幅的問題，接下來MappingQualityRankSumTest（MQRankSum）和ReadPOSRankSumTest（ReadPOSRankSum）的閾值設定原理模捂，我不打算再細說下去捶朵，思路和上面的4個是完全一樣的。都是通過比較VQSR之后的密度分布圖狂男，最后確定了硬過濾的閾值综看。
但請不要以為這只是適用于GATK得到的變異，實際上岖食，只要我們弄懂了這些指標選擇的原因和過濾的思路红碑，那么通過任何其他的變異檢測工具也是依舊可以適用的，區(qū)別就在于GATK幫我們把這些要用的指標算好了泡垃。
同樣地析珊，這些指標也不是一成不變的，可以根據(jù)實際的情況換成其他蔑穴，或者我們自己重新計算忠寻。
Ti/Tv處于合理的范圍
Ti/Tv的值是物種在與自然相互作用和演化過程中在基因組上留下來的一個統(tǒng)計標記，在物種中這個值具有一定的穩(wěn)定性存和。因此奕剃，一般來說，在完成了以上的質控之后捐腿，還會看一下這些變異位點Ti/Tv的值是多少纵朋，以此來進一步確定結果的可靠程度。
Ti（Transition）指的是嘌呤轉嘌呤叙量，或者嘧啶轉嘧啶的變異位點數(shù)目倡蝙，即A<->G或C<->T；
Tv（Transversion）指的則是嘌呤和嘧啶互轉的變異位點數(shù)目绞佩，即A<->C寺鸥，A<->T猪钮，G<->C和G<->T。（如下圖）
另外胆建，在哺乳動物基因組上C->T的轉換比較多烤低，這是因為基因組上的胞嘧啶C在甲基化的修飾下容易發(fā)生C->T的轉變。
說了這么多Ti/Tv的比值應該是多少才是正常的呢笆载？如果沒有選擇壓力的存在扑馁，Ti/Tv將等于0.5，因為從概率上講Tv將是Ti的兩倍凉驻。但現(xiàn)實當然不是這樣的腻要，比如對于人來說，全基因組正常的Ti/Tv在2.1左右涝登，而外顯子區(qū)域是3.0左右雄家，新發(fā)的變異（Novel variants）則在1.5左右。
最后多說一句胀滚，Ti/Tv是一個被動指標趟济，它是對最后質控結果的一個反應，我們是不能夠在一開始的時候使用這個值來進行變異過濾的咽笼。
indel
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar VariantFiltration -V am43_all_indel.vcf --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O am43_all_indel.filed.vcf
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar VariantFiltration -V 9930_am43_all_gy14_indel.vcf --filter-expression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -8.0 || QUAL < 30.0" --filter-name "PASS" -O 9930_am43_all_gy14_indel.filed.vcf
snp
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar VariantFiltration -V 9930_am43_all_gy14_snp.vcf --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || QUAL < 30.0" --filter-name "PASS" -O 9930_am43_all_gy14_snp.filed.vcf
QD,描述單位深度的變異值顷编，越大可信度越高。一般過濾掉<2的值剑刑。
FS媳纬，描述正負鏈特異性，差異性較大叛甫，說明測序或組裝的過程中不夠隨機层宫。FS越小越好杨伙。一般過掉掉>40（嚴格）或60（普通）
MQ 使用bwa-mem的話其监，正常值應該是60，描述某個位點測序reads的質量值的離散程度限匣。MQ< 40.0
MQRankSum < -12.5
SOR抖苦，也是表示正負鏈特異性，正常值在0-3米死，過濾掉>3的值锌历。
重新合并過濾后的SNP和Indel
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar MergeVcfs -I am43_all_snp.filed.vcf -I am43_all_indel.filed.vcf -O am43_all.filed.vcf
5轉換vcf為table格式
java -Xmx256G -jar /public/home/zzp/biosoft/gatk-4.1.9.0/gatk-package-4.1.9.0-local.jar VariantsToTable -R Gy14_genome_v2.fa -V am43_all.filed.vcf -F CHROM -F POS -F REF -F ALT -GF AD -GF DP -GF GQ -GF PL -O am43_all.filed.vcf.table
三、使用qtl-seqr包（可以和第四步同時進行）
#安裝使用QTL-seqr#安裝installed.packages('devtools')library('devtools')install_github("bmansfeld/QTLseqr")#load the packagelibrary("QTLseqr")#set sample and file name HighBulk<-"樣品名稱1"LowBulk<-"樣品名稱2"file<-"common.table"#choose which chromosomes will be included in the analysis chroms<-paste0Chroms<-paste0(rep("Chr",10),1:10)#Import SNP data from filedf<-importFromGATK(file=file,highBulk=HighBulk,lowBulk=LowBulk,chromList=Chroms)#Filter SNPs based on some criteriadf_filt<-filterSNPs(SNPset=df,refAlleleFreq=0.20,minTotalDepth=100,maxTotalDepth=400,minSampleDepth=40,minGQ=99)#Run G' analysisdf_filt<-runGprimeAnalysis(SNPset=df_filt,windowSize=1e6,outlierFilter="deltaSNP")#Run QTLseq analysisdf_filt<-runQTLseqAnalysis(SNPset=df_filt,windowSize=1e6,popStruc="F2",#根據(jù)群體結構調(diào)整F2或者RILbulkSize=c(25,25),#根據(jù)池中的樣本單株數(shù)量replications=10000,intervals=c(95,99))#PlotplotQTLStats(SNPset=df_filt,var="Gprime",plotThreshold=TRUE,q=0.01)plotQTLStats(SNPset=df_filt,var="deltaSNP",plotIntervals=TRUE)#export summary CSVgetQTLTable(SNPset=df_filt,alpha=0.01,export=TRUE,fileName="my_BSA_QTL.csv")
四. 使用snpEff處理gatk4輸出的vcf峦筒。（此步驟和第三步驟的輸出結果可以互相對比究西，兩者具有相同的功能）
解壓unzip snpEff_latest_core.zip
我的路徑是/home/chaim/bsa/snpEff/
配置玉米zm437版本的數(shù)據(jù)庫
在路徑/home/chaim/bsa/snpEff/snpEff/目錄下創(chuàng)建文件夾data,
cd /home/chaim/bsa/snpEff/snpEff/mkdir datacd datamkdir genomesmkdir zm437#在zm437目錄存放基因組注釋文件genes.gff,? 蛋白庫，protein.fa#在genomes目錄放置基因組參考序列 zm437.fa
注意上述的基因組注釋文件是gff3格式物喷。
修改snpEff.config的參數(shù)
添加如下內(nèi)容
#maize genome,version zm437 zm437.genome:maize
回到snpEFF目錄卤材，運行命令
java -jar snpEff.jar build -gff3 -v zm437
對vcf格式文件進行注釋：
bwa目錄存放著GATK4處理之后的文件common_filtration.vcf
在/home/chaim/bsa/bwa目錄執(zhí)行下面命令
java -Xmx128g -jar /home/chaim/bsa/snpEff/snpEff/snpEff.jar zm437 common_filtration.vcf > common.eff.vcf
會輸出三個文件遮斥，
snpEff_genes.txt
snpEff_summary.html
common.eff.vcf
如果想更改使用其他注釋文件，
刪除/home/chaim/bsa/snpEff/snpEff//data/zm437/snpEffectPredictor.bin該文件刪除即可扇丛。
重新從步驟1開始即可术吗。
2.1. bwa建立索引文件
