Огромный файл данных и запуск нескольких параметров и проблемы с памятью, тест Фишера

У меня есть код R, который я пытаюсь запустить на сервере. Но он останавливается посередине / зависает, вероятно, из-за ограничения памяти. Файлы данных огромны / массивны (один имеет 20 миллионов строк), и если вы посмотрите на двойной цикл for в коде, length(ratSplit) = 281 а также length(humanSplit) = 36, Данные имеют конкретные данные о генах человека и крысы, и у человека 36 повторов, а у крысы 281. Таким образом, цикл в основном состоит из 281*36 шагов. Что я хочу сделать, это обработать данные с помощью функции getGeneType и посмотрите, как различны / независимы выражения различных повторяющихся комбинаций. Используя тест Фишера. Данные rat_processed_7_25_FDR_05.out выглядят так:

2       Sptbn1  114201107       114200202       chr14|Sptbn1:114201107|Sptbn1:114200202|reg|-   2       Thymus_M_GSM1328751     reg
2       Ndufb7  35680273        35683909        chr19|Ndufb7:35680273|Ndufb7:35683909|reg|+     2       Thymus_M_GSM1328751     rev
2       Ndufb10 13906408        13906289        chr10|Ndufb10:13906408|Ndufb10:13906289|reg|-   2       Thymus_M_GSM1328751     reg
3       Cdc14b  1719665 1719190 chr17|Cdc14b:1719665|Cdc14b:1719190|reg|-       3       Thymus_M_GSM1328751     reg

и данные fetal_output_7_2.out имеет форму

SPTLC2  78018438        77987924        chr14|SPTLC2:78018438|SPTLC2:77987924|reg|-     11      Fetal_Brain_408_AGTCAA_L006_R1_report.txt       reg
EXOSC1  99202993        99201016        chr10|EXOSC1:99202993|EXOSC1:99201016|rev|-     5       Fetal_Brain_408_AGTCAA_L006_R1_report.txt       reg
SHMT2   57627893        57628016        chr12|SHMT2:57627893|SHMT2:57628016|reg|+       8       Fetal_Brain_408_AGTCAA_L006_R1_report.txt       reg
ZNF510  99538281        99537128        chr9|ZNF510:99538281|ZNF510:99537128|reg|-      8       Fetal_Brain_408_AGTCAA_L006_R1_report.txt       reg
PPFIBP1 27820253        27824363        chr12|PPFIBP1:27820253|PPFIBP1:27824363|reg|+   10      Fetal_Brain_408_AGTCAA_L006_R1_report.txt       reg

Теперь у меня есть несколько вопросов о том, как сделать это более эффективным. Я думаю, что когда я запускаю этот код, R занимает много памяти, что в конечном итоге вызывает проблемы. Мне интересно, есть ли способ сделать это более эффективно

Другая возможность - использование двойного цикла for. Поможет ли? В таком случае, как я должен применять саппли?

В конце хочу конвертировать result в CSV-файл. Я знаю, это немного ошеломляюще, чтобы поместить код, подобный этому. Но любая оптимизация / эффективное кодирование / программирование будет МНОГО! Мне действительно нужно запустить все, по крайней мере, один, чтобы получить данные в ближайшее время.

#this one compares reg vs rev
date()
ratRawData <- read.table("rat_processed_7_25_FDR_05.out",col.names = c("alignment", "ratGene", "start", "end", "chrom", "align", "ratReplicate", "RNAtype"), fill = TRUE)
humanRawData <- read.table("fetal_output_7_2.out", col.names = c("humanGene", "start", "end", "chrom", "alignment", "humanReplicate", "RNAtype"), fill = TRUE)
geneList <- read.table("geneList.txt", col.names = c("human", "rat"), sep = ',')
#keeping only information about gene, alignment number, replicate and RNAtype, discard other columns
ratRawData <- ratRawData[,c("ratGene", "ratReplicate", "alignment", "RNAtype")]
humanRawData <- humanRawData[, c( "humanGene", "humanReplicate", "alignment", "RNAtype")]

#function to capitalize
capitalize <- function(x){
  capital <- toupper(x) ## capitalize 
  paste0(capital)
}

#capitalizing the rna type naming for rat. So, reg ->REG, dup ->DUP, rev ->REV
#doing this to make data manipulation for making contingency table easier.
levels(ratRawData$RNAtype) <- capitalize(levels(ratRawData$RNAtype))

#spliting data in replicates
ratSplit <- split(ratRawData, ratRawData$ratReplicate)
humanSplit <- split(humanRawData, humanRawData$humanReplicate)

print("done splitting")
#HyRy :when some gene has only reg, rev , REG, REV
#HnRy : when some gene has only reg,REG,REV
#HyRn : add 1 when some gene has only reg,rev,REG
#HnRn : add 1 when some gene has only reg,REG

#function to be used to aggregate 
getGeneType <- function(types) {
  types <- as.character(types)
  if ('rev' %in% types) {
    return(ifelse(('REV' %in% types), 'HyRy', 'HyRn'))
  } 
  else {
    return(ifelse(('REV' %in% types), 'HnRy', 'HnRn'))
  }
}

#logical function to see whether x is integer(0) ..It's used the for loop bellow in case any one HmYn is equal to zero
is.integer0 <- function(x) {
  is.integer(x) && length(x) == 0L
}

result <- data.frame(humanReplicate = "human_replicate", ratReplicate = "rat_replicate", pvalue = "p-value", alternative = "alternative_hypothesis", 
                     Conf.int1 = "conf.int1", Conf.int2 ="conf.int2", oddratio = "Odd_Ratio")
for(i in 1:length(ratSplit)) {
  for(j in 1:length(humanSplit)) {
    ratReplicateName <- names(ratSplit[i])
    humanReplicateName <- names(humanSplit[j])

    #merging above two based on the one-to-one gene mapping as in geneList defined above.
    mergedHumanData <-merge(geneList,humanSplit[[j]], by.x = "human", by.y = "humanGene")
    mergedRatData <- merge(geneList, ratSplit[[i]], by.x = "rat", by.y = "ratGene")

    mergedHumanData <- mergedHumanData[,c(1,2,4,5)] #rearrange column
    mergedRatData <- mergedRatData[,c(2,1,4,5)]  #rearrange column
    mergedHumanRatData <- rbind(mergedHumanData,mergedRatData) #now the columns are "human", "rat", "alignment", "RNAtype"

    agg <- aggregate(RNAtype ~ human+rat, data= mergedHumanRatData, FUN=getGeneType) #agg to make HmYn form
    HmRnTable <- table(agg$RNAtype) #table of HmRn ie RNAtype in human and rat.

    #now assign these numbers to variables HmYn. Consider cases when some form of HmRy is not present in the table. That's why
    #is.integer0 function is used
    HyRy <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HyRy"]), 0, HmRnTable[names(HmRnTable) == "HyRy"][[1]])
    HnRn <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HnRn"]), 0, HmRnTable[names(HmRnTable) == "HnRn"][[1]])
    HyRn <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HyRn"]), 0, HmRnTable[names(HmRnTable) == "HyRn"][[1]])
    HnRy <- ifelse(is.integer0(HmRnTable[names(HmRnTable) == "HnRy"]), 0, HmRnTable[names(HmRnTable) == "HnRy"][[1]])

    contingencyTable <- matrix(c(HnRn,HnRy,HyRn,HyRy), nrow = 2)
    # contingencyTable: 
    # HnRn --|--HyRn
    # |------|-----|
    # HnRy --|-- HyRy
    #

    fisherTest <- fisher.test(contingencyTable)

    #make new line out of the result of fisherTest
    newLine <- data.frame(t(c(humanReplicate = humanReplicateName, ratReplicate = ratReplicateName, pvalue = fisherTest$p,
                              alternative = fisherTest$alternative, Conf.int1 = fisherTest$conf.int[1], Conf.int2 =fisherTest$conf.int[2], 
                              oddratio = fisherTest$estimate[[1]])))

    result <-rbind(result,newLine) #append newline to result
    if(j%%10 = 0) print(c(i,j))
  }
}

write.table(result, file = "compareRegAndRev.csv", row.names = FALSE, append = FALSE, col.names = TRUE, sep = ",")
Источник

1 ответ

Решение

Ссылаясь на принятый ответ на монитор использования памяти в R, объем памяти, используемый R можно отслеживать с gc(),

Если сценарию действительно не хватает памяти (что меня не удивит), самый простой способ решить эту проблему - переместить write.table() снаружи внутрь петли, чтобы заменить rbind(), Было бы просто необходимо создать новое имя файла для файла CSV, которое записывается из каждого вывода, например:

csvFileName <- sprintf("compareRegAndRev%03d_%03d.csv",i,j)

Если файлы CSV пишутся без заголовков, их можно объединить отдельно R (например, используя cat в Unix) и заголовок добавлен позже.

Хотя при таком подходе может быть успешно создан искомый CSV-файл, возможно, этот файл слишком велик для последующей обработки. Если это так, может быть предпочтительнее обрабатывать файлы CSV отдельно, а не объединять их вообще.

Источник
Другие вопросы по тегам