作者遍蟋，Evil Genius

我們來分享一下關(guān)于WES分析部分的內(nèi)容，其中內(nèi)容主要包括三部分

• 基礎(chǔ)部分发绢，包括WES數(shù)據(jù)的call SNP硬耍、INDEL、MSI边酒、fusion经柴、annovar、CNV等部分墩朦，我們不僅要會每個部分坯认，也要構(gòu)建全自動化腳本，實現(xiàn)一鍵式分析出結(jié)果氓涣，包括甄別配對樣本和單腫瘤樣本牛哺。

• 高級注釋部分，基礎(chǔ)部分當(dāng)然annovar已經(jīng)做了很多的注釋了劳吠，但是仍然還差得遠(yuǎn)引润，還需要更多的數(shù)據(jù)庫注釋，其中最重要的就是OncoKB的自動化注釋痒玩，還有有很多人工添加的用藥位點淳附，同樣的道理，要實現(xiàn)一鍵式自動化出分析結(jié)果的內(nèi)容蠢古，其中還有計算TMB等內(nèi)容奴曙。

• 臨檢報告的出具，這部分是WES核心內(nèi)容草讶，分析數(shù)據(jù)的解讀以及如何根據(jù)分析得到的數(shù)據(jù)進(jìn)行臨床用藥指導(dǎo)洽糟，同樣，我們需要代碼幫我們實現(xiàn)一鍵式出結(jié)果，對于敏感位點要報出warning坤溃，人工審核拍霜。

實際過程中接收到的樣本往往只有tumor樣本，其中血液樣本和組織樣本的過濾條件稍有不同薪介。

這一篇我們來實現(xiàn)第一部分沉御，其中分析用到的軟件大家可以自行查閱。我們從拿到數(shù)據(jù)開始昭灵。

第一步：數(shù)據(jù)質(zhì)控

普通質(zhì)控

fastp=fastp軟件路徑
fq1=tumor外顯子read1
fq2=tumor外顯子read2
$fastp -i fq1 -I fq2 -o tumor_clean.R1.fq.gz" -O \ 
tumor_clean.R2.fq.gz" \
 -h tumor_fastp.html \
-j tumor_fastp.json -3 -5  

實際分析過程則需要更多的過濾條件

fastp=fastp軟件路徑
$fastp -i normal_1.fastq.gz -o normal_1.depumi.fastq.gz \
-I normal_2.fastq.gz -O normal_2.depumi.fastq.gz \
        -U --umi_loc=per_read --umi_len=5 --umi_skip=4 --umi_prefix=UMI \
        -Q -L -A -3 -5 -h normal.html -j fastp/normal.json


$fastp -i tumor_1.fastq.gz -o tumor_1.depumi.fastq.gz -I tumor_2.fastq.gz \
-O tumor_2.depumi.fastq.gz \        
    -U --umi_loc=per_read --umi_len=5 --umi_skip=4 --umi_prefix=UMI \        
    -Q -L -A -3 -5 -h tumor.html -j tumor.json

第二步，數(shù)據(jù)比對

bwa=bwa路徑
samtools=samtools軟件路徑
$bwa mem -t 12 ucsc.hg19.fasta參考基因組路徑 \ 
tumor_clean.R1.fq.gz tumor_clean.R2.fq.gz  \
-R "@RG\tID:tumor\tLB:tumor\tPL:ILLUMINA\tSM:tumor\tPU:ILLUMINA" \
 | $samtools view -bS - -o tumor.raw.bam

這里面首先需要構(gòu)建好參考基因組伐谈，通常是hg19的烂完，還有就是-R的參數(shù)設(shè)置。示例是針對tumor樣本诵棵，對于normal也一樣的操作抠蚣。-t參數(shù)根據(jù)服務(wù)器性能自行調(diào)整。

第三步履澳，bam文件的sort和去重

gatk=gatk軟件路徑嘶窄，推薦最新版倒退一兩個版本

$samtools sort -@ 12 -m 10G tumor.raw.bam -o tumor.sort.bam

$gatk --java-options "-Xms5g -Xmx5g -Djava.io.tmpdir=臨時路徑" MarkDuplicates \
    -I tumor.sort.bam \
    -O tumor.markdup.bam \
    -M tumor.markdup_metrics.txt

第四步，BQSR,其中bed文件一般試劑盒會提供

target=bed文件路徑
bedtools=bedtools軟件路徑
DBSNP=dbSNP數(shù)據(jù)庫路徑
MILLS=千人計劃基礎(chǔ)SNP位點
INDELS=千人計劃基礎(chǔ)indel位點

cat $target | awk -v OFS="\t" -v gap=500 '{$2=$2-gap;$3=$3+gap;print $0}' | $bedtools merge -i - > ./targetplus_bed
$gatk BaseRecalibrator \
    -R 參考基因組路徑 \
    -I tumor.markdup.bam \
    -L  ./targetplus_bed\
    --known-sites $DBSNP \
    --known-sites $MILLS \
    --known-sites $INDELS \
    -O tumor.recal_data

$gatk ApplyBQSR \
    --bqsr-recal-file tumor.recal_data \
    -L  ./targetplus_bed\
    -R 參考基因組路徑 \
    -I tumor.markdup.bam \
    -O tumor.markdup.BQSR.bam

rm ./targetplus_bed -rf

需要注意的是分析必須配備bed文件距贷，否則后續(xù)的數(shù)據(jù)解讀會出問題柄冲。

如果有normal樣本需要再來一次上述分析。

第五步忠蝗，HaplotypeCaller现横，其中參數(shù)的設(shè)置需要大家研究一下

有normal樣本

$gatk --java-options "-Xms5g -Xmx5g -Djava.io.tmpdir=臨時路徑" HaplotypeCaller \
    -R 參考基因組路徑 \
    -I normal.markdup.BQSR.bam \
    --output normal.raw.vcf \
    -L $target \
    -stand-call-conf 30.0 \
    --max-mnp-distance 5 \
    -RF GoodCigarReadFilter \
    --dont-use-soft-clipped-bases \
    --minimum-mapping-quality 20 \
    --native-pair-hmm-threads 8

只有tumor樣本

$gatk --java-options "-Xms5g -Xmx5g -Djava.io.tmpdir=臨時路徑" HaplotypeCaller \
    -R 參考基因組路徑 \
    -I tumor.markdup.BQSR.bam \
    --output tumor.raw.vcf \
    -L $target \
    -stand-call-conf 30.0 \
    --max-mnp-distance 5 \
    -RF GoodCigarReadFilter \
    --dont-use-soft-clipped-bases \
    --minimum-mapping-quality 20 \
    --native-pair-hmm-threads 8

filter，這一步的參數(shù)選擇前面已經(jīng)介紹過阁最，大家需要思考的是為什么只有tumor的時候也需要對tumor數(shù)據(jù)進(jìn)行HaplotypeCaller分析

$gatk SelectVariants \
    --select-type-to-exclude INDEL \
    -V tumor.raw.vcf \
    -O tumor.snp.vcf 
$gatk SelectVariants \
    --select-type-to-include INDEL \
    -V tumor.raw.vcf \
    -O tumor.indel.vcf 
$gatk VariantFiltration \
    -R 參考基因組路徑 \
    --filter-expression "QD < 2.0 || ReadPosRankSum < -20.0 || FS > 200.0 || SOR > 10.0" \
    --filter-name "Standard" \
    --filter-expression "DP <= 5" \
    --filter-name "DP_5" \
    --filter-expression  "DP < 20" \
    --filter-name "DP_20" \
    -V tumor.indel.vcf \
    -O tumor.indel.filter.vcf
$gatk VariantFiltration \
    -R 參考基因組路徑 \
    --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || SOR > 3.0" \
    --filter-name "Standard" \
    --filter-expression "DP <= 5" \
    --filter-name "DP_5" \
    --filter-expression  "DP < 20" \
    --filter-name "DP_20" \
    -V tumor.snp.vcf \
    -O tumor.snp.filter.vcf

第六步戒祠，mutect，af_only_gnomad這個文件是什么大家自行查閱速种，這里就不過多介紹了姜盈。

在使用GATK4 Mutect2軟件分析候選體細(xì)胞變異時，最好使用Panel of Normals (PON)文件配阵，以提高變異分析準(zhǔn)確度馏颂。

GATK開發(fā)團(tuán)隊認(rèn)為，在生成PON文件時闸餐，正常對照樣本越多越好（至少40個）饱亮，但是其實一些小樣本（4-10個）也是可以的，有總比沒有好舍沙。

only tumor

af_only_gnomad=af_only_gnomad路徑
$gatk --java-options "-Xmx5g -Xms5g -Djava.io.tmpdir=臨時路徑" Mutect2 \
    -R 參考基因組路徑 \
    -I tumor.markdup.BQSR.bam \
    -tumor tumor \
    -L $target \
    --output tumor.raw.vcf \
    --germline-resource $af_only_gnomad \
    --af-of-alleles-not-in-resource 0.00003125 \
    --f1r2-tar-gz tumor.f1r2.tar.gz \
    --min-base-quality-score 20 \
    --max-mnp-distance 5 \
    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
    --native-pair-hmm-threads 8 \
    -G StandardMutectAnnotation \
    --dont-use-soft-clipped-bases

$gatk LearnReadOrientationModel \
    -I tumor.f1r2.tar.gz \
    -O tumor.read-orientation-model.tar.gz 
  $gatk GetPileupSummaries \
    -I tumor.markdup.BQSR.bam \
    -L $small_exac_common \
    -V $small_exac_common \
    -O tumor.getpileupsummaries.table 
  $gatk CalculateContamination \
    -I tumor.getpileupsummaries.table \
    -O tumor_calculatecontamination.table 
  $gatk FilterMutectCalls \
    -V tumor.raw.vcf \
    -R 參考基因組路徑 \
    --ob-priors tumor.read-orientation-model.tar.gz \
    -contamination-table tumor_calculatecontamination.table \
    -O tumor.raw.filtered.vcf

有配對樣本

$gatk --java-options "-Xmx5g -Xms5g -Djava.io.tmpdir=臨時路徑" Mutect2 \
    -R 參考基因組路徑 \
    -I tumor.markdup.BQSR.bam \
    -tumor tumor \
    -I normal.markdup.BQSR.bam \
    -normal normal \
    -L $target \
    --output tumor.raw.vcf \
    --germline-resource $af_only_gnomad \
    --af-of-alleles-not-in-resource 0.00003125 \
    --f1r2-tar-gz tumor.f1r2.tar.gz \
    --min-base-quality-score 20 \
    --max-mnp-distance 5 \
    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
    --native-pair-hmm-threads 8 \
    -G StandardMutectAnnotation \
    --dont-use-soft-clipped-bases \
    -RF NotSupplementaryAlignmentReadFilter

后面的filter過程跟第六步就是一致的近上。

當(dāng)然這樣的清洗還是不夠，還需要對數(shù)據(jù)進(jìn)行分析調(diào)整

Mutect2_Filter=python的mutect的清洗腳本
python $Mutect2_Filter tumor.raw.filtered.vcf 

第七步拂铡，annovar壹无，數(shù)據(jù)庫filter和注釋葱绒，需要提前下載好注釋的數(shù)據(jù)庫文件,注意這里對HaplotypeCaller和mutect的分析結(jié)果均要進(jìn)行該操作,mutect文件就是tumor.raw.filtered.vcf，下面是示例斗锭。

humandb=數(shù)據(jù)庫路徑地淀，就在annovar軟件下面，常見的寫在了protocol參數(shù)下面
table_annovar=annovar軟件注釋路徑
perl $table_annovar  \
    --buildver hg19 \
    --otherinfo \
    --nastring . tumor.snp.filter.vcf $humandb \
    -protocol refGene,dbnsfp31a_interpro,avsnp150,snp138,1000g2015aug_all,exac03,clinvar_20210501,intervar_20180118,cosmic95,dbnsfp30a \
    -operation g,f,f,f,f,f,f,f,f,f \
    --vcfinput --remove

分析得到的文件示例如下：

第八步岖是，MSI

msisensor2=msisensor2軟件路徑
msisensor_pro=msisensor_pro軟件路徑
microsatellites=微衛(wèi)星位點的基線文件（軟件msisensor2）
microsatellites_pro=微衛(wèi)星位點的基線文件（軟件msisensor_pro）
modelhg19=models_hg19_GRCh37
####有normal樣本
$msisensor msi \
    -d $microsatellites \
    -n normal.BQSR.bam \
    -t tumor.BQSR.bam \
    -o tumor.MSI
####only tumor
$msisensor_pro pro -d $microsatellites_pro -t tumor.markdup.BQSR.bam -o tumor.MSI

第九步帮毁，F(xiàn)usion，這里大家要考慮為什么用factera的時候需要對數(shù)據(jù)重新比對豺撑。

factera=factera.pl路徑
ref_bit=ucsc.hg19.2bit
exons_bed=外顯子的bed文件

#####重新比對
$bwa mem -Y -t 8 參考基因組路徑 tumor_clean.R1.fq.gz tumor_clean.R2.fq.gz  -R "@RG\tID:$prefix_t\tLB:tumor\tPL:ILLUMINA\tSM:tumor\tPU:ILLUMINA" | $samtools view -bS - -o tumor.raw.bam

$samtools sort -@ 12 -m 20G tumor.raw.bam -o tumor.sort.bam
perl $factera -F -o ./ tumor.sort.bam $exons_bed $ref_bit $target

如果有normal樣本的話需要再來一次這樣的分析烈疚。

第10步，CNV

cnvkit=cnvkit軟件路徑
access_hg19=access.hg19.bed文件
###有配對樣本
python3 $cnvkit batch \
    tumor.markdup.bam \
    --normal normal.markdup.bam \
    --targets $target \
    --annotate $refFlat \
    --fasta 參考基因組數(shù)據(jù)庫 \
    --access $access_hg19 \
    --output-reference tumor.reference.cnn \
    --output-dir ./ \
    --diagram --scatter

####只有tumor樣本
reference_cnn=CNV基線文件
python3 $cnvkit batch \
    tumor.markdup.bam \
    --reference $reference_cnn \
    --output-dir ./

第11步聪轿，bamdst爷肝，測序位點的突變內(nèi)容

bamdst=bamdst軟件路徑

$bamdst -p $target tumor.markdup.BQSR.bam -o ./

當(dāng)然我們實際分析過程不可能是這樣的一步一步操作，需要一鍵式出所有結(jié)果陆错，這就需要大家有更高的代碼能力灯抛，下面提供一個示例，用法是

只有tumor樣本

bash gatk.panel.sh sample.clean.R1.fastq.gz sample.clean.R2.fastq.gz sample(樣本名稱) 輸出路徑 -t bed文件 -c CNV基線文件

有配對樣本

bash gatk.panel.sh sample.clean.R1.fastq.gz sample.clean.R2.fastq.gz sample(樣本名稱) 輸出路徑 -t bed文件  -n normal_fq1 -N normal_fq2

核心腳本gatk.panel.sh