鉴于上次发表的"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)