2012-07-31 4 views
2

나는 ggplot2를 사용하여 4 개의 성장 곡선 그래프를 만들었습니다.성장 곡선의 최대 기울기 찾기

누구든지 시도하려는 경우 아래 코드는 그래프를 생성해야합니다.

각 선의 최대 경사 값을 4 개의 시간 지점에서 가져오고 싶습니다.

아무에게도이 문제를 어떻게 해결할 수 있습니까?

library(ggplot2) 
dat <- structure(list(TIME = c(0L, 2L, 4L, 6L, 8L, 10L, 12L, 14L, 16L, 
          18L, 20L, 22L, 24L, 26L, 28L, 30L, 0L, 2L, 4L, 6L, 8L, 10L, 12L, 
          14L, 16L, 18L, 20L, 22L, 24L, 26L, 28L, 30L, 0L, 2L, 4L, 6L, 
          8L, 10L, 12L, 14L, 16L, 18L, 20L, 22L, 24L, 26L, 28L, 30L, 0L, 
          2L, 4L, 6L, 8L, 10L, 12L, 14L, 16L, 18L, 20L, 22L, 24L, 26L, 
          28L, 30L), OD600 = c(0.2202, 0.2177, 0.2199, 0.2471, 0.2834, 
               0.357, 0.4734, 0.647, 0.898, 1.1959, 1.3765, 1.3978, 1.3948, 
               1.3928, 1.3961, 1.4018, 0.24, 0.2317, 0.2328, 0.2522, 0.2748, 
               0.3257, 0.4098, 0.5455, 0.7387, 0.9904, 1.2516, 1.3711, 1.3713, 
               1.3703, 1.3686, 1.3761, 0.2266, 0.2219, 0.2245, 0.2401, 0.2506, 
               0.2645, 0.3018, 0.3484, 0.4216, 0.5197, 0.666, 0.872, 1.1181, 
               1.2744, 1.3079, 1.2949, 0.2389, 0.2242, 0.2315, 0.2364, 0.2372, 
               0.2373, 0.2306, 0.2385, 0.236, 0.2331, 0.2379, 0.2334, 0.2336, 
               0.2339, 0.2389, 0.2349), MMS = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                       0, 0, 0, 0, 0, 0, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 
                       0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 0.005, 
                       0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 
                       0.01, 0.01, 0.01, 0.01, 0.01, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 
                       0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02)), .Names = c("TIME", 
                                          "OD600", "MMS"), class = "data.frame", row.names = c(NA, -64L 
                                                       )) 
graph = ggplot(data=dat, aes(x=TIME, y=OD600)) 
graph + geom_line(aes(colour=factor(MMS)), alpha=1) + 
opts(title="Log growth curves: change in cell density with increasing concentrations of MMS")+ 
scale_y_log10() 

많은 감사

이 같은
+0

큰 문제없이. 기본 솔루션과 ggplot 솔루션을 모두 얻을 수있어서 기쁩니다. 나는 태그 : [그라디언트]가 이미 그래픽 군중에 의해 촬영 되었기 때문에 [편집]에서 태그를 고소했다. –

답변

3

뭔가?

cbind(
    MMS = unique(dat$MMS), 
    do.call(
    rbind, 
    lapply(
     unique(dat$MMS), 
     function(x) { 
     tdat <- dat[dat$MMS == x, ] 
     response <- tdat$OD600 
     timepoints <- tdat$TIME 
     rise <- (response[4:length(response)] - response[1:(length(response) - 3)]) 
     run <- (timepoints[4:length(timepoints)] - timepoints[1:(length(timepoints) - 3)]) 
     slopes <- c(rep(NA, 3), rise/run) 
     return(
      list(
      max_slope = max(slopes, na.rm = T), 
      time = timepoints[which(slopes == max(slopes, na.rm = T)) - 3] 
     ) 
     ) 
     } 
    ) 
) 
) 

을 제공합니다 :

 MMS max_slope time 
[1,] 0  0.1215833 14 
[2,] 0.005 0.1176833 14 
[3,] 0.01 0.1014  20 
[4,] 0.02 0.002166667 2 
+0

이것은 잘 작동하는 것 같습니다. 고마워요. y 축은 대수이므로 출력 값은 플롯에 비해 약간 떨어져 보입니다. 데이터를 로그 OD600 값으로 변환하면 수동으로 계산 한 값과 매우 유사한 최대 기울기 값을 얻을 수 있습니다. – ds9000

4

당신은, @lockedoff의 솔루션은 괜찮 보간이 필요하지 않은 경우,하지만 당신은 당신이 모두 초기 농도 14 하시겠습니까?

더 나은 값을 얻으려면 기울기 시간, 즉 2 차 미분 값이 0 인 위치를 찾아야합니다. 실제 데이터에서는이 작업이 까다로울 수 있으므로 가능한 경우 파생 상품을 먼저 플롯해야합니다.

0.02의 농도는 절망적이며 이것이 내 실험 인 경우 실험실로 돌아가 실제로 이것이 0.02 또는 0.2인지 확인합니다. 그렇지 않다면, 당신은 매우 특이한 물질을 가지고 있습니다, 조심하십시오, 검토자가 좋은 설명없이 그것을 되돌려 보낼 것입니다.

미분을 계산하려면 predict.smooth.spline을 사용하고, 기울기 == 0 인 지점을 찾으려면 uniroot을 사용하십시오.

library(plyr) 
smoothingDf = 8 # Adujst this. Larger values-> Smoother curves 
# Check smoothing of second derivatives 
deriv2 = ddply(dat,.(MMS),function(x){ 
    data.frame(predict(smooth.spline(x$TIME,x$OD600,df=smoothingDf),0:max(x$TIME),2)) 
}) 
ggplot(data=deriv2, aes(x=x, y=y))+ geom_line(aes(colour=factor(MMS))) 
# No chance to get a good value for 0.02, remove it 
dat1 = dat[dat$MMS != 0.02,] 

ld50 = ddply(dat1,.(MMS),function(x){ 
    sp = smooth.spline(x$TIME, x$OD600, df=smoothingDf) 
    # Try to find a good initial range 
    app = predict(sp,min(x$TIME):max(x$TIME),2) 
    lower = app$x[which.max(app$y)] 
    upper = app$x[which.min(app$y)] 
    uniroot(function(t) predict(sp,t,2)$y ,lower=lower,upper=upper)$root 
}) 

결과 확인 보이지만, 0.02

MMS  V1 
1 0.000 16.23093 
2 0.005 17.43714 
3 0.010 22.29317 

Second derivatives. Note that 0.02 is not useful

+0

고맙습니다.이 같은 속도가 가장 높은 시간을 찾을 수있는 것은 매우 유용합니다. 다른 의견을 보내 주시면 감사하겠습니다. 0.02 % 농도에서이 효과를 보는 것이 너무 드문 이유가 궁금합니다. – ds9000

+0

나는 약리학에서 일하고 있는데, 2의 인자 내에서 그러한 가파른 억압을 보게되는 것은 매우 드문 일이며, 전에는 아주 순탄 한 전환이었습니다. 그것은 정확할 수 있지만 출판을 위해서는 0.015 정도의 중간 농도가 필요할 것입니다. 질문은 SO에 대한 답변으로 표시되어야합니다. –