GATK4的mutect2流程

生信技能树

发布于 2018-07-27 09:52:40

2.9K0

发布于 2018-07-27 09:52:40

本来以为肿瘤外显子教程分享完了，经粉丝提醒才发现原来是我在自己的生信菜鸟团博客连载完毕，却没有上传到微信公众号，给大家说一声抱歉，漏掉几个知识点。首先看看GATK4的mutect2和GATK3的相比有哪些改动，图片来源：https://gatkforums.broadinstitute.org/gatk/discussion/10911/differences-between-gatk3-mutect2-and-gatk4-mutect2

往期GATK4教程目录：

GATK4的gvcf流程

你以为的可能不是你以为的

新鲜出炉的GATK4培训教材全套PPT，赶快下载学习吧

曾老湿最新私已：GATK4实战教程

GATK4的CNV流程-hg38

官网教程

https://gatkforums.broadinstitute.org/gatk/discussion/9183/how-to-call-somatic-snvs-and-indels-using-mutect2 https://software.broadinstitute.org/gatk/documentation/article?id=9183

非常复杂，而且步骤繁多，如果只是想测试一下：

首先制作配置文件，如下；

oscc_01 /oscc/WES/alignment/OSCC_01_N_recal.bam /oscc/WES/alignment/OSCC_01_T_recal.bam
oscc_04 /oscc/WES/alignment/OSCC_04_N_recal.bam /oscc/WES/alignment/OSCC_04_T_recal.bam
oscc_06 /oscc/WES/alignment/OSCC_06_N_recal.bam /oscc/WES/alignment/OSCC_06_T_recal.bam
oscc_09 /oscc/WES/alignment/OSCC_09_N_recal.bam /oscc/WES/alignment/OSCC_09_T_recal.bam
oscc_10 /oscc/WES/alignment/OSCC_10_N_recal.bam /oscc/WES/alignment/OSCC_10_T_recal.bam
oscc_11 /oscc/WES/alignment/OSCC_11_N_recal.bam /oscc/WES/alignment/OSCC_11_T_recal.bam
oscc_13 /oscc/WES/alignment/OSCC_13_N_recal.bam /oscc/WES/alignment/OSCC_13_T_recal.bam
oscc_14 /oscc/WES/alignment/OSCC_14_N_recal.bam /oscc/WES/alignment/OSCC_14_T_recal.bam
oscc_15 /oscc/WES/alignment/OSCC_15_N_recal.bam /oscc/WES/alignment/OSCC_15_T_recal.bam
oscc_16 /oscc/WES/alignment/OSCC_16_N_recal.bam /oscc/WES/alignment/OSCC_16_T_recal.bam

需要根据以往的教程安装好GATK并且下载好配套文件。

然后运行下面的代码：

module load java/1.8.0_91
GENOME=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
GATK=/home/jianmingzeng/biosoft/GATK/gatk-4.0.3.0/gatk
DBSNP=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
reference=/home/jianmingzeng/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta


cat $config_file |while read id
do
    arr=($id)
    normal_bam=${arr[1]}
    tumor_bam=${arr[2]}
    sample=${arr[0]} 

start=$(date +%s.%N)
echo Mutect2 `date`
time $GATK  --java-options "-Xmx10G -Djava.io.tmpdir=./"  Mutect2 -R $reference \
-I $tumor_bam  -tumor $(basename "$tumor_bam" _recal.bam) \
-I $normal_bam -normal $(basename "$normal_bam" _recal.bam) \
-O ${sample}_mutect2.vcf
$GATK  FilterMutectCalls -V ${sample}_mutect2.vcf -O ${sample}_somatic.vcf
echo Mutect2 `date`
dur=$(echo "$(date +%s.%N) - $start" | bc)
printf "Execution time for Mutect2 : %.6f seconds" $dur
echo 

done

简单过滤

for i in *_somatic.vcf
do
j=$(basename "$i" _somatic.vcf )
echo $j
cat $i | perl -alne '{if(/^#/){print}else{next unless $F[6] eq ".";next if $F[0] =~/_/;print } }' > ${j}_filter.vcf
done

把vcf文件转为maf文件，需要参考我在生信菜鸟团前面的博客

用VEP对vcf格式的突变数据进行注释
把vcf文件转换为maf格式

cat config |while read id
do
    arr=($id)
    normal_bam=${arr[1]}
    tumor_bam=${arr[2]}
    sample=${arr[0]}

    perl ~/biosoft/vcf2maf/vcf2maf.pl --input-vcf ${sample}_filter.vcf   --output-maf ${sample}.maf  \
    --ref-fasta ~/.vep/homo_sapiens/86_GRCh38/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz \
    --tumor-id  $(basename "$tumor_bam" _recal.bam)  --normal-id  $(basename "$normal_bam" _recal.bam)  --ncbi-build GRCh38

done

得到的maf就可以用maftools去可视化啦！

提醒

GATK4目前主流分析选择的人不多，大部分公司或者科研院所仍然是以成熟版本的GATK4系列为流程！

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2018-06-23，如有侵权请联系 cloudcommunity@tencent.com 删除

大数据

数据分析

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

大数据

数据分析

登录后参与评论

0 条评论

热度