很容易谷歌搜索就可以找到最新版教程:https://software.broadinstitute.org/gatk/documentation/quickstart 我就不做文档搬运工了,感兴趣自己看原文,目录如下:
- Quick start for the impatient
- Requirements
- Get GATK
- Install it
- Test that it works
- Run GATK and Picard commands
- Grok the Best Practices
- Run pipelines
- Get help
我的规矩就是分享石锤,生信菜鸟团的一贯作风!
首先下载软件
最新版的GATK已经打包了PICARD软件,非常方便,它是一系列基因检测数据分析工具的合辑,想用好很不容易,如果只是想用起来,到还可以,下载命令如下,复制粘贴即可使用:
mkdir -p ~/biosoft/GATK/
cd ~/biosoft/GATK/
wget https://github.com/broadinstitute/gatk/releases/download/4.0.2.1/gatk-4.0.2.1.zip
unzip gatk-4.0.2.1.zip
~/biosoft/GATK/gatk-4.0.2.1/gatk --help
可以看到,所有的工具都被打包到一个命令里面了,进入文档可以查看所有命令详情:https://software.broadinstitute.org/gatk/documentation/tooldocs (有几百个命令哦)
同时,命令的参数也是有规律的:https://software.broadinstitute.org/gatk/documentation/article?id=11050
还是老规矩,生信菜鸟团不搬运英文文档,大家自己慢慢看。
值得注意的是 GATK Best Practices 也进行了更新:https://software.broadinstitute.org/gatk/best-practices/ 下面我就会对自己的数据走一波这些流程。
- Data Pre-processing
- Germline SNPs + Indels
- Somatic SNVs + Indels
- RNAseq SNPs + Indels
- Germline CNVs
- Somatic CNVs
虽然有最佳实践指导手册,但是不同的数据毕竟性质特征不一样,比如 Whole genomes. Exomes. Gene panels. RNAseq , 这些还是得靠自己把握。
下载必备文件
一般来说,可以选择最新版人类参考基因组了,如果你实在有其它需求,也可以自行选择其它版本。
mkdir -p ~/biosoft/GATK/resources/bundle/hg38
cd ~/biosoft/GATK/resources/bundle/hg38
## ftp://ftp.broadinstitute.org/bundle/hg38/
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/dbsnp_146.hg38.vcf.gz.tbi &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.fasta.fai &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/Homo_sapiens_assembly38.dict &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz &
nohup wget ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi &
mkdir bwa_index
cd bwa_index
~/biosoft/bwa/bwa-0.7.15/bwa index -a bwtsw \
-p gatk_hg38 ../Homo_sapiens_assembly38.fasta \
1>hg38.bwa_index.log 2>&1 &
下载全部成功后,应该如下:
|-- [1.8G] 1000G_phase1.snps.high_confidence.hg38.vcf.gz
|-- [2.0M] 1000G_phase1.snps.high_confidence.hg38.vcf.gz.tbi
|-- [4.0K] bwa_index
| |-- [ 20K] gatk_hg38.amb
| |-- [445K] gatk_hg38.ann
| |-- [3.0G] gatk_hg38.bwt
| |-- [767M] gatk_hg38.pac
| |-- [1.5G] gatk_hg38.sa
| |-- [6.2K] hg38.bwa_index.log
| |-- [ 0] index.129248.err
| |-- [ 0] index.129248.out
| `-- [ 566] run.sh
|-- [3.2G] dbsnp_146.hg38.vcf.gz
|-- [2.4M] dbsnp_146.hg38.vcf.gz.tbi
|-- [568K] Homo_sapiens_assembly38.dict
|-- [3.0G] Homo_sapiens_assembly38.fasta
|-- [157K] Homo_sapiens_assembly38.fasta.fai
|-- [ 20M] Mills_and_1000G_gold_standard.indels.hg38.vcf.gz
|-- [1.4M] Mills_and_1000G_gold_standard.indels.hg38.vcf.gz.tbi
一般来说下载大文件是需要进行md5检验的,但是我下载这么多次从来没有出过问题,为了简化教程, 就略过这一步咯。
首先测试GATK找germline变异流程
选取一个测试数据,是一个WES时间,介绍如下:
SRX1969919 Illumina HiSeq 2000 PAIRED Hybrid Selection GENOMIC CGU_OSCC_32_WXS_N
Homo sapiens SRP078156 SRR3943929 female
很明显,数据在SRA数据库里面,用下面的代码下载,并且转换
mkdir -p ~/biosoft/GATK/test
cd ~/biosoft/GATK/test
prefetch SRR3943929
# 这个下载速度还是蛮快的,就是有时候链接会无缘无故断开,需要重新下载
ls -lh ~/ncbi/public/sra/
## firtly sra2fastq
mkdir raw
dump='/home/jianmingzeng/biosoft/sratoolkit/sratoolkit.2.8.2-1-centos_linux64/bin/fastq-dump'
sample='CGU_OSCC_32_WXS_N'
$dump -A $sample -O raw --gzip --split-3 ~/ncbi/public/sra/$srr.sra
## 确保这个时候文件转换成功哦
## then qc for fastq
bin_trim_galore="trim_galore"
mkdir -p clean_fastq
$bin_trim_galore -q 25 --phred33 --length 50 -e 0.1 --stringency 3 --paired -o clean_fastq $fq1 $fq2
得到文件如下:
-rw-rw-r-- 1 jmzeng jmzeng 7.5G Mar 22 12:02 CGU_OSCC_32_WXS_N_1_val_1.fq.gz
-rw-rw-r-- 1 jmzeng jmzeng 7.7G Mar 22 13:28 CGU_OSCC_32_WXS_N_2_val_2.fq.gz
得到了clean的fastq 测序数据,接下来就可以走GATK流程啦,请注意看下面的代码,格式都修改成了GATK4要求的啦。
## the mapping and GATK best practice
## 需要设置好 $INDEX $fq1 $fq2
INDEX=/umac/biosoft/GATK/resources/bundle/hg38/bwa_index/gatk_hg38
fq1=clean_fastq/CGU_OSCC_32_WXS_N_1_val_1.fq.gz
fq2=clean_fastq/CGU_OSCC_32_WXS_N_2_val_2.fq.gz
sample='oscc'
## 还需要设置好软件地址
GENOME=/umac/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
GATK=/umac/biosoft/GATK/gatk-4.0.2.1/gatk
DBSNP=/umac/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
nohup ~/biosoft/bwa/bwa-0.7.15/bwa mem -t 5 -M -R "@RG\tID:$sample\tSM:$sample\tLB:WES\tPL:Illumina" $INDEX $fq1 $fq2 1>$sample.sam 2>/dev/null &
## 80G Mar 22 15:22 oscc.sam
nohup $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" SortSam -SO coordinate -I $sample.sam -O $sample.bam 1>log.sort 2>&1 &
samtools index $sample.bam
samtools flagstat $sample.bam > ${sample}.alignment.flagstat
samtools stats $sample.bam > ${sample}.alignment.stat
echo plot-bamstats -p ${sample}_QC ${sample}.alignment.stat
nohup $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" MarkDuplicates -I $sample.bam -O ${sample}_marked.bam -M $sample.metrics 1>log.mark 2>&1 &
nohup $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" FixMateInformation -I ${sample}_marked.bam -O ${sample}_marked_fixed.bam -SO coordinate 1>log.fix 2>&1 &
samtools index ${sample}_marked_fixed.bam
为了节省磁盘空间,最后只需要留学一个bam文件即可:
-rw-rw-r-- 1 jmzeng jmzeng 22G Mar 22 18:44 oscc.bam
-rw-rw-r-- 1 jmzeng jmzeng 4.5M Mar 22 21:28 oscc.bam.bai
-rw-rw-r-- 1 jmzeng jmzeng 2.8K Mar 22 22:54 oscc.metrics
-rw-rw-r-- 1 jmzeng jmzeng 23G Mar 22 22:54 oscc_marked.bam
-rw-rw-r-- 1 jmzeng jmzeng 23G Mar 23 00:21 oscc_marked_fixed.bam
-rw-rw-r-- 1 jmzeng jmzeng 4.5M Mar 23 00:37 oscc_marked_fixed.bam.bai
这3个bam文件的时间消耗如下:
[Thu Mar 22 18:44:53 CST 2018] picard.sam.SortSam done. Elapsed time: 66.31 minutes.
[Thu Mar 22 22:54:14 CST 2018] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 48.91 minutes.
[Fri Mar 23 00:21:05 CST 2018] picard.sam.FixMateInformation done. Elapsed time: 86.68 minutes.
值得注意的是 前面的比对其实耗时也只有一个多小时,但是sortsam耗时比较长,还有简简单单的 samtools index也会耗时十几分钟。
最后就可以直接找变异啦,使用GATK的HaplotypeCaller
命令,如下:
## 首先设置好软件地址
GENOME=/umac/biosoft/GATK/resources/bundle/hg38/Homo_sapiens_assembly38.fasta
GATK=/umac/biosoft/GATK/gatk-4.0.2.1/gatk
DBSNP=/umac/biosoft/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
sample='oscc'
nohup $GATK --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller \
-R $GENOME -I ${sample}_marked_fixed.bam --dbsnp $DBSNP -O ${sample}_raw.vcf 1>log.HC 2>&1 &
最后得到的vcf文件如下;
-rw-rw-r-- 1 jmzeng jmzeng 245M Mar 23 15:28 oscc_1.vcf (对 oscc.bam)
-rw-rw-r-- 1 jmzeng jmzeng 220M Mar 23 15:59 oscc_2.vcf (对 oscc_marked.bam)
-rw-rw-r-- 1 jmzeng jmzeng 220M Mar 23 14:34 oscc_raw.vcf (对 oscc_marked_fixed.bam)
至于这些VCF格式的变异记录该如何理解,以及GATK的每一个步骤是否仍然必须,我就不赘述啦。
请直接围观我的直播基因组活动:
参考资料:
https://gatkforums.broadinstitute.org/gatk/discussion/10502/gatk-4-0-will-be-released-jan-9-2018
https://software.broadinstitute.org/gatk/documentation/topic.php?name=tutorials
云上流程
如果有云服务器,其实可以直接运行GATK团队写好的流程:
居然还有阿里云专属GATK的帮助文档: