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

1. 创建基因组索引

bwa index genome.fa

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

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

解压，按read group分组。（RG（read group）信息非常重要，GATK需要通过RG来判断碱基测序质量。我的一个样品的测序数据可能会来自不同的title，不同的lane、flowcell，甚至不同的机器，这在重测序中比较常见。因此，我将一个fastq文件按RG分割，分别比对、添加RG信息。
比对，添加RG信息
用picard软件合并来自同一个材料的sam文件，并输出排序后的bam文件，建立索引。
用picard软件注释 Duplicates 信息。特殊的建库方式，如GBS、RAD-seq、RNA-seq数据可以忽略，因为它们或者用酶切建库，基因组打断的位置不随机，或者基于转录本建库，本身的表达丰度不同，都不适合去重。
用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.fqgunzip -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.tablewhile read rgdoaffix=`echo $rg|tr ":" "-"`#根据read group分隔文件grep -A 3 ${rg} ${id}_1.fq >${id}-${affix}_1.fqgrep -A 3 ${rg} ${id}_2.fq >${id}-${affix}_2.fqdone <${id}.RG.table#删除fastq文件，释放存储空间rm ${id}_1.fq ${id}_2.fq# bwa 比对#比对成功后，删除相应fastq文件while read rgdoaffix=`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}_*.fqdone <${id}.RG.table#用picard合并sam，排序，输出bam,indexls ${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 duplicatespicard -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" ];thenmkdir -p ${tmpdir}/${id}_helserm -r ${tmpdir}/${id}_h/*figatk --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

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

用GATK进行二代测序数据 SNP Calling 流程：（四）变异过滤
GATK推荐的最好的过滤方式是用 VQSR功能,它通过机器学习算法来判断SNP的优劣,因此至少需要两个已存在的 SNP 数据集,一个是经过验证的高质量 SNP 数据集作为真集(如 HapMap),还需 ...
【Sentieon】PacBio HiFi三代测序数据SNP/Indel加速分析
Sentieon软件在二代测序中SNP/Indel变异检测流程已非常成熟,并以其检测准确性高和检测速度快而广受业内人士认可.近日,Sentieon推出了DNAscope LongReads分析流程,深 ...
关于文献中二代测序数据下载（NCBI）的问题
关于文献中二代测序数据下载(NCBI)的问题现在二代测序用于生物学研究非常广泛,大部分文章的序列会上传到Sequence Read Archive(SRA)上,这东西也属于NCBI数据库中的吧,我理 ...
DNA甲基化测序数据的分析流程及相关软件总结
目前检测DNA甲基化的方法众多,主要可以分为以下几类(如表1所示): 图片来源(凡时财等,中国科学: 生命科学,2015) <更多精彩,可关注微信公众号:AIPuFuBio,和大型免费综合生物信 ...
二代测序数据统计分析中为什么是负二项分布？
本文转载自"universebiologygirl",已获授权 1. 转录组数据统计推断的难题在RNA-seq中进行两组间的差异分析是最正常不过的了我们在其它实验中同样会遇到类 ...
二代测序linux软件,二代测序数据分析软件包大全
二代测序数据分析软件包大全 Integrated solutions*CLCbio Genomics Workbench-de novoand reference assembly of Sanger ...
三代测序数据纠错的方法、装置和计算机可读存储介质与流程
三代测序数据纠错的方法.装置和计算机可读存储介质与流程文档序号:15616049发布日期:2018-10-09 21:24 导航: X技术> 最新专利>计算;推算;计数设备的制造及其应用 ...
经典：基因组测序数据从头拼接或组装算法的原理
基因组测序数据的拼接/组装 (图片来源:google) 每一个物种的参考基因组序列(reference genome)的产生都要先通过测序的方法,获得基因组的测序读段(reads),然后再进行从头拼接 ...
二代测序的原理和简介
二代测序的简介第二代测序(Next-generation sequencing,NGS)又称为高通量测序(High-throughput sequencing),是基于PCR和基因芯片发展而来的 ...

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

1. 创建基因组索引

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

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

最新文章

热门文章