PCR analysis

1. Indroducción

A. Generalidades

Como comentamos en clase, podemos dividir la técnica en dos grandes tipos (no son los únicos, pero por cuestiones didacticas los podemos agrupar así):

qRT-PCR: PCR en tiempo real cuantitativo.
endpoint PCR: PCR punto final.

La mayor ventaja del primero es que nos da datos cuantitativos de la concentración de nuestro amplificado de interes.

Amplificado: es un segmento de la secuencia que es el objetivo de la amplificación.

Estos datos cuantitavos puden ser relativos, si solo tomamos la intesidad de fluorescencia como parámetro, o pueden ser absolutos si tomamos una escala de referencia (una concentración de cDNA conocida a diferentes diluciones, lo que nos permite hacer una curva y aplicar una regresión lineal para saber la concentración de nuestra muestra en ng/mL).

La segunda técnica nos va servir especialmente para saber si está presente o no está presente nuestro gen de interes (del cual, seleccionamos una pequeña porción para amplificar). Un uso común para esta técnica es la genotipificación.

Genotipificar: es tomar el genoma de un animal (usualmente transgenico) y amplificar un gen de interes (usualemnte el transgen). Si el gen insertado está presente en las células del animal/bacteria/levadura entonces decimos que es positivo para el transgen (decimo que está en la célula y no insertado en el genoma, porque algunos genes cuando los transfectmos/infectamos en una celula, este puede quedar separa del genoma del animal y estar presente como un cromosoma extra o como un “cassete”)

2. Fundamentos del PCR

Cuando diseñamos nuestros experimentos, lo primero que tenemos que elegir es un par de primer adecuados para identificar nuestro gen. Estos primers que usualmente son secuencias de 20 bases nitrogenadas, preferentemente deben ser ricos en GC en la parte 3’ para obtener una unión mas firme entre el primer y el template (DNA/RNA objetivo del que se va a amplificar nuestra secuencia). Existen otras reglas o recomenciones al momento de diseñar el primer, pero no son el objetivo de este capítulo.

Una vez que definimos los primers adecuados y agregamos la DNApolimerasa con los nuecleotidos para comenzar el proceso de transcripción, es importante entender, que usualmente amplificamos DNA bicatenario. Y de cada hebra original obtenemos una hebra hija, es decir que comenzamos con 2 hebras, y terminamos el proceso con 4 hebras. A esto llamamos un ciclo. Para el siguiente ciclo obtenemos 8, y para el siguiente 16, y así sucesivamente. Podemos modelar este crecimiento de la siguiente manera:

$$ [Amplificado]_{total} = 2^n \

$$ Donde ‘n’ es el número de ciclo. Esto solo ocurre cuando la reacción es perfecta, es decir, tiene una eficiencia del 100%. Pero no es así siempre (de hecho es lo menos usual). La eficiencia puede variar por la calidad de nuestra muestra, con concentración de magnesio en la solución, la eficiencia de unión de nuestros primers, la presencia de contaminantes en nuestra solución como alcholes. Entonces nuestra eficiencia usualmente es menor 2. Basados en una muestra de referencia o usando diferentes diluciones de nuestra muestra problema podemos determinar la eficiencia de nuestra reacción en nuestras condiciones. Existen algunas funciones de R que nos permiten hacer estos calculos y que los veremos mas adelante.

Un detalle técnico que debemos tener en cuanta es que los equipos tienen limitaciones, en primer lugar tienen un umbral mínimo para detectar fluorescencia, por debajo de ese umbral, aunque haya amplificación y esté fluroesciendo, nuestro equipo será incapaz de recononcerlo (esto pasa en lor primeros ciclos de nuestro qRT-PCR). Por otro lado, cuando la fluroescencia es mucha, satura nuestro sensor y ya no puede registrar o darnos valores mas grandes. Entonces la señal que vemos se asemeja a una meseta. Además con mayor cantidad de amplificado la reacción va perdiendo progresivamente su eficiencia.

Veamos como se vería una gráfica de intensidad de flurescencia en un qRT-PCR

library(ggplot2)
fluore <- read.table("/Users/alfredocardenasrivera/Project/bioinformatics_2022/data/qRT-PCR/qPCR_intensities.txt")

ggplot(data = fluore, aes(x = Cycles, y = Fluorescence)) +
  geom_line() +
  geom_point() +
  labs(title = "PCR  fluorescencia", y = "Intensidad de fluroescencia (u.a.)", y = "Ciclos (# de ciclo)")

Como podrá ver los niveles de fluroescencia son indetectable hasta el ciclo número 16-17 (aunque la reacción de PCR se está dando y por eso vemos que la curva crece posteriormente).

Por otro lado a partir del ciclo 30 al 40, la curva comienza a cambiar y se termina siendo horizontal (cuando la reacción se satura o nuestro equipo se satura o ambos).

Como sabemos que la reacción de PCR se asemeja a una curva exponencial dada por $2^n$ entonces podemos tratar de linearizar la curva escalando los valores de fluorescencia por su logaritmo en base dos.

ggplot(data = fluore, aes(x = Cycles, y = Fluorescence)) +
  geom_line() +
  geom_point() +
  scale_y_continuous(trans = "log2") +
  labs(title = "PCR  fluorescencia", y = "Intensidad de fluroescencia (log2(u.a.)", y = "Ciclos (# de ciclo)")

## Warning in self$trans$transform(x): NaNs produced

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning in self$trans$transform(x): NaNs produced

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Removed 2 rows containing missing values (`geom_line()`).

## Warning: Removed 11 rows containing missing values (`geom_point()`).

Aqui es evidente que entre los ciclos 5 a 14, los valores no estan definidos en el logaritmo base 2 y por lo tanto no se grafican. Tambien, es evidente que la amplificación comienza en el ciclo 15, y es mas o menos lineal hasta el ciclo 21 o 22.

B. Archivos RDML

1. Importación de archivos RDML

Como mencionamos anteriormente existe dos formatos globales para subir y compartir los resultados de los ensayos de qRT-PCR. El formato RDML y el MIQE. El objetivo de ambos, es contener la información mínima necesaria de los parámetros de la muestra. Como, por ejemplo, la loacilización de cada muestra en los pozos, volumen utilizado de template, concentraciones, tipos de muestra, numero de ciclos, intensidad de fluorescencia por ciclo por pozo, etc.

Para entener un poco mejor este formato, vamos a cargar una ejemplo usando el paquete RDML y le paquete qpcR

Primero instalamos los paquetes

Luego cargamos las librerias en nuestro entorno de trabajo con la función library() y luego “construimos” el objeto RDML con la función RDML$new() del paquete RDML.

library(RDML)
library(qpcR)

## Loading required package: MASS

## Loading required package: minpack.lm

## Loading required package: rgl

## Loading required package: robustbase

## Loading required package: Matrix

ex.rdml <- RDML::RDML$new("/Users/alfredocardenasrivera/Project/bioinformatics_2022/data/qRT-PCR/qPCR_data_naive_18S_IL2.rdml")

## 
## Loading experiment: All Wells
##  run: Amp Step 4_SYBR

ex.rdml

##  dateMade: 2013-01-08T15:00:32.921+00:00
##  dateUpdated: 2015-01-15T13:54:13.750+02:00
##  id: [Bio-Rad Laboratories, Inc.]
##  experimenter: [admin]
##  documentation: []
##  dye: [SYBR]
##  sample: [koB 0 h, wt 0 h, koA 0 h, wt 4 h, koA 4 h, koB 4 h, wt 12 h, koA 12 h, koB 12 h, wt 24 h, koA 24 h, wt 48 h, koA 48 h]
##  target: [18s, IL-2]
##  thermalCyclingConditions: [IL_2.prcl]
##  experiment: [All Wells]

Con esta función acabamos de generar un objeto de la clase “RDML”

class(ex.rdml)

## [1] "RDML"         "rdmlBaseType" "R6"

Para ver una representación gráfica del experimento podemos usar la función AsDendrogram() que esta definida dentro del mimso objeto ** RDML**, y que lo podemos llamar usando el signo $

ex.rdml$AsDendrogram()

## 'dendrogram' with 1 branches and 4 members total, at height 5

A primera vista podemos entender que en este experiemnto se evaluaron la presencia de dos genes:

18S: que es la fracción 18S ribosomal, usualmente usada como control ya que su expresion es constitucional y puede funcionar como un buen housekeeping
IL-2: que es la interleucina 2, una interlucina relacionada con respuesta pro-inflamatoria.

Por otro lado, vemos que para cada uno se tomaron 39 muestras diferentes en 39 pozos distintos, unos para 18S y otros para IL-2. Las palabras adp y mdp hacen referencia a la cuantificacion de la fluorescencia para cada ciclo de amplificacion y la curva de fusión, respectivamente.

Para extraer los datos de este objeto podemos usar la función AsTable que esta definida, al igual que la anterior, dentro del objeto y que podemos llamarla de la siguiente manera

pcr.data <- ex.rdml$AsTable()
knitr::kable(head(pcr.data))

fdata.name	exp.id	run.id	react.id	position	sample	target	target.dyeId	sample.type	adp	mdp
A01_koB 0 h_unkn_18s	All Wells	Amp Step 4_SYBR	1	A01	koB 0 h	18s	SYBR	unkn	TRUE	TRUE
A02_koB 0 h_unkn_18s	All Wells	Amp Step 4_SYBR	2	A02	koB 0 h	18s	SYBR	unkn	TRUE	TRUE
A03_koB 0 h_unkn_18s	All Wells	Amp Step 4_SYBR	3	A03	koB 0 h	18s	SYBR	unkn	TRUE	TRUE
B01_wt 0 h_unkn_18s	All Wells	Amp Step 4_SYBR	13	B01	wt 0 h	18s	SYBR	unkn	TRUE	TRUE
B02_koA 0 h_unkn_18s	All Wells	Amp Step 4_SYBR	14	B02	koA 0 h	18s	SYBR	unkn	TRUE	TRUE
B03_wt 4 h_unkn_18s	All Wells	Amp Step 4_SYBR	15	B03	wt 4 h	18s	SYBR	unkn	TRUE	TRUE

str(pcr.data)

## Classes 'data.table' and 'data.frame':   78 obs. of  11 variables:
##  $ fdata.name  : chr  "A01_koB 0 h_unkn_18s" "A02_koB 0 h_unkn_18s" "A03_koB 0 h_unkn_18s" "B01_wt 0 h_unkn_18s" ...
##  $ exp.id      : chr  "All Wells" "All Wells" "All Wells" "All Wells" ...
##  $ run.id      : chr  "Amp Step 4_SYBR" "Amp Step 4_SYBR" "Amp Step 4_SYBR" "Amp Step 4_SYBR" ...
##  $ react.id    : int  1 2 3 13 14 15 16 17 18 19 ...
##  $ position    : chr  "A01" "A02" "A03" "B01" ...
##  $ sample      : chr  "koB 0 h" "koB 0 h" "koB 0 h" "wt 0 h" ...
##  $ target      : chr  "18s" "18s" "18s" "18s" ...
##  $ target.dyeId: chr  "SYBR" "SYBR" "SYBR" "SYBR" ...
##  $ sample.type : chr  "unkn" "unkn" "unkn" "unkn" ...
##  $ adp         : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  $ mdp         : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "sorted")= chr "fdata.name"

Podemos ver que es una tabla de 78 observaciones y 11 variables

fdata.name: nombre o ID de cada muestra.
exp.id : que es el ID del experimento, este experimento puede tener o ocupar varias placas de 96 pozos cada una, entonces este código se puede repetir varias veces en diferentes placas
run.id : define el ID de la corrida, esto sería como el equivalente de ID de la placa que estamos usando.
react.id: define el ID de cada muestra
position: es el codigo que determina la posición de la muestra en la placa de 96 pozos. Este codigo es alfanumerico, la letra define la fila y el número la columna.
sample : podemos considerarla como un factor para el lenguaje R, y nos define varias categorias “temporales” de nuestro experimento (luego veremos) por cada tipo de ratón.
target: hace referencia al gen que se esta amplificando, como vimos anteriormente son dos en este experimento: el control 18S y el gen objetivo IL-2.
target.dyeID: en este caso, vemos que se uso SYBR green para todas las reacciones.
sample.type: en este ejemplo no se definene, pero aqui podríamos poner si la muestra proviene de cerebro, piel, hígado, cultivo celular, etc.
adp: valor lógico para saber si contamos con los datos de flurescencia
mdp: valor lógico para saber si contamos con los datos de la curva de fusión.

Podemos explorar un poco mas la varible 6 sample para tener una mejor idea del diseño experimental.

table(pcr.data$sample)

## 
##  koA 0 h koA 12 h koA 24 h  koA 4 h koA 48 h  koB 0 h koB 12 h  koB 4 h 
##        6        6        6        6        6        6        6        6 
##   wt 0 h  wt 12 h  wt 24 h   wt 4 h  wt 48 h 
##        6        6        6        6        6

Como podemos ver, tenemos tres razas de ratón

wt: tipo silvestre
koA: knock out para un gen A (o variante A del KO)
koB: Knock out para un gen B (o variante B del mismo KO)

Además podemos ver que el experiment se realizó en distintos puntos temporales

A las 0 horas
A las 4 horas
A las 6 horas
A las 12 horas
A las 24 horas
A las 48 horas

Esto nos debería dar 18 combinaciones posibles, pero solo tenemos 13 combinaciones, eso quiere decir, que en este ensayo, va a faltarnos grupos para algunos puntos temporales o para algunas razas de ratón.

2. Extracción de los datos de fluorescencia

Para obtener los datos de fluorescencia en bruto, es necesario utilizar la función RDML::GetFData() que está definida en el paquete RDML

fluo.data <- RDML::GetFData(ex.rdml)
dim(fluo.data)

## [1] 40 79

knitr::kable(str(fluo.data))

## Classes 'data.table' and 'data.frame':   40 obs. of  79 variables:
##  $ cyc                   : num  1 2 3 4 5 6 7 8 9 10 ...
##  $ A01_koB 0 h_unkn_18s  : num  32.41 13.65 18.45 8.76 4.13 ...
##  $ A02_koB 0 h_unkn_18s  : num  28.2 21.74 12.16 6.87 3.78 ...
##  $ A03_koB 0 h_unkn_18s  : num  23.48 20.39 13.65 6.42 7.21 ...
##  $ B01_wt 0 h_unkn_18s   : num  23.87 16.93 10.98 7.48 3.14 ...
##  $ B02_koA 0 h_unkn_18s  : num  22.93 14.38 12.92 7.71 5.32 ...
##  $ B03_wt 4 h_unkn_18s   : num  15.09 11.62 7.74 10.76 8.18 ...
##  $ B04_koA 4 h_unkn_18s  : num  19.13 11.61 11.83 5.96 6.99 ...
##  $ B05_koB 4 h_unkn_18s  : num  17.3 10.03 12.82 5.69 2.66 ...
##  $ B06_wt 12 h_unkn_18s  : num  18.17 11.98 12.86 8.46 10.4 ...
##  $ B07_koA 12 h_unkn_18s : num  15 7.46 5.02 3.76 6.02 ...
##  $ B08_koB 12 h_unkn_18s : num  19.93 9.5 11.47 12.15 9.71 ...
##  $ B09_wt 24 h_unkn_18s  : num  18.72 7.45 13.15 10.42 1.96 ...
##  $ B10_koA 24 h_unkn_18s : num  5.12 3.68 9.18 3.55 9.02 ...
##  $ B11_wt 48 h_unkn_18s  : num  8.21 11.07 4.93 7.35 6.18 ...
##  $ B12_koA 48 h_unkn_18s : num  22.06 12.76 7.42 5.32 5.57 ...
##  $ C01_wt 0 h_unkn_18s   : num  28.42 17.98 11.8 6.81 5.1 ...
##  $ C02_koA 0 h_unkn_18s  : num  22.6 14.06 11.86 6.02 2.36 ...
##  $ C03_wt 4 h_unkn_18s   : num  12.53 13.26 11.46 3.24 6.53 ...
##  $ C04_koA 4 h_unkn_18s  : num  18.39 13.74 11.05 7.84 5.11 ...
##  $ C05_koB 4 h_unkn_18s  : num  14.14 5.61 12.37 7.36 5.03 ...
##  $ C06_wt 12 h_unkn_18s  : num  23.5 15.66 16.17 8.01 2.36 ...
##  $ C07_koA 12 h_unkn_18s : num  16.4 12.77 14.47 8.03 7.02 ...
##  $ C08_koB 12 h_unkn_18s : num  14.29 7.6 10.18 7.41 2.84 ...
##  $ C09_wt 24 h_unkn_18s  : num  12.78 6.67 7.93 9.1 6.95 ...
##  $ C10_koA 24 h_unkn_18s : num  2.17 8.66 6.27 5.59 3.86 ...
##  $ C11_wt 48 h_unkn_18s  : num  2.45 -2.37 6.71 7.22 3.17 ...
##  $ C12_koA 48 h_unkn_18s : num  16.44 17.48 12.71 4.34 7.63 ...
##  $ D01_wt 0 h_unkn_18s   : num  21.82 9.86 11.15 9.61 4.47 ...
##  $ D02_koA 0 h_unkn_18s  : num  7.19 8 8.04 6.32 7.68 ...
##  $ D03_wt 4 h_unkn_IL-2  : num  10.69 10.12 6.73 6.69 6.68 ...
##  $ D04_koA 4 h_unkn_IL-2 : num  7.75 4.11 6.35 5.85 7.56 ...
##  $ D05_koB 4 h_unkn_IL-2 : num  7.99 9.41 7.56 6.38 5.1 ...
##  $ D06_wt 12 h_unkn_18s  : num  11.83 10.27 14.83 6.56 3.83 ...
##  $ D07_koA 12 h_unkn_18s : num  16.9 11.2 12.2 10.5 4.1 ...
##  $ D08_koB 12 h_unkn_18s : num  17.7 13.4 12.06 8.1 7.11 ...
##  $ D09_wt 24 h_unkn_18s  : num  13 6.72 9.39 8.41 8.52 ...
##  $ D10_koA 24 h_unkn_18s : num  15.47 3.85 7.86 13.67 10.22 ...
##  $ D11_wt 48 h_unkn_18s  : num  1.38 -5.9 1.02 2.87 4.83 ...
##  $ D12_koA 48 h_unkn_18s : num  27.3 27.82 11.28 5.32 8.61 ...
##  $ E01_wt 0 h_unkn_IL-2  : num  20.88 11.26 10.29 8.17 6.8 ...
##  $ E02_koA 0 h_unkn_IL-2 : num  9.66 5.69 4.83 4.1 6.52 ...
##  $ E03_wt 4 h_unkn_18s   : num  15.19 11.15 14.27 9.28 10.63 ...
##  $ E04_koA 4 h_unkn_18s  : num  17.441 -0.163 7.512 8.548 2.369 ...
##  $ E05_koB 4 h_unkn_18s  : num  20.66 11.92 14.81 9.37 4.64 ...
##  $ E06_wt 12 h_unkn_IL-2 : num  14.05 8.07 12.97 8.18 4.56 ...
##  $ E07_koA 12 h_unkn_IL-2: num  3.35 5.38 7.5 1.53 5.49 ...
##  $ E08_koB 12 h_unkn_IL-2: num  1.84 2.32 5.47 8.12 8.08 ...
##  $ E09_wt 24 h_unkn_IL-2 : num  4.38 6.71 9.31 4.93 2.97 ...
##  $ E10_koA 24 h_unkn_IL-2: num  8.38 2.78 7.35 6.45 4.16 ...
##  $ E11_wt 48 h_unkn_IL-2 : num  -1.44 -2.17 1.27 2.76 9.87 ...
##  $ E12_koA 48 h_unkn_IL-2: num  14.17 10.59 7.14 7.11 3.69 ...
##  $ F01_wt 0 h_unkn_IL-2  : num  0.468 11.909 3.816 6.245 3.425 ...
##  $ F02_koA 0 h_unkn_IL-2 : num  12.06 5.54 8.53 3.34 3 ...
##  $ F03_wt 4 h_unkn_IL-2  : num  14.55 8.33 11.76 7.43 8.58 ...
##  $ F04_koA 4 h_unkn_IL-2 : num  7.2 3.9 8.79 3.79 7.26 ...
##  $ F05_koB 4 h_unkn_IL-2 : num  2.85 4.97 9.93 5.31 2.1 ...
##  $ F06_wt 12 h_unkn_IL-2 : num  12.71 4.73 12.33 6.71 7.39 ...
##  $ F07_koA 12 h_unkn_IL-2: num  8.13 7.19 4.61 7.15 4.46 ...
##  $ F08_koB 12 h_unkn_IL-2: num  7.18 8.08 7.82 5.6 4.78 ...
##  $ F09_wt 24 h_unkn_IL-2 : num  2.63 5.31 8.54 6.01 8.01 ...
##  $ F10_koA 24 h_unkn_IL-2: num  -3.163 0.996 2.677 4.037 7.759 ...
##  $ F11_wt 48 h_unkn_IL-2 : num  -1.46 -6.02 1.98 9.2 4.49 ...
##  $ F12_koA 48 h_unkn_IL-2: num  9.36 13.58 4.08 3.69 4.21 ...
##  $ G01_wt 0 h_unkn_IL-2  : num  12.79 6.63 5.89 5.66 5.62 ...
##  $ G02_koA 0 h_unkn_IL-2 : num  5.16 1.29 7.65 9.58 1.87 ...
##  $ G03_wt 4 h_unkn_IL-2  : num  16.46 9.75 7.24 11.48 3.75 ...
##  $ G04_koA 4 h_unkn_IL-2 : num  13.96 11.06 11 8.71 7.46 ...
##  $ G05_koB 4 h_unkn_IL-2 : num  7.52 4.95 8.62 3.52 5.69 ...
##  $ G06_wt 12 h_unkn_IL-2 : num  8.07 1.4 9.71 5.82 4.61 ...
##  $ G07_koA 12 h_unkn_IL-2: num  12.11 3.12 11.57 5.13 6.21 ...
##  $ G08_koB 12 h_unkn_IL-2: num  2.24 -0.57 8.9 3.35 5.34 ...
##  $ G09_wt 24 h_unkn_IL-2 : num  6.94 6.82 10.76 2.89 4.04 ...
##  $ G10_koA 24 h_unkn_IL-2: num  6.43 1.35 5.03 7.54 9.92 ...
##  $ G11_wt 48 h_unkn_IL-2 : num  -11.98 -11.27 1.23 2.87 4.49 ...
##  $ G12_koA 48 h_unkn_IL-2: num  3.09 7.84 7.75 2.02 6.55 ...
##  $ H01_koB 0 h_unkn_IL-2 : num  0.962 5.102 4.62 7.156 3.183 ...
##  $ H02_koB 0 h_unkn_IL-2 : num  10.87 5.38 9.65 3.65 4.05 ...
##  $ H03_koB 0 h_unkn_IL-2 : num  14.356 6.798 6.308 -0.264 0.486 ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "sorted")= chr "cyc"

El resultado es un data.frame/data.table que tiene 40 filas y 79 columnas. Las filas corresponden a los ciclos de amplificación del primero al cuatrigésimo. Por otro lado, de las columnas, la primera indica el número del ciclo y las 78 restantes corresponden a cada una de las muestras de los 78 pozos.

3. Calculo de Ct

Para poder comparar la expresión de un gen (en esta técnica de qRT-PCR), lo podemos hacer de dos formas basicamente:

De forma cuantitativa absoluta: para lo cual necesitaremos una curva de calibración del amplificado (el resultado es en ng/mL). Usualmente es usado para calcular la carga viral.
De forma cuantitativa relativa: en este caso se usa un gen de referencia y se expresa el nivel de expresión como “veces” en relación del gen de referencia, asi tenemos que el gen A se expresa 1.5 veces lo que se expresa el gen de referencia y el gen B se expresa 0.75 veces el gen de referencia. Entonces podemos concluir que el gen A se expresa dos veces mas que el gen B.

Antes de empezar con este analisis, es recomendable que Ud tenga en cuenta la diferencia entre una “Réplica técnica” y un Réplica experimental. En la primera, lo que se hace es tomar de una misma muestra dos o mas alicuatos y sometarlas al mismo proceso (en este caso el de la amplicación por PCR); al final lo que medimos es la variación entre las alicuotas que amplificamos y el error/diferencia entre ellas nos dice cuanta variación/precisión tenemmos en nuestra técnica. Por otro lado, una “réplica experimental” es tomar otro animal, cultivo celular o muestra, de un experimento nuevo sometido a las mismas condiciones para ver si nuestros resultados son reproducibles o están pasando de forma aleatoria entre nuestros experimentos. Entonces, mientras uno nos sirve para ver la precisión en nuestro proceso, el otro nos sirve para ver si nuestros hallazgos son respresentativos del proceso biológico que estamos estudiando.

Dicho esto, ahora vamos analizar un ejemplo en que se midio una misma muestra pero con diferentes diluciones.

db2 <- read.table("/Users/alfredocardenasrivera/Project/bioinformatics_2022/data/qRT-PCR/dilutions.txt") # cargamos los datos
# Exploracion rápida de la base de datos
dim(db2)                         # dimensiones de la base de datos

## [1] 40 11

knitr::kable(head(db2))          # primera 6 línas de la base de datos

Undiluted1	Undiluted2	Undiluted3	X5xs1	X5xs2	X5xs3	X25xs1	X25xs2	X25xs3	X125xs1	X125xs2
-12.943300	-36.872307	-54.8242492	17.339972	7.417517	24.362138	27.894176	17.453714	26.093649	23.307874	23.1790127
-7.926963	-31.064116	-27.5354159	1.818999	7.333289	16.392295	19.364755	12.065125	19.414690	15.465488	18.6837958
11.975101	-6.703292	-0.4937571	6.277274	8.482645	9.695062	13.222308	17.155351	12.259135	9.498506	14.4476463
9.039901	-2.146689	8.5439767	4.914373	9.514078	11.509139	7.617086	9.526847	10.776688	10.451743	8.0271782
4.933212	4.009910	3.2475677	6.808286	5.299639	6.871229	6.688790	7.089865	5.262589	2.158753	11.0289219
-1.008339	8.002981	8.4250450	5.383478	4.954822	4.208564	1.389350	4.824557	7.225718	4.454579	-0.7993775

Como vemos la base de datos está constituida por 40 filas (es decir, por 40 ciclos y 11 columnas). A diferencia de la base de datos anteriores,todas las columnas corresponden a muestras y no tenemos una columna para el número de ciclo. Las muestras correspoden a:

Las tres primeras a la muestra sin diluir
Las siguientes tres a las muestras diluidas 5 veces
Las siguientes tres a las muestras diluidas 25 veces
Las últimas dos a las muestras diluidas 125 veces

Ahora podemos graficar todas las muestras para ver como se ve la curva de intensidad de fluorescencia a lo largo de los 40 ciclos

matplot(db2,                             # plot especifico para matriz de datos
        t="l",                           # tipo de plot: lienal "l"
        col=c(rep(1:3,each=3),4,4),      # definimos loc colores de las lineas
        lty=1,                           # tipo de linea
        main="qRT-PCR Diluciones de una misma muestra",   # titulo de la grafica
        xlab="Ciclos (#)",                   # titulo del eje x
        ylab="Intensidad de Fluorescencia (a.u.)") # titulo del eje y
legend("topleft",                        # posición de la leyenda
       legend=c("Undiluted",             # nombres de las etiquetas en la leyenda 
                "5 x diluted", 
                "25 x diluted",
                "125 x diluted"),
       lty=1,                            # tipo de linea de la leyenda
       col=1:4)                          # definimos los colores

Aunque la forma gráfica ayuda mucho a entender como se está comportando la intensidad de flurescencia en cada muestra, podemo promediar las “réplicas técnicas” , para simplificar la gráfica y hacer más evidente los patrones (cabe aclarar, que partimos del supuesto que las “réplicas técnicas” tienen la misma concentración de template y por ende las curvas deberían ser la misma). Vamos a usar la función rowMeans() para calcular el promedio por linea (es decir para cada ciclo), de cada grupo de diluciones (es decir, cada tres columnas para las tres primeras diluciones y dos para la última)

no.dil <- rowMeans(db2[,1:3])            # promedio por ciclo par las primeras 3 columnas
dil5 <- rowMeans(db2[,4:6])              # promedio por ciclo par las siguientes 3 columnas
dil25 <- rowMeans(db2[,7:9])             # promedio por ciclo par las siguientes 3 columnas
dil125 <- rowMeans(db2[,10:11])          # promedio por ciclo par las últimas 2 columnas

db3 <- data.frame(no.dil = no.dil,
                  dil5 = dil5,
                  dil25 = dil25, 
                  dil125 = dil125)
# graficamos
matplot(db3,                             # plot especifico para matriz de datos
        t="l",                           # tipo de plot: lienal "l"
        col=1:4,                         # definimos loc colores de las lineas
        lty=1,                           # tipo de linea
        main="qRT-PCR Diluciones de una misma muestra",   # titulo de la grafica
        xlab="Ciclos (#)",                   # titulo del eje x
        ylab="Intensidad de Fluorescencia (a.u.)",  # titulo del eje y
        lwd = 2)                         # ancho de linea

# agregamos unas líneas extras para poder comparar más facilmente el nivel de fluorescencia en relación a los ciclos
abline(v=seq(0,40,1),lty=2, col="gray90")
abline(v=seq(0,40,5),lty=1, col="gray50")
abline(h=seq(0,8000,250),lty=2, col="gray90")
abline(h=seq(0,8000,1000),lty=1, col="gray50")
abline(v=27,lty=1, lwd=2, col="orange")

legend("topleft",                        # posición de la leyenda
       legend=c("Undiluted",             # nombres de las etiquetas en la leyenda 
                "5 x diluted", 
                "25 x diluted",
                "125 x diluted"),
       lty=1,                            # tipo de linea de la leyenda
       col=1:4)                          # definimos los colores

Como vemos en la figura, cuando una línea vertical (anaranjado) cruza las tres gráficas, podemos calcular a partir de las intensidades de flurorescencias que la relación de una con la otra es alrededor de 5 (como era de esperarse) Sin embargo, hacer este calculo de forma manual podria ser un poco tedioso y aumentar la variación o el sesgo si queremos comparar varios experimentos. Entonces se vió que se podría establecer un nivel de fluorescencia mímino, en el cual al curva empieza o es exponencial y que corte a todas las cruvas del experimento. A este nivel de fluorescencia se llamos Threshold, y a al ciclo en el cual el Threshold corta la curba de fluorescencia se llamó “cycle threshold” (Ct)

# graficamos
matplot(db3,                             # plot especifico para matriz de datos
        t="l",                           # tipo de plot: lienal "l"
        col=1:4,                         # definimos loc colores de las lineas
        lty=1,                           # tipo de linea
        main="qRT-PCR Diluciones de una misma muestra",   # titulo de la grafica
        xlab="Ciclos (#)",                   # titulo del eje x
        ylab="Intensidad de Fluorescencia (a.u.)",  # titulo del eje y
        lwd = 2)                         # ancho de linea

# agregamos unas líneas extras para poder comparar más facilmente el nivel de fluorescencia en relación a los ciclos
abline(v=seq(0,40,1),lty=2, col="gray90")
abline(v=seq(0,40,5),lty=1, col="gray50")
abline(h=seq(0,8000,250),lty=2, col="gray90")
abline(h=seq(0,8000,1000),lty=1, col="gray50")
abline(v=27,lty=1, lwd=2, col="orange")

abline(h=2000,lty=2, lwd = 2, col="magenta")
lines(c(25,25),c(2333,-300),lty=1,col="red")
lines(c(28,28),c(2130,-300),lty=1,col="green")
lines(c(31,31),c(2110,-300),lty=1,col="blue")
text(c(25,28,31),2500,
     labels=c("Ct=25","Ct=28","Ct=31"),
   srt=90,col=2:4)


legend("topleft",                        # posición de la leyenda
       legend=c("Undiluted",             # nombres de las etiquetas en la leyenda 
                "5 x diluted", 
                "25 x diluted",
                "125 x diluted"),
       lty=1,                            # tipo de linea de la leyenda
       col=1:4)                          # definimos los colores

Como sabemos, en cada ciclo duplicamos el numero de copias, por lo tanto el crecimiento de la curva es exponencial en base dos, es decir

\[ \frac{DNA_{dilB}}{DNA_{dilA}}=2^{Ct_B-Ct_A}= 2^{dCt} \]

Aunque estos calculos los podemos hacer manualmente, hay que cosiderar que la eficiencia de la reacción cambia y no siempre es “2” para todos los casos (de hecho, es el menos frecuente). R cuenta con varias paqueterías que nos permiten hacer estos calculos de forma sencilla. Vamos a utilizar la función qpcrfit del paquete qpcR para calcular las eficiencias de la reacción y así ajustar mejor nuestra curva teórica. Para esta función tenemos que darle como primer parámetro una base de datos cuya primera columna corresponda al número de cilos y, luego como segundo parámetro las columnas que corresponden a la fluorescencia.

library(qpcR)
 
no.dil <- qpcR::pcrfit(cbind(1:40,db2), fluo = 2:4)
dil5 <- qpcR::pcrfit(cbind(1:40, db2), fluo = 5:7)
dil25 <- qpcR::pcrfit(cbind(1:40, db2), fluo = 8:10)
dil125 <- qpcR::pcrfit(cbind(1:40, db2), fluo = 11:12)

Lo que hace esta función es aproximar o ajustar nuestros datos experimentales a una curva sigmoidea teórica. Aplicando un modelo no lineal, nos brinda el valor para cada uno de los coeficientes.

dil25

## Nonlinear regression model
##   model: Fluo ~ c + (d - c)/(1 + exp(b * (log(Cycles) - log(e))))
##    data: DATA
##       b       c       d       e 
##  -12.64  -19.71 7203.78   30.64 
##  residual sum-of-squares: 3073190
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 1.49e-08

summary(dil25)

## 
## Formula: Fluo ~ c + (d - c)/(1 + exp(b * (log(Cycles) - log(e))))
## 
## Parameters:
##     Estimate Std. Error t value Pr(>|t|)    
## b  -12.63796    0.36888 -34.261   <2e-16 ***
## c  -19.71107   19.99539  -0.986    0.326    
## d 7203.78434   76.01366  94.770   <2e-16 ***
## e   30.64374    0.08754 350.062   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 162.8 on 116 degrees of freedom
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 1.49e-08

Ahora podemos calcular la eficiencia para cada reacción usando la función efficeincy() del mismo paquete

qpcR::efficiency(dil25)

## [1] "pcrfit" "nls"

## $eff
## [1] 1.516266
## 
## $resVar
## [1] 26493.02
## 
## $AICc
## [1] 1568.981
## 
## $AIC
## [1] 1568.633
## 
## $Rsq
## [1] 0.9960621
## 
## $Rsq.ad
## [1] 0.9959602
## 
## $cpD1
## [1] 30.26
## 
## $cpD2
## [1] 27.07
## 
## $cpE
## [1] NA
## 
## $cpR
## [1] NA
## 
## $cpT
## [1] NA
## 
## $Cy0
## [1] NA
## 
## $cpCQ
## [1] NA
## 
## $cpMR
## [1] NA
## 
## $fluo
## [1] 1227.246
## 
## $init1
## [1] -19.71107
## 
## $init2
## [1] 0.01567985
## 
## $cf
## [1] NA

Este gráfico nos da mucha información:

El eje vertical de la izquierda vemos el nivel de fluorescencia.
El eje de coordenadas corresponde a los ciclos de la reacción.
El eje vertical de la derecha corresponde a la eficiencia de la reacción
Graficado con puntos no rellenos los datos experimentales, de los 3 conjuntos de datos
La linea negra continua corresponde al modelo logistico ajustado.
La línea curva azul corresponde a la eficiencia a lo largo de los ciclos.
La línea vertical verde corresponde a la primera derivada y corresponde al inicio de la fase exponencial, a la máxima eficiencia de la reacción, es decir al “threshold cycle” para estas reacciones.
La línea vertical roja corresponde a la segunda derivada y es la mitad de la parte exponencial de la curva, se considera como la eficiencia promedio de la reación.

A parte de los gráficos, la función efficienty() nos da un objeto lista con varios slots donde pondemos encontrar el valor de la primera y segunda derivada (ciclos), y los niveles de fluorescencia, entre otros varios. Un slot que vamos a utilizar es $fluo que nos permite saber el valor de fluorescencia “threshold” que vamos a utilizar para cruzar todas las graficas y obtener el valor del “cycle threshold”. Para hacer esto vamos a promediar todos los valores de $fluo de todas nuestras diluciones

eff.nodil   <- qpcR::efficiency(no.dil, plot = F)  # inactivamos la funcion de graficar para ejecutar estas lineas más rápidas/límpias 
eff.5dil    <- qpcR::efficiency(dil5, plot = F)
eff.25dil   <- qpcR::efficiency(dil25, plot = F)
eff.125dil  <- qpcR::efficiency(dil125, plot = F)

meanFluo <- mean(c(eff.nodil$fluo, eff.5dil$fluo, eff.25dil$fluo, eff.125dil$fluo))
sdFluo <- sd(c(eff.nodil$fluo, eff.5dil$fluo, eff.25dil$fluo, eff.125dil$fluo))

meanFluo

## [1] 1132.034

sdFluo

## [1] 116.065

Vemos que el nivel promedio de fluoresncia para el threshold es de 1132 +/- 116 unidades arbitrarias.

Ahora que ya hemos definido nuestro threshold, podemo volver a ejecutar nuestra función efficiency(), pero pre-estableciendo un threshold para que nos de el valor de “cycle threshold” que estamos buscando.

ct.nodil <- qpcR::efficiency(no.dil, plot = F, thresh = meanFluo)$cpT # para guardar solo el valor del cT
ct.dil5 <- qpcR::efficiency(dil5, plot = F, thresh = meanFluo)$cpT
ct.dil25 <- qpcR::efficiency(dil25, plot = F, thresh = meanFluo)$cpT
ct.dil125 <- qpcR::efficiency(dil125, plot = F, thresh = meanFluo)$cpT

ct.nodil

## [1] 21.58

ct.dil5

## [1] 23.31

ct.dil25

## [1] 26.87

ct.dil125

## [1] 29.7

\[ \frac{DNA_{sampleA}}{DNA_{sampleB}} = \frac{Efficiency_{sampleA}+Efficiency_{sampleB}}{2}^{Ct_A-Ct_B} \]

qRT-PCR

Cardenas-Rivera, A.

2023-05-18