GATK的Mutect2流程一直在变动,主要是GATK本身也更新频率有点高,所以基本上大家看到的教程很快就过时了,follow起来都是错误连连。现在这个教程的时间是:2020-09-22 (只能保证说未来半年内可能是ok的)
也就是说我搜索到了一个4小时前的教程,取代了之前的一个月前的教程,这,生活太苦了。
其它参考资料包括:
这样的资料是谷歌云数据文件(大约几百个G的磁盘空间消耗):
上来就摆这么多参考资料,可想而知最新最全的mutect2教程是多么的难写!
体细胞突变(somatic mutation)是指患者某些组织或者器官后天性地发生了体细胞变异,虽然它不会遗传给后代个体,却可以通过细胞分裂,遗传给子代细胞。体细胞突变对肿瘤的发生发展有关键性的作用,并且它也是制定肿瘤癌症靶向治疗措施的关键所在。
NGS使体细胞变异的检测更加全面,成本更低,在检测多种体细胞变异上具有很大的优势,但在使用过程中还存在着挑战:如**样品降解、覆盖度不足、遗传异质性和组织污染(杂质)等问题。为应对以上挑战,降低错误率,科学家采取了不同的算法和统计模型用于检测体细胞突变。目前最受欢迎的有Varscan、SomaticSniper、 Strelka 和MuTect2 **。这些软件大都是直接对肿瘤-正常样本的每个位点进行比较,对肿瘤样本中明显高于正常样本的次等位基因进行标记,作为体细胞变异,同时排除种系突变和杂合性丢失(LOH)情况。虽然这些软件具有较高的引用率,并在不断地更新,但仍存在不足:
下面是TCGA计划采取的软件:
大家可以去下载到TCGA计划的这4个软件输出的maf文件格式的somatic突变信息文件哦。
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) 流程
把所有的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 qc
GENOME=/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_bam
do
echo $normal_bam
if((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
fi
i=$((i+1))
done
一般来说,批量运行结束后输出如下的文件,如果你的测序数据很稳定,那么基本上每个normal的wes数据走完流程的文件大小是类似的:
71M Sep 15 16:30 N1092_mutect2_normal.vcf.gz
61M Sep 15 15:40 N1146_mutect2_normal.vcf.gz
50M Sep 15 15:41 N1254_mutect2_normal.vcf.gz
57M Sep 15 15:30 N1257_mutect2_normal.vcf.gz
54M Sep 15 15:18 N1268_mutect2_normal.vcf.gz
52M Sep 15 22:10 N1271_mutect2_normal.vcf.gz
55M Sep 15 20:09 N1281_mutect2_normal.vcf.gz
这个时候需要合并上述每个单独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
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 文件,供后续使用!
主要是 Mutect2 命令,以及 FilterMutectCalls 命令。
首先看Mutect2 命令的代码,前面步骤生成 pon.vcf.gz 文件这个时候就利用上来了,需要制作一个config文件,配合下面的脚本,主要是3列信息:
代码如下:
GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
reference=$GENOME
GATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk
cat $config_file |while read id
do
arr=($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。
真恶心,辣鸡软件!不用了···
后面我看了看教程,整理了就放弃了,浪费时间,浪费生命!真恶心,辣鸡软件!不用了···
真恶心,辣鸡软件!不用了···
参考: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=$GENOME
GATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk
cat $config_file |while read id
do
arr=($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
真恶心,辣鸡软件!不用了···
真恶心,辣鸡软件!不用了···
GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
reference=$GENOME
GATK=/home/jianmingzeng/biosoft/gatk/gatk-4.1.8.1/gatk
GATK=/home/jianmingzeng/biosoft/GATK/gatk-4.0.3.0/gatk
ls *_mutect2.vcf |while read id
do
sample=$(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.
真恶心,辣鸡软件!不用了···
虽然说,GATK是人类数据的首选工具,但是GATK的Mutect2流程并不是唯一的找somatic mutation选择。简单搜索几个关于的综述或者测评工具:
我们可以选择的 工具超级多,不要在一棵树上吊死!
我们后面就讲解提到方案,敬请期待哈!抵制辣鸡软件,人人有责!