mashengwei的个人博客分享 http://blog.sciencenet.cn/u/mashengwei

博文

GATK流程不用再分割小麦染色体

已有 5426 次阅读 2020-9-20 22:17 |系统分类:科研笔记| 小麦基因组

小麦GATK流程不支持小麦整条染色体,主要是指不支持bam文件和vcf文件的index。所以把index问题解决了,问题是不是就解决了呢?
GATK有些步骤会自动index,需要添加参数禁掉 --disable_bam_indexing。
GATK4支持csi的index,所以bam文件可以使用samtools index -c生成.csi格式的index文件。如果还是想要使用.bai格式的文件,建议安装和使用这个把版本的samtools(1.3.1)。
至于vcf的index,使用 bcftools index -t生成.tbi格式的index文件。bcftools的版本也是1.3.1。或者使用tabix,版本也是1.3.1。
这里推荐使用GATK4

remove duplicates

请先用samtools对bam文件进行index,生成csi或者bai格式都可以。
gatk MarkDuplicates -I test.filter.bam --REMOVE_DUPLICATES TRUE -M out.metrics -O test.rmdup.filter.bam
生成的test.rmdup.filter.bam文件也需要进行index。

HaplotypeCaller

gatk --java-options -Xmx8G HaplotypeCaller --emit-ref-confidence GVCF -OVI False -R ref.fa -I test.rmdup.filter.bam -O test.g.vcf.gz
为了加快运行速度,可以使用-L参数,我写了一个脚本来完成这个事。
python gatk4.py --process 8 --input test.rmdup.filter.bam --ref /data2/Fshare/FastaAndIndex/IWGSC_v1.0_STAR/IWGSC_v1.0.fasta --output test

gatk4.py内容如下

#!/usr/bin/ python
# -*- coding: utf-8 -*-
__author__ = 'shengwei ma'
__author_email__ = 'shengweima@icloud.com'
import argparse
import subprocess
from concurrent.futures import ProcessPoolExecutor


def gatk4(chrom):
    subprocess.run(['gatk','--java-options','-Xmx8G','HaplotypeCaller','-I',args.input,
                    '-O',args.output + '.' + chrom + '.g.vcf.gz','-R',args.ref,'--emit-ref-confidence','GVCF','-OVI','False''-L',chrom],shell=False)


parser = argparse.ArgumentParser(description='输入参数如下:')
parser.add_argument('--process''-P', type=int, help='进程数量,必要参数', required=True)
parser.add_argument('--input''-I', help='bam file,必要参数', required=True)
parser.add_argument('--ref''-R', help='reference genome,必要参数', required=True)
parser.add_argument('--output''-O', help='sample name,必要参数', required=True)
args = parser.parse_args()

if __name__ == '__main__':
    commom_wheat_chrom = ['chr1A','chr2A','chr3A','chr4A','chr5A','chr6A','chr7A',
                         'chr1B','chr2B','chr3B','chr4B','chr5B','chr6B','chr7B',
                         'chr1D','chr2D','chr3D','chr4D','chr5D','chr6D','chr7D','chrUn']
    chrom_gvcf = []
    for chrom in commom_wheat_chrom:
        chrom_gvcf.append( args.output + '.' + chrom + '.g.vcf.gz')
    try:
        with ProcessPoolExecutor(max_workers=args.process) as pool:
            future1 = pool.map(gatk4,commom_wheat_chrom)
            print(list(future1))
    except Exception as e: 
        print(e)

    subprocess.run('gatk GatherVcfs -RI TRUE -I ' + ' -I '.join(chrom_gvcf) + ' -O ' + args.output + '.g.vcf.gz',shell=True)
    subprocess.run('rm -fr ' + ' '.join(chrom_gvcf),shell=True)

最后生成的是test.g.vcf.gz 多个sample还需要使用gatk CombineGVCFs
gatk CombineGVCFs -OVI False -R ref.fa -O test.merge.g.vcf.gz -V test.parents1.g.vcf.gz -V test.parents2.g.vcf.gz -V test.pool1.g.vcf.gz -V test.parents2.g.vcf.gz

生成vcf文件。gatk GenotypeGVCFs -OVI False -R ref.fa -V test.merge.g.vcf.gz -O test.merge.vcf.gz

这个过程中生成的vcf.gz 或bam文件均需要使用开头所说的那样生成index文件。




http://blog.sciencenet.cn/blog-1094241-1251397.html

上一篇:综述|植物泛组学:方法、应用与进展
下一篇:浅谈转基因生物的安全性

0

该博文允许注册用户评论 请点击登录 评论 (1 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-11-30 02:42

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部