用GATK进行二代测序数据 SNP Calling 流程：（二）bwa比对和HaplotypeCaller 变异检测

1. 创建基因组索引

bwa index genome.fa

2. 查看read group信息，按read group分组, 比对、合并，生成gvcf

由于数据太多，无法存储过多的中间文件，因此写了一个脚本，边运行边删除中间文件，过程包括：

1. 解压，按read group分组。（RG（read group） 信息非常重要，GATK需要通过RG来判断碱基测序质量。我的一个样品的测序数据可能会来自不同的title，不同的lane、flowcell，甚至不同的机器，这在重测序中比较常见。因此，我将一个fastq文件按RG分割，分别比对、添加RG信息。
2. 比对，添加RG信息
3. 用picard软件合并来自同一个材料的sam文件，并输出排序后的bam文件，建立索引。
4. 用picard软件注释 Duplicates 信息。特殊的建库方式，如GBS、RAD-seq、RNA-seq数据可以忽略，因为它们或者用酶切建库，基因组打断的位置不随机，或者基于转录本建库，本身的表达丰度不同，都不适合去重。
5. 用GATK4 进行SNP Calling，因为我是用群体 call SNP的策略，所以这里产生的是GVCF文件。

Tips：

• gzip压缩文件和sra压缩文件大小相近。
• 一个fastq文件，建议给予其大小5倍的存储空间，以免存储空间不足，导致程序强行终止。
• 如果数据很大，如重测序数据，或者同时处理的数据很多，建议给bwa、picard、gatk建立临时文件夹，如果用默认的系统临时文件夹/tmp，可能会存储空间不足，导致程序强行终止。

以下是运行代码：

#参考基因组(绝对路径)
Ref=genome.fa

#bwa 参数
#线程
t=10
#索引
FM_index="genome.fa"

#picard 参数
#最大运行内存
p_Xmx=20g

#GATK4参数
# 临时文件夹（绝对路径）
tmpdir="tmpdir"
# bam 文件（绝对路径）
bam_dir="."
# java 运行内存
G_Xmx=20g


#对fastq文件按reads group分组
while read id
do
	#解压
	gunzip -c /database/public/reseq355_clean/${id}_1.fq.gz >${id}_1.fq
	gunzip -c /database/public/reseq355_clean/${id}_2.fq.gz >${id}_2.fq
	#提取read group信息
	grep "^@ST" ${id}_1.fq|cut -d ':' -f1-4|sort|uniq >${id}.RG.table

	while read rg
	do
		affix=`echo $rg|tr ":" "-"`
		#根据read group分隔文件
		grep -A 3 ${rg} ${id}_1.fq >${id}-${affix}_1.fq
		grep -A 3 ${rg} ${id}_2.fq >${id}-${affix}_2.fq
	done <${id}.RG.table
    #删除fastq文件，释放存储空间
	rm ${id}_1.fq ${id}_2.fq


	# bwa 比对
	#比对成功后，删除相应fastq文件
	while read rg
	do
		affix=`echo $rg|tr ":" "-"`
		bwa mem -t $t -M -R "@RG\tID:${affix}\tSM:${id}\tLB:${id}\tPL:illumina" $FM_index ${id}-${affix}_1.fq ${id}-${affix}_2.fq >${id}-${affix}.sam && rm ${id}-${affix}_*.fq
	done <${id}.RG.table

	
	
	#用picard合并sam，排序，输出bam,index
	ls ${id}-*.sam|xargs -I [] echo "I="[]|xargs -L 1000000 \
	picard -Xmx${p_Xmx} MergeSamFiles \
	O=${id}.merged.sorted.bam \
	SORT_ORDER=coordinate \
	CREATE_INDEX=true \
	VALIDATION_STRINGENCY=LENIENT \
	REFERENCE_SEQUENCE=$Ref	\
	TMP_DIR=${tmpdir} 2>${id%}.log \
	&&rm ${id}-*.sam
	
 # 用picard 进行Mark duplicates
	picard -Xmx${p_Xmx} MarkDuplicates \
	I=${id}.merged.sorted.bam \
	O=${id}.merged.sorted.dedup.bam \
	M=${id}.Marked_dup_metrics.txt \
	CREATE_INDEX=true \
	VALIDATION_STRINGENCY=LENIENT \
	REFERENCE_SEQUENCE=$Ref	\
	TMP_DIR=${tmpdir} 2>${id%}.log \
	&& rm ${id}.merged.sorted.bam*
	
#用gatk HaplotypeCaller进行SNP Calling	
#“-DF NotDuplicateReadFilter”：使用该参数则禁止对reads duplicate进行过滤，GBS或RNA-seq数据可能需要。
# 为gatk建立临时存储文件夹
	if [ ! -f "${tmpdir}/${id}_h" ];then
		mkdir -p ${tmpdir}/${id}_h
	else
		rm -r ${tmpdir}/${id}_h/*
	fi
	
	gatk --java-options "-Xmx${G_Xmx} -Djava.io.tmpdir=${tmpdir}/${id}_h" \
	HaplotypeCaller \
	-ERC GVCF \
	-R ${Ref} \
	-I ${id}.merged.sorted.dedup.bam \
	-O ${id}.gvcf.gz \
	1>${id}.gvcf.log 2>&1 \
	&& rm -r ${tmpdir}/${id}_h 
	

done < $1