Microarrays
Cardenas-Rivera, A.
2023-05-30
MICROARRAY
En español “Microarreglos”, fue una tecnología de vertiginoso ascenso en los mediados y finales de los 90’s, pero que fue rápidamente desplazadqa por la “Secuencienación de nueva Generación”. Todo el desarrollo de la técnologia de los microarreglos se vio pronto desfazado por este último, pese a todos los algortimos y otras fueciones que se desarrollaron
Principales empresas
Actualente, los proveedores mayoritarios de microaarreglos son tres:
- Affymetrix
- Illumina
- Agilent.
De estos, el más popular es Affymetric.
Cada marca de microarreglo tiene una paqueterîa propia y esclucisiva
que solo maneja la información de esa marca. Asi tenemos los paquetes
affy para Affymetrix, agilp para agilent y
lumi para Illumina.
Instalación de paqueteria
## Bioconductor version '3.16' is out-of-date; the current release version '3.17'
## is available with R version '4.3'; see https://bioconductor.org/install## Bioconductor version 3.16 (BiocManager 1.30.20), R 4.2.3 (2023-03-15)## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'affy'## Loading required package: BiocGenerics##
## Attaching package: 'BiocGenerics'## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
## table, tapply, union, unique, unsplit, which.max, which.min## Loading required package: Biobase## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.Importación de una base de datos
Vamos a trabajar on una base de datos de ejemplo hecha en la
plataforma de Affymetrix. Los archivos que brindan sus aparatos los
podemos obtener como archivos txt simples o como archivos comprrimidos
CEL. En caso de esta paque tería no necesitamos
descromprimir el archivo, la propia funcîon descomprime y lee el
archivo. Para importar los doce experimentos .CEL del
ejemplo vamos a utilizar la función affy::ReadAffy, y en el
parámetro celfile.path= vamos a colocar la raiz del folder
donde se encuentran nuestros archivos
dat <- affy::ReadAffy(celfile.path = "/Users/alfredocardenasrivera/Project/bioinformatics_2022/Micrroarray/GSE29797_subset")Una vez ya cargados, podemos revisar el iterior de este objeto con las funciones básicas que habitualmente usamos.
str(dat)## Formal class 'AffyBatch' [package "affy"] with 10 slots
## ..@ cdfName : chr "HT_MG-430_PM"
## ..@ nrow : Named int 744
## .. ..- attr(*, "names")= chr "Rows"
## ..@ ncol : Named int 744
## .. ..- attr(*, "names")= chr "Cols"
## ..@ assayData :<environment: 0x11ab9b200>
## ..@ phenoData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 1 obs. of 1 variable:
## .. .. .. ..$ labelDescription: chr "arbitrary numbering"
## .. .. ..@ data :'data.frame': 12 obs. of 1 variable:
## .. .. .. ..$ sample: int [1:12] 1 2 3 4 5 6 7 8 9 10 ...
## .. .. ..@ dimLabels : chr [1:2] "sampleNames" "sampleColumns"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## .. .. .. .. ..$ names: chr "AnnotatedDataFrame"
## ..@ featureData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 0 obs. of 1 variable:
## .. .. .. ..$ labelDescription: chr(0)
## .. .. ..@ data :'data.frame': 553536 obs. of 0 variables
## .. .. ..@ dimLabels : chr [1:2] "featureNames" "featureColumns"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## .. .. .. .. ..$ names: chr "AnnotatedDataFrame"
## ..@ experimentData :Formal class 'MIAME' [package "Biobase"] with 13 slots
## .. .. ..@ name : chr ""
## .. .. ..@ lab : chr ""
## .. .. ..@ contact : chr ""
## .. .. ..@ title : chr ""
## .. .. ..@ abstract : chr ""
## .. .. ..@ url : chr ""
## .. .. ..@ pubMedIds : chr ""
## .. .. ..@ samples : list()
## .. .. ..@ hybridizations : list()
## .. .. ..@ normControls : list()
## .. .. ..@ preprocessing :List of 2
## .. .. .. ..$ filenames : chr [1:12] "/Users/alfredocardenasrivera/Project/bioinformatics_2022/Micrroarray/GSE29797_subset/GSM738363_ykr264-htmg430pm.CEL.gz" "/Users/alfredocardenasrivera/Project/bioinformatics_2022/Micrroarray/GSE29797_subset/GSM738365_ykr265-htmg430pm.CEL.gz" "/Users/alfredocardenasrivera/Project/bioinformatics_2022/Micrroarray/GSE29797_subset/GSM738367_ykr267-htmg430pm.CEL.gz" "/Users/alfredocardenasrivera/Project/bioinformatics_2022/Micrroarray/GSE29797_subset/GSM738368_ykr268-htmg430pm.CEL.gz" ...
## .. .. .. ..$ affyversion: chr NA
## .. .. ..@ other :List of 1
## .. .. .. ..$ : chr ""
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 2
## .. .. .. .. .. ..$ : int [1:3] 1 0 0
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## .. .. .. .. ..$ names: chr [1:2] "MIAxE" "MIAME"
## ..@ annotation : chr "htmg430pm"
## ..@ protocolData :Formal class 'AnnotatedDataFrame' [package "Biobase"] with 4 slots
## .. .. ..@ varMetadata :'data.frame': 1 obs. of 1 variable:
## .. .. .. ..$ labelDescription: chr NA
## .. .. ..@ data :'data.frame': 12 obs. of 1 variable:
## .. .. .. ..$ ScanDate: chr [1:12] "2010-12-04T17:54:48Z" "2010-12-04T18:38:03Z" "2010-12-04T18:45:08Z" "2010-12-04T18:48:59Z" ...
## .. .. ..@ dimLabels : chr [1:2] "sampleNames" "sampleColumns"
## .. .. ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. .. .. ..@ .Data:List of 1
## .. .. .. .. .. ..$ : int [1:3] 1 1 0
## .. .. .. .. ..$ names: chr "AnnotatedDataFrame"
## ..@ .__classVersion__:Formal class 'Versions' [package "Biobase"] with 1 slot
## .. .. ..@ .Data:List of 4
## .. .. .. ..$ : int [1:3] 4 2 1
## .. .. .. ..$ : int [1:3] 2 58 0
## .. .. .. ..$ : int [1:3] 1 3 0
## .. .. .. ..$ : int [1:3] 1 2 0
## .. .. ..$ names: chr [1:4] "R" "Biobase" "eSet" "AffyBatch"summary(dat)## Length Class Mode
## 12 AffyBatch S4class(dat)## [1] "AffyBatch"
## attr(,"package")
## [1] "affy"Como vemos es un objeto S4, de la clase affy. Dentro de el encontramos 12 experientos separados (es decir, se corrieron 12 microarreglos de la misma marca y el mismo modelo)
Si nosotros llamamos al objeto, veremos que nos dice el tamaño/resolución del chip,el modelo del CHIP utilizado. Que en este caso es “HT_MG-430_PM”, número de genes cuantificados por CHIP es de 45141 y el número de muestras que se corrieron fue d e12.
Manejo de la base de datos “affy”
Para saber el número y nombre de cada una de las muestras en el
objeto podemos usar pData del mismo paquete.
Biobase::pData(dat)## sample
## GSM738363_ykr264-htmg430pm.CEL.gz 1
## GSM738365_ykr265-htmg430pm.CEL.gz 2
## GSM738367_ykr267-htmg430pm.CEL.gz 3
## GSM738368_ykr268-htmg430pm.CEL.gz 4
## GSM738370_ykr270-htmg430pm.CEL.gz 5
## GSM738372_ykr271-htmg430pm.CEL.gz 6
## GSM738373_ykr273-htmg430pm.CEL.gz 7
## GSM738375_ykr274-htmg430pm.CEL.gz 8
## GSM738377_ykr276-htmg430pm.CEL.gz 9
## GSM738378_ykr277-htmg430pm.CEL.gz 10
## GSM738379_ykr279-htmg430pm.CEL.gz 11
## GSM738380_ykr280-htmg430pm.CEL.gz 12Con esta información vamos a crear nuestro data.fame de la metadata de nuestro experimetno.
exp.des <- pData(dat)
exp.des$genotype <- factor(rep(c("WT","dTsc1"),each=6))
exp.des$Stimulation <- factor(rep(c("0h","4h"),6))
exp.des## sample genotype Stimulation
## GSM738363_ykr264-htmg430pm.CEL.gz 1 WT 0h
## GSM738365_ykr265-htmg430pm.CEL.gz 2 WT 4h
## GSM738367_ykr267-htmg430pm.CEL.gz 3 WT 0h
## GSM738368_ykr268-htmg430pm.CEL.gz 4 WT 4h
## GSM738370_ykr270-htmg430pm.CEL.gz 5 WT 0h
## GSM738372_ykr271-htmg430pm.CEL.gz 6 WT 4h
## GSM738373_ykr273-htmg430pm.CEL.gz 7 dTsc1 0h
## GSM738375_ykr274-htmg430pm.CEL.gz 8 dTsc1 4h
## GSM738377_ykr276-htmg430pm.CEL.gz 9 dTsc1 0h
## GSM738378_ykr277-htmg430pm.CEL.gz 10 dTsc1 4h
## GSM738379_ykr279-htmg430pm.CEL.gz 11 dTsc1 0h
## GSM738380_ykr280-htmg430pm.CEL.gz 12 dTsc1 4hpData(dat) <- exp.desimage(dat[,1])
hist(dat)## Warning: replacing previous import 'AnnotationDbi::tail' by 'utils::tail' when
## loading 'htmg430pmcdf'## Warning: replacing previous import 'AnnotationDbi::head' by 'utils::head' when
## loading 'htmg430pmcdf'## 
boxplot(exprs(dat))
boxplot(log2(exprs(dat)))
CONTROL DE CALIDAD
deg <- AffyRNAdeg(dat)
summaryAffyRNAdeg(deg)## GSM738363_ykr264-htmg430pm.CEL.gz GSM738365_ykr265-htmg430pm.CEL.gz
## slope 4.15e+00 3.81e+00
## pvalue 2.07e-14 2.89e-13
## GSM738367_ykr267-htmg430pm.CEL.gz GSM738368_ykr268-htmg430pm.CEL.gz
## slope 3.93e+00 3.77e+00
## pvalue 1.76e-13 5.49e-13
## GSM738370_ykr270-htmg430pm.CEL.gz GSM738372_ykr271-htmg430pm.CEL.gz
## slope 4.09e+00 3.83e+00
## pvalue 2.29e-14 3.10e-13
## GSM738373_ykr273-htmg430pm.CEL.gz GSM738375_ykr274-htmg430pm.CEL.gz
## slope 4.07e+00 4.94e+00
## pvalue 4.21e-14 2.40e-14
## GSM738377_ykr276-htmg430pm.CEL.gz GSM738378_ykr277-htmg430pm.CEL.gz
## slope 4.03e+00 3.80e+00
## pvalue 6.09e-14 2.77e-13
## GSM738379_ykr279-htmg430pm.CEL.gz GSM738380_ykr280-htmg430pm.CEL.gz
## slope 4.05e+00 3.93e+00
## pvalue 1.10e-13 1.92e-13plotAffyRNAdeg(deg)
Para saber el código o el id del Microarreglo podemos usar la función
cdfName()
cdfName(dat)## [1] "HT_MG-430_PM"O también llamar directamente al objeto AffyBatch
dat## AffyBatch object
## size of arrays=744x744 features (25 kb)
## cdf=HT_MG-430_PM (45141 affyids)
## number of samples=12
## number of genes=45141
## annotation=htmg430pm
## notes=Y la tercera forma es usando la función genérica
annotation()
annotation(dat)## [1] "htmg430pm"NOTA: Cuando correamos o cargamos los datos con la función
affyRead(), esta función tambien carga el paquete relacionado al microarreglo utilizado en los experimentos. En este caso htmg430. Si Ud. se fija en los paquetes instalados y cargados actualmente a su entorno, encontrará que “htmg430” esta descargado y activado.
NORMALIZACION DE LOS DATOS
Actualmente contamos con muchos métodos/algoritmos para realizar la normalización de los datos de microarreglos de diferentes ensayos. El objetivo principal de la normalización es compensar la diferencia de intensidad de las muestras producto del manejo de ellas en su procesamiento, de tal forma que las diferencias que se encuentren posteriormente sean por un factor biológico real y no solo producto de un artefacto técnico entre las muestras que crea una diferencia de intensidad aparente.
Los métodos de normalización mas utilizados son:
- MAS 5.0 : Es un método de normalización provisto por la empresa Affymetrix. Se caracteriza porque toma en cuenta las sondas con Miss Match (MM). Actualmente no es muy usado.
- RMA: El RMA o “Robust Multi-Array Average” utiliza los Perfect Match (PM) y hace una normalización por cuantiles.
- VSN: Nromalización por estabilización de la varianza. Este método trata de mantener constante la varianza entre todas las muestras.
- LOWESSS: Normalización por alisamiento local ponderado de la gráfica de dispersión. Es usado preferentemente con métodos que usan dos fluoróforos. Usa una gráfica en el cual uno de los ejes es igual a \(log(red/green)\) y el otro eje es \(log\sqrt{(Red*Green)}\). Este gráfico es llamada “M plot”
Para nuestro ejemplo, vamos a usar le método de normalizaci´øn RMA
dat.rma <- affy::rma(dat)## Background correcting
## Normalizing
## Calculating ExpressionPara ver la efeciencia de nuestra normalización, vamos a graficar los resultados antes y despuest de la normalización.
matexp<-exprs(dat.rma) # primero extraemos los datos de expresción
par(mfrow=c(1,2)) # programamos el área de graficado en dos mitades
boxplot(log2(exprs(dat)),main="Before normalization", # Boxplot de los dtos en bruto
ylab="log2(Intensity)")
boxplot(matexp,main="RMA normalized data") # Mismo datos normalizados 
par(mfrow=c(1,1))ANALISIS DE EXPRESI´ON
Identificación del agrupamiento de los datos
Para evaluar la expresión de los genes individuales, primero sería
bueno evaluar si los datos en cojunto se comienzan a agrupar. Es decir
si tenemos dos condiciones experimentales, esperariamos que nuestros
datos por cada microarreglo se agrupoen. Este análisis de PCA lo podemos
hacer con una función nativa o podemo usar una variante específica de
PCA que se aplica a estos microarreglos y que está definida en el
paquete affycortetools. Para utilizar esta función
necesitamos instalarlo el paquete affycoretools desde el
respositorio de Bioconductor, de la siguiente
manera.
## Bioconductor version 3.16 (BiocManager 1.30.20), R 4.2.3 (2023-03-15)## Warning: package(s) not installed when version(s) same as or greater than current; use
## `force = TRUE` to re-install: 'affycoretools'## ## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2## library(affycoretools)
par(mfrow=c(1,2))
plotPCA(matexp, # datos de expresión
groups = as.numeric(pData(dat)[,2]), # Agrupamiento de datos (genotype)
groupnames = levels(pData(dat)[,2])) # Nombre de los niveles
plotPCA(matexp, # datos de expresión
groups = as.numeric(pData(dat)[,3]), # Agrupamiento de datos (punto temporal)
groupnames = levels(pData(dat)[,3])) # Nombre de los niveles
par(mfrow=c(1,1))Como podemos ver las muestras se agrupan bastante bien por genotipo y por punto temporal de muestreo. En la parte superior de los gráficos se agrupan los genotipos dTsc1 y en la parte inferior los genotipos WT. A la izquierda se agrupan los grupos en tiempo cero y a la derecha se agrupa los grupos en tiempo 4h.
Construcción de las matrices para comparación
Ahora que ya sabemos que los datos se agrupan adecuadamente podemos
continuar con el análisis creando o “diseñando” las
matrices para ambas categorias. Primero tenemos que convertir las
columnas del data.frame pData en factores, y luego crear un
model.matrix para cada una de las variables (genotipo y
punto temporal).
head(pData(dat))## sample genotype Stimulation
## GSM738363_ykr264-htmg430pm.CEL.gz 1 WT 0h
## GSM738365_ykr265-htmg430pm.CEL.gz 2 WT 4h
## GSM738367_ykr267-htmg430pm.CEL.gz 3 WT 0h
## GSM738368_ykr268-htmg430pm.CEL.gz 4 WT 4h
## GSM738370_ykr270-htmg430pm.CEL.gz 5 WT 0h
## GSM738372_ykr271-htmg430pm.CEL.gz 6 WT 4h# diseño de la matriz con la columna "Genotipo"
genotype <- factor(pData(dat)[,2])
design.g <- model.matrix(~genotype)
knitr::kable(design.g)| (Intercept) | genotypeWT |
|---|---|
| 1 | 1 |
| 1 | 1 |
| 1 | 1 |
| 1 | 1 |
| 1 | 1 |
| 1 | 1 |
| 1 | 0 |
| 1 | 0 |
| 1 | 0 |
| 1 | 0 |
| 1 | 0 |
| 1 | 0 |
# diseño de la matriz con la columna "Punto temporal"
stimul <- factor(pData(dat)[,3])
design.t <- model.matrix(~stimul)
knitr::kable(design.t)| (Intercept) | stimul4h |
|---|---|
| 1 | 0 |
| 1 | 1 |
| 1 | 0 |
| 1 | 1 |
| 1 | 0 |
| 1 | 1 |
| 1 | 0 |
| 1 | 1 |
| 1 | 0 |
| 1 | 1 |
| 1 | 0 |
| 1 | 1 |
Ajuste a modelo lienal
Con los datos de expresión y con las matrices de agrupamiento, ahora podemos usar el paquete limma para hacer el ajuste de los datos a un modelo lineal. Vamos a comenzar el análisis dividiendo los datos por genotipo.
NOTA: El paquete limma no está actualizado para la última versión de R, entonces necesitamos activar el parámetro
force=TRUEcuando lo instalemos desde Bioconductor, usando la funciónBiocManager::install()
## Bioconductor version 3.16 (BiocManager 1.30.20), R 4.2.3 (2023-03-15)## Installing package(s) 'limma'##
## The downloaded binary packages are in
## /var/folders/70/4klcrc8j0nn3br1qq16nc_8w0000gn/T//RtmpmpSsEu/downloaded_packagesVamos a usar la función lmfit() del paquete
limma para generar un modelo lineal para la expresión
de cada gen. Esta función asume que estamos usando las mismas sondas
para todos los experimentos, es decir, estamos usando el mismo modelo de
microarreglo. Esta función acepta tanto datos de ensayos con un
fluoróforo como con dos fluoróforos. Los datos deben estar en escala
logarítmica.
library(limma)##
## Attaching package: 'limma'## The following object is masked from 'package:BiocGenerics':
##
## plotMAfit <- limma::lmFit(object = matexp, design = design.g)
fit## An object of class "MArrayLM"
## $coefficients
## (Intercept) genotypeWT
## 1415670_PM_at 8.280645 0.09075880
## 1415671_PM_at 9.256379 -0.64344298
## 1415672_PM_at 9.596247 -0.18578687
## 1415673_PM_at 7.793342 -0.89930999
## 1415674_PM_a_at 8.417905 0.07286682
## 45136 more rows ...
##
## $rank
## [1] 2
##
## $assign
## [1] 0 1
##
## $qr
## $qr
## (Intercept) genotypeWT
## 1 -3.4641016 -1.7320508
## 2 0.2886751 -1.7320508
## 3 0.2886751 0.2240092
## 4 0.2886751 0.2240092
## 5 0.2886751 0.2240092
## 7 more rows ...
##
## $qraux
## [1] 1.288675 1.224009
##
## $pivot
## [1] 1 2
##
## $tol
## [1] 1e-07
##
## $rank
## [1] 2
##
##
## $df.residual
## [1] 10 10 10 10 10
## 45136 more elements ...
##
## $sigma
## 1415670_PM_at 1415671_PM_at 1415672_PM_at 1415673_PM_at 1415674_PM_a_at
## 0.2384200 0.2007284 0.1439598 0.7019110 0.2357675
## 45136 more elements ...
##
## $cov.coefficients
## (Intercept) genotypeWT
## (Intercept) 0.1666667 -0.1666667
## genotypeWT -0.1666667 0.3333333
##
## $stdev.unscaled
## (Intercept) genotypeWT
## 1415670_PM_at 0.4082483 0.5773503
## 1415671_PM_at 0.4082483 0.5773503
## 1415672_PM_at 0.4082483 0.5773503
## 1415673_PM_at 0.4082483 0.5773503
## 1415674_PM_a_at 0.4082483 0.5773503
## 45136 more rows ...
##
## $pivot
## [1] 1 2
##
## $Amean
## 1415670_PM_at 1415671_PM_at 1415672_PM_at 1415673_PM_at 1415674_PM_a_at
## 8.326024 8.934657 9.503353 7.343687 8.454338
## 45136 more elements ...
##
## $method
## [1] "ls"
##
## $design
## (Intercept) genotypeWT
## 1 1 1
## 2 1 1
## 3 1 1
## 4 1 1
## 5 1 1
## 7 more rows ...Analisis de expresión diferencial
El resultado es un objeto complejo de la clase MArrayLM que contiene
varios slots. Para determinar si hay una diferencia significativa en la
expresión de los genes, usaremos la función eBayes del
paquete limma, y luego usaremos la función
topTable para mostrar los primeros resultados según su
valor de FDR
fit <- limma::eBayes(fit)
fit.genotype <- fit
res <- limma::topTable(fit, coef = 2, number = 50,adjust.method = "fdr")
knitr::kable(head(res))| logFC | AveExpr | t | P.Value | adj.P.Val | B | |
|---|---|---|---|---|---|---|
| 1419758_PM_at | -3.972031 | 6.200539 | -30.06466 | 0 | 0.0000000 | 14.348063 |
| 1419759_PM_at | -3.648735 | 6.532583 | -24.20075 | 0 | 0.0000002 | 13.384220 |
| 1434909_PM_at | -2.955792 | 5.339058 | -19.72618 | 0 | 0.0000018 | 12.226986 |
| 1434437_PM_x_at | -1.656676 | 5.756127 | -18.69385 | 0 | 0.0000026 | 11.881989 |
| 1449874_PM_at | -1.576391 | 4.107422 | -15.16105 | 0 | 0.0000242 | 10.384427 |
| 1448226_PM_at | -1.624405 | 4.320060 | -13.04231 | 0 | 0.0001153 | 9.173986 |
El resultado es un data.frame con 6 columnas, y 50 filas (en este
caso, así lo definimos en el parámetro number de la función
topTable). Ahora vamos a explicar con un poco de mas
detalle cada una de las columnas
- logFC: Es el logaritmo del “Fold change” de la expresión de cada gen. Entre más grande sea este número (positivo o negativamente) quiere decir que la expresión de ese gene es más diferente entre los grupos. El valor negativo indica respresción y el valor positivo indica sobreexpresión del gen.
- AveExpr: Es el promedio de expresión del gen para ambos grupos. Esto nos da idea de cuan abundantemente se expresa ese gen en comparación con el resto de genes.
- t: Este es el valor del estadístico de t-student para cada gen, comparando los dos grupos
- P.Value: Este es el valor de P significativa según el análisis de t-student
- adj.P.Val: Este es el valor ajustado de P, corrigiendo las comparaciones múltiples
- B: Es el valor de “Log-odds” y nos indica cuan “real” es la diferencia que encontramos. A valores mas altos, mayor grado de certeza.
Aunque el análisis ya está listo, lo que faltaría por hacer sería asignar el nombre de los genes, porque hasta el momento solo tenemos los códigos de los IDs para cada sonda (que sería equivalente a cada gen)
Anotación de los genes
Para asignar el nombre de los genes, según el código que tiene la sonda tenemos al menos tres formas de hacerlo:
- Descargar la base de datos de la plataforma/microarreglo. Usualmente estas se encuentran disponibles en el CRAN de R y el nombre del paquete es igual a: “nombre_del_microarreglo.db”. Así por ejemplo si la plataforma utilizada se llama: “Human Genome U133A 2.0”, el nombre sería: “hgu133a2.db”
- Descargarlo directamente de la página de GEO-dataset (NCBI). Quizá la forma más sencilla, todos los archivos subidos o almacenados en la base de datos GEO-dataset deben especificar las anotaciones de los genes de las plataformas utilizadas. Estos archivos son descargables como .txt y pueden ser cargados posteriormente a nuestro entorno de trabajo.
- BioMart. Este base de datos nos brinda los códigos de diferentes microarreglos para cada gen en su base de datos.
Para este ejercicio, vamos a suponer que descargamos las anotaciones de NCBI y tenemos ese archivo como un .txt en nuestra computadora. Entonces lo primero que tenemos que hacer es cargar el archivo
glp.annot <- read.table("/Users/alfredocardenasrivera/Project/bioinformatics_2022/Micrroarray/GPL11180_selected_annotation.txt")
head(glp.annot)## GB_ACC RefSeq.Transcript.ID Entrez.Gene Gene.Symbol
## 1415670_PM_at BC024686 NM_017477 /// NM_201244 54161 Copg
## 1415671_PM_at NM_013477 NM_013477 11972 Atp6v0d1
## 1415672_PM_at NM_020585 NM_001042484 /// NM_020585 57437 Golga7
## 1415673_PM_at NM_133900 NM_133900 100678 Psph
## 1415674_PM_a_at NM_021789 NM_021789 60409 Trappc4
## 1415675_PM_at BC008256 NM_010073 13481 Dpm2Ahora lo que temos que hacer es emparejar los códigos del
microarreglo con su respectivo gen. Para ello, vamos a llamar a los
nombres de las filas con la función row.names() de nuestros
resultados de 50 principales res y lo vamos a usar como
filtro selector de las filas que deseamos del objeto GLP
glp.annot que acabamos de cargar.
genotype.50 <- glp.annot[row.names(res),]
knitr::kable(head(genotype.50))| GB_ACC | RefSeq.Transcript.ID | Entrez.Gene | Gene.Symbol | |
|---|---|---|---|---|
| 1419758_PM_at | M30697 | NM_011076 | 18671 | Abcb1a |
| 1419759_PM_at | M30697 | NM_011076 | 18671 | Abcb1a |
| 1434909_PM_at | BF462770 | NM_027491 | 52187 | Rragd |
| 1434437_PM_x_at | AV301324 | NM_009104 | 20135 | Rrm2 |
| 1449874_PM_at | NM_016923 | NM_001159711 /// NM_016923 | 17087 | Ly96 |
| 1448226_PM_at | NM_009104 | NM_009104 | 20135 | Rrm2 |
Finalmente fusionaremos el data.frame resultante
genotype.50 con la base de datos de los 50 principales
res; para ello usaremos la función cbind()
genotype.50 <- cbind(res, genotype.50)
knitr::kable(head(genotype.50))| logFC | AveExpr | t | P.Value | adj.P.Val | B | GB_ACC | RefSeq.Transcript.ID | Entrez.Gene | Gene.Symbol | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1419758_PM_at | -3.972031 | 6.200539 | -30.06466 | 0 | 0.0000000 | 14.348063 | M30697 | NM_011076 | 18671 | Abcb1a |
| 1419759_PM_at | -3.648735 | 6.532583 | -24.20075 | 0 | 0.0000002 | 13.384220 | M30697 | NM_011076 | 18671 | Abcb1a |
| 1434909_PM_at | -2.955792 | 5.339058 | -19.72618 | 0 | 0.0000018 | 12.226986 | BF462770 | NM_027491 | 52187 | Rragd |
| 1434437_PM_x_at | -1.656676 | 5.756127 | -18.69385 | 0 | 0.0000026 | 11.881989 | AV301324 | NM_009104 | 20135 | Rrm2 |
| 1449874_PM_at | -1.576391 | 4.107422 | -15.16105 | 0 | 0.0000242 | 10.384427 | NM_016923 | NM_001159711 /// NM_016923 | 17087 | Ly96 |
| 1448226_PM_at | -1.624405 | 4.320060 | -13.04231 | 0 | 0.0001153 | 9.173986 | NM_009104 | NM_009104 | 20135 | Rrm2 |
En el paso de diseño de matrices, que hicimos anteriormente, tambien
realizamos una matriz para marcar o definir los grupos experimentales
según la fase del experimento: Al inicio del experimento (0 horas) y
luego del tratamiento (4h). Esta matriz la habíamos llamado
design.t, y usando la matriz de expresión
matexp, determinaremos la expresión diferencial de los
genes relacionando estos dos grupos con la función lmfit()
del paquete limma
res <- limma::lmFit(object = matexp,design = design.t)Nuevamente aplicamos la función eBayes() del mismo
paquete para determinar la significancia estadistica de los niveles de
expersión
fit.stim <- limma::eBayes(res)Y extraemos los 50 principales usando la función
topTable
res <- limma::topTable(fit.stim, coef = 2, adjust.method = "fdr", number = 50)
knitr::kable(head(res))| logFC | AveExpr | t | P.Value | adj.P.Val | B | |
|---|---|---|---|---|---|---|
| 1418322_PM_at | 5.451960 | 6.313120 | 62.71547 | 0 | 0 | 31.63014 |
| 1425822_PM_a_at | -4.981851 | 7.496269 | -44.71355 | 0 | 0 | 28.21168 |
| 1449037_PM_at | 6.061854 | 6.871326 | 42.99404 | 0 | 0 | 27.76799 |
| 1417038_PM_at | -3.884873 | 8.236228 | -40.30581 | 0 | 0 | 27.01878 |
| 1454893_PM_at | -2.367691 | 8.834408 | -39.51840 | 0 | 0 | 26.78542 |
| 1452621_PM_at | -2.834306 | 7.321527 | -38.26152 | 0 | 0 | 26.39887 |
Ahora que ya tenemos los 50 principales en la base de datos
res, usaremos la anotaciones del objeto GLP
glp.annot para emparejar los genes con los códigos de los
microarreglos.
stim.50 <- glp.annot[row.names(res),]
stim.50 <- cbind(res,stim.50)
knitr::kable(head(stim.50))| logFC | AveExpr | t | P.Value | adj.P.Val | B | GB_ACC | RefSeq.Transcript.ID | Entrez.Gene | Gene.Symbol | |
|---|---|---|---|---|---|---|---|---|---|---|
| 1418322_PM_at | 5.451960 | 6.313120 | 62.71547 | 0 | 0 | 31.63014 | AI467599 | NM_001110850 /// NM_001110851 /// NM_001110852 /// NM_001110853 /// NM_001110854 /// NM_001110855 /// NM_001110856 /// NM_001110857 /// NM_001110858 /// NM_001110859 /// NM_013498 | 12916 | Crem |
| 1425822_PM_a_at | -4.981851 | 7.496269 | -44.71355 | 0 | 0 | 28.21168 | AB015422 | NM_008052 | 14357 | Dtx1 |
| 1449037_PM_at | 6.061854 | 6.871326 | 42.99404 | 0 | 0 | 27.76799 | AI467599 | NM_001110850 /// NM_001110851 /// NM_001110852 /// NM_001110853 /// NM_001110854 /// NM_001110855 /// NM_001110856 /// NM_001110857 /// NM_001110858 /// NM_001110859 /// NM_013498 | 12916 | Crem |
| 1417038_PM_at | -3.884873 | 8.236228 | -40.30581 | 0 | 0 | 27.01878 | NM_017380 | NM_001113486 /// NM_001113487 /// NM_001113488 /// NM_017380 | 53860 | Sept9 |
| 1454893_PM_at | -2.367691 | 8.834408 | -39.51840 | 0 | 0 | 26.78542 | BB765852 | NM_001014995 | 68521 | Fam189b |
| 1452621_PM_at | -2.834306 | 7.321527 | -38.26152 | 0 | 0 | 26.39887 | AI324894 | NM_028281 | 72562 | Pcbd2 |
Ahora veremos si hay genes compartidos entre los diferencialmente expresados por el genotipo de los organismos y por los tratamientos. Para ello, vamos a usar un nuevo paquete llamado “ggven” para graficar diagramas de Venn. Este paqute, lo vamos a descargar
#if (!require(devtools)) install.packages("devtools")
#devtools::install_github("yanlinlin82/ggvenn")
library(ggvenn)## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following object is masked from 'package:Biobase':
##
## combine## The following objects are masked from 'package:BiocGenerics':
##
## combine, intersect, setdiff, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: grid## Loading required package: ggplot2ggvenn(data = list(Genotype = row.names(genotype.50), Stim = row.names(stim.50)),
fill_color = c("#0073C2FF", "#EFC000FF", "#868686FF"),
stroke_size = 0.5, set_name_size = 5
)
Como podemos ver en el diagrama de Venn, no hay genes que se compartan entre los dos grupos de genes diferencialmente expresados. Esto quiere decir que los genes que se expresan diferencialmente debido al tratamiento son diferentes a los relacionados a los genotipos de los animales.
Uno de los análisis más utilizados para analizar la expresión
diferencial de los genes son los mapas de calor. Primero seleccionaremos
todos los genes que tengan un nivel de expresión significativa mayor de
0.03 (usando el método de ajuste de “fdr” para comparaciones multiples).
Para ello, usaremos la función p.adjust(). Una vez
obtengamos la lista de genes con esta función, extraeremos los genes de
la base de datos matexp y la graficaremos con la función
heatmap()
idx <- p.adjust(p = fit.stim$p.value[,2],method ="fdr") < 0.03
matexp.stim <- matexp[idx,]
dim(matexp.stim)## [1] 14398 12sum(idx)## [1] 14398heatmap(matexp.stim)