首先創(chuàng)建一個項目目錄wes和存放每個步驟生成文件的子目錄raw, clean, align, genome, hg19_VCF
mkdir {raw,clean,align,genome,hg19_VCF}
此流程包含兩個raw外顯子測序文件：wes.1.fq.gz瘾境，wes.2.fq.gz颅湘；存放于raw目錄
原始數(shù)據(jù)質控，用fastp
# 項目主目錄下運行拗慨，在clean目錄下生成兩個過濾的fq文件
nohup fastp -i raw/wes.1.fq.gz -o clean/wes.1.clean.fq.gz -I raw/wes.2.fq.gz -O clean/wes.2.clean.fq.gz &
下載參考基因組（hg19）文件，存放于genome目錄沪袭，并建立索引：
for i in $(seq 1 22) X Y M;
do
nohup wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr${i}.fa.gz &
done
for i in $(seq 1 22) X Y M;
do cat chr${i}.fa >> hg19.fa;
done
nohup bwa index -a bwtsw -p hg19 hg19.fa &
比對
# clean目錄下運行此命令艇抠，在align目錄下生成sam文件
nohup bwa mem -t 4 -M -R "@RG\tID:lane1\tPL:illumina\tLB:library\tSM:wes" ../genome/hg19 wes.1.clean.fq.gz wes.2.clean.fq.gz >../align/wes.sam 2>wes.bwa.align.log &
-t，線程數(shù)讯榕；
-M , -M 將 shorter split hits 標記為次優(yōu)骤素，以兼容 Picard markDuplicates 軟件;
-R 接的是 Read Group的字符串信息，它是用來將比對的read進行分組的愚屁，這個信息對于我們后續(xù)對比對數(shù)據(jù)進行錯誤率分析和Mark duplicate時非常重要济竹。
(1) ID，這是Read Group的分組ID霎槐，一般設置為測序的lane ID
(2) PL送浊，指的是所用的測序平臺
(3) SM，樣本ID
(4) LB丘跌，測序文庫的名字
這些信息設置好之后袭景，在RG字符串中要用制表符（\t）將它們分開
samtools view -b -S wes.sam > wes.bam #sam2bam,便于存儲
samtools sort wes.bam -o wes.sorted.bam #有順序的排序唁桩，便于后面的操作
samtools flagstat wes.sorted.bam > wes.sorted.bam.flagstat #統(tǒng)計比對信息
1. QC pass的reads的數(shù)量為5002344 ，未通過QC的reads數(shù)量為0耸棒，意味著一共有5002344條reads荒澡；
2.標記為secondary的read
3. 嵌合體reads
4.PCR 重復引起的reads
5. 比對到參考基因組上的reads數(shù)量；
6. paired reads數(shù)據(jù)數(shù)量与殃；
7. read1的數(shù)量单山；
8. read2 的數(shù)量；
9. 正確地匹配到參考序列的reads數(shù)量奈籽；
10. 一對reads都比對到了參考序列上的數(shù)量饥侵，但是并不一定比對到同一條染色體上；
11. 一對reads中只有一條與參考序列相匹配的數(shù)量衣屏；
12. 一對reads比對到不同染色體的數(shù)量躏升；
13. 一對reads比對到不同染色體的且比對質量值大于5的數(shù)量。
外顯子區(qū)域覆蓋度
需要生成外顯子interval文件狼忱，生成這個文件的前提又需要dict文件和外顯子bed文件（此處用的是安捷倫外顯子bed文件膨疏，也可以去UCSC下載）。
gatk CreateSequenceDictionary -R hg19.fa -O hg19.dict
gatk BedToIntervalList -I S31285117_Regions.bed -O Exon.Interval.bed -SD ./genome/hg19.dict
gatk CollectHsMetrics -BI Exon.Interval.bed -TI Exon.Interval.bed -I wes.sorted.bam -O wes.stat.txt
標記PCR重復序列并建立索引
以前用picard標記重復序列钻弄，現(xiàn)在這個工具全部整合到gatk中了
nohup gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" MarkDuplicates -I wes.sorted.bam -O wes.sorted.MarkDuplicates.bam -M wes.sorted.bam.metrics > log.mark 2>&1 &
samtools view -f 1024 wes.sorted.MarkDuplicates.bam | less
samtools index wes.sorted.MarkDuplicates.bam
wes.sorted.bam.metrics有統(tǒng)計信息
wes.sorted.MarkDuplicates.bam創(chuàng)建索引文件佃却，他的作用能夠讓我們可以隨機訪問這個文件中的任意位置，而且后面的步驟也要求這個bam文件一定要有索引窘俺。
變異檢測
重新校正堿基質量值（BQSR）
變異檢測是一個極度依賴測序堿基質量值饲帅，因為這個質量值是衡量我們測序出來的這個堿基到底有多正確的重要指標。它來自于測序圖像數(shù)據(jù)的base calling瘤泪，因此灶泵，基本上是由測序儀和測序系統(tǒng)來決定的，計算出來的堿基質量值未必與真實結果統(tǒng)一对途。
BQSR（Base Quality Score Recalibration）這個步驟主要是通過機器學習的方法構建測序堿基的錯誤率模型赦邻，然后對這些堿基的質量值進行相應的調整。
這里包含了兩個步驟：
第一步实檀，BaseRecalibrator惶洲，這里計算出了所有需要進行重校正的read和特征值，然后把這些信息輸出為一份校準表文件（wes.recal_data.table）
第二步膳犹，ApplyBQSR恬吕，這一步利用第一步得到的校準表文件（wes.recal_data.table）重新調整原來BAM文件中的堿基質量值，并使用這個新的質量值重新輸出一份新的BAM文件须床。
首先下載變異注釋文件
wget -c -r -nd -np -k -L -p ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19
下載文件均為壓縮文件币呵，需要解壓才能使用
samtools faidx hg19.fa
GENOME=./genome/hg19.fa
hg19_VCF=./hg19_VCF/
gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" BaseRecalibrator \
-R $GENOME -I ./align/wes.sorted.MarkDuplicates.bam \
--known-sites $hg19_VCF/1000G_phase1.indels.hg19.sites.vcf \
--known-sites $hg19_VCF/Mills_and_1000G_gold_standard.indels.hg19.sites.vcf \
--known-sites $hg19_VCF/dbsnp_138.hg19.vcf \
-L S31285117_Regions.bed -O wes.recal_data.table
gatk --java-options "-Xmx10G -Djava.io.tmpdir=./" ApplyBQSR \
-R $GENOME -I ./align/wes.sorted.MarkDuplicates.bam \
-bqsr wes.recal_data.table -L S31285117_Regions.bed \
-O wes.sorted.MarkDuplicates.BQSR.bam
# gatk AnalyzeCovariates -bqsr wes.recal_data.table -plots wes.recal_data.table.plot
變異檢測
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結果將會越準確蚊俺，而如果樣本足夠多的話懈涛，在低測序深度的情況下也同樣可以獲得完整并且準確的結果。
第一種方法
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" HaplotypeCaller \
-R $GENOME -I wes.sorted.MarkDuplicates.BQSR.bam \
-D $hg19_VCF/dbsnp_138.hg19.vcf -L S31285117_Regions.bed \
-O wes.raw1.vcf
第二種方法
# 1.生成中間文件gvcf
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" HaplotypeCaller -R $GENOME \
--emit-ref-confidence GVCF -I wes.sorted.MarkDuplicates.BQSR.bam \
-D $hg19_VCF/dbsnp_138.hg19.vcf -L S31285117_Regions.bed -O wes.gvcf
# 2.通過gvcf檢測變異
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" GenotypeGVCFs \
-R $GENOME -V wes.gvcf -L S31285117_Regions.bed \
-O wes.raw.vcf
若有多個樣本的gvcf文件泳猬，運行gatk CombineGVCFs \
-V 1.gvcf –V 2.gvcf …… -O final.gvcf 批钠，再用final.gvcf運行第二步
變異質控與過濾
質控的含義和目的是指通過一定的標準，最大可能地剔除假陽性的結果得封，并盡可能地保留最多的正確數(shù)據(jù)埋心。
第一種方法 GATK VQSR，它通過機器學習的方法利用多個不同的數(shù)據(jù)特征訓練一個模型（高斯混合模型）對變異數(shù)據(jù)進行質控呛每，使用VQSR需要具備以下兩個條件：
第一，需要一個精心準備的已知變異集坡氯，它將作為訓練質控模型的真集晨横。比如，Hapmap箫柳、OMNI手形，1000G和dbsnp等這些國際性項目的數(shù)據(jù)，這些可以作為高質量的已知變異集悯恍。
第二库糠，要求新檢測的結果中有足夠多的變異，不然VQSR在進行模型訓練的時候會因為可用的變異位點數(shù)目不足而無法進行。
gatk --java-options "-Xmx8G -Djava.io.tmpdir=./" VariantRecalibrator -R $GENOME -V wes.raw.vcf \
-resource hapmap,known=false,training=true,truth=true,prior=15.0:$hg19_VCF/hapmap_3.3.hg19.sites.vcf \
-resource omini,known=false,training=true,truth=false,prior=12.0:$hg19_VCF/1000G_omni2.5.hg19.sites.vcf \
-resource 1000G,known=false,training=true,truth=false,prior=10.0:$hg19_VCF/1000G_phase1.snps.high_confidence.hg19.sites.vcf \
-resource dbsnp,known=true,training=false,truth=false,prior=6.0:$hg19_VCF/dbsnp_138.hg19.vcf \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -an DP -mode SNP \
-O wes.snps.recal.vcf --tranches-file wes.snps.tranches --rscript-file wes.snps.plots.R &
gatk ApplyRecalibration -V wes.raw.vcf -O wes.VQSR.vcf --recal-file wes.snps.recal.vcf \
--tranches-file wes.snps.tranches -mode SNP
此方法要求新檢測的結果中有足夠多的變異瞬欧，不然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)
通過過濾指標過濾
使用SelectVariants须蜗，選出SNP
gatk SelectVariants -select-type SNP -V wes.raw.vcf -O wes.snp.vcf
# 為SNP作過濾
gatk VariantFiltration -V wes.snp.vcf --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O wes.snp.filter.vcf
使用SelectVariants，選出Indel
gatk SelectVariants -select-type INDEL -V wes.raw.vcf -O wes.indel.vcf
# 為Indel作過濾
gatk VariantFiltration -V wes.indel.vcf --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O wes.indel.filter.vcf
突變注釋
ANNOVAR是一個perl編寫的命令行工具肿孵，能在安裝了perl解釋器的多種操作系統(tǒng)上執(zhí)行唠粥。允許多種輸入文件格式，包括最常被使用的VCF格式停做。輸出文件也有多種格式晤愧，包括注釋過的VCF文件、用tab或者逗號分隔的txt文件蛉腌，ANNOVAR能快速注釋遺傳變異并預測其功能官份。這個軟件需要edu郵箱注冊才能下載。
http://www.openbioinformatics.org/annovar/annovar_download_form.php
數(shù)據(jù)庫下載
perl annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/
# -buildver 表示version
# -downdb 下載數(shù)據(jù)庫的指令
# -webfrom annovar 從annovar提供的鏡像下載烙丛，不加此參數(shù)將尋找數(shù)據(jù)庫本身的源
# humandb/ 存放于humandb/目錄下
其它數(shù)據(jù)庫下載
perl annotate_variation.pl --buildver hg19 --downdb gwascatalog humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb ljb26_all --webfrom annovar humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb esp6500siv2_ea --webfrom annovar humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb esp6500siv2_all --webfrom annovar humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb 1000g2014oct humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb cytoBand humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb avsift -webfrom annovar humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb snp138 humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb genomicSuperDups humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb phastConsElements46way humandb/ &
perl annotate_variation.pl --buildver hg19 --downdb tfbs humandb/
三種注釋方式
Gene-based Annotation(基于基因的注釋)：基于基因的注釋（gene-based annotation）揭示variant與已知基因直接的關系以及對其產(chǎn)生的功能性影響舅巷。
Region-based Annotation（基于區(qū)域的注釋）：基于過濾的注釋精確匹配查詢變異與數(shù)據(jù)庫中的記錄：如果它們有相同的染色體，起始位置河咽，結束位置钠右，REF的等位基因和ALT的等位基因，才能認為匹配忘蟹§浚基于區(qū)域的注釋看起來更像一個區(qū)域的查詢（這個區(qū)域也可以是一個單一的位點），在一個數(shù)據(jù)庫中媚值，它不在乎位置的精確匹配狠毯，它不在乎核苷酸的識別∪烀ⅲ基于區(qū)域的注釋（region-based annotation）揭示variant與不同基因組特定段的關系嚼松。
Filter-based Annotation（基于過濾的注釋）：filter-based和region-based主要的區(qū)別是，filter-based針對mutation（核苷酸的變化）而region-based針對染色體上的位置。如在全基因組數(shù)據(jù)中的變異頻率献酗，可使用1000g2015aug寝受、kaviar_20150923等數(shù)據(jù)庫；在全外顯組數(shù)據(jù)中的變異頻率凌摄，可使用exac03羡蛾、esp6500siv2等；在孤立的或者低代表人群中的變異頻率锨亏，可使用ajews等數(shù)據(jù)庫痴怨。
用table_annovar.pl進行注釋，可一次性完成三種類型的注釋器予。
convert2annovar.pl # 將多種格式轉為.avinput的程序
convert2annovar.pl -format vcf4 wes.snp.filter.vcf >snp.avinput
格式如下：
chr1 69511 69511 A G hom 1736.03 53
chr1 877831 877831 T C hom 85.10 4
chr1 881627 881627 G A het 399.60 37
chr1 887801 887801 A G hom 1623.03 39
chr1 888639 888639 T C hom 3991.03 106
chr1 888659 888659 T C hom 4114.03 109
chr1 897325 897325 G C hom 3998.03 113
chr1 898852 898852 C T het 647.60 38
chr1 900505 900505 G C het 478.60 35
chr1 906272 906272 A C het 540.60 30
avinput文件由tab分割浪藻，最重要的地方為前5列，分別是:
- 染色體(Chromosome)
- 起始位置(Start)
- 結束位置(End)
- 參考等位基因(Reference Allele)
- 替代等位基因(Alternative Allele)
ANNOVAR主要也是利用前五列信息對數(shù)據(jù)庫進行比對乾翔，注釋變異爱葵。
SNP注釋：
用table_annovar.pl進行注釋，可一次性完成三種類型的注釋反浓。
table_annovar.pl snp.avinput $humandb -buildver hg19 -out snpanno \
-remove -protocol refGene,cytoBand,genomicSuperDups,esp6500siv2_all \
-operation g,r,r,f -nastring . -csvout
-buildver hg19 表示使用hg19版本
-out snpanno 表示輸出文件的前綴為snpanno
-remove 表示刪除注釋過程中的臨時文件
-protocol 表示注釋使用的數(shù)據(jù)庫萌丈，用逗號隔開，且要注意順序
-operation 表示對應順序的數(shù)據(jù)庫的類型（g代表gene-based雷则、r代表region-based辆雾、f代表filter-based），用逗號隔開月劈，注意順序
-nastring . 表示用點號替代缺省的值
-csvout 表示最后輸出.csv文件
Indel注釋同上
測試數(shù)據(jù)：鏈接：https://pan.baidu.com/s/16tX5GrCLBAov9k_GjYYOyA
提取碼：9pab
