将位置连接成基因组片段答案

【问题标题】：Concatenating positions into genomic segments将位置连接成基因组片段
【发布时间】：2015-07-28 10:15:41
【问题描述】：

我想连接所有相似度得分超过0.955 的行。 Abo 和 Bel 列分别代表上下行的相似度得分。在以下输入 df 我有 10 个基因组探针（NAME 列），它们仅连接在 4 个基因组片段中（dfout）。

df <- " NAME Abo  Bel Chr GD Position
 BovineHD0100009217 NA 1.0000000   1  0  31691781
 BovineHD0100009218 1.0000000 0.6185430   1  0  31695808
 BovineHD0100019600 0.6185430 0.9973510   1  0  69211537
 BovineHD0100019601 0.9973510 1.0000000   1  0  69213650
 BovineHD0100019602 1.0000000 1.0000000   1  0  69214650
 BovineHD0100019603 1.0000000 0.6600000   1  0  69217942
 BovineHD0100047112 0.6600000 1.0000000   1  0  93797691
 BovineHD0100026604 1.0000000 1.0000000   1  0  93815774
 BovineHD0100026605 1.0000000 0.4649007   1  0  93819471
 BovineHD0100029861 0.4649007 NA   1  0 105042452"
df <- read.table(text=df, header=T)

我的预期输出dfout：

dfout <- "Chr start end startp endp nprob 
           1  31691781 31695808 BovineHD0100009217 BovineHD0100009218 2
           1  69211537 69217942 BovineHD0100019600 BovineHD0100019603 4
           1  93797691 93819471 BovineHD0100047112 BovineHD0100026605 3
           1  105042452 105042452 BovineHD0100029861 BovineHD0100029861 1"
dfout <- read.table(text=dfout, header=T)

有什么想法吗？

【问题讨论】：

附带问题，您使用什么函数来获取探针的相似度分数？
给定种群中给定动物的后续探针之间的基因型相似性。推断具有非常相似基因型分布的探针的基因组片段（CNV 等位基因）的方法。

标签： r concatenation bioinformatics bioconductor genome

【解决方案1】：

我想不出任何使用基本数据框操作的漂亮解决方案，所以这是一个看起来很糟糕的解决方案：

首先，将stringsAsFactors添加到df创建：

df <- read.table(text=df, header=T, stringsAsFactors = FALSE)

start <- df$Position[1]
end <- integer()
output <- NULL
count <- 1
for (i in 1:(nrow(df)-1)) {
  if(df$Bel[i] < 0.955)  {
    end <- df$Position[i]
    output <- rbind(output, c(start, end, df$NAME[which(df$Position == start)], df$NAME[which(df$Position == end)], count))
    start <- df$Position[i+1]
    count <- 0
  } 
  count <- count + 1
}
end <- df$Position[nrow(df)]
output <- as.data.frame(rbind(output, c(start, end, df$NAME[which(df$Position == start)], df$NAME[which(df$Position == end)], count)))
colnames(output) <- c("start", "end", "startp", "endp", "nprob")

这里的基本思想是遍历行并检查是否应该将下一个添加到当前段 (Bel > 0.955) 或是否应该开始一个新段 (Bel <= 0.955)。当必须开始一个新序列时，定义 endrow，将相应的行添加到输出中，并定义新的起始段。 count 用于添加用于创建段的行数 (nprob)。

最后，在 for 循环之外添加最后一段，输出接收其列名并转换为数据帧。我没有使用Chr，因为 1. 它们都是平等的，2. 如果不是，您没有给出任何选择/总结它们的方法。

结果：

> output
      start       end             startp               endp nprob
1  31691781  31695808 BovineHD0100009217 BovineHD0100009218     2
2  69211537  69217942 BovineHD0100019600 BovineHD0100019603     4
3  93797691  93819471 BovineHD0100047112 BovineHD0100026605     3
4 105042452 105042452 BovineHD0100029861 BovineHD0100029861     1

我很确定您或其他人可以对此进行改进，以使其更短、更简洁。

【讨论】：

【解决方案2】：

这里是dplyr 版本。首先我们需要定义组，这就是mutate bit 所做的，然后是组内的简单summarise 函数。

library(dplyr)

df %>% 
  mutate(
   Abo955=ifelse(Abo<0.955,NA,Abo),
   myGroup=cumsum(is.na(Abo955)*1)) %>%
  group_by(myGroup) %>% 
  summarise(
    Chr=min(Chr),
    start=min(Position),
    end=max(Position),
    startp=first(NAME),
    lastp=last(NAME),
    nprob=n()) %>% 
  select(-myGroup)

【讨论】：

【解决方案3】：

此解决方案完全基于逻辑向量并与提供的示例一起使用。

正如 Molx 所说，让我们添加 stringsAsFactors=F

df <- read.table(text=df, header=T, stringAsFactors = F)

为了使逻辑评估起作用，我们将 NA 更改为 0s

df(is.na(df)) <- 0

现在，对于将要连接的连续行，让我们使用逻辑评估找到“开始”和“结束”行

starts <- df$Bel >= 0.955 &  df$Abo < 0.955
ends <- df$Bel < 0.955 &  df$Abo >= 0.955

有了这个我们已经可以构造一个data.frame连接需要连接的行

concatenated <- data.frame(Chr = df[starts, "Chr"], 
                            start = df[starts, "Position"], 
                            end = df[ends, "Position"],
                            startp = df[starts, "NAME"],
                            endp = df[ends, "NAME"],
                            nprob = c( diff (which(starts))[1]  ,diff (which(ends)))
                            )

我们还要构造一个data.frame，其中包含未连接的行，即那些不具有所需相似度分数的行，既不是上面也不是下面的行

notConcatenate <- df$Abo < 0.955 & df$Bel < 0.955

non_concatenated <- data.frame(Chr = df[notConcatenate, "Chr"], 
                            start = df[notConcatenate, "Position"], 
                            end = df[notConcatenate, "Position"],
                            startp = df[notConcatenate, "NAME"],
                            endp = df[notConcatenate, "NAME"],
                            nprob = 1
                            )

最后绑定两个data.frames

dfout <- rbind(concataneted,non_concatenated)

导致

> dfout
  Chr     start       end             startp               endp nprob
1   1  31691781  31695808 BovineHD0100009217 BovineHD0100009218     2
2   1  69211537  69217942 BovineHD0100019600 BovineHD0100019603     4
3   1  93797691  93819471 BovineHD0100047112 BovineHD0100026605     3
4   1 105042452 105042452 BovineHD0100029861 BovineHD0100029861     1

注意：此代码假定相关探针位于同一染色体内

干杯！

【讨论】：