Применяя sapply или другую функцию apply вместо вложенного цикла for для списков фреймов данных

Question

Применяя sapply или другую функцию apply вместо вложенного цикла for для списков фреймов данных

У меня есть два фрейма данных / списки данных HumanSplitandratSplit` и они имеют вид

> ratSplit$Kidney_F_GSM1328570
  ratGene        ratReplicate alignment RNAtype
1    Crot Kidney_F_GSM1328570         7     REV
2    Crot Kidney_F_GSM1328570        12     REV
3    Crot Kidney_F_GSM1328570         4     REV

а также

> humanSplit$Fetal_Brain_408_AGTCAA_L009_R1_report.txt
   humanGene                            humanReplicate alignment RNAtype
53     ZFP28 Fetal_Brain_408_AGTCAA_L009_R1_report.txt         5     reg
55     RC3H1 Fetal_Brain_408_AGTCAA_L009_R1_report.txt         9     reg
56     IFI27 Fetal_Brain_408_AGTCAA_L009_R1_report.txt         4     reg

И еще один файл, используемый ниже, это geneList вида:

ABAT,Abat
ABCA1,Abca1
ABCA12,Abca12
ABCA2,Abca2
ABCA3,Abca17
ABCA4,Abca4
ABCA5,Abca5

Теперь я хочу сделать точный тест Фишера между всеми сочетаниями пары элементов между ratSplit а также humanSplit после некоторых манипуляций с данными. И в конечном итоге хочу написать результат теста Фишера в csv файл. Прямо сейчас я делаю двойной цикл. Но мне интересно, как я могу использовать sapply или другие связанные вещи, чтобы сделать его более эффективным.

В настоящее время я делаю следующее: здесь я впервые делаю data.frameresult где я храню / добавляю всю информацию, полученную из теста Фишера, в паре на каждом шаге. Наконец, когда весь цикл завершен, я пишу окончательный result в csv файл. мое понимание заключается в использовании sapply Мне нужно преобразовать внутреннюю часть цикла в функцию, а затем вызвать sapply. Но я не уверен, что это лучший способ оптимизировать его. Любая помощь приветствуется

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)

    fisherTest <- fisher.test(contingencyTable) 
    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)
  }
}

write.table(result, file = "newData5.csv", row.names = FALSE, append = FALSE, col.names = TRUE, sep = ",")

0

r list for-loop apply sapply

Источник

user3587282 09 авг '14 в 22:44

1 ответ

Решение

Другие вопросы по тегам r list for-loop apply sapply

user3358272 10 авг '14 в 03:12 2014-08-10 03:12 · Accepted Answer · 2014-08-10 03:12

Это сложно проверить, так как нам не хватает geneList, но я предполагаю, что код работает правильно, так что вы просто хотите ускорить процесс. Вот несколько советов, чтобы помочь:

Не предопределять result как это. Установив в качестве первой записи каждого столбца строку имени столбца, вы гарантируете, что все последующие записи будут преобразованы в строки. Хотя это не обязательно кусает вас, так как вы все равно пишете в CSV, это плохая форма. Он укусит вас, если вы намереваетесь прочитать его обратно в R, поскольку первая строка будет использоваться в качестве имени столбца, но вторая строка будет строкой, поэтому последующие данные также будут строками. (Затем вам придется очистить свои собственные данные, расточительно.)
В конце вашего сценария вы звоните rbind, Это может сделать хорошо на некоторое время, но повторные звонки rbind в результате каждый раз копируется весь data.frame. При добавлении большего количества строк это приведет к значительному замедлению вашего кода. Это можно исправить одним из двух способов, перечисленных ниже.
Так как вы используете каждый из names(HmRnTable) == "HyRy" дважды, моя техника состоит в том, чтобы сначала сохранить это в вектор (или скаляр, если вы используете which(...)), а затем использовать эту переменную в подмножестве HmRnTable, Скорее всего, это немного ускорит процесс, но также может немного облегчить чтение кода. Фактически, вы можете сократить каждое из этих четырех назначений до чего-то вроде (не проверено):
```
idx <- is.integer0(HmRnTable[names(HmRnTable) == 'HyRy'])
HyRy <- HmRnTable[idx][[1]]
HyRy[idx] <- 0
## repeat for HyRn, HnRy, and HnRn
```
Я настоятельно призываю вас взять большую часть кода и поместить его в функцию, которая принимает два аргумента (i а также j) или четыре аргумента (list1, index1, list2, index2). (То, что вы делаете, зависит от того, сколько у вас OCD в отношении ссылок на переменную область.) Я предполагаю, что функция будет буквально возвращать результаты из fisher.test без массажа. Это значительно упростит тестирование при создании функции и позволит в дальнейшем включить его в этот скрипт. Я буду ссылаться на функцию как myfisher(i,j) ниже.
Я предполагаю, что вам нужно выполнить большое количество сравнений (поскольку каждая итерация не должна занимать так много времени). @Floo0 комментарий о outer может работать, как может expand.grid, В любом случае, вы сравниваете каждый элемент list1 с каждым элементом list2. Если вы начнете с:
```
(eg <- expand.grid(i = 1:length(ratSplit),
                   j = 1:length(humanSplit)))
##    i j
## 1  1 1
## 2  2 1
## 3  3 1
## 4  4 1
## ...
```
Это дает нам простой data.frame, на котором мы можем использовать apply, Честно говоря, мне нравится элегантность ddply в этом случае, поскольку он легко переходит из data.frame в data.frame, легко.
```
library(plyr)
ret <- ddply(eg, .(i, j), function(df) {
    with(myfisher(df$i, df$j),
         data.frame(pv = p.value, ci1 = conf.int[1], ci2 = conf.int[2],
                    alt = alternative))
})
```
Обратите внимание, что я явно не включил humanReplicateName или же ratReplicateName так как они могут быть добавлены до ddply (или после, ссылаясь на ret$ вместо eg$) с:
```
eg$ratReplicateName <- names(ratSplit[ eg$i ])
eg$humanReplicateName <- names(humanSplit[ eg$i ])
```
и имена будут волшебным образом в выводе тоже. Меньше иметь дело с внутри цикла.
До настоящего времени это приводило к одному data.frame, который вы затем можете сохранить в CSV.
Я сделаю еще одну рекомендацию, которая может быть излишней в зависимости от того, как долго это будет продолжаться. Время от времени я обнаруживал, что мои долгосрочные сценарии прерываются, возможно, потому, что мне пришлось что-то настраивать на лету; Я нашел ошибку; или если компьютер должен был перезагрузиться. Нет простого способа разрешить продолжение в середине потока, но я разработал несколько методов, чтобы смягчить это.
Вместо ddplyиспользовать d_ply который не принимает тип возвращаемого значения. Вместо того, чтобы возвращать (однорядный) data.frame, немедленно сохраните эту строку (с заголовками или без них, ваш вызов) в файл. Хотя приведенный ниже код не является действительно "атомарным" и, следовательно, могут возникать условия гонки, он достаточно устойчив, чтобы удовлетворить большинство наших потребностей:
```
library(plyr)
savedir <- './output'
ret <- ddply(eg, .(i, j), function(df) {
    fn <- file.path(savedir, sprintf('%s-%s.csv', df$i, df$j))
    if (! file.exists(fn)) {
        ret <- with(myfisher(df$i, df$j),
                    data.frame(pv = p.value, ci1 = conf.int[1], ci2 = conf.int[2],
                               alt = alternative))
        write.table(ret, file = fn, sep = ",", append = FALSE, 
                    row.names = FALSE, col.names = TRUE)
    }
})
```
Одним из достоинств этого является то, что если / когда оно прервано, самое большее, что вам нужно будет сделать, - это просмотреть все файлы в "./output/", удалить файлы с 0 байтами, а затем повторно запустить, и это будет только выполнить на недостающие файлы. О, и это становится лучше.
Распараллеливая это. Если вам нужно зайти так далеко (и некоторые функции не видят такого улучшения, как другие), вы можете вместо этого использовать несколько ядер в своей системе. Вы могли бы сделать что-то вроде этого:
```
library(parallel)
cl <- makeCluster(detectCores() - 1) # I like to keep one free
clusterEvalQ(cl, {
    load('myTwoLists.rda') # with ratSplit and humanSplit variables
    source('mycode.R') # with myfisher(i,j)
    ## any libraries you may want/need to add, if you get more advanced
})
eg <- expand.grid(i = 1:length(ratSplit),
                  j = 1:length(humanSplit))
eg$ratReplicateName <- names(ratSplit[ eg$i ])
eg$humanReplicateName <- names(humanSplit[ eg$i ])
ign <- parApply(cl, eg, 1, function(r) {
    i <- r[1] ; j <- r[2]
    fn <- file.path(savedir, sprintf('%s-%s.csv', i, j)
    if (! file.exists(fn)) {
        ret <- with(myfisher(i, j),
                    data.frame(ratName = r[3], humanName = r[4],
                               pv = p.value, ci1 = conf.int[1], ci2 = conf.int[2],
                               alt = alternative))
        write.table(ret, file = fn, sep = ",", append = FALSE, 
                    row.names = FALSE, col.names = TRUE)
    }
})
stopCluster(cl)
```
Обратите внимание, что мы больше не ссылаемся на строку как на data.frame. Я был бы рад, если бы код Хэдли изначально работал с параллелью. ~~(Насколько я знаю, это так, и я просто еще не нашел его!)~~ Редактировать: ни один из моих прошлых проектов не использовал parallel использовали ddplyи наоборот, поэтому я никогда не играл с ddply(..., .parallel = TRUE) который использует foreach,

Предостережение Emptor: этот код не был протестирован в этом контексте, хотя он почти скопирован / вставлен из рабочего кода. Из-за правок я могу быть один на один или пропустить парен. Надеюсь, поможет!