如何使用多条染色体计算基因之间的距离

Question

我想知道某些基因是否聚集在一起。现在，我已经有了一个基因列表，以及它们的起始和终止位置，并且我已经知道如何计算这些基因之间的距离。问题是我不知道如何考虑染色体的转换。

您无法测量 1 号染色体上的基因和 2 号染色体上的基因之间的距离。

我想过这样计算距离：基因 2 的起始位置 - 基因 1 的停止位置。然后，你就有了这些基因之间的距离。

但是我该如何解释：当你到达下一条染色体时，R 代码会抓取 2 号染色体上基因的起始位置，但会抓取 1 号染色体上基因的终止位置，这是不可能的（对于至少我的研究）。

所以我想知道如何在 R 中解释这一点。如果基因位于不同的染色体上，我只需要以某种方式跳过它们。

希望你们能帮帮我。

关于下面的代码：三个向量只是起点和终点位置的向量，以及染色体。它们都等长。 染色体是包含每个基因的染色体编号的向量

start_vector <- as.vector(sorted_coords$start_position)
end_vector <- as.vector(sorted_coords$end_position)
chromosomes <- as.vector(sorted_coords$chromosome_name)

chromosomes[is.na(chromosomes)] <- 24

count = 0
for(i in 1:length(chromosomes)){
  if(count != chromosomes[i]){
    start <- i - 1
    end <- i + 1

    start_vector <- start_vector[-start]
    end_vector <- end_vector[-end]
    count <- count + 1  
    }

}

我期望一个包含所有基因距离的向量，不包括位于不同染色体上的基因的距离。

Answer 1

library(tidyverse) # for all the tidyverse goodies
library(reshape2) # For the melt function

由于您没有提供可重现的示例，我冒昧地制作了自己的玩具数据框，如下所示。它只有2条染色体，但这种方法应该适用于任意数量的染色体和基因。

sorted_coords <- tibble(start_position = abs(rnorm(10)*10),
                        end_position = abs(rnorm(10)*10),
                        chromosome_name = c(rep(1,5),rep(2,5)))

编辑：OP 澄清说，他们想找到与基因的距离，而不是与其他所有基因的距离。后半部分的方法在底部，因为我觉得它很有趣。新的解决方案在这里：

sorted_coords %>% 
  group_by(chromosome_name) %>%  
  arrange(chromosome_name, start_position) %>% 
  mutate(distance = start_position - lag(end_position, n = 1, default = 0))

1. 我们按染色体分组，这样我们就不会在染色体之间进行任何错误的计算。

2. 我们按染色体名称排列，以便在最后进行排序。我们按起始位置排列，因此基因的顺序正确。

3. 我们按照建议计算距离。当前行的开始位置 - 上一行的结束位置。我们指定（尽管它是默认值）我们查看之前的行，并且如果之前没有行，则结束位置的值默认为 0。

---

旧答案

如果您想将每个基因与其他基因进行比较，最快的方法是创建一个矩阵。正如你所指定的，我们想将基因 1 的开头减去基因 2 的结尾。这对我来说感觉不对，但我已经有一段时间没有做 biochem 了:)。因为你想要一个单一的对列表，我们可以折叠它（melt 函数）。

下面的代码理解起来有点模糊，所以让我们分解一下。

sorted_coords %>% 
  group_by(chromosome_name) %>% 
  do( outer(.$start_position, .$end_position) %>% 
        melt() %>% 
        setNames(c("rows", "columns", "distance")))

1. 我们采用数据框并根据您的需要按每条染色体对其进行分组。
2. do 命令可以让我们进行复杂的操作。 group_by 命令可确保我们所做的任何事情对于每条染色体基本上都是独立的。
3. 外部函数为我们创建矩阵。 . 是我们传递的数据框（子集到特定染色体）。我们传递需要找出差异的两列。
4. melt 函数将矩阵转换为数据框，以便我们指定它用于计算差异的两个基因。它们以数字形式列出，您可以返回并进行比较。我建议使用安排来设置顺序，以便您可以轻松地引用它。
5. setNames 只是将列名设置为更可红色的名称。

这应该比为所有进程运行 for 循环要快得多。如果您提供更多信息，我可能会更多地清理答案。