2

我有一个代表细菌模型基因结构的文件。每行代表一个模型。行是一个固定长度的二进制字符串，其中存在基因（1 表示存在，0 表示不存在）。我的任务是比较每对模型的基因序列，并计算它们的相似程度，然后计算出一个相异矩阵。

一个文件总共有450个模型（行），有250个文件。我有一个工作代码，但是只为一个文件完成整个工作大约需要 1.6 小时。

#Sample Data    
Generation: 0
    [0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 1, 0, 0]
    [1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1]
    [1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0]
    [0, 1, 1, 0, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0, 0, 0]
    [0, 1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0]
    [1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0]

我的代码做了什么：

  1. 读取文件
  2. 将二进制字符串转换为数据框 Gene, Model_1, Model_2, Model_3, ... Model_450
  3. 运行嵌套的 for 循环来进行成对比较（仅矩阵的上半部分）——我取两个相应的列并将它们相加，然后计算总和为 2 的位置（意味着两个模型中都存在）
  4. 将数据写入文件
  5. 稍后创建矩阵

比较代码

generationFiles = list.files(pattern = "^Generation.*\\_\\d+.txt$")

start.time = Sys.time()

for(a in 1:length(generationFiles)){

  fname = generationFiles[a]

  geneData = read.table(generationFiles[a], sep = "\n", header = T, stringsAsFactors = F)

  geneCount = str_count(geneData[1,1],"[1|0]")

  geneDF <- data.frame(Gene = paste0("Gene_", c(1:geneCount)), stringsAsFactors = F)

  #convert the string into a data frame
      for(i in 1:nrow(geneData)){

    #remove the square brackets
    dataRow = substring(geneData[i,1], 2, nchar(geneData[i,1]) - 1)

    #removing white spaces
    dataRow = gsub(" ", "", dataRow, fixed = T)

    #splitting the string 
    dataRow = strsplit(dataRow, ",")

    #converting to numeric
    dataRow = as.numeric(unlist(dataRow))

    colName = paste("M_",i,sep = "")
    geneDF <- cbind(geneDF, dataRow)
    colnames(geneDF)[colnames(geneDF) == 'dataRow'] <- colName

    dataRow <- NULL
  }

  summaryDF <- data.frame(Model1 = character(), Model2 = character(), Common = integer(),
                          Uncommon = integer(), Absent = integer(), stringsAsFactors = F)

  modelNames = paste0("M_",c(1:450))

  secondaryLevel = modelNames

  fileName = paste0("D://BellosData//GC_3//Summary//",substr(fname, 1, nchar(fname) - 4),"_Summary.txt")

  for(x in 1:449){

    secondaryLevel = secondaryLevel[-1]

    for(y in 1:length(secondaryLevel)){

      result = geneDF[modelNames[x]] + geneDF[secondaryLevel[y]]

      summaryDF <- rbind(summaryDF, data.frame(Model1 = modelNames[x],
                                               Model2 = secondaryLevel[y],
                                               Common = sum(result == 2),
                                               Uncommon = sum(result == 1),
                                               Absent = sum(result == 0)))

    }


  }

  write.table(summaryDF, fileName, sep = ",", quote = F, row.names = F)
  geneDF <- NULL
  summaryDF <- NULL
  geneData <-NULL

}

转换为矩阵

maxNum = max(summaryDF$Common)
  normalizeData = summaryDF[,c(1:3)]
  normalizeData[c('Common')] <- lapply(normalizeData[c('Common')], function(x) 1 - x/maxNum)

  normalizeData[1:2] <- lapply(normalizeData[1:2], factor, levels=unique(unlist(normalizeData[1:2]))) 

  distMatrixN = xtabs(Common~Model1+Model2, data=normalizeData)

  distMatrixN = distMatrixN + t(distMatrixN)

有没有办法让这个过程运行得更快？有没有更有效的比较方法？

4

1 回答 1

3

这段代码应该更快。嵌套循环在 R 中是噩梦般的缓慢。像rbind-ing一次一行这样的操作也是 R 编程中最糟糕和最慢的想法之一。

生成 450 行，每行包含 20 个 0、1 的元素。

M = do.call(rbind, replicate(450, sample(0:1, 20, replace = T), simplify = F))

生成组合(450, 2) 行对数的列表

L = split(v<-t(utils::combn(450, 2)), seq(nrow(v))); rm(v)

应用您想要的任何比较功能。在这种情况下，每个行组合的相同位置的 1 的数量。如果要计算不同的指标，只需编写另一个函数（x），其中M[x[1],]第一行M[x[2],]是第二行。

O = lapply(L, function(x) sum(M[x[1],]&M[x[2],]))

代码在相当慢的 2.6 Ghz Sandy Bridge 上耗时约 4 秒

使用您的结果获取一个干净的 data.frame，三列：第 1 行、第 2 行、两行之间的度量

data.frame(row1 = sapply(L, `[`, 1),
           row2 = sapply(L, `[`, 2),
           similarity_metric = do.call(rbind, O))

老实说，我没有彻底梳理您的代码以准确复制您正在做的事情。如果这不是您要查找的内容（或无法修改以实现您要查找的内容），请发表评论。

于 2018-04-19T21:08:41.687 回答