最新最全的mutect2教程

GATK的Mutect2流程一直在变动,主要是GATK本身也更新频率有点高,所以基本上大家看到的教程很快就过时了,follow起来都是错误连连。现在这个教程的时间是:2020-09-22 (只能保证说未来半年内可能是ok的)

也就是说我搜索到了一个4小时前的教程,取代了之前的一个月前的教程,这,生活太苦了

  • https://gatk.broadinstitute.org/hc/en-us/articles/360035889791?id=11136
  • https://gatk.broadinstitute.org/hc/en-us/articles/360035531132

其它参考资料包括:

  • https://gatk.broadinstitute.org/hc/en-us/articles/360035531132--How-to-Call-somatic-mutations-using-GATK4-Mutect2
  • https://github.com/broadinstitute/gatk/blob/master/docs/mutect/mutect.pdf
  • http://rstudio-pubs-static.s3.amazonaws.com/491640_6c12de260e734538b01777c0cc4ff92b.html
  • https://gatk.broadinstitute.org/hc/en-us/articles/360035531852-Intervals-and-interval-lists

这样的资料是谷歌云数据文件(大约几百个G的磁盘空间消耗):

  • https://console.cloud.google.com/storage/browser/gatk-best-practices
  • https://console.cloud.google.com/storage/browser/genomics-public-data

上来就摆这么多参考资料,可想而知最新最全的mutect2教程是多么的难写!

背景知识

体细胞突变(somatic mutation)是指患者某些组织或者器官后天性地发生了体细胞变异,虽然它不会遗传给后代个体,却可以通过细胞分裂,遗传给子代细胞。体细胞突变对肿瘤的发生发展有关键性的作用,并且它也是制定肿瘤癌症靶向治疗措施的关键所在。

NGS使体细胞变异的检测更加全面,成本更低,在检测多种体细胞变异上具有很大的优势,但在使用过程中还存在着挑战:如**样品降解、覆盖度不足、遗传异质性和组织污染(杂质)等问题。为应对以上挑战,降低错误率,科学家采取了不同的算法和统计模型用于检测体细胞突变。目前最受欢迎的有Varscan、SomaticSniper、 Strelka 和MuTect2 **。这些软件大都是直接对肿瘤-正常样本的每个位点进行比较,对肿瘤样本中明显高于正常样本的次等位基因进行标记,作为体细胞变异,同时排除种系突变和杂合性丢失(LOH)情况。虽然这些软件具有较高的引用率,并在不断地更新,但仍存在不足:

  • a 、缺乏完整可靠的实验来评估检测结果;

  • b、 缺乏金标准,不能保证检测到的灵敏度和特异性最高;

  • c、 在实际应用中,各软件的相对优缺点在很大程度上是未知的。

下面是TCGA计划采取的软件:

  • MuSE
  • varscan
  • MuTect
  • SomaticSniper

大家可以去下载到TCGA计划的这4个软件输出的maf文件格式的somatic突变信息文件哦。

第一步:构建PoN

First, run gatk Mutect2 in tumor-only mode on each normal sample to call all detectable variants: Needs an interval list. 可以使用gatk 的 BedToIntervalList和PreprocessIntervals命令来根据bed制作interval list.,注意一下wgs和wes的差异即可。

最新版的GATK里面的CreateSomaticPanelOfNormals命令更新了,需要下面3个步骤完成 panel of normal (PoN) 流程

Step 1: Run Mutect2 in tumor-only mode for each normal sample:

把所有的normal的bam文件的路径整理好,然后批量运行,脚本如下:

## $1 for the config file, such as : normal_bam_for_pon.txt ## $2 and $3 for submit jobs.# for i in {0..5};do sbatch run_PoN-step1.sh normal_bam_for_pon.txt 6 $i;done# conda activate qcGENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta GATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk cat $1 |while read normal_bamdoecho $normal_bamif((i%$2==$3))then id=$(basename "$normal_bam" _recal.bam) $GATK Mutect2 \ -R ${GENOME} \ -I $normal_bam \ --max-mnp-distance 0 \ -O ${id}_mutect2_normal.vcf.gz \ 1> ${id}_mutect2_normal.log 2>&1
fii=$((i+1))done

一般来说,批量运行结束后输出如下的文件,如果你的测序数据很稳定,那么基本上每个normal的wes数据走完流程的文件大小是类似的:

71M Sep 15 16:30 N1092_mutect2_normal.vcf.gz61M Sep 15 15:40 N1146_mutect2_normal.vcf.gz50M Sep 15 15:41 N1254_mutect2_normal.vcf.gz57M Sep 15 15:30 N1257_mutect2_normal.vcf.gz54M Sep 15 15:18 N1268_mutect2_normal.vcf.gz52M Sep 15 22:10 N1271_mutect2_normal.vcf.gz55M Sep 15 20:09 N1281_mutect2_normal.vcf.gz

Step 2: Create a GenomicsDB from the normal Mutect2 calls:

这个时候需要合并上述每个单独normal的wes数据得到的vcf文件啦,需要指定wes的范围文件(bed格式),使用GenomicsDBImport命令进行合并操作,代码如下:

GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta GATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk interval=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/resources_broad_hg38_v0_wgs_calling_regions.hg38.interval_list # normal_bam_for_pon.txt${GATK} GenomicsDBImport \-R ${GENOME} \-L ${interval} \--merge-input-intervals TRUE \$(ls *.vcf.gz | awk '{print "-V "$0" "}') \--genomicsdb-workspace-path pon_db \1> GenomicsDBImport.log 2>&1

Step 3: Combine the normal calls using CreateSomaticPanelOfNormals:

GenomicsDBImport命令的合并输出的是临时文件夹pon_db,接下来仍然是需要使用衔接命令CreateSomaticPanelOfNormals来输出真正的pon.vcf.gz 文件。

这个时候还加入了somatic-hg38_af-only-gnomad文件,来源于GATK的官方谷歌云下载中心。

gnomad=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/somatic-hg38_af-only-gnomad.hg38.vcf ${GATK} CreateSomaticPanelOfNormals \-R ${GENOME} \--germline-resource ${gnomad} \-V gendb://pon_db \-O pon.vcf.gz \1> CreateSomaticPanelOfNormals.log 2>&1

In the third step, the AF-only gnomAD resource is the same public GATK resource used by Mutect2 (when calling tumor samples, but not in Step 1, above).

走这3个步骤,就是为了生成 pon.vcf.gz 文件,供后续使用!

第二步:call mutations in the tumor in somatic mode.

主要是 Mutect2 命令,以及 FilterMutectCalls 命令。

首先看Mutect2 命令的代码,前面步骤生成 pon.vcf.gz 文件这个时候就利用上来了,需要制作一个config文件,配合下面的脚本,主要是3列信息:

  • 第一列是肿瘤命名
  • 第二列是肿瘤病人的normal组织的bam文件地址
  • 第三列是肿瘤病人的肿瘤组织的bam文件地址。

代码如下:

GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta reference=$GENOMEGATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk cat $config_file |while read iddoarr=($id)normal_bam=${arr[1]}tumor_bam=${arr[2]}sample=${arr[0]} # 上面是解析配置文件config time $GATK Mutect2 -R $reference \ -I $tumor_bam -tumor $(basename "$tumor_bam" _recal.bam) \ -I $normal_bam -normal $(basename "$normal_bam" _recal.bam) \ -pon pon.vcf.gz \ -O ${sample}_mutect2.vcf
done ## end for $config_file

这样的话,每个样本大约可以拿到四五千左右的可能的somatic位点,需要进行一定程度的过滤操作。(去年我的教程过滤很简单,就是FilterMutectCalls 命令而已,现在它已经成为了辣鸡,不用它了。)

这个时候的过滤需要两个文件的配合,首先是read-orientation-model.tar.gz,其次是 contamination.table和segments.table。

真恶心,辣鸡软件!不用了···

后面我看了看教程,整理了就放弃了,浪费时间,浪费生命!真恶心,辣鸡软件!不用了···

制作contamination.table和segments.table

真恶心,辣鸡软件!不用了···

参考:https://github.com/broadinstitute/gatk/blob/master/docs/mutect/mutect.pdf

需要 GetPileupSummaries以及CalculateContamination命令:

  • https://gatk.broadinstitute.org/hc/en-us/articles/360037593451-GetPileupSummaries
  • https://gatk.broadinstitute.org/hc/en-us/articles/360036888972-CalculateContamination

同样的解析配置文件config,走下面的代码:

GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta reference=$GENOMEGATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk cat $config_file |while read iddoarr=($id)normal_bam=${arr[1]}tumor_bam=${arr[2]}sample=${arr[0]} # 上面是解析配置文件config time $GATK Mutect2 -R $reference \ -I $tumor_bam -tumor $(basename "$tumor_bam" _recal.bam) \ -I $normal_bam -normal $(basename "$normal_bam" _recal.bam) \ -pon pon.vcf.gz \ -O ${sample}_mutect2.vcf
done ## end for $config_file

制作read-orientation-model.tar.gz文件

真恶心,辣鸡软件!不用了···

执行 FilterMutectCalls 命令

真恶心,辣鸡软件!不用了···

GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta reference=$GENOMEGATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk GATK=/home/jianmingzeng/biosoft/GATK/gatk-4.0.3.0/gatkls *_mutect2.vcf |while read iddosample=$(basename "$id" _mutect2.vcf)$GATK FilterMutectCalls -R $reference \ -V $id \ -O ${sample}_filtered.vcf done

会报错,gatk官方论坛的意思是,在集群运行的过程中,会丢失后缀为.vcf.stats的文件,所以FilterMutectCalls 命令失败。如下:

A USER ERROR has occurred: Mutect stats table _mutect2.vcf.stats not found. When Mutect2 outputs a file calls.vcf it also creates a calls.vcf.stats file. Perhaps this file was not moved along with the vcf, or perhaps it was not delocalized from a virtual machine while running in the cloud.

真恶心,辣鸡软件!不用了···

我们我敢抛弃mutect2呢

虽然说,GATK是人类数据的首选工具,但是GATK的Mutect2流程并不是唯一的找somatic mutation选择。简单搜索几个关于的综述或者测评工具:

  • 2016年文章:Evaluation of Nine Somatic Variant Callers for Detection of Somatic Mutations in Exome and Targeted Deep Sequencing Data,链接是:https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0151664
  • 2018年文章:A review of somatic single nucleotide variant calling algorithms fornext-generation sequencing data,链接是https://www.sciencedirect.com/science/article/pii/S2001037017300946

我们可以选择的 工具超级多,不要在一棵树上吊死!

我们后面就讲解提到方案,敬请期待哈!抵制辣鸡软件,人人有责!

(0)

相关推荐