我们已经分析了 Somatic mutations，并进行了注释和可视化，接下来我们进行拷贝数变异的分析。

这里我们还是先从 GATK 的 somatic cnv 的最佳实践开始

拷贝数变异（copy number variations, CNVs）是属于基因组结构变异（structural variation, SV），是指 DNA 片 段长度在 1Kb-3Mb 的基因组结构变异。我们首先从 GATK 的 CNV 流程开始 CNV 的分析。

首先是流程图，先从 bam 文件，结合坐标文件计算每个外显子的 reads counts 数，然后 call segment，最后是画图：

11 拷贝数变异分析（GATK流程） - 图1

在下面这个链接中，给出了详细的教程：https://gatkforums.broadinstitute.org/gatk/discussion/11682#2
或者最近更新的教程，一样的：https://gatk.broadinstitute.org/hc/en-us/articles/360035531092

1. 外显子坐标的interval文件

interval 文件其实也是一个坐标文件，类似 bed，只不过 bed 文件的坐标是从 0 开始记录，而 interval 文件的坐标是从1开始记录。使用基因组 interval 文件可以定义软件分析的分辨率。如果是全基因组测序，interval 文件就用全基因组坐标的等间隔区间就好。对于外显子组的数据，我们使用捕获试剂盒的目标区域，理论上应该返回原文查找对应试剂盒，去搜索其捕获外显子的区域，一般是外显子侧翼上下游 250bp 以内。我懒得去查，就用软件的默认参数 250bp。

首先使用 BedToIntervalList 工具将 bed 转成 interval 格式（其实前面已经完成），然后用PreprocessIntervals工具获取 target 区间，即外显子侧翼上下游 250bp。

  1. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  2. ref=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
  3. bed=~/wes_cancer/data/hg38.exon.bed
  4. dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
  6. ## bed to intervals_list
  7. $GATK BedToIntervalList -I ${bed} -O ~/wes_cancer/data/hg38.exon.interval_list -SD ${dict}
  9. ## Preprocess Intervals
  10. $GATK  PreprocessIntervals \
  11. -L ~/wes_cancer/data/hg38.exon.interval_list \
  12. --sequence-dictionary ${dict} \
  13. --reference ${ref}  \
  14. --padding 250 \
  15. --bin-length 0 \
  16. --interval-merging-rule OVERLAPPING_ONLY \
  17. --output ~/wes_cancer/data/targets.preprocessed.interval.list

2. 获取样本的read counts

这里分为两小步进行：

  • 首先是获取所有样本的 read counts，用到了工具是 CollectReadCounts，其根据所提供的 interval 文件，对 bam 文件进行 reads 计数，可以简单理解为把 bam 文件转换成 interval 区间的 reads 数。最后会生成一个 HDF5 格式的文件，需要用第三方软件 HDFView 来查看，这里不做展示。这个文件记录了每个基因组interval 文件的 CONTIG，START，END 和原始 COUNT 值，并制成表格。
  • 然后构建正常样本的 CNV panel of normals，生成正常样本的 cnvponM.pon.hdf5 文件。对于外显子组捕获测序数据，捕获过程会引入一定的的噪音。因此后面需要降噪，该文件就是用于后面第 3 步 DenoiseReadCounts。

实际使用到的脚本如下：

  1. interval=~/wes_cancer/data/targets.preprocessed.interval.list
  2. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  3. ref=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
  5. cat config3 | while read id
  6. do
  7.     i=./5.gatk/${id}_bqsr.bam
  8.     echo ${i}
  9.     ## step1 : CollectReadCounts
  10.     time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"     CollectReadCounts \
  11.     -I ${i} \
  12.     -L ${interval} \
  13.     -R ${ref} \
  14.     --format HDF5  \
  15.     --interval-merging-rule OVERLAPPING_ONLY \
  16.     --output ./8.cnv/gatk/counts/${id}.clean_counts.hdf5
  18.     ## step2 : Generate a CNV panel of normals:cnvponM.pon.hdf5
  19.     $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./" CreateReadCountPanelOfNormals \
  20.     --minimum-interval-median-percentile 5.0 \
  21.     --output ./8.cnv/gatk/cnvponM.pon.hdf5 \
  22.     --input ./8.cnv/gatk/counts/case1_germline.clean_counts.hdf5 \
  23.     --input ./8.cnv/gatk/counts/case2_germline.clean_counts.hdf5 \
  24.     --input ./8.cnv/gatk/counts/case3_germline.clean_counts.hdf5 \
  25.     --input ./8.cnv/gatk/counts/case4_germline.clean_counts.hdf5 \
  26.     --input ./8.cnv/gatk/counts/case5_germline.clean_counts.hdf5 \
  27.     --input ./8.cnv/gatk/counts/case6_germline.clean_counts.hdf5 
  28. done

这里的 input 重复了 6 次，为的是想让初学者容易理解。如果追求代码整洁，可以把这 6 句改成

  1. $(for i in {1..6} ;do echo "--input ./8.cnv/gatk/counts/case${i}_germline.clean_counts.hdf5" ;done)

3. 降噪DenoiseReadCounts

该步骤用到的工具是 DenoiseReadCounts，主要是做了一个标准化和降噪，会生成两个文件 ${id}.clean.standardizedCR.tsv 和 ${id}.clean.denoisedCR.tsv，该工具会根据 PoN 的 counts 中位数对输入文件 ${id}.clean_counts.hdf5 进行一个标准化，包括 log2 转换。然后使用 PoN 的主成分进行标准化后的 copy ratios 降噪。实际上这里只需要对 tumor 样本进行降噪的，为了比较，我就把 tumor 和 normal 样本都分析了一遍，后面可以做个对比。

  1. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  2. cat config3 | while read id
  3. do
  4.     i=./8.cnv/gatk/counts/${id}.clean_counts.hdf5
  5.     $GATK  --java-options "-Xmx20g" DenoiseReadCounts \
  6.     -I ${i} \
  7.     --count-panel-of-normals ./8.cnv/gatk/cnvponM.pon.hdf5 \
  8.     --standardized-copy-ratios ./8.cnv/gatk/standardizedCR/${id}.clean.standardizedCR.tsv \
  9.     --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv
  10. done

4. 可视化降噪后的copy ratios

我们使用 PlotDenoisedCopyRatios 可视化标准化和去噪的 read counts。这些图可以直观地评估去噪的效果。

  1. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  2. dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
  4. cat config3 | while read id
  5. do
  6.     $GATK   --java-options "-Xmx20g" PlotDenoisedCopyRatios \
  7.     --standardized-copy-ratios     ./8.cnv/gatk/standardizedCR/${id}.clean.standardizedCR.tsv \
  8.     --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv \
  9.     --sequence-dictionary ${dict} \
  10.     --output ./8.cnv/gatk/cnv_plots \
  11.     --output-prefix ${id}
  12. done

对于每一个样本，如 case1_biorep_A_techrep，将生成 6 个文件：

  1. ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoised.png
  2. ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoisedLimit4.png  
  3. ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.deltaMAD.txt        
  4. ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.scaledDeltaMAD.txt
  5. ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.standardizedMAD.txt
  6. ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoisedMAD.txt

其中4个文件是文本文件，里面就是一个数字，记录几个拷贝数变化比值 copy ratio 的中位数绝对偏差（median absolute deviation, MAD）

  1. ## 标准化后的 copy ratios 的 MAD
  2. $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.standardizedMAD.txt
  3. 0.229
  4. ## 降噪后的 copy ratios 的 MAD
  5. $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.denoisedMAD.txt
  6. 0.231
  7. ## 标准化后的 MAD 和降噪后的 MAD 的差
  8. $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.deltaMAD.txt
  9. -0.002
  10. ## (降噪后的 MAD - 标准化后的 MAD ) / (标准化后的 MAD )
  11. $ cat ./8.cnv/gatk/cnv_plots/case1_biorep_A_techrep.scaledDeltaMAD.txt
  12. -0.01

另外两个文件是图片，表达同一个意思，标准化和降噪后的 copy ratios，底下还有中位数绝对偏差（median absolute deviation, MAD），只不过一张图片把 y 轴设置为 0 到 4。假如 tumor 样本的拷贝数没有发生变化，copy ratio 应该稳定在 1 附近。当然，要是发生了 CNV 事件，那应该就在 1 附近波动。

11 拷贝数变异分析（GATK流程） - 图2

5. 计算常见的germline mutation位点

这一步用到了 CollectAllelicCounts 工具，对输入的 bam 文件，根据指定的 interval 区间，进行 germline mutation 的检测（仅仅是 SNPs 位点，不包括 INDELs ），并计算该位点覆盖的 reads 数，即该位点的测序深度。值得注意的是，该工具一个默认参数是 MAPQ 值大于 20 的 reads 才会被纳入计数，最后生成一个 tsv 文件。

  1. GENOME=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
  2. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  3. interval=~/wes_cancer/data/targets.preprocessed.interval.list
  5. cat config3 | while read id
  6. do
  7.     i=./5.gatk/${id}_bqsr.bam
  8.     echo ${i}
  9.     time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  CollectAllelicCounts \
  10.     -I ${i} \
  11.     -L ${interval} \
  12.     -R ${GENOME} \
  13.     -O ./8.cnv/gatk/allelicCounts/${id}.allelicCounts.tsv
  14. done

生成的 tsv 文件主要内容如下（这里我过滤掉了第六列为 N 的位点）：

  1. $ less ./8.cnv/gatk/allelicCounts/case1_biorep_A_techrep.allelicCounts.tsv | grep -v ^@ | awk '{if($6 != "N") print $0}' |less
  2. CONTIG  POSITION        REF_COUNT       ALT_COUNT       REF_NUCLEOTIDE  ALT_NUCLEOTIDE
  3. chr1    925873  34      1       G       T
  4. chr1    925918  71      1       G       T
  5. chr1    925938  92      1       G       T
  6. chr1    925953  92      1       G       T
  7. chr1    925974  96      1       G       T
  8. chr1    925979  102     2       G       A

6. ModelSegments

在第 3 步我们拿到了标准化和降噪后的两个 tsv 文件，记录了某个区间的 LOG2_COPY_RATIO 值，内容大致如下：

  1. $ less ./8.cnv/gatk/denoisedCR/case1_biorep_A_techrep.clean.denoisedCR.tsv | grep -v ^@| less
  2. CONTIG  START   END     LOG2_COPY_RATIO
  3. chr1    925692  926262  -1.178123
  4. chr1    929905  930585  -0.569447
  5. chr1    930789  931338  -0.686712
  6. chr1    935522  936145  -0.404623
  7. chr1    938790  939201  -0.727205
  8. chr1    939202  939709  -0.882516

而第 5 步拿到的记录等位基因测序深度的 tsv 文件已经在上面展示过了。接下来第 6 步将利用这两个结果进行 call segment，需要注意的是输入文件要求 tumor match normal。不过好像也可以不输入第 5 步 CollectAllelicCounts 的结果，等有时间再比较一下两者的区别吧。这一步用到的工具是ModelSegments，它根据去噪后的第三步的 reads counts 值对 copy ratios 进行分割，并根据第五步的CollectAllelicCounts等位基因计数对分割片段进行分类。代码如下：

  1. GATK=/~/wec_cancer/biosoft/gatk-4.1.4.1/gatk
  2. cat config3  | while read id
  3. do
  4.     germline=${id:0:5}_germline
  5.     ## ModelSegments
  6.     $GATK   --java-options "-Xmx20g" ModelSegments \
  7.     --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv \
  8.     --allelic-counts ./8.cnv/gatk/allelicCounts/${id}.allelicCounts.tsv \
  9.     --normal-allelic-counts ./8.cnv/gatk/allelicCounts/${germline}.allelicCounts.tsv \
  10.     --output ./8.cnv/gatk/segments \
  11.     --output-prefix ${id}
  13. done

生成的文件有点多，每个样本生成 11 个文件。 param 文件包含用于 copy ratios（cr）和 allele fractions（af）的全局参数，而 seg 文件包含有关片段的数据。 具体说明可以查看这个链接：https://software.broadinstitute.org/gatk/documentation/tooldocs/4.1.4.0/org_broadinstitute_hellbender_tools_copynumber_ModelSegments.php

  1. ./8.cnv/gatk/segments/case1_biorep_A_techrep.af.igv.seg
  2. ./8.cnv/gatk/segments/case1_biorep_A_techrep.cr.igv.seg
  3. ./8.cnv/gatk/segments/case1_biorep_A_techrep.cr.seg
  4. ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelBegin.af.param
  5. ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelBegin.cr.param
  6. ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelBegin.seg
  7. ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelFinal.af.param
  8. ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelFinal.cr.param
  9. ./8.cnv/gatk/segments/case1_biorep_A_techrep.modelFinal.seg
  10. ./8.cnv/gatk/segments/case1_biorep_A_techrep.hets.normal.tsv
  11. ./8.cnv/gatk/segments/case1_biorep_A_techrep.hets.tsv

其不过上面拿到的文件中的 ${id}*.igv.seg 文件，可以直接载入到 IGV 中进行可视化，如：

（不知道为什么，case4_biorep_B_techrep和case4_techrep_2的CNV事件是碎片化的，而 case5 病人的 X 染色体直接缺失）

11 拷贝数变异分析（GATK流程） - 图3

7. CallCopyRatioSegments

这一步用来判断 copy ratio segments 是扩增、缺失、还是正常的可能性。对上一步拿到的${id}.cr.seg进行推断，得到${id}.clean.called.seg文件会增加一列 CALL，用 +、-、0分别表示扩增、缺失和正常。基本上MEAN_LOG2_COPY_RATIO大于 0.14 就是扩增，小于 -0.15 就是缺失，其他的为正常。

  1. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  2. cat config3 | while read id
  3. do
  4.     $GATK   --java-options "-Xmx20g" CallCopyRatioSegments \
  5.     -I ./8.cnv/gatk/segments/${id}.cr.seg \
  6.     -O ./8.cnv/gatk/segments/${id}.clean.called.seg
  7. done

8. 可视化CNV结果

通过上面的分析，我们拿到了最后建模的 copy ratios 和 allele fractions segment，接下来用一个工具进行可视化：PlotModeledSegments

  1. dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
  2. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  4. cat config3 | while read id
  5. do
  6.     $GATK   --java-options "-Xmx20g" PlotModeledSegments \
  7.     --denoised-copy-ratios ./8.cnv/gatk/denoisedCR/${id}.clean.denoisedCR.tsv \
  8.     --allelic-counts ./8.cnv/gatk/segments/${id}.hets.tsv \
  9.     --segments ./8.cnv/gatk/segments/${id}.modelFinal.seg \
  10.     --sequence-dictionary ${dict} \
  11.     --output ./8.cnv/gatk/cnv_plots \
  12.     --output-prefix ${id}.clean
  13. done

11 拷贝数变异分析（GATK流程） - 图4

整个gatk cnv流程

上面整个流程的代码，其实可以合并为一个脚本 gatk_cnv.sh：

  1. GENOME=~/wes_cancer/data/Homo_sapiens_assembly38.fasta
  2. dict=~/wes_cancer/data/Homo_sapiens_assembly38.dict
  3. INDEX=~/wes_cancer/data/bwa_index/gatk_hg38
  4. GATK=~/wes_cancer/biosoft/gatk-4.1.4.1/gatk
  5. DBSNP=~/wes_cancer/data/dbsnp_146.hg38.vcf.gz
  6. kgSNP=~/wes_cancer/data/1000G_phase1.snps.high_confidence.hg38.vcf.gz
  7. kgINDEL=~/wes_cancer/data/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
  8. interval=~/wes_cancer/data/targets.preprocessed.interval.list
  10. cd ~/wes_cancer/project/8.cnv/gatk
  11. ###################################
  12. #### 把bam文件转为外显子reads数 ######
  13. ###################################
  16. cat ~/wes_cancer/project/config3 | while read id
  17. do
  18.     i=~/wes_cancer/project/5.gatk/${id}_bqsr.bam
  19.     echo ${i}
  20.     ## step1 : CollectReadCounts
  21.     time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"     CollectReadCounts \
  22.     -I ${i} \
  23.     -L ${interval} \
  24.     -R ${GENOME} \
  25.     --format HDF5  \
  26.     --interval-merging-rule OVERLAPPING_ONLY \
  27.     --output ${id}.clean_counts.hdf5
  28.     ## step2 : CollectAllelicCounts
  29.     time $GATK  --java-options "-Xmx20G -Djava.io.tmpdir=./"  CollectAllelicCounts \
  30.     -I ${i} \
  31.     -L ${interval} \
  32.     -R ${GENOME} \
  33.     -O ${id}.allelicCounts.tsv
  34. done
  35. ### 注意这个CollectAllelicCounts步骤非常耗时，而且占空间
  37. mkdir allelicCounts
  38. mv *.allelicCounts.tsv ./allelicCounts
  39. mkdir counts
  40. mv *.clean_counts.hdf5  ./counts
  41. ##################################################
  42. # 接着合并所有的normal样本的数据创建 cnvponM.pon.hdf5 #
  43. ##################################################
  45. $GATK  --java-options "-Xmx20g" CreateReadCountPanelOfNormals \
  46. --minimum-interval-median-percentile 5.0 \
  47. --output cnvponM.pon.hdf5 \
  48. --input counts/case1_germline.clean_counts.hdf5 \
  49. --input counts/case2_germline.clean_counts.hdf5 \
  50. --input counts/case3_germline.clean_counts.hdf5 \
  51. --input counts/case4_germline.clean_counts.hdf5 \
  52. --input counts/case5_germline.clean_counts.hdf5 \
  53. --input counts/case6_germline.clean_counts.hdf5
  56. ############################################
  57. ############# 最后走真正的CNV流程 #############
  58. ############################################
  60. cat config | while read id
  61. do
  62.     i=./counts/${id}.clean_counts.hdf5
  63.     $GATK  --java-options "-Xmx20g" DenoiseReadCounts \
  64.     -I $i \
  65.     --count-panel-of-normals cnvponM.pon.hdf5 \
  66.     --standardized-copy-ratios ${id}.clean.standardizedCR.tsv \
  67.     --denoised-copy-ratios ${id}.clean.denoisedCR.tsv
  68. done
  70. mkdir denoisedCR standardizedCR segments cnv_plots
  71. mv *denoisedCR.tsv ./denoisedCR
  72. mv *standardizedCR.tsv ./standardizedCR
  74. cat config | while read id
  75. do
  76.     i=./denoisedCR/${id}.clean.denoisedCR.tsv
  77.     ## ModelSegments的时候有两个策略，是否利用CollectAllelicCounts的结果
  78.     $GATK   --java-options "-Xmx20g" ModelSegments \
  79.     --denoised-copy-ratios $i \
  80.     --output segments \
  81.     --output-prefix ${id}
  82.     ## 如果要利用CollectAllelicCounts的结果就需要增加两个参数，这里就不讲解了。
  84.     $GATK   --java-options "-Xmx20g" CallCopyRatioSegments \
  85.     -I segments/${id}.cr.seg \
  86.     -O segments/${id}.clean.called.seg
  89.     ## 这里面有两个绘图函数，PlotDenoisedCopyRatios 和 PlotModeledSegments ，可以选择性运行。
  91.     $GATK   --java-options "-Xmx20g" PlotDenoisedCopyRatios \
  92.     --standardized-copy-ratios     ./standardizedCR/${id}.clean.standardizedCR.tsv \
  93.     --denoised-copy-ratios $i \
  94.     --sequence-dictionary ${dict} \
  95.     --output cnv_plots \
  96.     --output-prefix ${id}
  98.     $GATK   --java-options "-Xmx20g" PlotModeledSegments \
  99.     --denoised-copy-ratios $i \
  100.     --segments segments/${id}.modelFinal.seg \
  101.     --sequence-dictionary ${dict} \
  102.     --output cnv_plots \
  103.     --output-prefix ${id}.clean
  104. done

对于每一个样本，就会拿到拷贝数变异的结果，如：

11 拷贝数变异分析（GATK流程） - 图5

11 拷贝数变异分析（GATK流程） - 图6

11 拷贝数变异分析（GATK流程） - 图7