如何在多列中使用相同编号的连续行的频率答案

【问题标题】：How to frequency of consecutive rows with the same number for several columns如何在多列中使用相同编号的连续行的频率
【发布时间】：2015-12-05 13:17:48
【问题描述】：

我有一个数据集如下：

structure(list(chr = c(1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 1, 
1, 0, 0, 1, 1, 1, 1), leftPos = c(240000, 1080000, 1200000, 1320000, 
1440000, 1800000, 2400000, 2520000, 3120000, 3360000, 3480000, 
3600000, 3720000, 4200000, 4560000, 4920000, 5040000, 5160000, 
5280000, 6e+06), chr.1 = c(1, 0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 1, 
1, 0, 0, 1, 1, 1, 1, 1), leftPos.1 = c(240000, 1080000, 1200000, 
1320000, 1440000, 1800000, 2400000, 2520000, 3120000, 3360000, 
3480000, 3600000, 3720000, 4200000, 4560000, 4920000, 5040000, 
5160000, 5280000, 6e+06), ASample = c(0, 
0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0), Sample1 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 1, 0, 1), Sample2 = c(0, 
1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1, 1), Sample3 = c(0, 
1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1), Sample4 = c(0, 
0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1), Sample5 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1), Sample6 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1), Sample7 = c(0, 
0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 1), Sample8 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1), Sample9 = c(0, 
0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1), Sample10 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1), Sample11 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1), Sample12 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1), Sample13 = c(0, 
0, 1, 0, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0), Sample14 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 1), Sample15 = c(0, 
1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 1)), .Names = c("chr", 
"leftPos", "chr.1", "leftPos.1", "Sample1", 
"Sample2", 
"Sample3", "Sample4", 
"Sample5", "Sample6", 
"Sample7", "Sample8", 
"Sample9", "Sample10", 
"Sample11", "Sample12", 
"Sample13", "Sample14", 
"Sample15"), row.names = c(NA, 
20L), class = "data.frame")

我需要计算每列有多个相同的 1 或 -1 的行数

我希望能够计算每列的连续行数，按 chr 分组，染色体中有三个连续的 1 或 -1（称为chr 的列）。

理想的输出是这样的（不是从上面的 dput 数据中获取的）

chr numberOfConsecutive1s FreqSample1  FreqSample2  FreqSample3 etc
1          2                3           2               14
1          3                5           2               2
1          4                5           0               6
1          5                4           3               5
1          6                3           0               3
1          7                7           5               7
1          8                5           0               2
1          9                54          2               6
1          10               34          77              7
2          2                6           4               2
2          3                23          34              34
2          4                5           37              2
2          5                55          24              22
2          6                2           0               11
2          7                3           14              5
2          8                2           5               77
2          9                5           23              34
2          10               5           11              34
3          1                32          0               2

到目前为止，我已经尝试了以下方法，它只是将非连续的 1 转换为 0，所以我只剩下连续的 1。我不知道如何根据所需的输出来计算它们。

dx<-DAT_list2res
f0 <- function( colNr, dx )
{
  col <- dx[,colNr]
  n1 <- which( col == 1 )            # The `1`-rows.
  d0 <- which( diff(col) == 0 )      # Consecutive entries are equal.
  dc0 <- which( diff(dx[,1]) == 0 )  # Same chromosome.
  m <- intersect( n1-1, intersect( d0, dc0 ) )
  return ( setdiff( 1:nrow(dx), union(m,m+1) ) )
}
g <- function( dx )
{
  for ( i in 3:ncol(dx) ) { dx[f0(i,dx),i] <- 0 }  
  return ( dx )
}
dx<-g(dx)

编辑

我也按照 bramtayl 的建议尝试了这个：

result = 
  consecFreq %>%
  select(-chr) %>%
  gather(variable, chr,  5:190) %>%
  group_by(variable) %>%
  mutate(ID = 
           chr %>%
           lag %>%
           `!=`(chr) %>%
           plyr::mapvalues(NA, FALSE) %>%
           cumsum) %>%
  count(variable, chr, ID) %>%
  rename(numberOfConsecutive1s = n) %>%
  count(variable, chr, numberOfConsecutive1s) %>%
  spread(variable, n, fill = 0)

但它给了我一个“索引超出范围”的错误。如果我忽略扩展线，我也会得到一个奇怪的输出，所以我不确定这是不是答案

【问题讨论】：

标签： r

【解决方案1】：

修订

根据澄清，此方法对每个染色体使用rle 函数来查找连续1 或-1 的运行，然后table 计算每个值的运行次数。这会为没有特定值计数的样本提供NA，因此如果有帮助，代码的最后一行会将NA's 转换为0's。最后，您的structure 输入似乎存在问题，因为structure 的.Names 部分缺少Cytospongex10_SLX.9395.FastSeqK.fq.gz.res。这会导致所有列名被移动，最后一个列名是NA，这可能会导致执行问题。

下面的代码将正确的名称分配给输入数据（在data.framedf），然后如上所述计算频率。

    colnames(data) <- c("chr", 
                    "leftPos", "chr.1", "leftPos.1", "Cytospongex10_SLX.9395.FastSeqK.fq.gz.res", "Sample1", 
                    "Sample2", 
                    "Sample3", "Sample4", 
                    "Sample5", "Sample6", 
                    "Sample7", "Sample8", 
                    "Sample9", "Sample10", 
                    "Sample11", "Sample12", 
                    "Sample13", "Sample14", 
                    "Sample15")

 chr_labels <- sort(unique(data$chr))
 sampl_freqs <- data.frame(chr=1,  numberOfConsecutive1s=1, count=0)

for( sampl in colnames(data)[-(1:5)]) {
  freqs <- data.frame()
  for( chr in chr_labels )  {
     runs  <-  rle(data[data$chr == chr,sampl]) 
     freqs_chr <- data.frame(chr=chr, table(runs$length[runs$values %in% c(-1,1)], dnn = "numberOfConsecutive1s") )
     freqs <- rbind(freqs, freqs_chr)
   }
  sampl_freqs <- merge.data.frame(sampl_freqs, freqs, by = c("chr","numberOfConsecutive1s"), all=TRUE)
  colnames(sampl_freqs) <- c(head(colnames(sampl_freqs),-1),paste("Freq",sampl,sep=""))
}
# clean up from sampl_freqs definition
 sampl_freqs <- sampl_freqs[,-3]
 #  To convert NA's to 0
 sampl_freqs <- data.frame(sampl_freqs[,1:2], sapply(sampl_freqs[,-(1:2)], function(x) ifelse(is.na(x), 0, x)))

与上面类似，但使用dplyr

library(reshape2)
library(dplyr)

df <- melt(data[,-(2:5)], id.vars="chr",  variable.name="sample")
sampl_freqs <- df %>% group_by(sample, chr )   %>%
 do(data.frame(unclass(rle(.$value))) %>%
      filter(values %in% c(-1,1)) ) %>%
 group_by(sample, chr, lengths) %>%
  summarize(Freq = n() ) %>%
 dcast( chr + lengths ~ sample, value.var = "Freq" ) 
sampl_freqs <- with(sampl_freqs,data.frame( chr, numberOfConsecutive1s = lengths , 
                                            sapply(sampl_freqs[,-(1:2)], function(x) ifelse(is.na(x), 0, x))))

【讨论】：

不错，但染色体编号似乎不存在，由于某种原因，它们给了我-1和0和1的混合，而不是实际的染色体编号1到22。等一下.. .
不，实际的染色体编号似乎不是我整个数据集中的染色体 - 它给了我 -1,0 和 1。在我的数据集中，chr 是第一列，所以我不确定这是怎么回事正在提出这些价值观
数据框中的第一列，标签 chr，仅包含 0"s 和 1's 。这些是染色体编号吗？另外，标准 3 个连续 1 是否具有相同的 chr 列号？如果所以，我在你的数据中没有看到任何这样的例子。你能指出哪些样本集包含这样的例子吗？谢谢。
是标签 chr 是染色体编号，我很欣赏它们不是 1 到 22，尽管它仍然不能解释为什么我在结果中的该列中得到 -1 0 和 1。标准是在每个样本的 chr 内连续行数为 1 或 -1 的频率。数据中有两个和三个连续 1（但没有 -1）的示例。
您的第二个答案效果很好。正是我想要的。非常感谢

【解决方案2】：

我想你想要这样的东西：

library(dplyr)
library(tidyr)

min_chunk_length = 1

result = 
  data %>%
  rename(chromosome = chr) %>%
  select(chromosome, Sample1:Sample15) %>%
  gather(sample, value, Sample1:Sample15) %>%
  group_by(chromosome, sample) %>%
  mutate(non_zero = value %in% c(1, -1),
         chunk_ID = 
           non_zero %>%
           lag %>%
           `!=`(non_zero) %>%
           plyr::mapvalues(NA, FALSE) %>%
           cumsum) %>%
  filter(non_zero = TRUE) %>%
  group_by(chromosome, sample, chunk_ID) %>%
  mutate(length_of_chunk = n()) %>%
  filter(length_of_chunk > min_chunk_length) %>%
  count(chromosome, sample) %>%
  spread(sample, n, fill = 0)

【讨论】：

OK 看起来很有希望，尽管我遇到的一个问题是我的样本并不总是连续命名，并且可以有随机名称。要将其合并到收集语句中，我该怎么做。我尝试了 st
好的，所以我只输入了列号，它到达了扩展语句，然后我得到了“索引超出范围”错误。知道为什么吗？到目前为止，它似乎运行良好
我不知道我的代码对你的示例数据有用
刚刚用发布的数据和发布的代码再次尝试。我仍然得到同样的错误。我的 dplyr 是 0.3.0.2 版本，tidyr 是 0.1
我有 dplyr 0.4.3 和 tidyr 0.3.1！我添加了一个 select 语句，因此您可以选择要对其进行分析的任何列。