RepeatMasker是一个通过已有数据库预测重复序列的软件,可以筛选DNA序列中的散在重复序列和低复杂序列,是重复序列注释的重要软件。
我们想对RepeatMasker预测的结果文件进行重复序列的合并,也就是去除染色体之间的overlap区域同时将基因间距小于50个bp的也同样视为overlap,我们应该如何用python处理并生成新的预测结果?
我们这里只需要结果文件的前三列,可以使用awk命令获取
awk '{for(i = 1; i <= 3; i++) printf("%s ", $i); printf("\n")}' result.txt > pretreatment.txt #result.txt为结果文件,pretreatment.txt为预处理结果文件
with open ('pretreatment.txt','r')as f: for i in f.readlines(): if i.strip() == '//': continue c = i.strip().split('\t') b.append(c[0]) a.append((c[0],int(c[1]),int(c[2]))) print ("全部染色体数量: "+str(len(a)))
c = [i for i in b_set if b.count(i) == 1] for i in a: if i[0] not in c: continue a.remove(i) result.append((i[0],int(i[1]),int(i[2]))) print ("去重后染色体数量: "+str(len(a)))
a.sort(key = lambda x : (x[0], x[1], x[2])) #按照第一列,第二列,第三列分别排降升序
q = '' start = 0 end = 0 tem1 = [] tem2 = [] gap = 50 for i in a: if i[0] != q: if tem1: if tem1 not in tem2: tem2.append(tem1) tem1 = [] q = I[0] start = int(i[1]) end = int(i[2]) continue if int(i[1]) < end or int(i[1]) - end < gap: if int(i[2]) > end: end = int(i[2]) continue else: continue tem1.append([q,start,end]) start = int(i[1]) end = int(i[2])
with open ('new_result.txt','w')as f: for i in tem2: for o in I: print (o[0],o[1],o[2],file=f) for i in result: print (i[0],i[1],i[2],file=f)
a = [] b = [] with open ('pretreatment.txt','r')as f: for i in f.readlines(): if i.strip() == '//': continue c = i.strip().split('\t') b.append(c[0]) a.append((c[0],int(c[1]),int(c[2]))) print ("全部染色体数量: "+str(len(a))) b_set = set(b) result = [] c = [i for i in b_set if b.count(i) == 1] for i in a: if i[0] not in c: continue a.remove(i) result.append((i[0],int(i[1]),int(i[2]))) print ("去重后染色体数量: "+str(len(a))) a.sort(key = lambda x : (x[0], x[1], x[2])) q = '' start = 0 end = 0 tem1 = [] tem2 = [] gap = 50 for i in a: if i[0] != q: if tem1: if tem1 not in tem2: tem2.append(tem1) tem1 = [] q = I[0] start = int(i[1]) end = int(i[2]) continue if int(i[1]) < end or int(i[1]) - end < gap: if int(i[2]) > end: end = int(i[2]) continue else: continue tem1.append([q,start,end]) start = int(i[1]) end = int(i[2]) with open ('new_result.txt','w')as f: for i in tem2: for o in I: print (o[0],o[1],o[2],file=f) for i in result: print (i[0],i[1],i[2],file=f)