Невозможно воспроизвести R data.table::foverlaps в Python

Я использую data.table::foverlaps в контексте проблемы интервала перекрытия геномики. Недавно я попытался найти эквивалент foverlaps в Python, потому что было бы гораздо лучше использовать только один язык, а не комбинировать Python и R каждый раз, когда мне приходится копаться в результатах анализа. Конечно, я не первый, кто задает вопрос о поиске эквивалента R foverlaps в python, применимого в Python pandas. Это наиболее актуальный пост, который я нашел на SO:

2015 Объединить кадры данных pandas, где одно значение находится между двумя другими

2016 Эквивалент R foverlaps в Python

2017 Как объединить два кадра данных, для которых значения столбца находятся в определенном диапазоне?

2018 Как воспроизвести тот же вывод foverlaps в R со слиянием панд в python?

Дело в том, что я вообще не специалист по 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)

И вот как я попытался воспроизвести это на питоне:

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 будет слишком большой.

Заключить :

  • Мой первый вопрос: где я ошибся в своем коде?
  • Есть ли подход, который был бы более подогнанным и близким к foverlaps? по скорости вычислений?

Спасибо за чтение, я надеюсь, что этот пост подходит для SO.

Анализ настроения постов в Twitter с помощью Python, Tweepy и Flair
Анализ настроения постов в Twitter с помощью Python, Tweepy и Flair
Анализ настроения текстовых сообщений может быть настолько сложным или простым, насколько вы его сделаете. Как и в любом ML-проекте, вы можете выбрать...
7 лайфхаков для начинающих Python-программистов
7 лайфхаков для начинающих Python-программистов
В этой статье мы расскажем о хитростях и советах по Python, которые должны быть известны разработчику Python.
Установка Apache Cassandra на Mac OS
Установка Apache Cassandra на Mac OS
Это краткое руководство по установке Apache Cassandra.
Сертификатная программа "Кванты Python": Бэктестер ансамблевых методов на основе ООП
Сертификатная программа "Кванты Python": Бэктестер ансамблевых методов на основе ООП
В одном из недавних постов я рассказал о том, как я использую навыки количественных исследований, которые я совершенствую в рамках программы TPQ...
Создание персонального файлового хранилища
Создание персонального файлового хранилища
Вы когда-нибудь хотели поделиться с кем-то файлом, но он содержал конфиденциальную информацию? Многие думают, что электронная почта безопасна, но это...
Создание приборной панели для анализа данных на GCP - часть I
Создание приборной панели для анализа данных на GCP - часть I
Недавно я столкнулся с интересной бизнес-задачей - визуализацией сбоев в цепочке поставок лекарств, которую могут просматривать врачи и...
2
0
203
2
Перейти к ответу Данный вопрос помечен как решенный

Ответы 2

По сути, логика слияния и интервала отличается между R и Python.

р

Согласно документам foverlaps, вы используете тип Любые по умолчанию, который выполняется в следующих условиях:

Let [a,b] and [c,d] be intervals in x and y with a<=b and c<=d.
...
For type="any", as long as c<=b and d>=a, they overlap.

Кроме того, вы присоединяетесь к другим столбцам ключей. В целом вы навязываете следующую логику (переведенную на столбцы 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

питон

Вы выполняете интервальный запрос 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

Заметные отличия Python от R: 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 и создавать индексы в полях соединения: цвет, прядь, Начало и остановка.

Приятно! Это помогло с некоторыми изменениями. Это делает вещи более ясными. Я чувствую, что должен копать больше этих проблем слияния.

Paul Endymion 23.05.2019 19:32

Рад слышать! Да, я думал, что потребуются некоторые корректировки. Удачного кодирования!

Parfait 23.05.2019 19:45

Индексы SQL не похожи на setkey, в data.table есть индексы (setindex). setkey - это конкретный clustered index, не уверен, что sqlite имеет кластеризованный индекс.

jangorecki 24.05.2019 16:01

Спасибо @jangorecki. Отредактировал строку подобия. Обратите внимание: некоторые СУБД поддерживают кластерные индексы, но зависят от поставщика. OP использует SQLite в качестве поставщика, а его WITHOUT ROWID использует индекс кластера в качестве первичного ключа.

Parfait 24.05.2019 16:12
Ответ принят как подходящий

Чтобы сделать перекрытие интервалов в Python, просто используйте пиранжи:

import pyranges as pr

c1 = """Chromosome Start End Gene
1 10 20 blo
1 45 46 bla"""

c2 = """Chromosome Start End Gene
1 10 35 bip
1 25 50 P53
1 40 10000 boop"""


gr1, gr2 = pr.from_string(c1), pr.from_string(c2)

j = gr1.join(gr2)
# +--------------+-----------+-----------+------------+-----------+-----------+------------+
# |   Chromosome |     Start |       End | Gene       |   Start_b |     End_b | Gene_b     |
# |   (category) |   (int32) |   (int32) | (object)   |   (int32) |   (int32) | (object)   |
# |--------------+-----------+-----------+------------+-----------+-----------+------------|
# |            1 |        10 |        20 | blo        |        10 |        35 | bip        |
# |            1 |        45 |        46 | bla        |        25 |        50 | P53        |
# |            1 |        45 |        46 | bla        |        40 |     10000 | boop       |
# +--------------+-----------+-----------+------------+-----------+-----------+------------+
# Unstranded PyRanges object has 3 rows and 7 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

Спасибо за ваш ответ, это выглядит здорово! Вы написали эту библиотеку? Я вижу эталон внизу страницы, случайно ли вы не делали его с foverlaps data.table?

Paul Endymion 24.04.2020 14:17

Я не знал, я даже не знал о существовании data.tables :) Но pyranges легко справляется с большими наборами данных. Да, я написал :)

The Unfun Cat 24.04.2020 15:45

Другие вопросы по теме