外显子（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 -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过程跟第六步就是一致的。

当然这样的清洗还是不够，还需要对数据进行分析调整

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

第七步，annovar，数据库filter和注释，需要提前下载好注释的数据库文件,注意这里对HaplotypeCaller和mutect的分析结果均要进行该操作,mutect文件就是tumor.raw.filtered.vcf，下面是示例。

humandb=数据库路径，就在annovar软件下面，常见的写在了protocol参数下面
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=微卫星位点的基线文件（软件msisensor2）
microsatellites_pro=微卫星位点的基线文件（软件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

第九步，Fusion，这里大家要考虑为什么用factera的时候需要对数据重新比对。

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 参考基因组数据库 \
    --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，测序位点的突变内容

bamdst=bamdst软件路径

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

当然我们实际分析过程不可能是这样的一步一步操作，需要一键式出所有结果，这就需要大家有更高的代码能力，下面提供一个示例，用法是

只有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