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

博文

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

已有 1073 次阅读 2021-6-11 16:36 |系统分类:科研笔记

鉴于上次发表的"GATK流程不用再分割小麦染色体",有些内容说的不是很清楚, 且科学网博客发表之后无法更改,遂再次更新一次。

对小麦来说,染色体较大,bam index不能使用.bai,所以要使用.csi。GATK支持.csi格式。

在运行过程中,HaplotypeCaller这一步,需要加入参数 -OVI False,即禁止gatk自动生成index文件。GatherVcfsCloud 这一步同样需要设置 -OVI False。SortVcf 这一步需要设置 --CREATE_INDEX False。

该流程输入文件是fastq.gtz,输出文件为bwa mapping的结果文件以及gvcf文件。其中mapping结果文件,为节省存储空间以cram格式保存。

另外为了节省空间的需要,fastq均使用了gtz格式压缩,所以bwa-mem也使用了能读gtz格式的版本。具体请参见https://github.com/Genetalks/gtz#bwa-mem2。如果fastq文件使用普通的gz格式压缩,那就使用正常的bwa-mem就可以。

关于genome intervals文件,格式是“chrom:star-pos”。一般是在gap的位置进行分割。根据个人的实际情况,可将整个小麦基因组切割成N个区间。

gatk --versionUsing 
GATK jar /public1/home/sca2417/software/gatk-4.2.0.0/gatk-package-4.2.0.0-local.jar
Running:    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /public1/home/sca2417/software/gatk-4.2.0.0/gatk-package-4.2.0.0-local.jar --version
The Genome Analysis Toolkit (GATK) v4.2.0.0
HTSJDK Version: 2.24.0
Picard Version: 2.25.0


samblaster --version
samblaster: Version 0.1.26


samtools --version
samtools 1.12
Using htslib 1.12
Copyright (C) 2021 Genome Research Ltd.

具体的运行脚本如下。注意,请根据自己的服务器的内存、缓存和CPU等调整相应的参数。

#!/usr/bin/env python
# -*- coding: utf-8 -*-
__author__ = 'shengwei ma'
__author_email__ = 'shengweima@icloud.com'

import argparse

import subprocess

from concurrent.futures import ProcessPoolExecutor
from os.path import isdir, isfile, join, exists

def gatk4(chrom):
    subprocess.run('gatk --java-options "-Xmx20G -Djava.io.tmpdir=./" HaplotypeCaller -I ' + args.output + '.mrkdup.sorted.bam -OVI False --minimum-mapping-quality 20 --read-filter ProperlyPairedReadFilter -O ' + args.output + '.' + chrom.replace(':','.') + '.g.vcf.gz -R ' + args.ref + ' --emit-ref-confidence GVCF -L ' + chrom,shell=True)


parser = argparse.ArgumentParser(description='输入参数如下:')
parser.add_argument('--input1', '-1', help='*_1_fastq.gtz file,必要参数', required=True)
parser.add_argument('--input2', '-2', help='*_2_fastq.gtz file,必要参数', required=True)
parser.add_argument('--ref', '-R', help='reference genome,必要参数', required=True)
parser.add_argument('--output', '-O', help='sample name,必要参数', required=True)
parser.add_argument('--process', '-t', type =int, help='process number,必要参数', required=True)
parser.add_argument('--intervals', '-L', help='gatk4 interval file,必要参数', required=True)
args = parser.parse_args()

if __name__ == '__main__':
    commom_wheat_chrom = []
    with open(args.intervals, 'r') as f:
        for line in f:
            commom_wheat_chrom.append(line.strip())

    if exists(args.output +'g.sorted.vcf.gz'):
        pass
    else:
        # bwa mapping and filter amd mrdup duplicated. bwa-mem and samtools should have the same number thread (-t and -@).
        if exists(args.output + '.mrkdup.sorted.bam.csi'):
            pass
        else:
            subprocess.run('bwa-mem2-gtz mem -K 100000000 -Y -t 64 -R ' + "'@RG\\tID:" + args.output + "\\tSM:" + args.output + "\\tPL:ILLUMINA'" + ' -M '+'-V ' + args.ref + ' ' + args.input1 + ' ' + args.input2 + '| samblaster | samtools sort -m 500M -@ 64 -o ' +  args.output + '.mrkdup.sorted.bam' + ' - ',shell=True)

            # bam index
            subprocess.run('samtools index -@ 64 -c ' + args.output + '.mrkdup.sorted.bam',shell=True)


            chrom_gvcf = []
            for chrom in commom_wheat_chrom:
                chrom_gvcf.append(args.output + '.' + chrom.replace(':','.') + '.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)
    
        # combine gvcf file
            subprocess.run('gatk --java-options "-Xmx60G -Djava.io.tmpdir=./" GatherVcfsCloud -OVI False -I ' + ' -I '.join(chrom_gvcf) + ' -O ' + args.output + '.g.vcf.gz',shell=True)
        # sort gvcf file
            subprocess.run('gatk --java-options "-Xmx60G -Djava.io.tmpdir=./" SortVcf --CREATE_INDEX False -I ' +  args.output + '.g.vcf.gz' + ' -O ' + args.output + '.g.sorted.vcf.gz',shell=True)
        #remove chr.g.vcf.gz
            subprocess.run('rm -fr ' + ' '.join(chrom_gvcf),shell=True)
        # get cram file
            subprocess.run('samtools view -@ 64 -C -T ' + args.ref + ' -o ' + args.output + '.mrkdup.sorted.cram ' + args.output + '.mrkdup.sorted.bam',shell=True)
        # remmove bam file
            subprocess.run('rm -fr *.bam*',shell=True)





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

上一篇:NCBI上传fastq.gz以及bam文件有文件大小限制
下一篇:翻译:植物泛基因组综述(Plant pan-genomes are the new reference)

0

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

数据加载中...

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

GMT+8, 2021-10-20 15:34

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部