本篇帖子的分享內(nèi)容主要也是和前面大基因組拆分相呼應(yīng)的，利用前面拆分的基因組文件，對變異檢測結(jié)果文件中拆分的染色體部分進(jìn)行合并，并輸出最終結(jié)果。

work_path=工作目錄

#前期拆分基因組時(shí)生成的文件，也可手動(dòng)做成
path/to/contrast.id
## 格式如下：
#拆分后染色體號    舊染色體號   要更改為的新染色體號 原有坐標(biāo)-現(xiàn)有坐標(biāo)    拆分后對應(yīng)染色體長度
Chr01_S1       Chr01  new_01  0       551383505
Chr01_S2       Chr01  new_01  91897250        551383505
Chr01_S3       Chr01  new_01  183794500       551383505
Chr01_S4       Chr01  new_01  275691750       551383505
Chr01_S5       Chr01  new_01  367589000       551383505

#要更改的vcf,壓縮非壓縮均可
input.vcf(.gz)

#輸出vcf文件，gz格式的
out.vcf.gz

python add_vcf.py -p ${work_path} \
                -c path/to/contrast.id \
                -v input.vcf(.gz) -o out.cvf.gz

腳本內(nèi)容如下：

import argparse
import gzip
import re

parser = argparse.ArgumentParser()

parser.add_argument('-p', "--path",dest = "p", default=".", #metavar=", defining metavar is causing an error for some reason
                    help = "work path"
                    )
parser.add_argument('-c', "--contrast_id",dest = "c", default="contrast.id", #metavar=", defining metavar is causing an error for some reason
                    help = "contrast.id"
                    )
parser.add_argument('-v', "--vcf",dest = "v", default="input.vcf.gz", #metavar=", defining metavar is causing an error for some reason
                    help = "old vcf"
                    )
parser.add_argument('-o', "--out_file",dest = "o", default = "input_add.vcf.gz", #metavar=", defining metavar is causing an error for some reason
                    help = "outfile"
                    )

args = parser.parse_args()

contarst = args.p + "/" + args.c
work_dir = args.p
outfile = args.p + "/" + args.o
vcffile = args.p + "/" + args.v

name_dic = {}
with open(contarst,'r')as f1:
    for l in f1:
        name_dic[l.strip().split("\t")[0]] = l.strip().split("\t")[1:]

change_id = []
with gzip.open(outfile , 'wb')as outvcf:
    if 'gz' in vcffile:
        with gzip.open(vcffile , 'rb')as vcf:
            for l in vcf:
                line = l.decode()
                if line[0] == '#':
                    if line[:10] == "##contig=<":
                        Id = re.search("ID=.*,{1}",line).group()[3:-1].split(",")[0]
                        length = re.search("length=.*,{1}",line).group()[7:-1]
                        if name_dic[Id][1] != Id and name_dic[Id][1] not in change_id:
                            outline = "##contig=<ID=" + name_dic[Id][1] +",length=" + name_dic[Id][-1] + ",assembly=unknown>\n"
                            outvcf.write(outline.encode())
                            change_id.append(name_dic[Id][1])
                        elif name_dic[Id][1] == Id and name_dic[Id][1] not in change_id:
                            outvcf.write(l)
                    else:
                        outvcf.write(l)
                else:
                    lst = line.strip().split("\t")
                    Id = lst[0]
                    new_pos = str(int(lst[1]) + int(name_dic[Id][-2]))
                    outline = name_dic[Id][1] + "\t" + new_pos + "\t" + "\t".join(lst[2:]) + "\n"
                    outvcf.write(outline.encode())
    else:
        with open(vcffile , 'r')as vcf:
            for line in vcf:
                if line[0] == '#':
                    if line[:10] == "##contig=<":
                        Id = re.search("ID=.*,{1}",line).group()[3:-1].split(",")[0]
                        length = re.search("length=.*,{1}",line).group()[7:-1]
                        if name_dic[Id][1] != Id and name_dic[Id][1] not in change_id:
                            outline = "##contig=<ID=" + name_dic[Id][1] +",length=" + name_dic[Id][-1] + ",assembly=unknown>\n"
                            outvcf.write(outline.encode())
                            change_id.append(name_dic[Id][1])
                        elif name_dic[Id][1] == Id and name_dic[Id][1] not in change_id:
                            outvcf.write(line.encode())
                    else:
                        outvcf.write(line.encode())
                else:
                    lst = line.strip().split("\t")
                    Id = lst[0]
                    new_pos = str(int(lst[1]) + int(name_dic[Id][-2]))
                    outline = name_dic[Id][1] + "\t" + new_pos + "\t" + "\t".join(lst[2:]) + "\n"
                    outvcf.write(outline.encode())

往期內(nèi)容

大基因組染色體拆分

整理不易，給個(gè)贊再走唄~~