如何比较值范围并合并两个pandas数据框（或转移值）

Question

在以下数据中：

data01 =

contig  start    end    haplotype_block 
2   5207    5867    1856
2   155667    155670    2816
2   67910    68022  2
2   68464    68483  3
2   525    775  132
2   118938    119559    1157

data02 =

contig    start   last    feature gene_id gene_name   transcript_id
2   5262    5496    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   5579    5750    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   5856    6032    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   6115    6198    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   916 1201    exon    scaffold_200001.1   NA  scaffold_200001.1
2   614 789 exon    scaffold_200001.1   NA  scaffold_200001.1
2   171 435 exon    scaffold_200001.1   NA  scaffold_200001.1
2   2677    2806    exon    scaffold_200002.1   NA  scaffold_200002.1
2   2899    3125    exon    scaffold_200002.1   NA  scaffold_200002.1

问题：

• 我想比较这两个数据框的范围（起始-结束）。
• 如果范围重叠，我想要将data02中的gene_id和gene_name值转移到data01的新列中。

我尝试了以下方法（使用pandas）：

data01['gene_id'] = ""
data01['gene_name'] = ""

data01['gene_id'] = data01['gene_id'].\
apply(lambda x: data02['gene_id']\
        if range(data01['start'], data01['end'])\
           <= range(data02['start'], data02['last']) else 'NA')

如何改进这段代码？我目前使用的是pandas，但如果使用字典能更好地解决问题，我也愿意尝试。但请解释一下过程，我希望学习而不仅仅是得到答案。

谢谢，

期望输出：

contig  start    end    haplotype_block    gene_id    gene_name
2   5207    5867    1856    scaffold_200003.1,scaffold_200003.1,scaffold_200003.1   CP5,CP5,CP5

# the gene_id and gene_name are repeated 3 times because three intervals (i.e 5262-5496, 5579-5750, 5856-6032) from data02 overlap(or touch) the interval ranges from data01 (5207-5867)

# So, whenever there is overlap of the ranges between two dataframe, copy the gene_id and gene_name.

# and simply NA on gene_id and gene_name for non overlapping ranges

2   155667    155670    2816    NA    NA
2   67910    68022  2    NA    NA
2   68464    68483  3    NA    NA
2   525    775  132    scaffold_200001.1   NA
2   118938    119559    1157    NA    NA

- everestial007

---

我对目标感到困惑。如果您能详细说明输出的样子，我很乐意再看一遍。 - piRSquared

@piRSquared：刚刚新增了所需输出。希望能够理解。解决问题的更好方法是通过将数据按照“conting”和“start”进行排序（以防止大规模的for循环），然后先对两个数据框进行处理，我已经完成了这一步；但是上面的data01和data02还没有这样做。 - everestial007

Answer 1

我知道你在使用Python，但是你的问题可以很容易地使用经典的生物信息学工具“bedtools intersect”解决：http://bedtools.readthedocs.io/en/latest/content/tools/intersect.html 你的两个输入文件都遵循标准的BED格式：http://bedtools.readthedocs.io/en/latest/content/general-usage.html Bedtools intersect 提供了高级逻辑来确定两个区域之间的交集或重叠。我认为它也可以直接操作bgzipped输入。

- Gordon Bean

---

我刚刚重新审视了bedtools。对于我的目的，我最好使用Python和Pandas（或者字典）。原因是我正在做的事情是大型管道的一部分，引入bedtools只会使事情变得更加复杂而不是简单。你可能有其他想法/意见。如果有的话，我将不胜感激。谢谢。 - everestial007

Bedtools是为了解决你的问题而开发的 - 这是生物信息学中常见的问题。算法本身可能有些微妙，特别是如果你希望它高效。我建议在重新实现算法之前先尝试一下。 ;) - Gordon Bean

对我有帮助的唯一方法是 bedtools intersect，但这将使我丢失 data01 中的数据，或者 bedtools merge，但我只想更新 data01 中的一列值而不是合并所有信息。我不确定。我会花些时间在 bedtools 上，如果不行的话，我对 Python 很熟悉。哈哈。 - everestial007

不用担心 bedtools merge 产生的额外列，Bedtools 的速度非常快。你可以将结果加载到 pandas 中，并选择感兴趣的列。 - Gordon Bean

让我们在聊天中继续这个讨论。 - everestial007

Answer 2

你应该在Python中使用区间树函数，它们非常高效且友好，我尝试了类似的东西，遇到了一些问题，但后来解决了，这是我写的代码， 使用区间树查找重叠区域 你可以在此代码基础上进行扩展。

- sbradbio

---

谢谢，我会研究一下。 - everestial007

Answer 3

如果您想要比bedtools更快的工具，或者想要使用Python科学栈的本地工具，请尝试pyranges。

import pyranges as pr

c1 = """Chromosome  Start    End    haplotype_block
2   5207    5867    1856
2   155667    155670    2816
2   67910    68022  2
2   68464    68483  3
2   525    775  132
2   118938    119559    1157"""

c2 = """Chromosome Start End  feature gene_id gene_name   transcript_id
2   5262    5496    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   5579    5750    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   5856    6032    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   6115    6198    exon    scaffold_200003.1   CP5 scaffold_200003.1
2   916 1201    exon    scaffold_200001.1   NA  scaffold_200001.1
2   614 789 exon    scaffold_200001.1   NA  scaffold_200001.1
2   171 435 exon    scaffold_200001.1   NA  scaffold_200001.1
2   2677    2806    exon    scaffold_200002.1   NA  scaffold_200002.1
2   2899    3125    exon    scaffold_200002.1   NA  scaffold_200002.1"""

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

j = gr1.join(gr2).sort()

print(j)
# +--------------+-----------+-----------+-------------------+-----------+-----------+------------+-------------------+-------------+-------------------+
# |   Chromosome |     Start |       End |   haplotype_block |   Start_b |     End_b | feature    | gene_id           | gene_name   | transcript_id     |
# |   (category) |   (int32) |   (int32) |           (int64) |   (int32) |   (int32) | (object)   | (object)          | (object)    | (object)          |
# |--------------+-----------+-----------+-------------------+-----------+-----------+------------+-------------------+-------------+-------------------|
# |            2 |       525 |       775 |               132 |       614 |       789 | exon       | scaffold_200001.1 | nan         | scaffold_200001.1 |
# |            2 |      5207 |      5867 |              1856 |      5262 |      5496 | exon       | scaffold_200003.1 | CP5         | scaffold_200003.1 |
# |            2 |      5207 |      5867 |              1856 |      5579 |      5750 | exon       | scaffold_200003.1 | CP5         | scaffold_200003.1 |
# |            2 |      5207 |      5867 |              1856 |      5856 |      6032 | exon       | scaffold_200003.1 | CP5         | scaffold_200003.1 |
# +--------------+-----------+-----------+-------------------+-----------+-----------+------------+-------------------+-------------+-------------------+
# Unstranded PyRanges object has 4 rows and 10 columns from 1 chromosomes.
# For printing, the PyRanges was sorted on Chromosome.

print(j.df)
#   Chromosome  Start   End  haplotype_block  Start_b  End_b feature            gene_id gene_name      transcript_id
# 0          2    525   775              132      614    789    exon  scaffold_200001.1       NaN  scaffold_200001.1
# 1          2   5207  5867             1856     5262   5496    exon  scaffold_200003.1       CP5  scaffold_200003.1
# 2          2   5207  5867             1856     5579   5750    exon  scaffold_200003.1       CP5  scaffold_200003.1
# 3          2   5207  5867             1856     5856   6032    exon  scaffold_200003.1       CP5  scaffold_200003.1

- The Unfun Cat

Answer 4

s1 = data01.start.values
e1 = data01.end.values
s2 = data02.start.values
e2 = data02['last'].values

overlap = (
    (s1[:, None] <= s2) & (e1[:, None] >= s2)
) | (
    (s1[:, None] <= e2) & (e1[:, None] >= e2)
)

g = data02.gene_id.values
n = data02.gene_name.values

i, j = np.where(overlap)
idx_map = {i_: data01.index[i_] for i_ in pd.unique(i)}

def make_series(m):
    s = pd.Series(m[j]).fillna('').groupby(i).agg(','.join)
    return s.rename_axis(idx_map).replace('', np.nan)

data01.assign(
    gene_id=make_series(g),
    gene_name=make_series(n),
)

- piRSquared