用GATK进行二代测序数据 SNP Calling 流程：（三）GenomicsDBImport 的多样本变异检测

1. Genomics Database

对于群体数据来说，多样本同时时行 SNP Calling 的准确度要优于单个样本的 SNP Calling.
GATK3 的多样本 SNP Calling 功能是 CombineGVCFs，GATK4 新出了 GenomicsDBImport功能，官网建议它适合1000个样本以上的 SNP Calling，但是它的另一个优点是可扩展性，即随时可以向 database 里添加新的材料，以扩大群体数量，而不用对旧的数据再从头操作一次。

Tips：GenomicsDBImport 会按 intervals 进行变异检测，即每计算一个 interval 都会将所有的 GVCF 文件打开和检索一次，当样本量很大或intervals 很多时，打开和关闭 GVCF 文件会消耗大量的时间。因此，对基因组组装不好的物种，如存在大量 Scarffolds，建议仅对染色体进行 SNP Calling，或者将部分 Scarffolds 合并后再进行 SNP Calling。

gatk=gatk
# 必需的文件，具体格式参考GATK
interval_list=chr.list
# Java内存限制
Xmx=100g
Xms=10g
# genomicsdb，必须是未创建的或者空的文件夹
my_database=Genomicsdb
# 每批次读入多少样本，默认 0, 表示一次性全部读入。注意，一次读入过多样本会发生错误，建议每次读入不多于100个;
batch_size=1
# gvcf文件列表，每行格式为“gvcf文件名	路径”，具体格式参考GATK
map_file=all.map
# 每个batch的并行线程
reader_threads=1
#同时读入的最大intervals数量,值越高，速度越快，但消耗更多的内存，且需要有权限创建足够多的文件.
num=26
#缓存文件
tmp=tmpdir2


#第一次建库用的参数：
# --genomicsdb-workspace-path $my_database \
#向已存在的库添加新数据的参数：
#--genomicsdb-update-workspace-path $my_database \

$gatk --java-options "-Xmx$Xmx -Xms$Xms" \
GenomicsDBImport \
--genomicsdb-workspace-path $my_database \
--intervals $interval_list \
--batch-size  $batch_size \
--sample-name-map $map_file \
--reader-threads $reader_threads \
--max-num-intervals-to-import-in-parallel $num \
--tmp-dir ${tmp}

2. 分染色体SNP Calling

#gatk4.0版

#gatk路径
gatk=gatk
#参考基因组
reference=genome.fa
#GenomicsImport产生的数据库
mydatabase=Genomicsdb
#输出文件名称
output=raw.vcf
#intervals文件
intervals=chr.list
#运行内存
Xmx=20g
Xms=1g
#缓存文件夹
tmp=tmpdir

cat $intervals|while read id
do
	if [ ! -f "${tmp}/chr_${id}" ];then

		mkdir -p ${tmp}/chr_${id}
	else
		rm -r ${tmpdir}/chr_${id}/*
	fi
nohup $gatk --java-options "-Xmx$Xmx -Xms$Xms" GenotypeGVCFs \
-R $reference \
-V gendb://$mydatabase \
-L $id \
-O $output \
--tmp-dir ${tmp}/chr_${id} \
>chr_${id}.log 2>&1 &

done <$1

3. 合并VCF

nohup picard MergeVcfs I=chr_vcf.list O=raw.vcf.gz TMP_DIR=tmpdir/ >mergevcfs.log 2>&1 &