外显子（WES）panel分析基础篇

作者，Evil Genius

我们来分享一下关于WES分析部分的内容，其中内容主要包括三部分

• 基础部分，包括WES数据的call SNP、INDEL、MSI、fusion、annovar、CNV等部分，我们不仅要会每个部分，也要构建全自动化脚本，实现一键式分析出结果，包括甄别配对样本和单肿瘤样本。

• 高级注释部分，基础部分当然annovar已经做了很多的注释了，但是仍然还差得远，还需要更多的数据库注释，其中最重要的就是OncoKB的自动化注释，还有有很多人工添加的用药位点，同样的道理，要实现一键式自动化出分析结果的内容，其中还有计算TMB等内容。

• 临检报告的出具，这部分是WES核心内容，分析数据的解读以及如何根据分析得到的数据进行临床用药指导，同样，我们需要代码帮我们实现一键式出结果，对于敏感位点要报出warning，人工审核。

实际过程中接收到的样本往往只有tumor样本，其中血液样本和组织样本的过滤条件稍有不同。

这一篇我们来实现第一部分，其中分析用到的软件大家可以自行查阅。我们从拿到数据开始。

第一步：数据质控

普通质控

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

第二步，数据比对

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

这里面首先需要构建好参考基因组，通常是hg19的，还有就是-R的参数设置。示例是针对tumor样本，对于normal也一样的操作。-t参数根据服务器性能自行调整。

第三步，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数据库路径
MILLS=千人计划基础SNP位点
INDELS=千人计划基础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文件，否则后续的数据解读会出问题。

如果有normal样本需要再来一次上述分析。

第五步，HaplotypeCaller，其中参数的设置需要大家研究一下

有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，这一步的参数选择前面已经介绍过，大家需要思考的是为什么只有tumor的时候也需要对tumor数据进行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软件分析候选体细胞变异时，最好使用Panel of Normals (PON)文件，以提高变异分析准确度。

GATK开发团队认为，在生成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 -Djav