查看原文
其他

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

曾建明 biobabble 2020-02-03

很容易谷歌搜索就可以找到最新版教程:https://software.broadinstitute.org/gatk/documentation/quickstart,我就不做文档搬运工了,感兴趣自己看原文,目录如下:

  1. Quick start for the impatient

  2. Requirements

  3. Get GATK

  4. Install it

  5. Test that it works

  6. Run GATK and Picard commands

  7. Grok the Best Practices

  8. Run pipelines

  9. 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讲

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存