时间戳

首页 > 解决方案 > 无法在 Python 中重现 R data.table::foverlaps

python - 无法在 Python 中重现 R data.table::foverlaps

问题描述

我在重叠基因组区间问题data.table::foverlaps的背景下使用。我最近开始尝试在 Python 中找到等效的 foverlaps，因为每次我必须挖掘分析输出时，只使用一种语言而不是结合 Python 和 R 会更好。当然，我不是第一个提出在 Python 中找到适用于 Python 熊猫的 R foverlaps 的等价物的问题。这些是我在 SO 上找到的最相关的帖子：

2015合并熊猫数据框，其中一个值位于其他两个值之间

2016 R foverlaps 等效于 Python

2017如何连接列值在一定范围内的两个数据框？

2018如何通过 python 中的 pandas 合并重现 R 中 foverlaps 的相同输出？

问题是我根本不是 Python 专家。所以我选择了对我来说最相关/最容易理解的答案，那个sqlite3。

这就是我在 R 中的做法：

library(data.table)

interv1 <- cbind(seq(from = 3, to = 40, by = 4),seq(from = 5, to = 50, by = 5), c(rep("blue",5), rep("red", 5)), rep("+",10))
interv2 <- cbind(seq(from = 3, to = 40, by = 4),seq(from = 5, to = 50, by = 5), c(rep("blue",5), rep("red", 5)), rep("-",10))
interv  <- rbind(interv1, interv2)
interv <- data.table(interv)
colnames(interv) <- c('start', 'stop', 'color', 'strand')
interv$start <- as.integer(interv$start)
interv$stop <- as.integer(interv$stop)
interv$stop <- interv$stop -1
interv$cov <- runif(n=nrow(interv), min = 10, max = 200)

to_match <- data.table(cbind(rep(seq(from = 4, to = 43, by = 4),2), rep(c(rep("blue", 5), rep("red", 5)), 2), c(rep("-", 10), rep("+", 10))))
colnames(to_match) <- c('start', 'color', 'strand')
to_match$stop <-  to_match$start 
to_match$start <- as.integer(to_match$start)
to_match$stop <- as.integer(to_match$stop)

setkey(interv, color, strand, start, stop)
setkey(to_match, color, strand, start, stop)

overlapping_df <- foverlaps(to_match,interv)

#write.csv(x = interv, file = "Documents/script/SO/wig_foverlaps_test.txt", row.names = F)
#write.csv(x = to_match, file = "Documents/script/SO/cov_foverlaps_test.txt", row.names = F)

这就是我试图在 python 中重现它的方式：

import pandas as pd
import sqlite3

cov_table = pd.DataFrame(pd.read_csv('SO/cov_foverlaps_test.txt', skiprows = [0], header=None))
cov_table.columns = ['start', 'stop', 'chrm', 'strand', 'cov']
cov_table.stop = cov_table.stop - 1


wig_file = pd.DataFrame(pd.read_csv('SO/wig_foverlaps_test.txt', header=None, skiprows = [0]))
wig_file.columns = ['i_start', 'chrm', 'i_strand', 'i_stop']

cov_cols = ['start','stop','chrm','strand','cov']
fract_cols = ['i_start','i_stop','chrm','i_strand']

cov_table = cov_table.reindex(columns = cov_cols)
wig_file = wig_file.reindex(columns = fract_cols)

cov_table.start = pd.to_numeric(cov_table['start'])
cov_table.stop = pd.to_numeric(cov_table['stop'])

wig_file.i_start = pd.to_numeric(wig_file['i_start'])
wig_file.i_stop = pd.to_numeric(wig_file['i_stop'])



conn = sqlite3.connect(':memory:')

cov_table.to_sql('cov_table', conn, index=False)
wig_file.to_sql('wig', conn, index=False)

qry = '''
    select  
        start PresTermStart,
        stop PresTermEnd,
        cov RightCov,
        i_start pos,
        strand Strand
    from
        cov_table join wig on
        i_start between start and stop and 
        cov_table.strand = wig.i_strand
     '''

test = pd.read_sql_query(qry, conn)

无论我如何更改代码，我总是在输出（测试）中发现一些小的差异，在这个例子中，我在 python 结果表中缺少两行，其中的值应该在范围内并且是相等的到范围的末尾：

缺线：

> 19   24  141.306318     24      +
> 
> 19   24  122.923700     24      -

最后，我担心如果我找到正确的方法来做它sqlite3，计算时间差异data.table::foverlaps会太大。

总结：

感谢您的阅读，我希望这篇文章适合 SO。

标签： pythonrpandasdata.table

解决方案

本质上，R 和 Python 之间的合并逻辑和区间逻辑不同。

R

根据foverlaps文档，您使用的是默认值，任何在以下条件下运行的类型：

令 [a,b] 和 [c,d] 为 x 和 y 中的区间，其中 a<=b 且 c<=d。
...
对于 type="any"，只要 c<=b 和 d>=a，它们就会重叠。

此外，您加入其他列的键。总而言之，您正在强加以下逻辑（转换为 SQLite 列以进行比较）：

foverlaps(to_match, interv) --> foverlaps(cov_table, wig)

  1. wig.i_start <= cov_table.stop (i.e., c <= b)
  2. wig.i_stop >= cov_table.start (i.e., d >= a)
  3. wig.color == cov_table.color
  4. wig.strand == cov_table.strand

Python

您正在运行一个INNER JOIN强制执行以下逻辑的 + 间隔查询：

  1. wig.i_start >= cov_table.start (i.e., i_start between start and stop)
  2. wig.i_start <= cov_table.stop (i.e., i_start between start and stop)
  3. wig.strand == cov_table.strand

与 R 相比，Python 的显着差异：wig.i_stop从未使用过；wig.i_chrm（或颜色）从未使用过；并且wig.i_start被调节了两次。

要解决此问题，请考虑以下未经测试的 SQL 调整，以希望实现 R 结果。顺便说一句，SQL 中的最佳实践是为JOIN子句中的所有列（甚至SELECT）别名：

select  
   cov_table.start as PresTermStart,
   cov_table.stop as PresTermEnd,
   cov_table.cov as RightCov,
   wig.i_start as pos,
   wig.strand as Strand
from
   cov_table 
join wig 
    on cov_table.color = wig.i_chrm
   and cov_table.strand = wig.i_strand
   and wig.i_start <= cov_table.stop 
   and wig.i_stop  >= cov_table.start

为了获得更好的性能，请考虑使用持久（非内存）SQLite 数据库并在连接字段上创建索引： color、strand、start和stop。

推荐阅读