GATK(全称The Genome Analysis Toolkit)是Broad Institute开发的用于二代重测序数据分析的一款软件,里面包含了很多有用的工具。
需要的文件
参考基因组:reference.fasta
测序文件:bam 格式比对结果文件
所需软件: samtools 、gatk 、 Picard
一、创建gatk所需的索引文件
GATK需要reference序列是经过index的,而且需要两个index文件,一个是后缀名为.fai的,另外一个是后缀名称为.dict的,缺少这些文件,或者两个文件中的内容不一致都可能导致程序报错。
samtools faidx reference.fasta #生成后缀名为.fai
gatk CreateSequenceDictionary-R reference.fasta -O reference.dict #生成后缀名为.dict
二、 增加GATk 要求的read group的格式
GATK要求输入的bam文件包含Read groups,如果没有就会报错。
Read group是@RG开始,包括以下几个部分:
ID= Read group identifier
每一个Read group独有的ID;
PU= Platform Unit
PL= Platform/technology used to produce the read
测序使用的平台: ILLUMINA, SOLID, LS454, HELICOS and PACBIO。
LB= DNA preparation library identifier
对一个read group的reads进行重复序列标记时,需要使用LB来区分reads来自那条lane;有时候,同一个库可能在不同的lane上完成测序;为了加以区分,同一个或不同库只要是在不同的lane产生的reads都要单独给一个ID。
SM= Sample
reads属于的样品名,自由设定
可以在BWA比对时增加read group:
bwa mem -R '@RG\tID:group\tLB:library\tPL:illumina\tPU:unit1\tSM:676R' ~/ref/reference.fasta read1.fq read2.fq > bulk.sam
或者使用Picard增加:
Picard: 它是目前最著名的组学研究中心-Broad研究所开发的一款强大的NGS数据处理工具,功能方面和Samtools有些重叠,但更多的是互补,它是由java编写的,我们直接下载最新的.jar包就行了。
下载链接:wget https://github.com/broadinstitute/picard/releases/download/2.25.5/picard.jar
安装
java -jar picard.jar
执行命令:
java -jar ~/biosoft/picard.jar AddOrReplaceReadGroups I=bulk.bam O=bulk.add.bam RGID=4 RGLB=library1 RGPL=illumina RGPU=unit1 SORT_ORDER=coordinate RGSM=M05
三、去除PCR重复序列
gatk MarkDuplicates -I bulk.add.bam -O bulk.marked.bam -M bu.metrics 1>log.mark 2>&1
四、VCF输出
#先生成gvcf格式文件 gvcf可记录所有位点的变异情况
gatk HaplotypeCaller -R ~/ref/reference.fasta -I bulk.marked.bam -O output.g.vcf.gz -ERC GVCF
#然后在从gvcf提取变异情况
gatk GenotypeGVCFs -R ~/ref/reference.fasta -V output.g.vcf.gz -O output.vcf.gz
或者直接生成vcf文件
gatk HaplotypeCaller -R ~/ref/reference.fasta -I bulk.marked.bam -O out.vcf
五、执行筛选
gatk VariantFiltration \
-V out.vcf \
-filter "QD < 2.0" --filter-name "QD2" \
-filter "QUAL < 30.0" --filter-name "QUAL30" \
-filter "FS > 60.0" --filter-name "FS60" \
-filter "MQ < 40.0" --filter-name "MQ40" \
-filter "MQRankSum < -12.5" --filter-name "MQRankSum-12.5" \
-filter "ReadPosRankSum < -8.0" --filter-name "ReadPosRankSum-8" \
-O out_prefix.vcf
六、合并VCF文件
第一种方法:
gatk -T CombineVariants -V file1.vcf.gz -V file2.vcf.gz -o merge.vcf.gz -R ref.fa
第二种方法:
bcftools merge file1.vcf.gz file2.vcf.g -o merge_bcftools.vcf
网友评论