終於看到了一個完整的mutect2使用指令碼
終於看到了一個完整的mutect2使用指令碼
因為嫌麻煩,所以一直使用的是簡化版mutect2流程,其實就一個命令:
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
但是很明顯,這其實不符合官網教程的理念
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
但是官網教程的確太繁瑣,裡面涉及到6個步驟,而我只是執行了最後一個!
因為種種原因,沒能抽出時間細緻的探索mutect2的用法,但是無意中搜尋到指令碼一個: https://figshare.com/articles/scripts_sh/4542397
可以說是非常良心了:
#!/bin/sh
##GATK bundle download
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/ucsc.hg19.fasta.gz -O /path/to/GATK/bundle/ucsc.hg19.fasta.gz
wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg19/dbsnp_138.hg19.vcf.gz -O /path/to/GATK/bundle/dbsnp_138.hg19.vcf.gz
wget https://ndownloader.figshare.com/files/7354246 -O /path/to/GATK/bundle/TP53.sorted.bed
wget https://ndownloader.figshare.com/files/7354213 -O /path/to/GATK/bundle/CosmicCodingMuts.chr.sort.head.vcf
##ANNOVAR database files download
export PATH=$PATH:/path/to/annovar
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar refGene humandb/
annotate_variation.pl -buildver hg19 -downdb cytoBand humandb/
annotate_variation.pl -buildver hg19 -downdb genomicSuperDups humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar esp6500siv2_all humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar 1000g2015aug humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar snp138 humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar dbnsfp30a humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar cosmic70 humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar exac03 humandb/
annotate_variation.pl -buildver hg19 -downdb -webfrom annovar clinvar_20160302 humandb/
#define the reference file path
GATK=/path/to/GATK/GenomeAnalysisTK.jar
REF=/path/to/GATK/bundle/ucsc.hg19.fasta
DBSNP=/path/to/GATK/bundle/dbsnp_138.hg19.vcf
COSMIC=/path/to/GATK/bundle/CosmicCodingMuts.chr.sort.head.vcf
BED=/path/to/GATK/bundle/TP53.sorted.bed
#create a Panel of Normals (PoN) vcf file from the 4 low-grade tumour samples
for pathandfile in /path/to/EGA/normal/*.clean.recal.TP53.bam ; do
basewithpath=${pathandfile%.clean.recal.TP53.*}
basenopath=$(basename $basewithpath)
java -jar $GATK \
-T MuTect2 \
-R $REF \
-I:tumor $(echo $basewithpath).clean.recal.TP53.bam \
--dbsnp $DBSNP \
--cosmic $COSMIC \
--artifact_detection_mode \
-L $BED \
-o $(echo $basewithpath).clean.recal.TP53.normal.vcf
done
java -jar $GATK \
-T CombineVariants \
-R $REF \
-V /path/to/EGA/normal/GB544-10_S18.clean.recal.TP53.normal.vcf -V /path/to/EGA/normal/GB624-11_S81.clean.recal.TP53.normal.vcf -V /path/to/EGA/normal/GB730-12_S41.clean.recal.TP53.normal.vcf -V /path/to/EGA/normal/GB909-13_S90.clean.recal.TP53.normal.vcf \
-minN 2 \
--setKey "null" \
--filteredAreUncalled \
--filteredrecordsmergetype KEEP_IF_ANY_UNFILTERED \
-L $BED \
-o /path/to/EGA/normal/MuTect2_PON.vcf
#call somatic variants
for pathandfile in /path/to/EGA/tumor/*.clean.recal.TP53.bam ; do
basewithpath=${pathandfile%.clean.recal.TP53.*}
basenopath=$(basename $basewithpath)
java -jar $GATK \
-T MuTect2 \
-R $REF \
--dbsnp $DBSNP \
--cosmic $COSMIC \
-dt NONE \
--input_file:tumor $(echo $basewithpath).clean.recal.TP53.bam \
--intervals $BED \
-PON /path/to/EGA/normal/MuTect2_PON.vcf \
-o $(echo $basewithpath).clean.recal.TP53.vcf
done
#ANNOVAR annotation
mkdir /path/to/EGA/tumor/ANNOVAR
cp /path/to/EGA/tumor/*.vcf /path/to/EGA/tumor/ANNOVAR
##convert file format
for pathandfile in /path/to/EGA/tumor/ANNOVAR/*.vcf ; do
filename=${pathandfile%.*}
convert2annovar.pl --format vcf4 --includeinfo --withzyg $pathandfile > $(echo $filename).avinput
done
##add sample information
for pathandfile in /path/to/EGA/tumor/ANNOVAR/*.avinput ; do
filename=$(basename $pathandfile)
mainfilename=${filename%.clean.recal.TP53.avinput}
pathandmain=${pathandfile%.*}
awk -v var1="$mainfilename" '{OFS = "\t" ; print $0, var1}' $pathandfile > $(echo $pathandmain).avinputs
done
##merge file
cat /path/to/EGA/tumor/ANNOVAR/*.avinputs > /path/to/EGA/tumor/ANNOVAR/TP53.Tonly.avinput
##annotate
table_annovar.pl /path/to/EGA/tumor/ANNOVAR/TP53.Tonly.avinput /path/to/annovar/humandb/ -buildver hg19 -out /path/to/EGA/tumor/ANNOVAR/TP53.Tonly -remove -protocol refGene,cytoBand,genomicSuperDups,esp6500siv2_all,1000g2015aug_all,snp138,dbnsfp30a,cosmic70,exac03,clinvar_20160302 -operation g,r,r,f,f,f,f,f,f,f --otherinfo
雖然這個流程是基於 hg19 參考基因組的,但是很容易就能改寫為 hg38版本的!
還有一個問題是他那個流程的專案背景是,得到的腫瘤樣品是沒有配對normal的,而是 create a Panel of Normals (PoN) vcf file from the 4 low-grade tumour samples
而且,這個流程是基於 GATK3成熟版本的,並不是GATK4哦。
如果大家有GATK4的MUTECT2經驗,可以交流一下哈。
相關文章
- [Python] 用python做一個遊戲輔助指令碼,完整思路Python遊戲指令碼
- 一個用於生成大量mac地址的python指令碼MacPython指令碼
- 使用Java實現一個JS指令碼引擎JavaJS指令碼
- 分享一個基於jQuery的鎖定表格行列的js指令碼。jQueryJS指令碼
- 建立批次AD域使用者的指令碼可以使用 PowerShell 來實現。以下是一個簡單的示例指令碼,用於批次建立使用者:指令碼
- 一個分詞指令碼分詞指令碼
- 開源一個Flutter編寫的完整終端模擬器Flutter
- 如何使用Tampermonkey開發並使用一個瀏覽器指令碼瀏覽器指令碼
- 京東2020雙十二活動終於來啦,指令碼助你領年終指令碼
- 一個方便 LeetCode 複習的指令碼LeetCode指令碼
- 新增多個使用者的shell指令碼指令碼
- 30個關於Shell指令碼的經典案例(中)指令碼
- 30個關於Shell指令碼的經典案例(上)指令碼
- 30個關於Shell指令碼的經典案例(下)指令碼
- 向大家分享一個shell指令碼的坑指令碼
- 一個快速檢視trace的小指令碼指令碼
- 共享一個iptables的shell指令碼檔案指令碼
- 如何在Linux中使用Shell指令碼終止使用者會話?Linux指令碼會話
- 一個centos初始化指令碼CentOS指令碼
- 將一個Python指令碼做成一個Windows服務Python指令碼Windows
- 如何在Windows上使用Git建立一個可執行指令碼?WindowsGit指令碼
- PHP 使用檔案鎖 避免同時執行一個指令碼PHP指令碼
- Mac 終端執行 shell 指令碼Mac指令碼
- 一文詳解基於指令碼的攻擊指令碼
- 分享工作中常用的一個Git指令碼Git指令碼
- 分享一個提高運維效率的 Python 指令碼運維Python指令碼
- 一個能夠生成 Markdown 表格的 Bash 指令碼指令碼
- 油猴外掛及油猴指令碼使用手冊(完整版)指令碼
- 使用代理池用py完整的爬取一個網站(尾部有github原始碼)網站Github原始碼
- 一個網站故障排查的、程式碼更新的簡便指令碼網站指令碼
- 擼一個 iOS 重簽名指令碼iOS指令碼
- 基於python編寫一個簡單的多執行緒埠掃描指令碼Python執行緒指令碼
- 獲取sql完整指令碼,get_fulltext.shSQL指令碼
- 終於,幫開發寫了一個bug
- 一個“指令碼執行夯死”問題的分析指令碼
- Shell:如何寫一個多選選單的指令碼指令碼
- 用於自動監控磁碟使用情況的 Shell 指令碼指令碼
- 幾個Linux命令及指令碼使用中的奇淫巧技Linux指令碼