Как я могу написать свою собственную функцию для автоматического вычисления между клиническими признаками и генами в R - PullRequest
0 голосов
/ 24 февраля 2020

Я пытаюсь написать свою собственную функцию в R, роль которой заключается в автоматическом вычислении корреляции между генами и интересующими клиническими признаками. Это мои строки кода:

#Empty data.frame
cc1 <- data.frame(Estimate=paste("Site", 1:35), P.value="")
estimates = numeric(35)
pvalues = numeric(35)

#compute correlation between clinical feature and genes
computeCC = function(x)
{
  if (x = ""){
  for (i in 1:35) {
      cc<-cor.test(cor[,i], cor[,x],    
                   method = "spearman")
      estimates[i] = cc$estimate
      pvalues[i] = cc$p.value
      cc1$Estimate <- estimates
      cc1$P.value <- pvalues
      rownames(cc1) = colnames(cor)[1:35]}}}

, в которых cor - это фрейм данных, включающий 1904 пациента и 38 столбцов (35 генов + «лимфа», «npi» и «stage»); «лимфа», «npi» и «стадия» являются названиями столбцов в cor и являются тремя клиническими признаками, т. е. числом положительных лимфатических узлов, Nottingham prognosti c index и стадией опухоли соответственно.

I ' Я хочу написать функцию, чтобы, когда я напишу что-то вроде:

computeCC(lymph)

, он покажет мне коэффициент корреляции и значение p между числом лимфатических узлов и каждым из 35 генов.

Точно так же, когда я пишу: compute CC (stage)

Это покажет мне коэффициент корреляции и p-значение между стадией опухоли и каждым из 35 генов.

Но сейчас, Я столкнулся с проблемой:

Ошибка: неожиданно '}' в: "cc1 $ P.value <- имена значений pvalues ​​(cc1) = colnames (cor) [1:35]}} "</p>

Это мои воспроизводимые данные:

cor <-  structure(list(NCOR1 = c(0.6488, 0.3312, -0.3336, 0.2663, -1.3986), ZFP36L1 = c(-1.4278, -1.9684, -1.4047, -1.1984, 0.397), SMAD4 = c(-0.5692, -2.5897, -1.4175, -2.2613, 0.6804), CDKN1B = c(-0.9829, -1.7246, -1.1409, -1.5033, -0.8475), CDH1 = c(-0.1387, 1.5924, -0.7637, 1.2737, 0.5298), PIK3R1 = c(0.2649, -0.2267, -0.6875, -0.8364, 1.3622), BRCA2 = c(0.6442, 1.2209, -0.6712, -1.0785, -0.296), KMT2C = c(-0.8759, -0.327, -0.0154, -0.7076, -0.0817), KRAS = c(0.5975, -0.0729, 0.0069, -1.3664, -0.9904), MUC16 = c(0.4375, -0.7318, -0.5569, -0.8224, -0.3882), CBFB = c(-0.9757, 0.9849, -0.9263, -1.7691, -0.7777), MAP2K4 = c(0.385, -0.6192, -1.5389, -0.1092, -2.4083), AHNAK2 = c(0.69, 0.2453, -0.0492, -1.0581, -0.2553), BAP1 = c(0.0535, -3.1571, 1.8473, -0.2338, -0.9715), ERBB2 = c(0.6171,4.4808, -0.643, 0.496, 1.1611), TP53 = c(-0.065, 1.3605, -0.0393, 1.6328, -0.3413), MAP3K1 = c(-1.241, -0.6619, -1.4874, -2.1246, 2.2862), ERBB3 = c(0.7237, -0.1072, -0.2926, -1.1115,0.5288), PTEN = c(-0.4454, -1.2554, -0.9175, -0.6936, -0.0996
), PIK3CA = c(-1.9252, -2.2674, -0.0451, -0.6883, -1.0361
), GPS2 = c(0.489, -0.363, 0.1914, -0.1519, 0.237), SF3B1 = c(1.0353, 
1.0428, 0.1198, -0.1978, 1.3932), AGTR2 = c(0.395, 1.7066, 
0.2963, 0.5277, 0.5876), SYNE1 = c(0.1814, -0.8717, -0.3925, 
-0.6181, 0.2515), GATA3 = c(0.727, -0.1693, 0.1266, 0.2467, 
0.7005), AKT1 = c(0.7579, 1.9675, -1.0293, -1.1985, -1.902
), FOXO3 = c(-0.1501, 0.0589, -0.3752, -0.4585, -0.8405), 
ARID1A = c(0.7732, -0.695, 0.0034, -0.9322, 0.5824), RB1 = c(-0.135, 
-0.6994, 0.487, 1.7919, 0.9048), CDKN2A = c(0.0647, 0.1072, 
-0.3117, -0.2668, -0.6555), MEN1 = c(-0.5376, 2.164, 1.2287, 
0.5037, 0.7852), NF1 = c(-0.5943, -0.2639, -0.8211, 0.2209, 
1.5184), TBX3 = c(-0.765, -0.2696, 0.1784, 0.6917, 0.3603
), CHEK2 = c(-0.5534, 1.8462, -0.8928, 0.7362, -0.3503), 
RUNX1 = c(-0.8007, -1.9473, 0.6226, -0.6965, 0.1434), lymph = c(1, 
5, 8, 1, 0), npi = c(4.036, 6.032, 6.03, 5.042, 3.046), stage = c(2, 
2, 3, 2, 2)), row.names = c("MB-0362", "MB-0346", "MB-0386", "MB-0574", "MB-0503"), class = "data.frame")

Может кто-нибудь предложить мне идею? Заранее спасибо.

Ответы [ 3 ]

1 голос
/ 24 февраля 2020

Попробуйте это. Ваш код имеет ряд проблем. Я минимально изменил его, чтобы (я думаю) получить то, что вы хотели.

computeCC = function(data, var) # Pass the data and a variable
{
  var <- eval(substitute(var), data, parent.frame()) # This may confuse you

  for (i in 1:35) {
    cc <- cor.test(data[,i], var, method = "spearman")
    estimates[i] = cc$estimate
    pvalues[i] = cc$p.value
  }

# These belong outside the loop
  cc1$Estimate <- estimates
  cc1$P.value <- pvalues
  rownames(cc1) = colnames(cor)[1:35]
  cc1
}

Затем вы вызываете его и сохраняете результаты:

cc2 <- computeCC(cor, lymph)

И затем посмотрите на результаты в cc2.

Существуют другие изменения, которые могут быть сделаны для улучшения кода, но по одному шагу за раз.


Данные, предоставленные вами:

cc1 <- data.frame(Estimate=paste("Site", 1:35), P.value="")
estimates = numeric(35)
pvalues = numeric(35)
1 голос
/ 24 февраля 2020

Я пытался сохранить структуру функции (которая, безусловно, может быть улучшена) до того, что вы пытались написать:

computeCC = function(x)
{
    cc1 <- data.frame(name=paste("Site", 1:35),Estimate=NA ,P.value=NA)
    estimates = numeric(35)
    pvalues = numeric(35)
    for (i in c(1:35)) {
      cc=cor.test(cor[,i],cor[,x])  # mention the method explicitly if you want to
      cc1$Estimate[i]=cc$estimate
      cc1$P.value[i]=cc$p.value

    }
     rownames(cc1) = colnames(cor)[1:35]
     cc1
  }


cc1=computeCC("lymph")
cc2=computeCC("npi")

> cc1
           name     Estimate   P.value
NCOR1    Site 1  0.123842113 0.8427233
ZFP36L1  Site 2 -0.568357574 0.3174529
SMAD4    Site 3 -0.480725033 0.4123901
CDKN1B   Site 4 -0.330612236 0.5868509
CDH1     Site 5 -0.339943380 0.5756579
PIK3R1   Site 6 -0.568101683 0.3177210
BRCA2    Site 7  0.065458611 0.9167151
KMT2C    Site 8  0.521021412 0.3679883
KRAS     Site 9  0.406528583 0.4970250
MUC16   Site 10 -0.338787942 0.5770417
CBFB    Site 11  0.375126137 0.5338257
MAP2K4  Site 12 -0.148551132 0.8115568
AHNAK2  Site 13  0.198285499 0.7491993
BAP1    Site 14  0.251796831 0.6828230
ERBB2   Site 15  0.001412623 0.9982014
TP53    Site 16  0.033297857 0.9576117
MAP3K1  Site 17 -0.380753401 0.5271922
ERBB3   Site 18 -0.251814676 0.6828010
PTEN    Site 19 -0.755059691 0.1400511
PIK3CA  Site 20  0.290714428 0.6351329
GPS2    Site 21 -0.252464062 0.6820009
SF3B1   Site 22 -0.343552477 0.5713393
AGTR2   Site 23  0.165631150 0.7900801
SYNE1   Site 24 -0.536445779 0.3513193
GATA3   Site 25 -0.719213282 0.1708848
AKT1    Site 26  0.248607276 0.6867549
FOXO3   Site 27  0.430492324 0.4693150
ARID1A  Site 28 -0.273943198 0.6556177
RB1     Site 29 -0.384189957 0.5231494
CDKN2A  Site 30  0.242973181 0.6937084
MEN1    Site 31  0.609644017 0.2749784
NF1     Site 32 -0.670039102 0.2159080
TBX3    Site 33 -0.075477911 0.9039899
CHEK2   Site 34 -0.005660548 0.9927928
RUNX1   Site 35  0.133233884 0.8308646
> cc2
           name       Estimate     P.value
NCOR1    Site 1  0.48465617704 0.408006568
ZFP36L1  Site 2 -0.82629260852 0.084607455
SMAD4    Site 3 -0.87281518266 0.053397753
CDKN1B   Site 4 -0.75022777439 0.144101826
CDH1     Site 5  0.08194436836 0.895782074
PIK3R1   Site 6 -0.84986938818 0.068235062
BRCA2    Site 7  0.09159854348 0.883536407
KMT2C    Site 8  0.14929428036 0.810621134
KRAS     Site 9  0.22613893825 0.714544197
MUC16   Site 10 -0.52100073451 0.368010756
CBFB    Site 11  0.35773390679 0.554429546
MAP2K4  Site 12  0.24169010029 0.695293378
AHNAK2  Site 13 -0.02307876193 0.970617816
BAP1    Site 14  0.00857051167 0.989087819
ERBB2   Site 15  0.21025246956 0.734283873
TP53    Site 16  0.54283256599 0.344473129
MAP3K1  Site 17 -0.68263449805 0.204095599
ERBB3   Site 18 -0.58993798568 0.295053985
PTEN    Site 19 -0.96008945668 0.009513672
PIK3CA  Site 20  0.10502883527 0.866519400
GPS2    Site 21 -0.60150437258 0.283225721
SF3B1   Site 22 -0.55953971804 0.326724402
AGTR2   Site 23  0.38051911497 0.527468083
SYNE1   Site 24 -0.87191165633 0.053960137
GATA3   Site 25 -0.91961553798 0.027026199
AKT1    Site 26  0.44418297389 0.453639104
FOXO3   Site 27  0.65557871536 0.229694002
ARID1A  Site 28 -0.68430594947 0.202542090
RB1     Site 29 -0.28148382493 0.646394380
CDKN2A  Site 30  0.50935264641 0.380721870
MEN1    Site 31  0.61847858859 0.266100528
NF1     Site 32 -0.72417687328 0.166510176
TBX3    Site 33 -0.00003856786 0.999950894
CHEK2   Site 34  0.40475255738 0.499091916
RUNX1   Site 35 -0.26198218509 0.670289901
>
0 голосов
/ 24 февраля 2020

Как уже упоминали другие, без минимального воспроизводимого примера будет трудно сказать, что пошло не так ... но вот мои пять центов:

if (x = "") Имеет две проблемы : Используйте == для проверки на равенство, и если я правильно понимаю ваше объяснение, x является вектором, и это выражение будет только первым элементом x. или, может быть, x должно быть именем столбца, но тогда предложение if не имеет смысла ... MRE действительно поможет вам помочь!

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...