脚本使用方法及参数介绍
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
