基因组 Tumor-only mode – Mutect2/CNVkit

Posted on by KevinZhou

Mutect2

官方脚本模板：

  gatk Mutect2 \
  -R reference.fa \
  -I sample.bam \
  --germline-resource af-only-gnomad.vcf.gz \
  --panel-of-normals pon.vcf.gz \
  -O single_sample.vcf.gz

构建批量代码

ls ../align/*_bqsr.bam | perl -ne 'chomp; if($_ =~ /([^\/]+)_bqsr\.bam/) { print "gatk --java-options \"-Xmx20G -Djava.io.tmpdir=./\" Mutect2 -R /data02/zhangmengmeng/database/hg38/hg38.fa -I ../align/$1_bqsr.bam -L /data02/zhangmengmeng/database/gatk_resource_bundle/hg38/AgilentV6_GRCh38_ex_region.sort.bed --germline-resource /data02/zhangmengmeng/database/gatk_resource_bundle/hg38/somatic-hg38_af-only-gnomad.hg38.vcf.gz --panel-of-normals /data02/zhangmengmeng/database/gatk_resource_bundle/hg38/somatic-hg38_1000g_pon.hg38.vcf.gz --f1r2-tar-gz $1.f1r2.tar.gz --O $1.raw.vcf && gatk LearnReadOrientationModel -I $1.f1r2.tar.gz -O $1.read-orientation-model.tar.gz && gatk GetPileupSummaries -I ../align/$1_bqsr.bam -L /data02/zhangmengmeng/database/gatk_resource_bundle/hg38/AgilentV6_GRCh38_ex_region.sort.bed -V /data02/zhangmengmeng/database/gatk_resource_bundle/hg38/somatic-hg38_af-only-gnomad.hg38.vcf.gz -O $1.pileups.table && gatk CalculateContamination -I $1.pileups.table -O $1.CalculateContamination.table && gatk FilterMutectCalls -R /data02/zhangmengmeng/database/hg38/hg38.fa -V $1.raw.vcf --contamination-table $1.CalculateContamination.table --ob-priors $1.read-orientation-model.tar.gz -O $1.somatic.vcf && gatk FilterAlignmentArtifacts -R /data02/zhangmengmeng/database/hg38/hg38.fa --bwa-mem-index-image /data02/zhangmengmeng/database/hg38/hg38.fa.img -V $1.somatic.vcf -O $1.somatic.FAA.vcf -I ../align/$1_bqsr.bam && awk \x27/^#/ {print \$0; next} \$7==\x22PASS\x22  {print \$0}\x27 $1.somatic.FAA.vcf > ./$1.somatic.passed.vcf  && echo $1 ok\n" }' > mutect2_TumorOnly.sh

annovar 注释

ls *.somatic.passed.vcf | perl -ne 'chomp; my $name = $1 if ($_ =~ /([^\/]+)\.somatic\.passed\.vcf/); print "table_annovar.pl $name.somatic.passed.vcf /home/zhoukaiwen/database/Annovar/humandb_hg38 -buildver hg38 -out $name.somatic.passed -remove -protocol refGene,cytoBand,avsnp151,exac03,gnomad41_exome,EAS.sites.2015_08,cosmic102_unique -operation g,r,f,f,f,f,f -nastring . -vcfinput \n"'>annovar.sh

CNVkit

构建正常对照

# 初步划分bin，给每个bin增加基因名称注释，有利于结果展示
cnvkit.py target ../GRCh38/genecode_GRCh38.p14.genome.interval.bed --annotate refFlat_hg38.txt --split -o genecode_GRCh38.p14.genome.interval_cnvkit_target.bed

# 根据样本测序深度自动划分bin(可选)
cnvkit.py autobin Tumor.bam -t /home/zhoukaiwen/database/cnvkit/genecode_GRCh38.p14.genome.interval_cnvkit_target.bed # 会生成my_targets.target.bed和my_targets.antitarget.bed
cnvkit.py reference -o Tumor_FlatReference.cnn -f hg38.fa -t Tumor_targets.target.bed

# 直接生成FlatReference.cnn
cnvkit.py reference -o FlatReference.cnn -f ../GRCh38/genecode_GRCh38.p14.genome.fa -t genecode_GRCh38.p14.genome_cnvkit.target.bed

运行CNVkit Tumor-only mode

cnvkit.py batch Tumor.bam -r FlatReference.cnn