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


猜你喜欢
- 一. 图片懒加载的目的大型网站如常用的淘宝,京东等页面,需要展示大量的商品图片信息,如果打开网页时让所有图片一次性加载完成,需要处理很多次网
- export default和export的区别export主要用于对外输出本模块变量的接口,一个文件就可以被理解为一个模块。export就
- 前端使用uniapp开发项目完成后,需要将页面打包,生成H5的静态文件,部署在服务器上。这样通过服务器链接地址,直接可以在手机上点开来访问。
- 一、数据集爬取现在的深度学习对数据集量的需求越来越大了,也有了许多现成的数据集可供大家查找下载,但是如果你只是想要做一下深度学习的实例以此熟
- 1.安装numpy进入python安装目录的lib下的site-packages文件夹下打开cmd输入pip install numpy下载
- 本文中介绍了一个MySQL的存储过程,其中涉及Cursor的使用,示例如下:CREATE PROCEDURE `justifyGroupNu
- 1在配置文件中添加skip-grant-tables后重启mysql,然后直接登录[root@tyjs09 ~]# cat /etc/my.
- 这个代码不是很完善,能实现基本的功能;另外有个问题,就是divOpenWin层的定位问题:发现如果其属性设置成display:none,那么
- 让我们来看一些例子:--获取表的count信息select count(*) from T with(nolock)--获取特定值的coun
- 指针的操作在Go语言中,指针是一种非常重要的类型,可以用来传递变量的地址而不是变量本身。定义指针在Go语言中,使用*运算符来定义指针。例如,
- 本文实例为大家分享了python实现图书管理系统的具体代码,供大家参考,具体内容如下添加新书查询借阅二次添加新书(读取已有的.xls并修改)
- 本文实例讲述了python执行等待程序直到第二天零点的方法。分享给大家供大家参考。具体分析如下:如果需要通过python每天凌晨定时执行执行
- llama Index是什么《零开始带你入门人工智能系列》第一篇:还用什么chatpdf,让llama Index 帮你训练pdf。Llam
- 首先,运行 Python 解释器,导入 re 模块并编译一个 RE:#!python Python 2.2.2 (#1, Feb 10 20
- 简单演示import matplotlib.pyplot as pltimport numpy as np# 从[-1,1]中等距去50个数
- 写在前面 众所周知python拥有众多的第三方库,据不完全统计python有1w多个第三方库(为什么是不完全统计,因为我也记不清了☺),
- auto-vue-fileauto create .vue file by shell command通过终端自动创建vue文件前言:1:
- Sql Server 中一个非常强大的日期格式化函数: 获得当前系统时间,GETDATE(): 2008年01月08日 星期二 14:59
- 引言语音端点检测最早应用于电话传输和检测系统当中,用于通信信道的时间分配,提高传输线路的利用效率.端点检测属于语音处理系统的前端操作,在语音
- 1、说明tqdm是一个方便且易于扩展的Python进度条,可以在python执行长循环时在命令行界面实时地显示一个进度提示信息,包括执行进度