肿瘤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的可视化结果

• 后记

0

背景

数据来源于这篇文章：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，批量跑的，所以代码中很多输入参数都用变量代替。

1

Raw fastq QC

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

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

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

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

__________________________________________

SRR6269851 样本的质量

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

1. trim_galore --length 50 \

1.    --stringency 5 \

      -q 25 \

2.    -e 0.1 \

3.    --paired \

4.    --phred33 -o $output_dir \

5.    $read_1 $read_2

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

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

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

__________________________________________

SRR6269851 样本trim后的质量

2

Alignment

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

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

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

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

1. bwa mem -t 8 -M -R "@RG\tID:${RGID}\tPL:${PL}\tLB:${SM}\tSM:${SM}" \

2.    ${ref_genome} \

3.    ${read_1} \

4.    ${read_2} | samtools sort -@ 4 -m 10G -o \

5.    ${output_dir}/sm_${SM}.hg38.sort.bam

7. $gatk MarkDuplicates \

8.        -I ${output_dir}/sm_${SM}.hg38.sort.bam \

9.        -M ${output_dir}/sm_${SM}.markdup.metrics.txt \

10.        -O ${output_dir}/sm_${SM}.markdup.bam \

11.        && rm ${output_dir}/sm_${SM}.hg38.sort.bam

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

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

3

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。

1. $gatk BaseRecalibrator \

2.    -R ${GATK_bundle}/Homo_sapiens_assembly38.fasta \

3.    -I ${output_dir}/sm_${SM}.markdup.bam ${Interval_list} \

4.    --known-sites ${GATK_bundle}/Homo_sapiens_assembly38.known_indels.vcf.gz \

5.    --known-sites ${GATK_bundle}/Homo_sapiens_assembly38.dbsnp138.vcf\

6.    --known-sites ${GATK_bundle}/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \

7.    -O ${output_dir}/sm_${SM}.BQSR.table

9. $gatk ApplyBQSR \

10.    --bqsr-recal-file ${output_dir}/sm_${SM}.BQSR.table \

11.    -R ${GATK_bundle}/Homo_sapiens_assembly38.fasta \

12.    -I ${output_dir}/sm_${SM}.markdup.bam ${Interval_list} \

13.    -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 就做了这么几件事：（说得好像很简单）

1. 找出高变异的活跃区域

2. 组装出可能的基因型

3. pairHMM确定每条read可能是那种基因型

4. 最后确定sample的genotype

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

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

1. $gatk HaplotypeCaller \

2.    --native-pair-hmm-threads 4 \

3.    -R ${ref_genome} \

4.    ${Interval_list} \

5.    -I ${output_dir}/sm_${SM}.ApplyBQSR.bam \

6.    -O ${output_dir}/gvcf/${SM}.g.vcf.gz \

7.    -ERC GVCF

1. samples=(${gvcf_dir}/N*.g.vcf.gz)# 给每个normal sample前加了 N 方便匹配

2. # 多个--variant参数输入可以通过shell元组解决

3. $gatk CombineGVCFs \

4.        -R ${ref_genome} \

5.        ${samples[@]/#/--variant } \

6.        -O ${gvcf_dir}/cohort_${num}_g.vcf.gz

7. $gatk GenotypeGVCFs \

8.    -R ${ref_genome} \

9.    -V ${gvcf_dir}/cohort_${num}_g.vcf.gz \

10.    -O ${vcf_output}/unfilter_cohort_${num}.vcf.gz

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

3.4 VQSR

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

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

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

1. $gatk VariantRecalibrator \

2.    -R ${ref_genome} \

3.    -V ${vcf_output}/unfilter_cohort_${num}.vcf.gz \

4.    --resource hapmap,known=false,training=true,truth=true,prior=15.0:${GATK_bundle}/hapmap_3.3.hg38.vcf.gz \

5.    --resource omni,known=false,training=true,truth=false,prior=12.0:${GATK_bundle}/1000G_omni2.5.hg38.vcf.gz \

6.    --resource 1000G,known=false,training=true,truth=false,prior=10.0:${GATK_bundle}/1000G_phase1.snps.high_confidence.hg38.vcf.gz \

7.    --resource dbsnp,known=true,training=false,truth=false,prior=2.0:${GATK_bundle}/Homo_sapiens_assembly38.dbsnp138.vcf \

8.    --resource mills,known=false,training=true,truth=true,prior=15.0:${GATK_bundle}/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \

9.    -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR \

10.    -mode BOTH \

11.    -O ${work_dir}/log/${num}_SNP_INDEL.VQSR.recal \

12.    --tranches-file ${work_dir}/log/${num}_SNP_INDEL.VQSR.tranches \

13.    --rscript-file ${work_dir}/log/${num}_SNP_SNP_INDEL.VQSR.plots.R

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

4

Somatic Variants Calling

4.1 Create PoN

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

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

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

1. $gatk Mutect2 \

2.    -R ${ref_genome} \

3.    -I ${work_dir}/output/sm_N${sm}.ApplyBQSR.bam \

4.    -tumor ${sm} \

5.    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \

6.    -L ${Interval} \

7.    -O ${PoN_dir}/sample/${sm}.vcf.gz

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

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

1. samples=(${PoN_dir}/sample/*.vcf.gz)

2. $gatk CreateSomaticPanelOfNormals \

3.    ${samples[@]/#/-vcfs } \

4.    -O ${PoN_dir}/PoN.vcf.gz

4.2 Mutect2 calling

1. $gatk Mutect2 \

2.    -R ${ref_genome} \

3.    -I ${work_dir}/output/sm_INV_${sm}.ApplyBQSR.bam \

4.    -I ${work_dir}/output/sm_N_${sm}.ApplyBQSR.bam \

5.    -tumor ${tumorsampletag} \

6.    -normal ${tumorsampletag} \

7.    -pon ${PoN_dir}/PoN.vcf.gz \

8.    --germline-resource ${germline_resource} \

9.    --af-of-alleles-not-in-resource 0.0000025 \

10.    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter \

11.    -L ${Interval}  \

12.    -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间的污染。

1. $gatk GetPileupSummaries \

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

      -V ${resource} \

1.    -L ${Interval} \

2.    -O ${filter_dir}/table/${sm}.getpileupsummaries.table

4. $gatk CalculateContamination \

5.    -I ${filter_dir}/table/${sm}.getpileupsummaries.table \

6.    -O ${filter_dir}/table/${sm}.calculatecontamination.table

7. sm1=${sm%%_*}.${sm##*_}

8. $gatk FilterMutectCalls \

9.    -V ${work_dir}/somatic/vcf/${sm1}.somatic_unfilter.vcf.gz \

10.    --contamination-table ${filter_dir}/table/${sm}.calculatecontamination.table \

11.    -O ${filter_dir}/vcf/${sm}.filter.vcf.gz

过滤完的VCF中会有对应的 PASS 标记，可以提取行。

1. awk -F '\t' '{if($0 ~ /\#/) print; else if($7 == "PASS") print}'

5

annovar注释

annovar是个perl脚本；用学校教育邮箱即可以注册获取下载链接

这里只注释到gene上，没有注释到SNP信息

1. ./table_annovar.pl ${sm} ./hg38_humandb/ \

2.    -buildver hg38 -out $n \

3.    -remove -protocol refGene,cytoBand \

4.    -operation g,r -nastring . -vcfinput

产生三种文件： 

• annovar的avinput文件

• 注释后的vcf

• 注释后的结果txt

其中txt和vcf内容一致。

最后R中用 maftools可视化。 maftools 的具体使用可参考往期 ：

肿瘤突变数据可视化神器-maftools

6

somatic的可视化结果

与之前处理过的肝癌和膀胱癌症相比，发现乳腺癌每个患者的突变数量要少，如果算TMB应该也不高。

这里看到 patient 1 、8 的突变数目最少。

乳腺癌中最多的是C>A突变，其次C>T突变，不同癌症的情况也会不同。（这样的突变模式是研究的热点，我在直播我的基因组活动中多次强调）

【直播】我的基因组46:SNV突变(96种)频谱的制作

SNV突变(96种)频谱的制作

【直播】我的基因组 45：SNV突变(6种)频谱的制作

除开1、8号patient，同一个患者的INV和DCIS样本都存在一致的高频突变，从这点出发也可以旁证下肿瘤的亚克隆进化。

这里是只展示突变率前30的突变，所以容易给人错觉：patient 8 没有突变。但其实只是patient 8 的突变恰好不是高频突变而已。

频突变基因的词云。不得不说maftools真是便利。

最后还有不同基因的突变事件是否显著共发或者互斥。

一趟下来，拿到了germline的vcf和somatic的vcf，也对somatic VCF做了注释和可视化。

当然后续还有CNV没做。而且这只是一次基础流程的实践，偏向于简略地记录个人步骤和浅簿理解。

中途有缺漏谬误，请大家多多指正。

7

后记

接触GATK已经有2个多月了。越来越觉得一套分析流程，要跑下来没报错很容易，难点在于理解为什么，以及针对自己的数据贴身修改。还有的是要保证结果的可信度。

最担心是程序虽然没报错，但拿到不可靠的结果。

比如参数不合理，数据输入混淆，参考注释版本不同、未标化当成以及标准化等等，这些地方弄错，程序没有报错，最后自己也没发现。

另外越来越觉得 数理统计、编程能力和对生物学问题的敏感性 是始终绕不开的核心。

果然别人的经验是不会骗你的，至于各类分析流程，应该看成随着任务不同而不断换用的工具。

混好生信和湿实验一样，得问题出发、思路先行，再结合工具方法回答问题。

8