肿瘤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

6. $gatk MarkDuplicates \

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

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

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

10.        && 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

8. $gatk ApplyBQSR \

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

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

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

12.    -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

3. $gatk CalculateContamination \

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

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

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

7. $gatk FilterMutectCalls \

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

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

10.    -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 的具体使用可参考往期 ：

6

somatic的可视化结果

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

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

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

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

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

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

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

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

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

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

7

后记

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

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

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

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

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

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

8

点评及学习班通知

关于肿瘤外显子项目数据分析实战的方方面面我都发布过零零散散的流程，而且遍布多个公众号平台及论坛、博客。这次安排任务给了十多个朋友，认真做的大概就6个，最后坚持下来，勉强走完这个肿瘤WES项目的就南医大的吴同学，也算是不容易，虽然中间很多质控过程是缺失的，比如查看病人样本是否混淆，正常对照是否混入肿瘤组织，各个质控步骤的QC可视化及统计表格，后期多款somatic突变工具的比较，再者是CNV相关分析以及公共数据库结合解释生物学现象都没有。

但毕竟是从头到尾完成了，是值得鼓励的！

最后就是广告时间啦，肿瘤外显子数据分析实战学习班，2天课程，学费1.2万每人！

有需求的朋友直接微信联系我本人即可：看招学徒  （这个里面是微信）