肿瘤WES项目实战演练分享及学习班通知 / 四六文摘

日前在朋友圈发布了非常好的一个肿瘤WES项目测序文件的质控谍照，如下：

很多人留言感兴趣这个是哪家公司测序的，但事实上，这是一个公共数据，来自于2018年发表的单细胞DNA测序数据文章，我只是纯粹感兴趣就找到了众多问过我 GATK问题的粉丝，嘱咐他们试一下这个数据集的肿瘤wes分析，也因此结识了南方医优秀的小伙伴，在大家的热情呼唤下邀请他加入了我们VIP群：
回顾问过Jimmy的8个问题

下面是他实战该肿瘤WES数据分析项目的一点心得体会和记录，如果你坚持看到最后，请留意我们的肿瘤WES学习班通知哦！

目

录

Raw fastq QC(质控)
Alignment(比对)
Germline variants calling(胚系突变流程)

Create exon interval bed文件
BQSR
HaplotypeCaller & Joint genotype
VQSR

Somatic Variants Calling(体细胞突变流程)

Create PoN
Mutect2 calling

annovar注释
somatic的可视化结果
后记

背景

数据来源于这篇文章：Multiclonal Invasion in Breast Tumors Identified by Topographic Single Cell Sequencing。

来自10个乳腺癌患者的normal、INV、DCIS组织的exon测序结果，30个样本。

当然文章里还有其他很多数据，例如45万X的超高深度测序数据，果然又是彻底不用愁经费的团队。

我知道GATK很慢，但没想到下数据更慢。（充分说明了外网的重要性，或者高级工具aspera）

从12月12日开始下数据，学校小水管润物细无声，到了12月25日圣诞节才可以开始分析。

以下是处理过程的简要记录，包括主要步骤和参数，以及个人的浅薄理解。

如有错漏，还请大佬指正。

说明：因为30个samples，批量跑的，所以代码中很多输入参数都用变量代替。

Raw fastq QC

拿到sra，转完 fastq 先做了质控。

样本量有点多，跑完 fastqc 再用 multiQC 统一看质量。

原始数据应该是有过处理的。因为长度不统一，长度分布从80b到100b。

总的来说：还有illumina接头残留。大部分数据尾部质量一般，需要去接头，然后过滤下质量。

__________________________________________

SRR6269851 样本的质量

用 trim_galore 去接头以及过滤低质量碱基， trim_galore 利用 Cutadapt 完成去接头和低质量碱基的工作，所以环境变量里需要有cutadapt。

trim_galore --length 50 \

--stringency 5 \

-q 25 \
-e 0.1 \
--paired \
--phred33 -o $output_dir \
$read_1 $read_2

实际使用下来发现adapt去得很干净，但是质量改善效果一般。

这和trim的原理有关—— 碱基score减去阈值，然后从末端往回累加，直到累加之和大于0，这时候对应的碱基之后都会被剪掉。

对于不是很糟糕的数据，一般不会剪掉很多。不过越靠近末端剪得越多，所以尾部剪得比较干净。

__________________________________________

SRR6269851 样本trim后的质量

Alignment

接下来用 BWA mem把fastq map到参考基因组 hg38 版本。

比对结果直接通过管道传给samtools处理，节省 I／O 时间。

再用 MarkDuplicates 去重复，picard集成在GATK里了（实际上是标记没有删掉，也改成可以rm模式）。（当然，这里也可以使用高级工具，sambamba）

然后把原来的bam删掉。节省磁盘空间。

bwa mem -t 8 -M -R "@RG\tID:${RGID}\tPL:${PL}\tLB:${SM}\tSM:${SM}" \
${ref_genome} \
${read_1} \
${read_2} | samtools sort -@ 4 -m 10G -o \
${output_dir}/sm_${SM}.hg38.sort.bam
$gatk MarkDuplicates \
-I ${output_dir}/sm_${SM}.hg38.sort.bam \
-M ${output_dir}/sm_${SM}.markdup.metrics.txt \
-O ${output_dir}/sm_${SM}.markdup.bam \
&& rm ${output_dir}/sm_${SM}.hg38.sort.bam

习惯给每个sample前面加上前缀，这样方便以后的匹配处理。例如 sm_ 代表sample， T、 N 代表tumor、normal。

如果每个样本有多个bam这时候可以merge成一个bam本文件。

Germline variants calling

3.1 Create exon interval bed文件

由于是外显子测序，把关注区域聚焦到 wes region。

一开始我是从GTF提取全部exon位点的，想法是外显子测序嘛，那就全部exon区域都看。

后来想想对于那些 lncRNA 的 exon 上的 SNP，不是说没用，而是一般实验室最后关注点都是落在基因上的 variants ，那其实只提取 CDS region 也可以了。

其实我问了实验的同学，她们也没有具体说明到底外显子测序是真的测全部exon，还是只测到CDS中的外显子。她们只是说通过探针捕获目标区域再测序。不过的确是，大家都忙，没有时间和精力每个细节都亲自去研究。(其实应该是看WES芯片的bed文件)

3.2 BQSR

按照GATK的说法，测序机器对质量会有系统性的误差。

例如：假设机器在连续call了两个AA碱基后，后面的一个碱基质量score会偏高10%。那我们就可以反过来对AA后一个碱基score降低10%点。（这个是为了方便理解，不是真的这么简单）

而 variant calling 需要考虑 base quality，所以GATK推荐对碱基质量进行校正。

另外校正碱基质量不是改善碱基质量，不是让 low score -> high score 而是要让 incorrect score -> correct score。

$gatk BaseRecalibrator \
-R ${GATK_bundle}/Homo_sapiens_assembly38.fasta \
-I ${output_dir}/sm_${SM}.markdup.bam ${Interval_list} \
--known-sites ${GATK_bundle}/Homo_sapiens_assembly38.known_indels.vcf.gz \
--known-sites ${GATK_bundle}/Homo_sapiens_assembly38.dbsnp138.vcf\
--known-sites ${GATK_bundle}/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-O ${output_dir}/sm_${SM}.BQSR.table
$gatk ApplyBQSR \
--bqsr-recal-file ${output_dir}/sm_${SM}.BQSR.table \
-R ${GATK_bundle}/Homo_sapiens_assembly38.fasta \
-I ${output_dir}/sm_${SM}.markdup.bam ${Interval_list} \
-O ${output_dir}/sm_${SM}.ApplyBQSR.bam

先计算哪些需要更正，然后再调整，之后生成的bam 文件会比原bam更大，这是因为

Base recalibration adds insertion and deletion tags, so that will increase the size of your bam file. --GATK team

对于是否有必要做这一步呢，个人是持不算积极的态度的，如果时间紧张我会跳过，因为感觉这个位点85%是A和98%是A都是一个估算，而45%是A和55%是A对我来说又同样是不可信。

时间上：平均10G+ 的 clean bam文件，先第一步 BaseRecalibrator 用了 1+hr，之后 ApplyBQSR 用了 1+hr，

3.3 HaplotypeCaller & joint genotype

先生成gvcf，再联合分析方便增量分析：

aplotypecaller 就做了这么几件事：（说得好像很简单）

找出高变异的活跃区域
组装出可能的基因型
pairHMM确定每条read可能是那种基因型
最后确定sample的genotype

HP过程很耗时间，而且GATK4.0里面没有提供多线程的设置。取消了 -nct 参数设置

但其中pairHMM这一步是可以进行多线程加速的，指定 --native-pair-hmm-threads

$gatk HaplotypeCaller \
--native-pair-hmm-threads 4 \
-R ${ref_genome} \
${Interval_list} \
-I ${output_dir}/sm_${SM}.ApplyBQSR.bam \
-O ${output_dir}/gvcf/${SM}.g.vcf.gz \
-ERC GVCF

samples=(${gvcf_dir}/N*.g.vcf.gz)# 给每个normal sample前加了 N 方便匹配
# 多个--variant参数输入可以通过shell元组解决
$gatk CombineGVCFs \
-R ${ref_genome} \
${samples[@]/#/--variant } \
-O ${gvcf_dir}/cohort_${num}_g.vcf.gz
$gatk GenotypeGVCFs \
-R ${ref_genome} \
-V ${gvcf_dir}/cohort_${num}_g.vcf.gz \
-O ${vcf_output}/unfilter_cohort_${num}.vcf.gz

这部分用的时间：13+ 个小时

3.4 VQSR

上步得到的VCF需要进行过滤，降低假阳性。

VQSR要提供很多已知数据信息，以及繁琐的 -an 参数。

SNP 和 INDEL 可以分开过滤，也可以一起过滤。如 -mode BOTH

$gatk VariantRecalibrator \
-R ${ref_genome} \
-V ${vcf_output}/unfilter_cohort_${num}.vcf.gz \
--resource hapmap,known=false,training=true,truth=true,prior=15.0:${GATK_bundle}/hapmap_3.3.hg38.vcf.gz \
--resource omni,known=false,training=true,truth=false,prior=12.0:${GATK_bundle}/1000G_omni2.5.hg38.vcf.gz \
--resource 1000G,known=false,training=true,truth=false,prior=10.0:${GATK_bundle}/1000G_phase1.snps.high_confidence.hg38.vcf.gz \
--resource dbsnp,known=true,training=false,truth=false,prior=2.0:${GATK_bundle}/Homo_sapiens_assembly38.dbsnp138.vcf \
--resource mills,known=false,training=true,truth=true,prior=15.0:${GATK_bundle}/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
-an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \
-mode BOTH \
-O ${work_dir}/log/${num}_SNP_INDEL.VQSR.recal \
--tranches-file ${work_dir}/log/${num}_SNP_INDEL.VQSR.tranches \
--rscript-file ${work_dir}/log/${num}_SNP_SNP_INDEL.VQSR.plots.R

最后会用R画结果图，所以需要提取装好 R 和 ggplot2

Somatic Variants Calling

4.1 Create PoN

对 somatic variants calling 朴素的理解是：从检测到的variants里尽可能地去掉germline 、common variants。

在进行somatic variants前，需要先用mutect的tumor模式创建 panel of normal。

mutect其实包含了HP过程在里头，所以也是耗时间大户。

$gatk Mutect2 \
-R ${ref_genome} \
-I ${work_dir}/output/sm_N${sm}.ApplyBQSR.bam \
-tumor ${sm} \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
-L ${Interval} \
-O ${PoN_dir}/sample/${sm}.vcf.gz

disable-read-filter 选项避免把没有properly mapped的reads丢掉。

创建Panel of normal，同样利用shell元组解决多个输入的问题。

samples=(${PoN_dir}/sample/*.vcf.gz)
$gatk CreateSomaticPanelOfNormals \
${samples[@]/#/-vcfs } \
-O ${PoN_dir}/PoN.vcf.gz

4.2 Mutect2 calling

$gatk Mutect2 \
-R ${ref_genome} \
-I ${work_dir}/output/sm_INV_${sm}.ApplyBQSR.bam \
-I ${work_dir}/output/sm_N_${sm}.ApplyBQSR.bam \
-tumor ${tumorsampletag} \
-normal ${tumorsampletag} \
-pon ${PoN_dir}/PoN.vcf.gz \
--germline-resource ${germline_resource} \
--af-of-alleles-not-in-resource 0.0000025 \
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \
-L ${Interval} \
-O ${somatic_vcf}/vcf/INV.${sm}.somatic_unfilter.vcf.gz \

${tumorsampletag} tumorsampletag 分别对应normalsample和tumorsample bam文件header里的sample记录信息。

${germline_resource} 来自 gnomAD resource，包含了～200K exomes测序的 allele-sepcific population frequency信息。作为参考的population germline AF。

${Interval} 是WES interval 的bed文件。

和germline过程一样，需要过滤。

首先 GetPileupSummaries 收集每个突变位点上ref 和 alt 的支持reads 数量。

calculatecontamination 估计样本cross contaminated的程度。这里的样本交叉是指不同sample间的污染不是normal和tumor间的污染。

$gatk GetPileupSummaries \
-I ${work_dir}/output/sm_${sm}/*.ApplyBQSR.bam \

-V ${resource} \

-L ${Interval} \
-O ${filter_dir}/table/${sm}.getpileupsummaries.table
$gatk CalculateContamination \
-I ${filter_dir}/table/${sm}.getpileupsummaries.table \
-O ${filter_dir}/table/${sm}.calculatecontamination.table
sm1=${sm%%_*}.${sm##*_}
$gatk FilterMutectCalls \
-V ${work_dir}/somatic/vcf/${sm1}.somatic_unfilter.vcf.gz \
--contamination-table ${filter_dir}/table/${sm}.calculatecontamination.table \
-O ${filter_dir}/vcf/${sm}.filter.vcf.gz