python合并RepeatMasker预测结果中染色体的overlap区域

前言

RepeatMasker是一个通过已有数据库预测重复序列的软件，可以筛选DNA序列中的散在重复序列和低复杂序列，是重复序列注释的重要软件。

问题

我们想对RepeatMasker预测的结果文件进行重复序列的合并，也就是去除染色体之间的overlap区域同时将基因间距小于50个bp的也同样视为overlap，我们应该如何用python处理并生成新的预测结果？

思路

• 首先需要对文件进行预处理提取出需要处理的列，'//'可以忽略

• 对相同染色体序列按照升序进行归并排序

• 分别取相应染色体按照滑动窗口的思路进行双指针比对，注意gap=50

1. 预处理

我们这里只需要结果文件的前三列，可以使用awk命令获取

awk '{for(i = 1; i <= 3; i++)
        printf("％s ", $i);
        printf("\n")}' result.txt >  pretreatment.txt  
        #result.txt为结果文件，pretreatment.txt为预处理结果文件

2. 将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)))

3.去重+归并排序

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]))
#按照第一列，第二列，第三列分别排降升序

4.开始比对，gap=50

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

5.将new_result.txt作为输出文件，生成结果

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)

6. 完整代码

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)

来源：https://www.jianshu.com/p/42faa3a0228e