awk 查找重叠

2024-02-21

我有一个包含列的文件,如下所示。

Group   Start        End
chr1    117132092    118875009
chr1    117027758    119458215
chr1    103756473    104864582
chr1    105093795    106219211
chr1    103354114    104747251
chr1    102741437    105235140
chr1    100090254    101094139
chr1    100426977    101614730
chr2    86644663     87767193
chr2    82473711     83636545
chr2    83896702     85079032
chr2    83876122     85091910
chr2    82943211     84350917
chr3    89410051     90485635
chr3    89405753     90485635
chr3    86491492     87593215
chr3    82507157     83738004
chr3    85059618     86362254

我想找到每个组中这些坐标之间的重叠(按 chr1,chr2,chr3 ..分组)。

必须检查起始坐标和结束坐标是否与同一组中的其他坐标至少有 50% 的重叠。如果至少有 50% 重叠,则必须在第 3 列和第 4 列中报告新的开始和结束坐标(这是重叠区域的范围)。如果它们不重叠,则必须在第 3 列和第 4 列中报告原始开始和结束位置。

为了更清楚地说明,我们取前两行

                 117132092..........118875009
         117027758...........................119458215

由于它们彼此至少重叠 50%,因此重叠范围在输出中报告为新的开始和新的结束。并且第 3 行和第 4 行不与其他行重叠,因此原始坐标在第 3 列和第 4 列中被报告为新的起点和新的结束。同样,由于第 5 行和第 6 行彼此有 50% 的重叠,因此它们的范围被报告为新的第 3 列和第 4 列中的开始和新结束。 这是预期的输出:

Group   Start     End         NewStart   NewEnd   
chr1 117132092 118875009  117027758   119458215
chr1 117027758 119458215  117027758   119458215
chr1 103756473 104864582  103354114   104864582
chr1 105093795 106219211  105093795   106219211
chr1 103354114 104747251  102741437   105235140
chr1 102741437 105235140  102741437   105235140
chr1 100090254 101094139  100090254   101614730
chr1 100426977 101614730  100090254   101614730
chr2 86644663 87767193    86644663    87767193
chr2 82473711 83636545    82473711    83636545 
chr2 83896702 85079032    83876122    85091910
chr2 83876122 85091910    83876122    85091910
chr2 82943211 84350917    82943211    84350917
chr3 89410051 90485635    89405753    90485635
chr3 89405753 90485635    89405753    90485635
chr3 86491492 87593215    86491492    87593215
chr3 82507157 83738004    82507157    83738004
chr3 85059618 86362254    85059618    86362254

我已经用 R 编程语言实现了这一点,但是原始文件太大并且需要很长时间才能运行。有人可以帮助在 awk 中做到这一点吗?


使用 Gnu Awk 版本 4,您可以尝试:

gawk -f a.awk file file

where a.awk is:

NR==FNR {
    if (FNR>1) {
        a[$1][++i]=$2
        b[$1][i]=$3
    }
    next
}
FNR==1 {
    fmt="%-7s%-10s%-10s%-10s%-10s\n"
    printf fmt,"Group","Start","End","NewStart","NewEnd" 
}
FNR>1{
    $4=$2; $5=$3
    n=checkInside($1,$2,$3)
    if (n>0) {
        ff=0; x=$2; y=$3
        for (i=1; i<=n; i++) {
            ar=a[$1][R[i]]; br=b[$1][R[i]];
            getIntersect($2,$3,ar,br)
            getLargest($2,$3,ar,br)
            ovl=((i2-i1)/($3-$2))*100;
            ovr=((i2-i1)/(br-ar))*100;
            if (ovl>50 && ovr>50) {
                if (r1<x) x=r1
                if (r2>y) y=r2
                ff=1
            }
        }
        if (ff) {
            $4=x; $5=y
        }
    }
    printf fmt,$1,$2,$3,$4,$5
}

function getLargest(x1,y1,x2,y2) {
    r1=(x1<=x2)?x1:x2
    r2=(y1>=y2)?y1:y2
}

function getIntersect(x1,y1,x2,y2) {
    if (x1>=x2 && x1<=y2) {
        i1=x1;
    } else {
        i1=x2;
    }
    i2=(y1<=y2)?y1:y2
}

function checkInside(g,x,y,i,j,x1,y1) {
    R["x"]=0
    for (i in a[g]) {
        x1=a[g][i]; y1=b[g][i];
        if ((x>=x1 && x<=y1) || (y>=x1 && y<=y1)) {
            if (!(x==x1 && y==y1))
                R[++j]=i
        }
    }
    return j
}

Output:

Group  Start     End       NewStart  NewEnd    
chr1   117132092 118875009 117027758 119458215 
chr1   117027758 119458215 117027758 119458215 
chr1   103756473 104864582 103354114 104864582 
chr1   105093795 106219211 105093795 106219211 
chr1   103354114 104747251 102741437 105235140 
chr1   102741437 105235140 102741437 105235140 
chr1   100090254 101094139 100090254 101614730 
chr1   100426977 101614730 100090254 101614730 
chr2   86644663  87767193  86644663  87767193  
chr2   82473711  83636545  82473711  83636545  
chr2   83896702  85079032  83876122  85091910  
chr2   83876122  85091910  83876122  85091910  
chr2   82943211  84350917  82943211  84350917  
chr3   89410051  90485635  89405753  90485635  
chr3   89405753  90485635  89405753  90485635  
chr3   86491492  87593215  86491492  87593215  
chr3   82507157  83738004  82507157  83738004  
chr3   85059618  86362254  85059618  86362254  
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

awk 查找重叠 的相关文章

随机推荐

  • 实体框架 ObjectContext 分享 - 优缺点

    在我的项目中 我使用实体框架 4 0 作为 ORM 将数据保存在 SQL Server 中 我的项目是应用程序的功能区 主窗体中有网格视图和导航树 其顶部有功能区面板 我的应用程序基本上是一个 CRUD UI 几乎没有业务逻辑 第一次使用
  • 位域元素的默认值

    在 C 11 中可以做 struct S int i 42 如果忘记初始化成员i它 默认初始化为 42 我刚刚尝试过 位域为 struct S int i 42 5 我正在得到 错误 预期为 在 标记之前 位域成员是否存在此功能 如果存在
  • JNI 的用处[关闭]

    就目前情况而言 这个问题不太适合我们的问答形式 我们希望答案得到事实 参考资料或专业知识的支持 但这个问题可能会引发辩论 争论 民意调查或扩展讨论 如果您觉得这个问题可以改进并可能重新开放 访问帮助中心 help reopen questi
  • 获取文本字段中最常用的 10 个单词

    我有一个包含数千个文档的索引 每个文档都有一个全文字段 我想搜索所有这些字段并获取最常出现的 10 个最常见的单词 如果可能的话 我还想要一种在 Kibana 上可视化它的方法 实现此目的的最常见方法是使用keyword datatype
  • Drupal CCK 的复选框

    我是 Drupal 的新人 到目前为止很喜欢 我正在创建 CCK 自定义内容类型 我需要以复选框格式制作便利设施列表 所以我做了 文件类型 Text 小部件类型 复选框 单选按钮 和允许值列表 onsite dining 现场用餐 Meet
  • 如何在启动时运行命令?

    我试图弄清楚如何在启动时运行命令 就像我将其输入控制台一样 我在 Raspberry Pi 上使用 Rasbian 但我认为这个问题对于 Debian 来说通常是相同的 我尝试运行的命令是 sudo screen mono server e
  • Rails 6 无法连接到 AWS Elastic Beanstalk 预置的 RDS。 Unix 域套接字“/var/run/postgresql/.s.PGSQL.5432”

    我在尝试向 Elastic Beanstalk 启动示例 Rails 6 应用程序时遇到了非常困难的情况 对于上下文 我遵循这些说明 将 RDS 添加到 Ruby 应用程序 https docs aws amazon com elastic
  • Python 中的线程队列挂起

    我正在尝试通过队列使解析器成为多线程 它似乎有效 但我的队列挂起 如果有人能告诉我如何解决这个问题 我将不胜感激 因为我很少编写多线程代码 此代码从 Q 中读取 from silk import import json import dat
  • 当父视图的任何子视图重绘时,如何触发父视图的重绘?

    背景 我写了一个基于 Android 的自定义视图LinearLayout我称之为ReflectingLayout 这个想法相当简单 在声明的任何子视图下面渲染反射效果ReflectingLayout e g 我通过覆盖来实现这一目标dis
  • C# NHibernate 简单问题

    我在用着NHibernate驱动存储库 Fluent映射并尝试使用Linq to NHibernate 但对于像这样的一些简单查询 Retrieve
  • 解析 HTML 内容时阻止 etree 解析 HTML 实体

    有没有办法阻止etree在解析HTML内容时解析HTML实体 html etree HTML amp html find body text 这给了我 但我想得到 本身 您始终可以对数据进行预处理 后处理 在输入 HTML 解析器之前将 替
  • 如何设置 10 点到 19 点每两小时执行一次 Cron 作业

    1 个月前我对此有一个疑问 那是1个小时的间隔 我得到了确切的答案 以下是旧问题的链接 如何设置从上午 9 00 到下午 6 00 周一至周五 每隔一小时执行一次 Cron 作业 https stackoverflow com questi
  • Python - Tkinter 文本大小未调整大小

    我正在尝试使用 Tkinter 制作一个可以调整大小的窗口 并且效果很好 但我希望字体大小也能按比例缩放 输入框的大小完美调整 但文本大小保持不变 我也可以更改输入文字大小吗 我该怎么做呢 先感谢您 这是到目前为止我的代码片段 defini
  • 压缩 XML 指标。

    我有一个客户端服务器应用程序 它通过 TCP IP 将 XML 从客户端发送到服务器 然后广播到其他客户端 我如何知道 XML 的最小大小可以通过压缩 XML 而不是通过常规流发送来保证性能改进 有什么好的指标或例子吗 Xml 通常压缩得很
  • 为什么有各种 JPEG 扩展名?

    在开发下载器时 我遇到了以下使用 Python 的情况mimetypes guess extension功能 In 2 mimetypes guess extension image jpeg strict False Out 2 jpe
  • 当 wifi 打开时,仅通过 Android 手机上的移动数据发送数据

    即使 wifi 已打开并连接到互联网 是否也可以通过移动数据以编程方式路由请求 我的应用程序需要调用运营商提供的服务 该服务只能通过移动数据使用 而且我认为要求关闭 wifi 不太方便 看一眼https developer android
  • 多对多关系的3表之间的SQL查询

    我有三张桌子 friends locations friend location friend location是一个连接表 允许多对多关系friends and locations 所以表格看起来像这样 Friends ID Name 1
  • WooCommerce:管理员手动创建订单时需要挂钩

    我的网站之一使用 WooCommerce 客户有时希望在订单管理中手动创建订单 WooCommerce gt 订单 gt 添加订单 当他们单击该页面上的 保存订单 时 我需要对订单进行一些额外的处理 有可用的钩子吗 我浏览了 WooComm
  • Microsoft 聊天机器人 (Node.js) 是否在单个 LUIS.AI 应用程序中支持多种语言?

    我有一个使用 Node js 在 Microsoft 机器人框架中构建的聊天机器人 并将该机器人与名为 LUIS AI 智能的 NLP 框架集成 以根据用户的意图和实体处理用户对话 在这里 我需要这个机器人在单个 LUIS 应用程序中支持多
  • awk 查找重叠

    我有一个包含列的文件 如下所示 Group Start End chr1 117132092 118875009 chr1 117027758 119458215 chr1 103756473 104864582 chr1 10509379