脚本使用方法及参数介绍

ref=要拆分基因组
out_path=结果输出路径
rename=染色体重命名对应表：old_id\tnew_id,不做更改两列一致即可。比如把CHR01改为Chr01,文件为CHR01\tChr01
split_len=拆分染色体长度的最小值,只拆分长度大于该指定长度的染色体，默认是500000000
split_num=拆分份数，将长度大于指定长度的染色体拆分成指定份数，默认是5份，即拆分结果为，chr_S1,chr_S2,chr_S3,chr_S4,chr_S5,chr_S6
-c -o 参数对应都是输出文件

#命令行
python split_genome.py -r ${ref} -p ${out_path} \
                       -new_name ${rename} \
                       -l ${split_len} \
                       -n ${split_num} \
                       -c ${out_path}/contrast.id -o genome/genome.split.fa >${out_path}/log

##参数
usage: split_genome.py [-h] -r R [-p P] [-new_name N] [-c C] [-n N] [-l L] [-o O]

optional arguments:
  -h, --help            show this help message and exit
  -r R, --ref R         full path of reference
  -p P, --path P        work path
  -new_name N, --new_name N
                        information of new_name:old_name new_name
  -c C, --contrast_id C
                        information of new_name:old_name new_name
  -n N, --split_num N   split num
  -l L, --split_len L   max length of chr
  -o O, --out_file O    outfile

脚本内容如下

import os
import argparse

parser = argparse.ArgumentParser()


parser.add_argument('-r', "--ref",dest = "r", required= True, 
                    help = "full path of reference"
                    )
parser.add_argument('-p', "--path",dest = "p", default=".",
                    help = "work path"
                    )
parser.add_argument('-new_name', "--new_name",dest = "N", default="rename.txt", 
                    help = "information of new_name:old_name\tnew_name"
                    )
parser.add_argument('-c', "--contrast_id",dest = "c", default="contrast.id", 
                    help = "information of new_name:old_name\tnew_name"
                    )
parser.add_argument('-n', "--split_num",dest = "n",type=int, default = 5, 
                    help = "split num"
                    )
parser.add_argument('-l', "--split_len",dest = "l", type=int,default = 500000000, 
                    help = "max length of chr"
                    )
parser.add_argument('-o', "--out_file",dest = "o", default = "genome/genome.split.fa", 
                    help = "outfile"
                    )

args = parser.parse_args()


ref = args.r
work_dir = args.p
outfile = args.o
name_change = args.N
split_len = args.l
split_num = args.n
con =  args.c

#将要拆分基因组软链接到工作目录下，并用samtools建立索引
if "/" in outfile:
    os.system("mkdir -p " + work_dir + "/" + outfile.split('/')[0])
    os.system("ln -sf " + ref + " " +  work_dir + "/" + outfile.split('/')[0] + "/")
    os.system("samtools faidx " + work_dir + "/" + outfile.split('/')[0] + "/" + ref.split("/")[-1] )
    fai = work_dir + "/" + outfile.split('/')[0] + "/" + ref.split("/")[-1] + ".fai"
else:
    os.system("ln -sf " + ref + " " +  work_dir + "/")
    os.system("samtools faidx " + work_dir + "/" + ref.split("/")[-1])
    fai = work_dir + "/" + ref.split("/")[-1] + ".fai"

name_dic = {}
with open(name_change,'r')as f1:
    for l in f1:
        if len(l.split("\t")) >=2:
            name_dic[l.strip().split("\t")[0]] = l.strip().split("\t")[1]
        else:
            name_dic[l.strip().split("\t")[0]] = l.strip().split("\t")[0]

def pick_seq(chr,start,end,ref,new_chr,length):
    os.system("echo -e \"" + chr + "\\t" + str(start) + "\\t" + str(end) + "\" >> test.bed")
    os.system("echo -e \"" + chr + "\\t" + str(start) + "\\t" + str(end) + "\" > mid.bed")
    seq_lst = os.popen("/work1/Software/bedtools/bin/bedtools getfasta -fi "+ ref + " -bed mid.bed").readlines()
    seq_lst[0] = ">" + new_chr + " " + str(length) + " " + chr + "\n"
    con_lst = [new_chr,chr,name_dic[chr],str(length)]
    return seq_lst,con_lst

fai_dic = {}
contrast_dic = {}
with open(fai,'r')as f,open(outfile,'w')as outf,open(con,'w')as log:
    for line in f:
        lst = line.strip().split("\t")
        Id = lst[0]
        n = 0
        if int(lst[1]) >= int(split_len):
            print(lst[0])
            for n in range(split_num):
                add_len = int(int(lst[1]) / int(split_num))
                n = n + 1
                new_id = name_dic[Id] + "_S"+ str(n)
                if n == 1:
                    contrast_dic[new_id] = 0
                    seq,contrast_dic[new_id]=pick_seq(Id,0,add_len,ref,new_id,0)
                elif n < split_num:
                    seq,contrast_dic[new_id]=pick_seq(Id,add_len * (n-1),add_len * n,ref,new_id,add_len * (n-1))
                else:
                    seq,contrast_dic[new_id]=pick_seq(Id,add_len * (n-1),int(lst[1]),ref,new_id,add_len * (n-1))
                log.write("\t".join(contrast_dic[new_id])+ "\t" + lst[1] +"\n")
                outf.write("".join(seq))
        else:
            new_id = name_dic[Id]
            seq,contrast_dic[new_id]=pick_seq(Id,0,int(lst[1]),ref,new_id,0)
            log.write("\t".join(contrast_dic[new_id]) + "\t" + lst[1] +"\n")
            outf.write("".join(seq))

os.system("rm test.bed mid.bed")

拆分结果

#输出拆分结果默认为${out_path}/genome/genome.split.fa
格式如下：
>Chr01_S1 0 Chr01
AAAAAACCCTAAAAACCCTAA
>Chr01_S2 91897250 Chr01
TGCTAAAGAGGGGCTGATTGGGTCCTGGCCCACC
>Chr01_S3 183794500 Chr01
ACATAGCTGAGATGGCCAACAACAATCAAGG
格式解读：
奇数行（如第一行）
第一列：拆分后染色体ID
第二列：原有坐标-现有坐标(后续也是根据这个数字将位点坐标恢复，比如Chr01_S2的第一个碱基对应的物理位置为原来Chr01的91897250 + 1)
第三列：对应旧的染色体号
偶数行（如第二行）是序列

#同时也会输出一个id对应文件,后续也可以根据这个文件进行合并
#拆分后染色体号    旧染色体号   要更改为的新染色体号 原有坐标-现有坐标    拆分后对应染色体长度
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