Mutect2使用

2022-05-16  本文已影响0人  日月其除

首先参考一下gatk官网:
https://gatk.broadinstitute.org/hc/en-us/articles/360047232772--Notebook-Intro-to-using-Mutect2-for-somatic-data
https://gatk.broadinstitute.org/hc/en-us/articles/360037593851-Mutect2
以及建明老师的解说:http://www.360doc.com/content/21/0714/12/76149697_986501266.shtml
第一步:make your own PoN
分为三步:

  1. Run Mutect2 in tumor-only mode on each normal BAM individually
    只把normal的处理后的bam文件生成normal.vcf
gatk Mutect2 -R reference.fasta -I normal1.bam --max-mnp-distance 0 -O normal1.vcf.gz 
gatk Mutect2 -R reference.fasta -I normal2.bam --max-mnp-distance 0 -O normal2.vcf.gz 
...
gatk Mutect2 -R reference.fasta -I normal40.bam --max-mnp-distance 0 -O normal40.vcf.gz 
  1. Create a GenomicsDB from the normal Mutect2 calls
gatk GenomicsDBImport -R reference.fasta -L intervals.interval_list \
  --genomicsdb-workspace-path pon_db \
  -V normal1.vcf.gz \
  -V normal2.vcf.gz \
  ...
  -V normal40.vcf.gz

其中-L后面这文件,我看到有的人说WGS可以不用,WES需要测序时提供的位置文件转化。-L指定你需要操作的目标区域。
我是用的全基因组测序,我使用全基因组参考基因组fa文件生成这个文件,参考 https://www.jianshu.com/p/3982ea467976
根据这个网址可以得到每个染色体的位置信息
生成这个list文件后,我又继续使用 gatk PreprocessIntervals 参考:https://gatk.broadinstitute.org/hc/en-us/articles/360037226712-PreprocessIntervals
把文件切分为1000bp为bin的文件

To generate consecutive bins of 1000 bases from the reference (e.g., for whole genome sequencing analyses):
gatk PreprocessIntervals \
          -R reference.fa \
          --bin-length 1000 \
          --padding 0 \
          -O preprocessed_intervals.interval_list

这个网址里面也有Interval list比较好的解说:https://zhuanlan.zhihu.com/p/66034477

  1. and then Combine the normal calls using CreateSomaticPanelOfNormals.
gatk CreateSomaticPanelOfNormals -R reference.fasta \
  --germline-resource af-only-gnomad.vcf.gz \
  -V gendb://pon_db \
  -O pon.vcf.gz

这里的这个af-only-gnomad.vcf.gz文件来源于gatk bundle
那怎么获得呢,参考:Resource bundle – GATK (broadinstitute.org)
可以通过谷歌云网盘以及ftp下载,我没有云账号,因此我使用ftp下载

wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/Mutect2/af-only-gnomad.hg38.vcf.gz.tbi

可以参考如何使用lftp通过服务器直接链接gatk的下载界面
2 下载GATK需要的参考基因组文件 - 简书 (jianshu.com)

最后就是call varations

gatk Mutect2 \
     -R reference.fa \
     -I tumor.bam \
     -I normal.bam \
     -normal normal_sample_name \
     --germline-resource af-only-gnomad.vcf.gz \
     --panel-of-normals pon.vcf.gz \
     -O somatic.vcf.gz

怎么得到这里的normal_sample_name,可以使用gatk的getsamplename,参考网站:https://gatk.broadinstitute.org/hc/en-us/articles/360037224912-GetSampleName

 gatk GetSampleName \
     -I input.bam \
     -O sample_name.txt

后续过滤还有很多步操作,后面补充!这个确实不太好用!

上一篇 下一篇

猜你喜欢

热点阅读