GATK4.0和全基因组数据分析实践(上)
前言
一系列WGS文章 没有提供实际可用的数据来实战完成具体的流程,不能得到直观的体会如何入门生物信息学WGS系列第四节 更好地理解整个WGS数据的分析和处理过程
项目目录结构
清晰的目录结构是管理众多项目的有效途径,
./201802_wgs_practice/
├── bin
├── input
└── output
使用E.coli K12完成比对和变异检测
数据比对 变异检测 两个
下载E.coli K12的参考基因组序列
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
$ wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/005/845/GCF_000005845.2_ASM584v2/GCF_000005845.2_ASM584v2_genomic.fna.gz
为了接下来表达上的清晰和操作上的方便,我们使用bgzip将这个序列文件进行解压并把名字重命名为E.coli_K12_MG1655.fa,这样就一目了然了。
gzip -dc GCF_000005845.2_ASM584v2_genomic.fna.gz > E.coli_K12_MG1655.fa
>NC_000913.3 Escherichia coli str. K-12 substr. MG1655, complete genome
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGCTTCTGAACTG
GTTACCTGCCGTGAGTAAATTAAAATTTTATTGACTTAGGTCACTAAATACTTTAACCAATATAGGCATAGCGCACAGAC
$ /Tools/common/bin/samtools faidx E.coli_K12_MG1655.fa
$ samtools faidx E.coli_K12_MG1655.fa NC_000913.3:1000000-1000200
>NC_000913.3:1000000-1000200
GTGTCAGCTTTCGTGGTGTGCAGCTGGCGTCAGATGACAACATGCTGCCAGACAGCCTGA
AAGGGTTTGCGCCTGTGGTGCGTGGTATCGCCAAAAGCAATGCCCAGATAACGATTAAGC
AAAATGGTTACACCATTTACCAAACTTATGTATCGCCTGGTGCTTTTGAAATTAGTGATC
TCTATTCCACGTCGTCGAGCG
这个小技巧在特定的时候非常实用
下载E.coli K12的测序数据
基因组参考序列
https://www.ncbi.nlm.nih.gov/sra/?term=SRR1770413$ wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR177/SRR1770413/SRR1770413.sra
注意,下载下来的这个SRA文件虽然只有一份,但是里面其实存了read1和read2的测序数据,我们要把它解出来,转换成为我们所需的fastq格式
$ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413.sra
SRR1770413_1.fastq
SRR1770413_2.fastq
$ /Tools/sratoolkit/2.8.2/bin/fastq-dump --split-files SRR1770413
$ /Tools/common/bin/bgzip -f SRR1770413_1.fastq
$ /Tools/common/bin/bgzip -f SRR1770413_1.fastq
./201802_wgs_practice/
├── bin
├── input
│ ├── E.coli
│ ├── fasta
│ │ ├── E.coli_K12_MG1655.fa
│ │ ├── E.coli_K12_MG1655.fa.fai
│ │ └── work.log.sh
│ └── fastq
│ ├── SRR1770413_1.fastq.gz
│ ├── SRR1770413_2.fastq.gz
│ └── work.log.sh
├── output
└── work.log.sh
work.log.sh记录了我在该目录下的所有重要操作
质控
比对
$ /Tools/common/bin/bwa index E.coli_K12_MG1655.fa
[bwa_index] Pack FASTA... 0.06 sec
[bwa_index] Construct BWT for the packed sequence...
[bwa_index] 1.78 seconds elapse.
[bwa_index] Update BWT... 0.04 sec
[bwa_index] Pack forward-only FASTA... 0.03 sec
[bwa_index] Construct SA from BWT and Occ... 0.63 sec
[main] Version: 0.7.16a-r1181
[main] CMD: /Tools/common/bin/bwa index E.coli_K12_MG1655.fa
[main] Real time: 3.266 sec; CPU: 2.550 sec
E.coli_K12_MG1655.fa.amb
E.coli_K12_MG1655.fa.ann
E.coli_K12_MG1655.fa.bwt
E.coli_K12_MG1655.fa.pac
E.coli_K12_MG1655.fa.sa
现在我们使用bwa完成比对,用samtools完成BAM格式转换、排序并标记PCR重复序列。步骤分解如下:
#1 比对
time /Tools/common/bin/bwa mem -t 4 -R '@RG\tID:foo\tPL:illumina\tSM:E.coli_K12' /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_1.fastq.gz /Project/201802_wgs_practice/input/E.coli/fastq/SRR1770413_2.fastq.gz | /Tools/common/bin/samtools view -Sb - > /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** bwa mapping done **"#2 排序
time /Tools/common/bin/samtools sort -@ 4 -m 4G -O bam -o /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam && echo "** BAM sort done"rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam
#3 标记PCR重复
time /Tools/common/bin/gatk/4.0.1.2/gatk MarkDuplicates -I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam -O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam -M /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup_metrics.txt && echo "** markdup done **"#4 删除不必要文件(可选)
rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.bam
rm -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.bam
#5 创建比对索引文件
time /Tools/common/bin/samtools index /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam && echo "** index done **"
-R 设置Read Group信息,虽然我在以前的文章中已经反复强调过它的重要性,但这里还是再说一次,它是read数据的组别标识,并且其中的ID,PL和SM信息在正式的项目中是不能缺少的(如果样本包含多个测序文库的话,LB信息也不要省略),另外由于考虑到与GATK的兼容关系,PL(测序平台)信息不能随意指定,必须是:ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这12个中的一个。
变异检测
$ /Tools/common/bin/gatk/4.0.1.2/gatk CreateSequenceDictionary -R E.coli_K12_MG1655.fa -O E.coli_K12_MG1655.dict && echo "** dict done **"
#1 生成中间文件gvcf
time /Tools/common/bin/gatk/4.0.1.2/gatk HaplotypeCaller \
-R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
--emit-ref-confidence GVCF \
-I /Project/201802_wgs_practice/output/E.coli/E_coli_K12.sorted.markdup.bam \
-O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf && echo "** gvcf done **"
#2 通过gvcf检测变异
time /Tools/common/bin/gatk/4.0.1.2/gatk GenotypeGVCFs \
-R /Project/201802_wgs_practice/input/E.coli/fasta/E.coli_K12_MG1655.fa \
-V /Project/201802_wgs_practice/output/E.coli/E_coli_K12.g.vcf \
-O /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf && echo "** vcf done **"
变异检测不是一个样本的事情,有越多的同类样本放在一起joint calling结果将会越准确,而如果样本足够多的话,在低测序深度的情况下也同样可以获得完整并且准确的结果
#1 压缩
time /Tools/common/bin/bgzip -f /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf
#2 构建tabix索引
time /Tools/common/bin/tabix -p vcf /Project/201802_wgs_practice/output/E.coli/E_coli_K12.vcf.gz
bgzip压缩完成之后,原来的VCF文件会被自动删除。
为了保持一致,现在再看一下完成到这里之后我们的目录长什么样了,供大家对照。
./201802_wgs_practice/
├── bin
│ ├── bwa_and_markdup.sh
│ └── gatk.sh
├── input
│ └── E.coli
│ ├── fasta
│ │ ├── E.coli_K12_MG1655.dict
│ │ ├── E.coli_K12_MG1655.fa
│ │ ├── E.coli_K12_MG1655.fa.amb
│ │ ├── E.coli_K12_MG1655.fa.ann
│ │ ├── E.coli_K12_MG1655.fa.bwt
│ │ ├── E.coli_K12_MG1655.fa.fai
│ │ ├── E.coli_K12_MG1655.fa.pac
│ │ ├── E.coli_K12_MG1655.fa.sa
│ │ └── work.log.sh
│ └── fastq
│ ├── SRR1770413_1.fastq.gz
│ ├── SRR1770413_2.fastq.gz
│ └── work.log.sh
├── output
│ └── E.coli
│ ├── E_coli_K12.g.vcf
│ ├── E_coli_K12.g.vcf.idx
│ ├── E_coli_K12.sorted.markdup.bam
│ ├── E_coli_K12.sorted.markdup.bam.bai
│ ├── E_coli_K12.sorted.markdup_metrics.txt
│ ├── E_coli_K12.vcf
│ └── E_coli_K12.vcf.idx
└── work.log.sh
小结
在我不断回答公众号后台一个一个的问题之后,我强烈意识到需要有一个地方,来把和朋友们共同讨论的有价值内容汇集起来。于是我在知识星球上开通了一个圈子,名字是:解螺旋技术交流圈,这是与读者们的私人朋友圈,圈子中的牛人已经越来越多了,它是付费的,这是知识星球上第一个真正与基因组学和生物信息学强相关的圈子。