python合并RepeatMasker预测结果中染色体的overlap区域
作者:生信工具箱 发布时间:2021-12-28 21:29:35
标签: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
0
投稿
猜你喜欢
- 淘宝商城的 detail 页面“产品详情”部分是商家自定义区块,曾出现这样一个问题:推荐:css行高:line-height属性详解 <
- 1.概述pyecharts 是百度开源的,适用于数据可视化的工具,配置灵活,展示图表相对美观,顺滑。2.安装python3环境下的安装:pi
- 最近项目中需要Python的打包,看到网上也没有很详细的资料,于是做了一些示例程序。研究了一下,Python如何在Windows和Linux
- 作为收费应用方面的数据库管理员(DBA),公司首席信息官(CIO)经常邀请我与Sarbanes-Oxley审查员开会讨 * 司数据的安全与整合
- 1. 实验说明问题要求:针对静态单赋值(SSA)形式的函数中间代码输入,输出函数返回值的范围实现思路: 基本根据 2013年在CGO会议上提
- 1.定时器Timer定时器源码实现,和自定义一个线程方式一样,都是继承Thread类,重写了run()方法,只是实现的功能是延时执行一个函数
- 人一旦习惯了某些东西就很难去改,以及各种各样的原因,新的浏览器越来越多,而老的总淘汰不了。增长总是快于消亡导致了浏览器兼容是成了谈不完的话题
- 问题描述101/100,想要的结果是1,但是算出来的结果是1.01解决方案101 // 100结果图补充:python2 与 3:一般除法、
- 在客户端,Get方式在通过URL提交数据,数据在URL中可以看到;POST方式,数据放置在HTML HEADER内提交。GET方式提交的数据
- unittest是python的一个单元测试框架关于断言它是用于对一个确定结果和预测结果的一种判断,如果结果正确无任何返回效果,如果结果错误
- python语言本身没有提供const,但实际开发中经常会遇到需要使用const的情形,由于语言本身没有这种支出,因此需要使用一些技巧来实现
- 本文实例讲述了python中global用法。分享给大家供大家参考。具体分析如下:1、global---将变量定义为全局变量。可以通过定义为
- 本文实例讲述了Python使用matplotlib实现在坐标系中画一个矩形的方法。分享给大家供大家参考。具体实现方法如下:import ma
- 实现功能QuestType 1->查询语句, 2->更新语句, 3->删除语句, 4->插入语句<
- 适配器模式Adapter Pattern是什么适配器模式是一种结构型模式,它可以将一个类的接口转换成客户端所期望的接口,从而使原本不兼容的类
- 前言激活函数在机器学习中常常用在神经网络隐含层节点与神经网络的输出层节点上,激活函数的作用是赋予神经网络更多的非线性因素,如果不用激励函数,
- str.join即sequence – 要连接的元素序列。返回通过指定字符连接序列中元素后生成的新字符串。n =
- 原来的题目设想为界面视觉效果的统一性,但是“统一”这个词似乎有点敏感,怕触动萌点无数,而我也无意去设定一个什么什么的统一性来侃侃而谈,极为专
- 如何在ADO中使用SQL函数?代码见下:<%Set conn1 = Server.CreateObjec
- 一、概念介绍Thread 是threading模块中最重要的类之一,可以使用它来创建线程。有两种方式来创建线程:一种是通过继承Thread类