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

查看原文

其他

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

Original 曾建明 biobabble 2020-02-03

很容易谷歌搜索就可以找到最新版教程：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的帮助文档：

https://help.aliyun.com/document_detail/60414.html?spm=5176.11065259.1996646101.searchclickresult.14c51b0bIlpQgI

投稿系列

欢迎大家来投稿:)

猜你喜欢： 生信基础知识100讲

常德悲剧：让谴责无差别杀戮之声更加响亮一点

魏加宁：日本之所以能走出“大衰退”，靠的是不断改革，而不是所谓“积极的财政政策”

国内突然宣布，将大规模调整经济布局！又一朝阳行业悄然崛起，新一轮机会出现了！

这家头部券商资管总经理"迎新"

泪目！8死17伤！江苏一职校持刀伤人案，背后隐情令人心惊！