利用gatk进行变异检测的n步走战略之二
生信人R语言学习必备
立刻拥有一个Rstudio账号
开启升级模式吧
（56线程，256G内存，个人存储1T）
hello！小果又来啦！
前情回顾：上期我们利用gatk生成了索引文件并对sam文件进行了排序，折起我们接着往下走。
一个小tip：索引文件最好要和参考基因组的fasta文放到同一个目录哦！不然后面会出错的，小果的血泪教训！
好了我们继续吧，现在我们需要利用GATK的MarkDuplicates工具去除PCR重复，并将结果存储到指定的输出文件中。
for sample in `cat /media/desk16/iyun003/GBS_test/ERR_SRA/sample.txt`; do #对sample.txt中的每一个样本循环操作echo ${sample} #输出当前样本的名称gatk --java-options "-Xmx16g -Djava.io.tmpdir=./tmp" MarkDuplicates #调用GATK的MarkDuplicates工具-I /media/desk16/iyun003/GBS_test/ERR_SRA/tair_sam/${sample}.sort.bam #输入文件为排序后的BAM文件-M /media/desk16/iyun003/GBS_test/ERR_SRA/tair_sam/sort/${sample}.metrics --CREATE_INDEX #指定输出重复信息的文件以及输出文件的创建-O /media/desk16/iyun003/GBS_test/ERR_SRA/tair_sam/sort/${sample}.sort.MarkDup.bam #输出文件名done
其中，${sample}代表sample.txt文件中每一行所存储的样本名称。
-I选项用于指定输入的BAM文件，-M选项指定输出重复信息的文件，–CREATE_INDEX选项指定输出文件的同时也创建其对应的索引文件，-O选项指定输出文件名。
–java-options “-Xmx16g -Djava.io.tmpdir=./tmp”选项则是在运行GATK时为Java虚拟机指定参数，其中-Xmx16g指定了最大堆内存为16GB，-Djava.io.tmpdir=./tmp指定了临时文件的目录为当前工作目录下的tmp文件夹。
接下来我们对每个样本进行HaplotypeCaller变异检测，并生成GVCF格式的文件，注意替换自己的文件路径：
REF=/media/desk16/iyun003/download/ninanjie_data/GCF_000001735.4_TAIR10.1_genomic.fnafor sample in `cat /media/desk16/iyun003/GBS_test/ERR_SRA/sample.txt`;doecho "$sample"gatk --java-options "-Xmx4g" HaplotypeCaller-R "${REF}"-I "/media/desk16/iyun003/GBS_test/ERR_SRA/tair_sam/sort/${sample}.sort.MarkDup.bam"-O "/media/desk16/iyun003/GBS_test/ERR_SRA/gvcf/${sample}.g.vcf.gz"-ERC GVCFdone
其中：REF 指定参考基因组文件路径
sample 指定样本名称，通过读取 sample.txt 文件获取gatk --java-options "-Xmx4g" HaplotypeCaller 运行 GATK 的 HaplotypeCaller 变异检测程序，--java-options "-Xmx4g" 指定程序运行时最大占用4GB内存-R "${REF}" 指定参考基因组文件-I "/media/desk16/iyun003/GBS_test/ERR_SRA/tair_sam/sort/${sample}.sort.MarkDup.bam" 指定输入BAM文件路径-O "/media/desk16/iyun003/GBS_test/ERR_SRA/gvcf/${sample}.g.vcf.gz" 指定输出文件路径及文件名-ERC GVCF 指定输出为GVCF格式这一步运行时间较久，可以通过查看运行日志来了解运行情况哦，小果把它放到脚本里后台运行，就不用担心意外中断了。nohup bash gvcf.sh > gvcf.log &
查看日志cat gvcf.log
小果这边放一下运行过程中的结果，可以看到其他结果还没有出来。
好啦，今天的内容暂时就到这里了，我们下期继续！
欢迎使用：云生信 – 学生物信息学 (biocloudservice.com)
如果想用服务器可以联系微信：18502195490（快来联系我们使用吧！）
扫码加小果
领取生信大礼包
点击“阅读原文”立刻拥有
↓↓↓