腳本使用方法及參數(shù)介紹

ref=要拆分基因組
out_path=結(jié)果輸出路徑
rename=染色體重命名對應表：old_id\tnew_id,不做更改兩列一致即可。比如把CHR01改為Chr01,文件為CHR01\tChr01
split_len=拆分染色體長度的最小值,只拆分長度大于該指定長度的染色體，默認是500000000
split_num=拆分份數(shù)，將長度大于指定長度的染色體拆分成指定份數(shù)，默認是5份，即拆分結(jié)果為，chr_S1,chr_S2,chr_S3,chr_S4,chr_S5,chr_S6
-c -o 參數(shù)對應都是輸出文件

#命令行
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

##參數(shù)
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

腳本內(nèi)容如下

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")

拆分結(jié)果

#輸出拆分結(jié)果默認為${out_path}/genome/genome.split.fa
格式如下：
>Chr01_S1 0 Chr01
AAAAAACCCTAAAAACCCTAA
>Chr01_S2 91897250 Chr01
TGCTAAAGAGGGGCTGATTGGGTCCTGGCCCACC
>Chr01_S3 183794500 Chr01
ACATAGCTGAGATGGCCAACAACAATCAAGG
格式解讀：
奇數(shù)行（如第一行）
第一列：拆分后染色體ID
第二列：原有坐標-現(xiàn)有坐標(后續(xù)也是根據(jù)這個數(shù)字將位點坐標恢復，比如Chr01_S2的第一個堿基對應的物理位置為原來Chr01的91897250 + 1)
第三列：對應舊的染色體號
偶數(shù)行（如第二行）是序列

#同時也會輸出一個id對應文件,后續(xù)也可以根據(jù)這個文件進行合并
#拆分后染色體號    舊染色體號   要更改為的新染色體號 原有坐標-現(xiàn)有坐標    拆分后對應染色體長度
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