也就是说我搜索到了一个4小时前的教程，取代了之前的一个月前的教程，这，生活太苦了。
其它参考资料包括：
这样的资料是谷歌云数据文件（大约几百个G的磁盘空间消耗）：
上来就摆这么多参考资料，可想而知最新最全的mutect2教程是多么的难写！ 背景知识
NGS使体细胞变异的检测更加全面，成本更低，在检测多种体细胞变异上具有很大的优势，但在使用过程中还存在着挑战：如**样品降解、覆盖度不足、遗传异质性和组织污染（杂质)等问题。为应对以上挑战，降低错误率，科学家采取了不同的算法和统计模型用于检测体细胞突变。目前最受欢迎的有Varscan、SomaticSniper、 Strelka 和MuTect2 **。这些软件大都是直接对肿瘤-正常样本的每个位点进行比较，对肿瘤样本中明显高于正常样本的次等位基因进行标记，作为体细胞变异，同时排除种系突变和杂合性丢失（LOH）情况。虽然这些软件具有较高的引用率，并在不断地更新，但仍存在不足：
下面是TCGA计划采取的软件：
大家可以去下载到TCGA计划的这4个软件输出的maf文件格式的somatic突变信息文件哦。 第一步：构建PoNFirst, 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>&1fii=$((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.gzStep 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>&1Step 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>&1In 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列信息：
代码如下： 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.vcfdone ## 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命令：
同样的解析配置文件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.vcfdone ## 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选择。简单搜索几个关于的综述或者测评工具：
我们可以选择的 工具超级多，不要在一棵树上吊死！ 我们后面就讲解提到方案，敬请期待哈！抵制辣鸡软件，人人有责！ |
|
|
