DESeq2_Intro

INTRODUCCIÓN a DEseq2 en R/Rstudio

[iBio] http://www.ibio.cl/2019/08/06/workshop-ii-ibio-bioinformatica iBio

1. Introducción

Existen diversos paquetes en R para realizar análisis de expresión diferencial. El más clásico es limma el cual está diseñado para expresión de microarreglos usando modelos lineales y que en un principio fue utilizado también en RNAseq. Con el tiempo la experiencia demostró que los modelos lineales no se ajustan a los resultados de RNAseq por lo cual se desarrolló un amplio debate sobre cuál es el mejor modelo para este tipo de experimentos contratando con la realidad. Esto dió origen a distintos paquetes de R donde hoy en día los más comunes son edgeR, cummeRbund y DESeq2.

En el siguiente link encontrarán distintas explicaciones sobre las diferencias de los métodos https://www.biostars.org/p/284775/

En este manual utilizaré el paquete DESeq2, el cual tiene un desarrollo más reciente y se está posicionando rápidamente como un estándar en el área.

Toda la información detallada sobre este paquete la pueden encontrar en (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)

2. Carga de datos en R

Iniciamos el trabajo instalando los paquetes requeridos (en caso que no estén instalados)

if (!requireNamespace(“BiocManager”, quietly = TRUE))
install.packages(“BiocManager”)
BiocManager::install(“DESeq2”)

Cargamos en memoria los paquetes requeridos (Se debe ejecutar cada vez que se inicia RStudio)

library("DESeq2")

## Loading required package: S4Vectors

## Loading required package: stats4

## Loading required package: BiocGenerics

## Loading required package: parallel

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, 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,
##     which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following object is masked from 'package:base':
## 
##     expand.grid

## Loading required package: IRanges

## 
## Attaching package: 'IRanges'

## The following object is masked from 'package:grDevices':
## 
##     windows

## Loading required package: GenomicRanges

## Loading required package: GenomeInfoDb

## Loading required package: SummarizedExperiment

## 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")'.

## Loading required package: DelayedArray

## Loading required package: matrixStats

## Warning: package 'matrixStats' was built under R version 3.6.1

## 
## Attaching package: 'matrixStats'

## The following objects are masked from 'package:Biobase':
## 
##     anyMissing, rowMedians

## Loading required package: BiocParallel

## Warning: package 'BiocParallel' was built under R version 3.6.1

## 
## Attaching package: 'DelayedArray'

## The following objects are masked from 'package:matrixStats':
## 
##     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges

## The following objects are masked from 'package:base':
## 
##     aperm, apply, rowsum

Seteamos el directorio de trabajo en el lugar donde tenemos el archivo de cuantificación

setwd("~/iBio/workshop2019/Sesion5_RNAseq/cyverse/At")

DESeq2 requiere el conteo de reads “crudo”, no RPKM, no FPKM ni ninguna variante.

The DESeq2 model internally corrects for library size, so transformed or normalized values such as counts scaled by library size should not be used as input.

Utilizaremos el archivo que obtuvimos al final del procedimiento de mapeo con “salmon”… ya sea usando los archivos que poseen todos los reads o los archivos “disminuídos” a 1 millón de reads (1M).

Leemos el archivo de cuantificación (At.1M.separated.quants.txt) y seteamos un nuevo nombre para las columnas

counts<-data.frame(read.table("At.1M.separated.quants.txt",row.names = 1,header = T,sep = "\t"))
colnames(counts)<-c("C1","C2","C3","T1","T2","T3")
head(counts,2)

##              C1  C2  C3  T1  T2 T3
## AT1G51370.2   0   0   1   3   5  5
## AT1G50920.1 134 124 117 104 111 99

Creamos un vector para describir a qué condición corresponde cada muestra

condition <- factor(rep(c("control", "treated"), each = 3))
condition

## [1] control control control treated treated treated
## Levels: control treated

En caso que nuestras muestras pertenezcan a un diseño multifactorial, podemos generar nuevos vectores que contengan dicha información.Por ejemplo, si algunas de las muestras pertenecen distintos genotipos y queremos realizar el análisis DEG ponderando por genotipos podemos agregar la info con el siguiente comando

En el siguiente link encontrarán más información al respecto (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs)

Ojo que es sólo un ejemplo pues las muestras de prueba pertenen al mismo genotipo

genotype <- factor(c(rep("col0", each = 3),rep("col1", each = 3)))
genotype

## [1] col0 col0 col0 col1 col1 col1
## Levels: col0 col1

Consolidamos la información de nuestras muestras en la variable targets.

targets <- data.frame(condition,genotype)
targets

##   condition genotype
## 1   control     col0
## 2   control     col0
## 3   control     col0
## 4   treated     col1
## 5   treated     col1
## 6   treated     col1

Leemos la tabla y aplicamos la función DESeq

dds=DESeqDataSetFromMatrix(counts,targets,design=~condition)
dds <- DESeq(dds)

## estimating size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## -- note: fitType='parametric', but the dispersion trend was not well captured by the
##    function: y = a/x + b, and a local regression fit was automatically substituted.
##    specify fitType='local' or 'mean' to avoid this message next time.

## final dispersion estimates

## fitting model and testing

Revisamos que los counts crudos ahora son parte del objeto dds

head(counts(dds))

##              C1  C2  C3  T1  T2 T3
## AT1G51370.2   0   0   1   3   5  5
## AT1G50920.1 134 124 117 104 111 99
## AT1G36960.1   0   0   0   0   0  0
## AT1G44020.1   2   2   0   3   2  1
## AT1G15970.1  13   9  10  17  12 13
## AT1G73440.1   8   8   7  11   5  3

Además tenemos counts normalizados

head(counts(dds,normalized=T))

##                     C1         C2         C3        T1         T2
## AT1G51370.2   0.000000   0.000000   1.015264  2.852986   5.059133
## AT1G50920.1 126.428561 118.104858 118.785920 98.903506 112.312745
## AT1G36960.1   0.000000   0.000000   0.000000  0.000000   0.000000
## AT1G44020.1   1.886993   1.904917   0.000000  2.852986   2.023653
## AT1G15970.1  12.265457   8.572127  10.152643 16.166919  12.141918
## AT1G73440.1   7.547974   7.619668   7.106850 10.460948   5.059133
##                     T3
## AT1G51370.2   5.110874
## AT1G50920.1 101.195301
## AT1G36960.1   0.000000
## AT1G44020.1   1.022175
## AT1G15970.1  13.288272
## AT1G73440.1   3.066524

Si por alguna razón desea grabar los valores de expresión crudos o los valores normalizados puede usar los siguientes

write.csv(counts(dds), file = "quants_normalized.csv")
write.csv(counts(dds,normalized=T), file = "quants_normalized.csv")

2.1 Visualización de datos

2.1.1 Boxplot

Podemos observar los perfiles de expresión de cada librería pre y post normalización

par(mfrow=c(1,2))
boxplot(log2(counts(dds, normalized=FALSE)+1), main="Raw counts", col=rep(c("lightblue","lightgreen"), each=3), ylim = c(0, 15),las=2)
boxplot(log2(counts(dds, normalized=TRUE)+1), main="Normalized counts", col=rep(c("lightblue","lightgreen"), each=3), ylim = c(0, 15),las=2)

2.1.2 PCAplot

pendiente

2.1.3 Heatmap

pendiente

3. Comparaciones entre muestras

DESeq automáticamente realizará las posibles comparaciones entre las librerías según el diseño que le pasamos como argumento. En este caso el diseño que le pedí fue “design=~condition”, y como sólo existen dos condiciones (control y treated), sólo hará esa comparación. En caso que hubieran “condition” tuviera 3 condiciones (control, treated1 y treated2), va a realizar la comparación “treated1_vs_control” y “treated2_vs_control” además de un resultado “común” que denomina “Intercept”.

Para revisar las comparaciones posibles ocupamos el siguiente comando

resultsNames(dds)

## [1] "Intercept"                    "condition_treated_vs_control"

Como era de esperar, la única comparación disponible es “condition_treated_vs_control”

Como la variable “dds” contiene diferentes comparaciones, es necesario extraer la comparación de interés para proceder con los siguientes análisis. Guardaré en la variable “res” el resultado de la comparación “treated_vs_control”

res <- results(dds, name="condition_treated_vs_control")

Resumen de datos del objeto “res”

summary(res)

## 
## out of 24811 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 684, 2.8%
## LFC < 0 (down)     : 483, 1.9%
## outliers [1]       : 58, 0.23%
## low counts [2]     : 7092, 29%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Primeras líneas del objeto “res”

head(res)

## log2 fold change (MLE): condition treated vs control 
## Wald test p-value: condition treated vs control 
## DataFrame with 6 rows and 6 columns
##                     baseMean     log2FoldChange             lfcSE
##                    <numeric>          <numeric>         <numeric>
## AT1G51370.2 2.33970941669419   3.64019320886484  1.57207668466735
## AT1G50920.1 112.621815203343 -0.219027281001241 0.146200722103832
## AT1G36960.1                0                 NA                NA
## AT1G44020.1 1.61512068170417  0.627809053853954   1.4267505610525
## AT1G15970.1 12.0978894067029  0.426384099590713 0.471639948169337
## AT1G73440.1 6.81018277238528 -0.252147570834862 0.716498225603596
##                           stat             pvalue              padj
##                      <numeric>          <numeric>         <numeric>
## AT1G51370.2   2.31553157957757 0.0205838583757382                NA
## AT1G50920.1  -1.49812721749546  0.134100200854798   0.6404390609239
## AT1G36960.1                 NA                 NA                NA
## AT1G44020.1  0.440027199562357   0.65991740761844                NA
## AT1G15970.1  0.904045768908498  0.365971132749654  0.89704787042123
## AT1G73440.1 -0.351916532134392  0.724900860457775 0.997056853062986

3.1 Visualización de datos: Volcano

Una forma de visualizar los resultados de una comparación tipo “tratamiento_vs_control” para todos los genes es graficar el valor p ajustado (padj) vs el nivel de cambio (Fold of Change [FC]). Este tipo de gráfico es conocido como “Volcano” y permite visualizar los rangos entre los cuales se mueven dichas variables para estimar posibles puntos de corte. (https://en.wikipedia.org/wiki/Volcano_plot_(statistics))

Vamos a realizar un gráfico “volcano” con los datos de “res”

Primero seteamos algunos puntos de corte típicos para el gráfico

Iré paso a paso para que puedan entender cómo se configura este gráfico y puedan realiza cambios si lo requieren

Punto de corte máximo de padj (los genes que cumplan aparecerán como puntos rojos)

alpha=0.05

Punto de mínimo corte para el fold of change (los genes que cumplan aparecerán como puntos rojos)

FC=4

Punto de corte para el eje del padj (eje Y)

padjlim=NULL

Punto de corte para el eje del FC (ele X)

logFClim=7

Ajusto los valores padj=“0” a un valor muy cercano a cero. Esto previene problemas de escala.

res$padj[which(res$padj==0)] <- .Machine$double.xmin

Calcular el -log10 del padj para una mejor representación gráfica

log10pval <- -log10(res$padj)

Seteo el límite del “eje y” en la representación 99% de los datos a menos que haya seteado otro valor la valiable “padjlim”

if (is.null(padjlim)){
  ylim <- c(0,1) * quantile(log10pval, probs=0.99, na.rm=TRUE)
} else{
  ylim <- c(0,1) * -log10(padjlim)
}

Seteo el límite del “eje x” revisando si ha sido seteado en la variable logFClim

if (is.null(logFClim)){
  xlim <- NULL
  #xlim <- c(-1,1) * quantile(res$log2FoldChange, probs=0.999, na.rm=TRUE)
}else{
  xlim <- c(-1,1)*logFClim
}

Finalmente, creamos el gráfico

plot(res$log2FoldChange, pmin(ylim[2], log10pval), ylim=ylim, xlim=xlim, las=1, cex=0.45,
xlab=expression(log[2]~fold~change), ylab=expression(-log[10]~adjusted~P~value),
col=ifelse(res$padj <= alpha & (res$log2FoldChange <= -log2(FC) | res$log2FoldChange >= log2(FC)), "red", "black"), pch=ifelse(log10pval >= ylim[2], 2, 20),
main=paste0("Volcano plot - ",gsub(".*: ","",res@elementMetadata@listData[["description"]][2])))
abline(h=-log10(alpha), col="black", lty=4, lwd=2.0)
abline(v=0, col="black", lty=3, lwd=1.0)
abline(v=-log2(FC), col="black", lty=4, lwd=2.0)
abline(v=log2(FC), col="black", lty=4, lwd=2.0)

Los puntos rojos representan genes que “cumplen” con las condiciones de corte escogidas.

Las líneas representan los cuadrantes “límite” según las condiciones escogidas.

Los triángulos representan genes que poseen un pvalue por sobre el punto de corte escogido.

3.2 Selección de genes según puntos de corte

Vamos a revisar el número de genes que cumplen con distintos puntos de corte

3.2.1 Punto de corte sólo de pvalue < 0.05

* genes con pvalue <0.05

sum(res$padj<0.05, na.rm = T)

## [1] 897

* genes con pvalue < 0.05 y nivel de cambio mayor que 0 ( log2(X) > 0 )

sum(res$padj<0.05 & res$log2FoldChange>0, na.rm = T)

## [1] 542

* genes con pvalue < 0.05 y nivel de cambio menor que 0 ( log2(X) < 0 )

sum(res$padj<0.05 & res$log2FoldChange<0, na.rm = T)

## [1] 355

3.2.2 Punto de corte de pvalue < 0.05 y FC > 2

* genes con pvalue < 0.05 y nivel de cambio absoluto mayor que ( |log2(X)| > 1 )

sum(res$padj<0.05 & abs(res$log2FoldChange)>=0, na.rm = T)

## [1] 897

* genes con pvalue < 0.05 y nivel de cambio mayor que 2 ( log2(X) > 1 )

sum(res$padj<0.05 & res$log2FoldChange>=1, na.rm = T)

## [1] 409

* genes con pvalue < 0.05 y nivel de cambio meno que 2 ( log2(X) < -1 )

sum(res$padj<0.05 & res$log2FoldChange<=-1, na.rm = T)

## [1] 228

3.2.3 Selección de IDs

En cualquiera de los ejemplos anteriores, para obtener el listado de genes modificamos el comando utilizado de la siguiente manera

* Ejemplo para genes up

up_reg <- row.names(subset(res, res$padj < 0.05 & res$log2FoldChange >= 1))
head(up_reg)

## [1] "AT1G29020.1" "AT1G65610.1" "AT1G17380.1" "AT1G65690.1" "AT1G73220.1"
## [6] "AT1G17420.1"

* Ejemplo para genes down

down_reg <- row.names(subset(res, res$padj < 0.05 & res$log2FoldChange <=- 1))
head(down_reg)

## [1] "AT1G74770.1" "AT1G70850.1" "AT1G74940.1" "AT1G73600.2" "AT1G70890.1"
## [6] "AT1G69160.1"

Y con el siguiente comando podemos guardar el listado de IDs para usarlo, por ejemplo, en Cytoscape

write.table( up_reg, file = "up_reg.txt", col.names = FALSE, row.names = FALSE, quote = FALSE)
write.table( down_reg, file = "down_reg.txt", col.names = FALSE, row.names = FALSE, quote = FALSE)

FIN… por ahora

DESeq2_Intro

J.Maldonado

20-08-2019

INTRODUCCIÓN a DEseq2 en R/Rstudio

1. Introducción

En el siguiente link encontrarán distintas explicaciones sobre las diferencias de los métodos https://www.biostars.org/p/284775/

En este manual utilizaré el paquete DESeq2, el cual tiene un desarrollo más reciente y se está posicionando rápidamente como un estándar en el área.

Toda la información detallada sobre este paquete la pueden encontrar en (http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)

2. Carga de datos en R

Iniciamos el trabajo instalando los paquetes requeridos (en caso que no estén instalados)

Cargamos en memoria los paquetes requeridos (Se debe ejecutar cada vez que se inicia RStudio)

Seteamos el directorio de trabajo en el lugar donde tenemos el archivo de cuantificación

DESeq2 requiere el conteo de reads “crudo”, no RPKM, no FPKM ni ninguna variante.

Utilizaremos el archivo que obtuvimos al final del procedimiento de mapeo con “salmon”… ya sea usando los archivos que poseen todos los reads o los archivos “disminuídos” a 1 millón de reads (1M).

Leemos el archivo de cuantificación (At.1M.separated.quants.txt) y seteamos un nuevo nombre para las columnas

Creamos un vector para describir a qué condición corresponde cada muestra

Ojo que es sólo un ejemplo pues las muestras de prueba pertenen al mismo genotipo

Consolidamos la información de nuestras muestras en la variable targets.

Leemos la tabla y aplicamos la función DESeq

Revisamos que los counts crudos ahora son parte del objeto dds

Además tenemos counts normalizados

Si por alguna razón desea grabar los valores de expresión crudos o los valores normalizados puede usar los siguientes

2.1 Visualización de datos

2.1.1 Boxplot

Podemos observar los perfiles de expresión de cada librería pre y post normalización

2.1.2 PCAplot

2.1.3 Heatmap

3. Comparaciones entre muestras

Para revisar las comparaciones posibles ocupamos el siguiente comando

Como era de esperar, la única comparación disponible es “condition_treated_vs_control”

Como la variable “dds” contiene diferentes comparaciones, es necesario extraer la comparación de interés para proceder con los siguientes análisis. Guardaré en la variable “res” el resultado de la comparación “treated_vs_control”

Resumen de datos del objeto “res”

Primeras líneas del objeto “res”

3.1 Visualización de datos: Volcano

Vamos a realizar un gráfico “volcano” con los datos de “res”

Primero seteamos algunos puntos de corte típicos para el gráfico

Iré paso a paso para que puedan entender cómo se configura este gráfico y puedan realiza cambios si lo requieren

Punto de corte máximo de padj (los genes que cumplan aparecerán como puntos rojos)

Punto de mínimo corte para el fold of change (los genes que cumplan aparecerán como puntos rojos)

Punto de corte para el eje del padj (eje Y)

Punto de corte para el eje del FC (ele X)

Ajusto los valores padj=“0” a un valor muy cercano a cero. Esto previene problemas de escala.

Calcular el -log10 del padj para una mejor representación gráfica

Seteo el límite del “eje y” en la representación 99% de los datos a menos que haya seteado otro valor la valiable “padjlim”

Seteo el límite del “eje x” revisando si ha sido seteado en la variable logFClim

Finalmente, creamos el gráfico

Los puntos rojos representan genes que “cumplen” con las condiciones de corte escogidas.

Las líneas representan los cuadrantes “límite” según las condiciones escogidas.

Los triángulos representan genes que poseen un pvalue por sobre el punto de corte escogido.

3.2 Selección de genes según puntos de corte

Vamos a revisar el número de genes que cumplen con distintos puntos de corte

3.2.1 Punto de corte sólo de pvalue < 0.05

* genes con pvalue <0.05

* genes con pvalue < 0.05 y nivel de cambio mayor que 0 ( log2(X) > 0 )

* genes con pvalue < 0.05 y nivel de cambio menor que 0 ( log2(X) < 0 )

3.2.2 Punto de corte de pvalue < 0.05 y FC > 2

* genes con pvalue < 0.05 y nivel de cambio absoluto mayor que ( |log2(X)| > 1 )

* genes con pvalue < 0.05 y nivel de cambio mayor que 2 ( log2(X) > 1 )

* genes con pvalue < 0.05 y nivel de cambio meno que 2 ( log2(X) < -1 )

3.2.3 Selección de IDs

En cualquiera de los ejemplos anteriores, para obtener el listado de genes modificamos el comando utilizado de la siguiente manera

* Ejemplo para genes up

* Ejemplo para genes down

Y con el siguiente comando podemos guardar el listado de IDs para usarlo, por ejemplo, en Cytoscape