if(!require(pacman)) install.packages("pacman", dependencies=TRUE)
## Loading required package: pacman
## Warning: package 'pacman' was built under R version 4.2.3
library("pacman")
p_load("vroom",
       "ggplot2",
       "ggrepel",
       "dplyr",
       "tidyverse")
Datos <- vroom(file = "https://raw.githubusercontent.com/ManuelLaraMVZ/resultados_PCR_practica/refs/heads/main/Cts.csv")
## Rows: 29 Columns: 5
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): Well, Grupo, Practica, Fluor
## dbl (1): Cq
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Datos
## # A tibble: 29 × 5
##    Well  Grupo Practica Fluor    Cq
##    <chr> <chr> <chr>    <chr> <dbl>
##  1 A01   G1    Relativa SYBR   30.8
##  2 B01   G1    Relativa SYBR   41  
##  3 C01   G1    Relativa SYBR   26.1
##  4 D01   G1    Relativa SYBR   41  
##  5 E01   G1    Relativa SYBR   41  
##  6 F01   G1    Relativa SYBR   22.3
##  7 A01   G2    Relativa SYBR   30.7
##  8 B01   G2    Relativa SYBR   20.3
##  9 C01   G2    Relativa SYBR   18.7
## 10 D01   G2    Relativa SYBR   20.5
## # … with 19 more rows

Seleccion de datos de la curva

datos_rt_pcr <- Datos %>% 
  filter(Grupo == "Curva") %>% 
  select(1,5) %>%
  mutate(Well = as.numeric(Well),
         Curva = Well,
         logConc = log10(Curva)) %>% 
  select(3,4,2)

datos_rt_pcr
## # A tibble: 5 × 3
##    Curva logConc    Cq
##    <dbl>   <dbl> <dbl>
## 1 100          2  11.2
## 2  10          1  15.2
## 3   1          0  17.6
## 4   0.1       -1  20.1
## 5   0.01      -2  23.9
grafica <- ggplot(mapping = aes(x = datos_rt_pcr$logConc, y = datos_rt_pcr$Cq)) +
  geom_point(color = "#900C3F", size = 5) +
  theme_classic()

grafica

REalizar cálculos para el ajuste lineal

# Ajustar modelo de regresión lineal
modelo_pcr <- lm(Cq ~ logConc, data = datos_rt_pcr)

modelo_pcr
## 
## Call:
## lm(formula = Cq ~ logConc, data = datos_rt_pcr)
## 
## Coefficients:
## (Intercept)      logConc  
##      17.608       -3.024
# Obtener coeficientes de la ecuación
coeficientes <- coef(modelo_pcr)
pendiente <- round(coeficientes[2], 3)
intercepto <- round(coeficientes[1], 3)
coeficientes
## (Intercept)     logConc 
##   17.608002   -3.024013
pendiente
## logConc 
##  -3.024
intercepto
## (Intercept) 
##      17.608
# Crear puntos de predicción para la línea ajustada
predicciones <- data.frame(logConc = seq(min(datos_rt_pcr$logConc), max(datos_rt_pcr$logConc), length.out = 100))
predicciones$Cq <- predict(modelo_pcr, newdata = predicciones)
predicciones
##         logConc       Cq
## 1   -2.00000000 23.65603
## 2   -1.95959596 23.53385
## 3   -1.91919192 23.41166
## 4   -1.87878788 23.28948
## 5   -1.83838384 23.16730
## 6   -1.79797980 23.04512
## 7   -1.75757576 22.92293
## 8   -1.71717172 22.80075
## 9   -1.67676768 22.67857
## 10  -1.63636364 22.55639
## 11  -1.59595960 22.43420
## 12  -1.55555556 22.31202
## 13  -1.51515152 22.18984
## 14  -1.47474747 22.06766
## 15  -1.43434343 21.94548
## 16  -1.39393939 21.82329
## 17  -1.35353535 21.70111
## 18  -1.31313131 21.57893
## 19  -1.27272727 21.45675
## 20  -1.23232323 21.33456
## 21  -1.19191919 21.21238
## 22  -1.15151515 21.09020
## 23  -1.11111111 20.96802
## 24  -1.07070707 20.84583
## 25  -1.03030303 20.72365
## 26  -0.98989899 20.60147
## 27  -0.94949495 20.47929
## 28  -0.90909091 20.35710
## 29  -0.86868687 20.23492
## 30  -0.82828283 20.11274
## 31  -0.78787879 19.99056
## 32  -0.74747475 19.86838
## 33  -0.70707071 19.74619
## 34  -0.66666667 19.62401
## 35  -0.62626263 19.50183
## 36  -0.58585859 19.37965
## 37  -0.54545455 19.25746
## 38  -0.50505051 19.13528
## 39  -0.46464646 19.01310
## 40  -0.42424242 18.89092
## 41  -0.38383838 18.76873
## 42  -0.34343434 18.64655
## 43  -0.30303030 18.52437
## 44  -0.26262626 18.40219
## 45  -0.22222222 18.28000
## 46  -0.18181818 18.15782
## 47  -0.14141414 18.03564
## 48  -0.10101010 17.91346
## 49  -0.06060606 17.79128
## 50  -0.02020202 17.66909
## 51   0.02020202 17.54691
## 52   0.06060606 17.42473
## 53   0.10101010 17.30255
## 54   0.14141414 17.18036
## 55   0.18181818 17.05818
## 56   0.22222222 16.93600
## 57   0.26262626 16.81382
## 58   0.30303030 16.69163
## 59   0.34343434 16.56945
## 60   0.38383838 16.44727
## 61   0.42424242 16.32509
## 62   0.46464646 16.20290
## 63   0.50505051 16.08072
## 64   0.54545455 15.95854
## 65   0.58585859 15.83636
## 66   0.62626263 15.71418
## 67   0.66666667 15.59199
## 68   0.70707071 15.46981
## 69   0.74747475 15.34763
## 70   0.78787879 15.22545
## 71   0.82828283 15.10326
## 72   0.86868687 14.98108
## 73   0.90909091 14.85890
## 74   0.94949495 14.73672
## 75   0.98989899 14.61453
## 76   1.03030303 14.49235
## 77   1.07070707 14.37017
## 78   1.11111111 14.24799
## 79   1.15151515 14.12580
## 80   1.19191919 14.00362
## 81   1.23232323 13.88144
## 82   1.27272727 13.75926
## 83   1.31313131 13.63708
## 84   1.35353535 13.51489
## 85   1.39393939 13.39271
## 86   1.43434343 13.27053
## 87   1.47474747 13.14835
## 88   1.51515152 13.02616
## 89   1.55555556 12.90398
## 90   1.59595960 12.78180
## 91   1.63636364 12.65962
## 92   1.67676768 12.53743
## 93   1.71717172 12.41525
## 94   1.75757576 12.29307
## 95   1.79797980 12.17089
## 96   1.83838384 12.04870
## 97   1.87878788 11.92652
## 98   1.91919192 11.80434
## 99   1.95959596 11.68216
## 100  2.00000000 11.55998

Obtener el valor de ajute (R^2)

# Obtener resumen del modelo
resumen_modelo <- summary(modelo_pcr)
resumen_modelo
## 
## Call:
## lm(formula = Cq ~ logConc, data = datos_rt_pcr)
## 
## Residuals:
##         1         2         3         4         5 
## -0.328525  0.611422 -0.000825 -0.518515  0.236443 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  17.6080     0.2319   75.94 5.03e-06 ***
## logConc      -3.0240     0.1640  -18.44 0.000348 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5185 on 3 degrees of freedom
## Multiple R-squared:  0.9913, Adjusted R-squared:  0.9883 
## F-statistic: 340.1 on 1 and 3 DF,  p-value: 0.0003478

Extraer la ecuación de la recta y el valor de R^2

# Extraer el valor de R²
R2 <- resumen_modelo$r.squared
# Crear la ecuación en formato de texto

ecuacion_texto <- paste0("y = ", pendiente, "x + ", intercepto, "\nR² = ", round(R2, 4))

# Mostrar el valor de R²
cat("y = ", pendiente, "x + ", intercepto, "\nR² = ", round(R2, 4),"\n")
## y =  -3.024 x +  17.608 
## R² =  0.9913

Graficación

grafica2 <- grafica +  # Puntos de datos reales
  geom_smooth(data = predicciones, aes(x = logConc, y = Cq), color = "#FFC300", size = 1.5) +  # Ajuste lineal
  labs(
    title = "Curva Estándar de RT-PCR en Tiempo Real",
    subtitle = "Simulación de datos de RT-PCR",
    caption = "Diseñó: Manuel Lara",
    x = "Log10(Concentración de ADN [pg/µL])",
    y = "Cq (Ct)") +
  theme_classic(base_size = 15) +
  scale_y_continuous(labels = scales::number_format(accuracy = 1),) +  # Redondeo de 3 dígitos en Y
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    axis.title.x = element_text(face = "bold"),
    axis.title.y = element_text(face = "bold")
  ) +
  #geom_smooth(method = "lm", se = FALSE, color="#FFC300") +  # Línea de regresión
  annotate("text", 
           x = max(datos_rt_pcr$logConc) -1.5,  # Posición hacia la derecha
           y = min(datos_rt_pcr$Cq) + 12,         # Ajustar para que no se sobreponga
           label = ecuacion_texto, 
           size = 5, color = "#581845", fontface = "bold", hjust = 0)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
grafica2
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Análisis de las muestras:

  1. Determinar los valores del equipo y del profesor
Datos
## # A tibble: 29 × 5
##    Well  Grupo Practica Fluor    Cq
##    <chr> <chr> <chr>    <chr> <dbl>
##  1 A01   G1    Relativa SYBR   30.8
##  2 B01   G1    Relativa SYBR   41  
##  3 C01   G1    Relativa SYBR   26.1
##  4 D01   G1    Relativa SYBR   41  
##  5 E01   G1    Relativa SYBR   41  
##  6 F01   G1    Relativa SYBR   22.3
##  7 A01   G2    Relativa SYBR   30.7
##  8 B01   G2    Relativa SYBR   20.3
##  9 C01   G2    Relativa SYBR   18.7
## 10 D01   G2    Relativa SYBR   20.5
## # … with 19 more rows
Analisis <- Datos %>% 
  filter(Practica == "Absoluta",
         Well == "Profesor", 
         Grupo == "G1") %>% 
  select(1,2,5)
Analisis
## # A tibble: 1 × 3
##   Well     Grupo    Cq
##   <chr>    <chr> <dbl>
## 1 Profesor G1     18.6

Colocar los Ct obtenidos del profesor

# Ct (Cq) observado
Ct_Control <- Analisis$Cq
Ct_Control
## [1] 18.57611
# Calcular log10(Concentración)
logConc_calculado_Cx <- (Ct_Control - intercepto) / pendiente
logConc_calculado_Cx
## (Intercept) 
##  -0.3201429
# Obtener concentración real
Conc_calculada_Cx <- 10^logConc_calculado_Cx
Conc_calculada_Cx
## (Intercept) 
##   0.4784726
# Crear un data frame con el nuevo punto
nuevo_punto_Cx <- data.frame(logConc = logConc_calculado_Cx, Cq = Ct_Control)
nuevo_punto_Cx
##                logConc       Cq
## (Intercept) -0.3201429 18.57611
# Mostrar resultado
cat("La concentración estimada para la muestra del profesor es:", round(Conc_calculada_Cx, 4), "pg/µL\n")
## La concentración estimada para la muestra del profesor es: 0.4785 pg/µL

Graficamos el punto en la recta

interpolacion <- grafica2+
  geom_point(data = nuevo_punto_Cx, 
             aes(x = logConc, y = Cq), 
             color = "black", 
             size = 6, shape = 20) +  # Nuevo punto
  geom_segment(aes(x = min(datos_rt_pcr$logConc), 
                   xend = logConc_calculado_Cx, 
                   y = Ct_Control, 
                   yend = Ct_Control), 
               linetype = "dotted", 
               color = "#2f2d2c", 
               size = 1.5) +  # Línea horizontal hasta el punto
  geom_segment(aes(x = logConc_calculado_Cx, xend = logConc_calculado_Cx, 
                   y = min(datos_rt_pcr$Cq), yend = Ct_Control),
               linetype = "dotted", 
               color = "#2f2d2c", 
               size = 1.5)+
   annotate("text", 
           x = logConc_calculado_Cx, 
           y = Ct_Control - 1,  
           label = paste0("  [Prof]= ", 
                          round(Conc_calculada_Cx, 4), " pg/µL"), 
           color = "black", size = 5, 
           fontface = "bold", 
           hjust = -0.007, vjust = -1.7)

interpolacion
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Datos
## # A tibble: 29 × 5
##    Well  Grupo Practica Fluor    Cq
##    <chr> <chr> <chr>    <chr> <dbl>
##  1 A01   G1    Relativa SYBR   30.8
##  2 B01   G1    Relativa SYBR   41  
##  3 C01   G1    Relativa SYBR   26.1
##  4 D01   G1    Relativa SYBR   41  
##  5 E01   G1    Relativa SYBR   41  
##  6 F01   G1    Relativa SYBR   22.3
##  7 A01   G2    Relativa SYBR   30.7
##  8 B01   G2    Relativa SYBR   20.3
##  9 C01   G2    Relativa SYBR   18.7
## 10 D01   G2    Relativa SYBR   20.5
## # … with 19 more rows
Analisis2 <- Datos %>% 
  filter(Practica == "Absoluta",
         Well == "Profesor", 
         Grupo == "G2") %>% 
  select(1,2,5)
Analisis2
## # A tibble: 1 × 3
##   Well     Grupo    Cq
##   <chr>    <chr> <dbl>
## 1 Profesor G2     20.9

Colocar los Ct obtenidos del equipo

# Ct (Cq) observado
Ct_equipo<- Analisis2$Cq
Ct_equipo
## [1] 20.90737
# Calcular log10(Concentración)
logConc_calculado_equipo <- (Ct_equipo - intercepto) / pendiente
logConc_calculado_equipo
## (Intercept) 
##   -1.091061
# Obtener concentración real
Conc_calculada_equipo <- 10^logConc_calculado_equipo
Conc_calculada_equipo
## (Intercept) 
##  0.08108467
# Crear un data frame con el nuevo punto
nuevo_punto_equipo <- data.frame(logConc = logConc_calculado_equipo, Cq = Ct_equipo)
nuevo_punto_equipo
##               logConc       Cq
## (Intercept) -1.091061 20.90737
# Mostrar resultado
cat("La concentración estimada para la muestra del equipo es:", round(Conc_calculada_equipo, 4), "pg/µL\n")
## La concentración estimada para la muestra del equipo es: 0.0811 pg/µL

Interpolación en la gráfica

interpolacion2 <- interpolacion+
    geom_point(data = nuevo_punto_equipo, 
             aes(x = logConc, y = Cq), 
             color = "black", 
             size = 6, shape = 20) +  # Nuevo punto en verde
  geom_segment(aes(x = min(datos_rt_pcr$logConc), 
                   xend = logConc_calculado_equipo, 
                   y = Ct_equipo, 
                   yend = Ct_equipo), 
               linetype = "dotted", 
               color = "#2f2d2c", 
               size = 1.5) +  # Línea horizontal hasta el punto
  geom_segment(aes(x = logConc_calculado_equipo, xend = logConc_calculado_equipo, 
                   y = min(datos_rt_pcr$Cq), yend = Ct_equipo),
               linetype = "dotted", 
               color = "#2f2d2c", 
               size = 1.5)+
   annotate("text", 
           x = logConc_calculado_equipo, 
           y = Ct_equipo - 1,  
           label = paste0("  [Equipo]= ", 
                          round(Conc_calculada_equipo, 4), " pg/µL"), 
           color = "black", size = 5, 
           fontface = "bold", 
           hjust = -0.007, vjust = -1.7)
interpolacion2
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Comparación con respecto al profesor

Comparacion <- data.frame(Condicion=c("Profesor", "Equipo"),
           Concentracion= c(Conc_calculada_Cx, Conc_calculada_equipo))
Comparacion
##   Condicion Concentracion
## 1  Profesor    0.47847260
## 2    Equipo    0.08108467
# Calcular el máximo
max_conc <- max(Comparacion$Concentracion, na.rm = TRUE)

# Calcular incremento, asegurando que no sea 0
break_by <- round(max_conc / 10, 1)
break_by <- ifelse(break_by == 0, 0.1, break_by)

# Crear la gráfica
Barras <- ggplot(Comparacion, aes(x = Condicion, y = Concentracion, fill = Condicion)) +
  geom_bar(stat = "identity", color = "black", width = 0.9) +
  labs(title = "Comparación de Concentraciones",
       x = "Condición",
       y = "Concentración (pg/µL)") +
  theme_classic(base_size = 14) +
  scale_fill_manual(values = c("Profesor" = "#f8c471", "Equipo" = "#2980b9")) +
  scale_y_continuous(limits = c(0, max_conc+ 0.1), 
                     breaks = seq(0, max_conc + 1, by = 0.1)) +
  theme(
    axis.title = element_text(face = "bold", size = 16),
    axis.text = element_text(face = "bold", size = 14),
    axis.line = element_line(size = 1.5, color = "black"),
    legend.position = "none")+
  geom_text(aes(label = paste0(round(Concentracion, 2), " pg/µL")),  
            vjust = -0.5, size = 5) 
## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
Barras

Porcentaje <- Comparacion %>%
  mutate(Porcentaje = (Concentracion ) *100/ Concentracion[Condicion == "Profesor"]) %>% 
  select(1,3)
  
Porcentaje
##   Condicion Porcentaje
## 1  Profesor  100.00000
## 2    Equipo   16.94656
# Calcular el máximo valor
max_porcentaje <- max(Porcentaje$Porcentaje, na.rm = TRUE)

# Calcular incremento, asegurando que no sea 0
break_by <- round(max_porcentaje / 10, 0)
break_by <- ifelse(break_by == 0, 1, break_by)  # Asegura que by no sea 0

# Crear la gráfica
Barras2 <- ggplot(Porcentaje, aes(x = Condicion, y = Porcentaje, fill = Condicion)) +
  geom_bar(stat = "identity", color = "black", width = 0.9) +
  labs(
    title = "Comparación de Concentraciones",
    x = "Condición",
    y = "Tasa de cambio (u. a.)"
  ) +
  theme_classic(base_size = 14) +
  scale_fill_manual(values = c("Profesor" = "#f8c471", 
                               "Equipo" = "#2980b9")) +
  scale_y_continuous(limits = c(0, max_porcentaje +20), 
                     breaks = seq(0, max_porcentaje + 20,  by = 10)) +
  theme(axis.title = element_text(face = "bold", size = 16),
        axis.text = element_text(face = "bold", size = 14),
        axis.line = element_line(size = 1.5, color = "black"),
        legend.position = "none") +
  geom_text(aes(label = paste0(round(Porcentaje, 2), " %")),  # <- Aquí está la etiqueta con unidad
            vjust = -0.5, size = 5) 

# Mostrar la gráfica
Barras2