将位置连接到基因组片段中

Question

我想连接相似度得分超过 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)

有什么想法吗？

Answer 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. 如果他们不平等，你就没有给 choose/summarize 他们任何方式。

结果：

> 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

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

Answer 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)

Answer 3

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

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

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

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

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

现在，对于将要连接的连续行，让我们使用逻辑评估找到 "start" 和 "end" 行

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

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

干杯！

将位置连接到基因组片段中

Concatenating positions into genomic segments

r

concatenation

bioinformatics

bioconductor

genome