I. Objetivos

  • Conocer el flujo de trabajo empleado en el análisis bioinformático de datos obtenidos por RNA-seq
  • Explicar el análisis de expresión diferencial empleando DESeq2
  • Analizar un conjunto de datos de expresión empleando un script de DESEq2

II. Requisitos

Esta práctica está basada en el flujo de trabajo o pipeline de Michael Love por lo que requerirémos de:

Datos

Los datos fueron descargados de un proyecto público titulado NRF2 signalling associated with radioresistance in colorectal cancer models de la base de datos GEO. El número de acceso del trabajo es GSE136011. Los datos que requeriremos para la práctico serán:

  • En formato FastQ para iniciar la prácitca desde el inicio

  • Las cuentas crudas generadas por Salmon ubicadas en este vínculo

  • Los objetos de R para iniciar el análisis de expresión diferencial. Se ubican en el siguiente vínculo

Programas

Shell

Descarga e instala los programas empleando el gestor de paquetes de python Conda. Sigue las instrucciones específicas para tu sistema operativo desde la opción Regular installation. Una vez que hayas descargado e instalado correctamente conda, prueba el siguiente comando:

conda list

Instala los siguientes programas:

  • FastQC
  • MultiQC
  • Cutadapt
  • Salmon

Corriendo el comando:

conda install -c  bioconda <paquete>

R

Los paquetes de R que vas a emplear en la práctica se descargarán desde los repositiorios de Bioconductor y CRAN. Para ello asegurate que tienes instalada la última versión de R y la IDE de RStudio. Una vez configurado correctamente R, procede a instalar la última versión de Bioconductor corriendo el siguiente comando:

install.packages("BiocManager")
BiocManager::install(version = 3.13)

Instalado Bioconductor, descarga e instala los siguientes paquetes:

  • tximeta
  • DESeq2
  • PCAtools
  • apeglm
  • Glimma
  • clusterProfiler
  • enrichplot
  • biomaRt
  • fgsea

Corriendo el siguiente comando:

BiocManager::install("paquete")

Los paquetes que descargaremos desde CRAN son:

  • tidyverse
  • RColorBrewer
  • pheatmap
  • here

Corriendo el comando:

install.packages("paquete")

III. Análisis bioinformático de los datos de secuenciación

Para el análisis bioinformático vamos a seguir este pipeline:

  1. Control de calidad de las secuencias crudas
  2. Filtrado y limpieza de las secuencias En los siguientes vínculos encontrarán los manuales de Trimmomatic y Cutadapt.
  3. Alineamineto Para una revisión detallada del algoritmo de Salmon, pueden revisar el artículo publicado por sus desarrolladores (dobinSTARUltrafastUniversal2013?). Asimismo la versión más reciente del manual.

Si decides realizar el ejercicio desde el inicio, empieza desde este punto. De lo contrario salta a la sección 4.2

1. Control de calidad de las secuencias crudas

  • Primer paso a realizar en los flujos de trabajo para el análisis de expresión diferencial. A este proceso también se le conoce como Quality Control o QC analysis

  • Permite verificar la calidad de las lecturas a distintos niveles:

    • Asignación de bases
    • Secuencias duplicadas
    • Presencia de contaminación específica o inespecífica

Pasos para correr FastQC

  1. Abre una ventana de la terminal y ubicate en el directorio bin/ del presente repositorio

  2. Activa tu ambiente de conda corriendo:

conda activate <ruta al directorio miniconda2 o minicoda3>
  1. Comprueba que el programa FastQC se encuentra correctamente instalado:
fastqc -h
  1. Analiza la calidad de las secuencias crudas:
fastqc ../data/*.fastq.gz -o ../Output/FastQC/

Para concatenar o juntar los resultados de todas las muestras analizadas por FastQC, usemos MultiQC

  1. Corrobora que MultiQC se encuentra correctamente instalado:
multiqc -h
  1. Ejecuta el programa para analizar los archivos con terminación .zip generados por FastQC
multiqc ../Output/FastQC/*.zip -o ../Output/MultiQC/

Explora y revisa los archivos generados por FastQC y MultiQC

QC

Importante revisar si los problemas detectados por FastQC en la calidad de las secuencias se presentan en una, varias o todas las muestras porque de esto dependerán las correcciones a las lecturas crudas para mejorar su calidad.

2. Limpieza y corrección de la calidad de las lecturas

  • Proceso que depende de los resultados del QC
  • Permite:
    • Filtrar las lecturas para seleccionar aquellas que tengan mayor calidad
    • Recortar lecturas para preservar los fragmentos de mayor calidad
    • Recortar lecturas para eliminar secuencias contaminantes (adaptadores)

Pasos para correr Cutadapt:

  1. De manera similar al análisis anterior, asegura que te encuentras en el directorio bin/, tienes activo el ambiente de conda y cutadapt se encuentra correctamente instalado:
cutadapt -h
  1. Ejecuta cutadapt para una sola muestra:
cutadapt -a AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT ##Para el read forward
-A AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT ##Para el read reverse
-o ../Output/Cutadapt/samp1_clean_R1.fastq.gz ##Archivo de salida para las lecturas forward
-p ../Output/Cutadapt/samp1_clean_R2.fastq.gz ##Archivo de salida para las lecturas reverse
../data/samp1_R1.fastq.gz ##Ruta hacia las lecturas crudas forward
../data/samp1_R2.fastq.gz ##Ruta hacia las lecturas crudas reverse

Para ejecutar cutadapt en todas las muestras abre el archivo titulado 2.Read_cleaning en el directorio bin/ y encotrarás el código anidado en un bucle o for loop. Desde la terminal corre:

bash 2.Read_cleaning

Al finalizar el filtrado y recorte o triming de la secuencias crudas, vuelve a analizar la calidad con fastQC y MultiQC de las lecturas limpias.

3. Alineamiento, mapeo y estimación de la abundancia de las lecturas

  • Proceso necesario para conocer la región del genoma o transcriptoma a partir de la cual se generaron las lecturas
  • Permita anotar las lecturas y categorizarlas como: mRNA, lncRNA, ncRNA, miRNA, …
  • Estimación del número de lecturas a nivel de gen, isoforma, exón, …

Aunque Salmon es un pseudoalineador y no requiere de mayor uso de recursos computacionales comparado con STAR, HISAT2 o Tophat2, es recomendable que el proceso lo realices en un cluster para que lo paralelices.

Pasos para alinear con salmon: 1. Generar el índice del transcriptoma (humano) empleando el archivo fasta descargado desde la página de GENECODE (sección de Fasta Files).

  1. Correr el siguiente código para generar el índice:
salmon index -t ../reference/Homo_sapiens.GRCh38.cdna.all.fa ##Ruta hacia el archivo fasta del transcriptoma de referencia
-i ../reference/hsa_v38_gencode ##Ruta para almacenar los archivos generados del índice
  1. Una vez generado el índice, realizar el pseudoalineamiento para una muestra:
salmon quant -i ../reference/hsa_v38_gencode ##Ruta hacia el índice
-l A ##Opción para permitir a salmon deducir de forma automática las características de la libreria
-1 ../Output/Cutadapt/samp1_clean_R1.fastq.gz ##Ruta a las lecturas forward
-2 ../Output/Cutadapt/samp1_clean_R2.fastq.gz ##Ruta a las lecturas reverse
-p 8 ##Número de cores para paralelizar el proceso
--validateMappings -o ##Opción  para permitir a salmon validar los archivos bam generados de forma temporal
../Output/salmon_quants/samp1_quant ##Ruta para guardar las cuentas crudas generadas por salmon

Si deseas realizar este proceso para todas las muestras, corre el siguiente código:

bash 3.Pseudoalignment

4. Análisis de expresión diferencial

4.1 Importe de archivos a R

Al finalizar el proceso de estimación de la abundancia, Salmon nos devuelve los siguientes archivos por muestra:

ls ../results/salmon_quants/
## metadata.txt
## samp10_quant
## samp11_quant
## samp12_quant
## samp1_quant
## samp2_quant
## samp3_quant
## samp4_quant
## samp5_quant
## samp6_quant
## samp7_quant
## samp8_quant
## samp9_quant

Se requieren importar los datos generados quant.sf de cada muestra a R. Para ello vamos a utilizar la librería tximeta [Cita]. Esta librería importa archivos de cuentas:

  • Generadas por salmon
  • Preserva metadatos asociados al proceso de pseudoalineamiento (versión del transcriptoma de referencia, código empleado, versión del pseudoalineador, …)

Para importar los datos a R vamos a requerir de una tabla de metadatos en formato de texto plano

list.files(here("results/salmon_quants/"))
##  [1] "metadata.txt" "samp1_quant"  "samp10_quant" "samp11_quant" "samp12_quant"
##  [6] "samp2_quant"  "samp3_quant"  "samp4_quant"  "samp5_quant"  "samp6_quant" 
## [11] "samp7_quant"  "samp8_quant"  "samp9_quant"
  1. Cargar a R la tabla de metadatos:
coldata <- read.table(here("results/salmon_quants/metadata.txt"), he = T, stringsAsFactors = T, sep = "\t")

rownames(coldata) <- coldata$unique_id

head(coldata)
##                                                                                   file_name
## siRNA_NRF2_1        SRR10004274_GSM4039376_NRF2_siRNA_rep_1_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_control_1 SRR10004275_GSM4039377_NRF2_siRNA_cont_rep1_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_NRF2_2        SRR10004276_GSM4039378_NRF2_siRNA_rep_2_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_control_2 SRR10004277_GSM4039379_NRF2_siRNA_cont_rep2_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_NRF2_3        SRR10004278_GSM4039380_NRF2_siRNA_rep_3_Homo_sapiens_RNA-Seq_1.fastq.gz
## siRNA_control_3 SRR10004279_GSM4039381_NRF2_siRNA_cont_rep3_Homo_sapiens_RNA-Seq_1.fastq.gz
##                     Treatment Replicate       unique_id   key Geno
## siRNA_NRF2_1       siRNA_NRF2         1    siRNA_NRF2_1 samp1   WT
## siRNA_control_1 siRNA_control         1 siRNA_control_1 samp2   WT
## siRNA_NRF2_2       siRNA_NRF2         2    siRNA_NRF2_2 samp3   WT
## siRNA_control_2 siRNA_control         2 siRNA_control_2 samp4   WT
## siRNA_NRF2_3       siRNA_NRF2         3    siRNA_NRF2_3 samp5   WT
## siRNA_control_3 siRNA_control         3 siRNA_control_3 samp6   WT
  1. Crear una nueva columna a la tabla que incluya la ruta hacia los archivos de las cuentas crudas generados por salmon:
coldata <- dplyr::mutate(coldata, files = file.path(here("results/salmon_quants"), paste0(coldata$key, "_quant"), "quant.sf"))
  1. Verificar que las rutas hacia los archivos de las cuentas fue correctamente creada:
coldata <- mutate(coldata, file_exist = file.exists(coldata$files))

coldata[, 2:8]
##                     Treatment Replicate       unique_id    key Geno
## siRNA_NRF2_1       siRNA_NRF2         1    siRNA_NRF2_1  samp1   WT
## siRNA_control_1 siRNA_control         1 siRNA_control_1  samp2   WT
## siRNA_NRF2_2       siRNA_NRF2         2    siRNA_NRF2_2  samp3   WT
## siRNA_control_2 siRNA_control         2 siRNA_control_2  samp4   WT
## siRNA_NRF2_3       siRNA_NRF2         3    siRNA_NRF2_3  samp5   WT
## siRNA_control_3 siRNA_control         3 siRNA_control_3  samp6   WT
## CT_1                       CT         1            CT_1  samp7   WT
## mut_NRF2_1           mut_NRF2         1      mut_NRF2_1  samp8  Mut
## CT_2                       CT         2            CT_2  samp9   WT
## mut_NRF2_2           mut_NRF2         2      mut_NRF2_2 samp10  Mut
## CT_3                       CT         3            CT_3 samp11   WT
## mut_NRF2_3           mut_NRF2         3      mut_NRF2_3 samp12  Mut
##                                                                                                       files
## siRNA_NRF2_1     /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp1_quant/quant.sf
## siRNA_control_1  /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp2_quant/quant.sf
## siRNA_NRF2_2     /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp3_quant/quant.sf
## siRNA_control_2  /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp4_quant/quant.sf
## siRNA_NRF2_3     /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp5_quant/quant.sf
## siRNA_control_3  /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp6_quant/quant.sf
## CT_1             /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp7_quant/quant.sf
## mut_NRF2_1       /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp8_quant/quant.sf
## CT_2             /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp9_quant/quant.sf
## mut_NRF2_2      /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp10_quant/quant.sf
## CT_3            /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp11_quant/quant.sf
## mut_NRF2_3      /Users/rodolfochavez/RNAseq/Tutorial_of_RNA_seq/results/salmon_quants/samp12_quant/quant.sf
##                 file_exist
## siRNA_NRF2_1          TRUE
## siRNA_control_1       TRUE
## siRNA_NRF2_2          TRUE
## siRNA_control_2       TRUE
## siRNA_NRF2_3          TRUE
## siRNA_control_3       TRUE
## CT_1                  TRUE
## mut_NRF2_1            TRUE
## CT_2                  TRUE
## mut_NRF2_2            TRUE
## CT_3                  TRUE
## mut_NRF2_3            TRUE
  1. tximeta requiere que exista una columna en la tabla de metadatos llamada names y que incluya el nombre de los directorios en donde se encuentra almacenadas las cuentas por cada muestra (samp1_quant)
coldata <- mutate(coldata, names = paste0(coldata$key, "_quant"))
  1. Si todas las rutas fueron creadas correctamente, importa los datos usando tximeta
se <- tximeta(coldata)
## importing quantifications
## reading in files with read_tsv
## 1 2 3 4 5 6 7 8 9 10 11 12
## Warning: DEPRECATION: As of BiocFileCache (>1.15.1), default caching location has changed.
##   Problematic cache: /Users/rodolfochavez/Library/Caches/BiocFileCache
##   See https://www.bioconductor.org/packages/devel/bioc/vignettes/BiocFileCache/inst/doc/BiocFileCache.html#default-caching-location-update
## found matching transcriptome:
## [ Ensembl - Homo sapiens - release 103 ]
## Warning: DEPRECATION: As of BiocFileCache (>1.15.1), default caching location has changed.
##   Problematic cache: /Users/rodolfochavez/Library/Caches/BiocFileCache
##   See https://www.bioconductor.org/packages/devel/bioc/vignettes/BiocFileCache/inst/doc/BiocFileCache.html#default-caching-location-update
## loading existing EnsDb created: 2021-03-09 20:52:43
## Warning: DEPRECATION: As of BiocFileCache (>1.15.1), default caching location has changed.
##   Problematic cache: /Users/rodolfochavez/Library/Caches/BiocFileCache
##   See https://www.bioconductor.org/packages/devel/bioc/vignettes/BiocFileCache/inst/doc/BiocFileCache.html#default-caching-location-update
## loading existing transcript ranges created: 2021-03-09 20:52:48

Ya que tximeta descarga información relacionada con el transcriptoma de referencia, este último paso requiere de una conexión estable a internet de alta velocidad por lo que puede llegar a tardar el proceso.

  1. El objeto se tiene clase de SummarizedExperiment y las cuentas están reportadas a nivel de transcrito (isoforma). Para obtener las cuentas a nivel de gen:
gse <- summarizeToGene(se)
## Warning: DEPRECATION: As of BiocFileCache (>1.15.1), default caching location has changed.
##   Problematic cache: /Users/rodolfochavez/Library/Caches/BiocFileCache
##   See https://www.bioconductor.org/packages/devel/bioc/vignettes/BiocFileCache/inst/doc/BiocFileCache.html#default-caching-location-update
## loading existing EnsDb created: 2021-03-09 20:52:43
## obtaining transcript-to-gene mapping from database
## loading existing gene ranges created: 2021-03-09 20:53:45
## summarizing abundance
## summarizing counts
## summarizing length

Explora el objeto gse (SummarizedExperiment)

##Acceder a la  tabla de metadatos
colData(gse)

##Acceder al objeto GRanges con la información de los rangos y anotación de cada gen
head(rowRanges(gse))

##Acceder a la matriz de cuantas crudas
head(assay(gse))

4.2 Cargando los archivos desde el objeto .rdata

Si has decidido comenzar con el ejercicio desde este punto, solo necesitas cargar el objeto SummarizedExperimentObj.rdata ubicado en el folder llamado data/

##Carga el archivo usando la función load
load("../data/SummarizedExperimentObj.rdata")

Verifica que en el ambiente de R se cargaron correctamente los siguientes objetos:

  • dds
  • gse
  • se

Exploración de los datos

Previamente al análisis de expresión diferencial, es necesario explorar los datos:

  • Conocer la estructura del set de datos
  • Detectar consistencia entre réplicas
  • Visualizar comportamientos atípicos (outliers)
  • Reportar las variables que separan o agrupan los datos

El análisis de componentes principales (PCA) nos permita reducir la dimensionalidad de nuestro set de datos y analizarlos de manera más intuitiva.

Los pasos para explorar el juego de datos son los siguientes:

  1. Crear un objeto DESeq (SummarizedExperiment) a partir del objeto gse
dds <- DESeqDataSet(se = gse, ##Objeto SummarizedExperiment que contiene la matriz de cuentas crudas
                    design = ~ Treatment) ##Variable con los grupos experimentales a contrastar. Forma parte de las columnas de coldata
## using counts and average transcript lengths from tximeta
  1. Visualizar el grupo de referencia o control de la variable seleccionada
##El primer nivel o elemento es el grupo control  ("CT")
levels(dds$Treatment)
## [1] "CT"            "mut_NRF2"      "siRNA_control" "siRNA_NRF2"
  1. Dependiendo de la pregunta a responder, recodifica el grupo de referencia:
dds$Treatment <- relevel(dds$Treatment, "siRNA_control")

levels(dds$Treatment)
## [1] "siRNA_control" "CT"            "mut_NRF2"      "siRNA_NRF2"
  1. Ya que los genes que muestran bajo número de cuentas (abundancia) no nos interesan en el análisis, selecciona aquellos que muestran mayor abundancia usando un punto de corte:
##Recupera los genes que muestran más de 1 cuenta en al menos 3 muestras
keep <- rowSums(counts(dds) >= 1) >=3

##Visualiza cuántos genes pasaron el filtro de abundancia
table(keep)
## keep
## FALSE  TRUE 
## 16846 21299
##Corta los genes que pasaron el filtro de abundancia del objeto dds
dds <- dds[keep, ]
  1. Para explorar los datos por medio del PCA necesitamos normalizar las cuentas crudas para que los datos sean más comparables entre sí. Empleamos la transformación regularized-logarithm transformation recomendada para juegos de datos con una n < 30:
rld <-  rlog(dds, blind = F)
## using 'avgTxLength' from assays(dds), correcting for library size
  1. Realiza el PCA empleando el paquete PCAtools:
##Indicarle a  la función pca que los datos deben de ser escalados para que sean más comparables entre sí
pca <- pca(assay(rld), metadata = colData(rld), scale = T)

##Graficar el resultado del PCA mostrando las etiquetas de acuerdo al nombr + tratamiento + réplica y coloreando los puntos de acuerdo al tipo de tratamiento
biplot(pca, lab = rownames(colData(rld)), colby = "Treatment")

  • ¿Cuál es la estructura del juego de datos?
  • ¿Cuál fue la variable que causó dicho separamiento?
  • ¿Existen datos outlier?
  • ¿La consistencia entre las réplicas es adecuada?

Si deseas generar una gráfica interactiva de un análisis de reducción de dimensiones, utiliza la funcion glimmaMDS() del paquete Glimma la cual realiza un escalamiento multidimensional del set de datos.

Análisis de expresión diferencial

El análisis de expresión diferencial en DESeq2 se encuentra completamente automatizado. Al correr la función DESeq(), el programa:

  • Emplea el método de las medianas de los ratios para calcular los factores de escalación (size-factors) y normalizar las cuentas crudas.
  • Estima la dispersión de cada gen (réplicas biológicas cobran importancia)
  • Ajusta los datos a un modelo binomial negativo linearizado (GLM)
  • Realiza la prueba de Wald para comparar los coeficientes asociados a cada gen (log2FC y padj)

El análisis de expresión diferencial en DESeq2 se lleva a cabo de la siguiente forma:

  1. Corre la función DESeq() para normalizar, estimar la dispersión y ajustar los datos al modelo binomial negativo.
dds  <- DESeq(dds)
## estimating size factors
## using 'avgTxLength' from assays(dds), correcting for library size
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
  1. Realiza la prueba estadística para obtener los genes diferencialmente expresados:
##Usemos la función results con los parámetros por default
res <- results(dds)

##Visualiza el objeto que contiene los resultados
res
## log2 fold change (MLE): Treatment siRNA NRF2 vs siRNA control 
## Wald test p-value: Treatment siRNA NRF2 vs siRNA control 
## DataFrame with 21299 rows and 6 columns
##                    baseMean log2FoldChange     lfcSE       stat      pvalue
##                   <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## ENSG00000000003 1497.757504      0.0982849 0.0902643   1.088857 2.76217e-01
## ENSG00000000419 1789.280526     -0.0107935 0.0702472  -0.153650 8.77885e-01
## ENSG00000000457  304.687311     -0.0243559 0.1151762  -0.211466 8.32523e-01
## ENSG00000000460  467.602306      0.5142010 0.1159843   4.433368 9.27721e-06
## ENSG00000000938    0.618937      1.3390032 2.8076874   0.476906 6.33429e-01
## ...                     ...            ...       ...        ...         ...
## ENSG00000288690    8.913454      0.0326943  0.663455  0.0492789    0.960697
## ENSG00000288694    0.258397      0.8179394  4.365460  0.1873661    0.851374
## ENSG00000288695   11.159112     -0.3769746  0.701013 -0.5377567    0.590745
## ENSG00000288698   23.403296      0.5109625  1.060520  0.4818036    0.629946
## ENSG00000288699    0.869679      1.5093784  3.007207  0.5019203    0.615724
##                        padj
##                   <numeric>
## ENSG00000000003 0.476825323
## ENSG00000000419 0.948127083
## ENSG00000000457 0.924197635
## ENSG00000000460 0.000121719
## ENSG00000000938          NA
## ...                     ...
## ENSG00000288690    0.980632
## ENSG00000288694          NA
## ENSG00000288695    0.768387
## ENSG00000288698    0.795601
## ENSG00000288699          NA
##Obten el número de genes sobre- y sub-expresados de acuerdo a la prueba estadística
summary(res)
## 
## out of 21299 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2271, 11%
## LFC < 0 (down)     : 2660, 12%
## outliers [1]       : 26, 0.12%
## low counts [2]     : 5778, 27%
## (mean count < 8)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Aspectos importantes al correr la función results() con los parámetros por default:

  • El contraste entre grupos experimentales de la variable tratamiento se realizó entre el factor re-codificado como control (siRNA_cotnrol) y el último factor de la variable (siRNA_NRF2), es decir, siRNA_NRF2 vs siRNA_control
  • El punto de corte de expresión (log2FoldChange) fue de 0
  • El punto de corte estadístico fue de padj < 0.1

…En otras palabras, estás recuperando genes que se expresan más de una o menos de una vez con respecto al control y aceptando 10 resultados falsos positivos por cada 100 pruebas…

¿Podemos modificar el contraste y los threshold de expresión y estadísticos?… ¡Sí!

Para modificar las variables a contrastar:

##Corre la función results names para conocer los contrastes disponibles
resultsNames(dds)
## [1] "Intercept"                            
## [2] "Treatment_CT_vs_siRNA_control"        
## [3] "Treatment_mut_NRF2_vs_siRNA_control"  
## [4] "Treatment_siRNA_NRF2_vs_siRNA_control"
##Elige el contraste de tu preferencia y especificalo en el argumento contrast de la función  res
results(dds, contrast = c("Treatment", "mut_NRF2", "siRNA_control"))
## log2 fold change (MLE): Treatment mut_NRF2 vs siRNA_control 
## Wald test p-value: Treatment mut NRF2 vs siRNA control 
## DataFrame with 21299 rows and 6 columns
##                    baseMean log2FoldChange     lfcSE      stat      pvalue
##                   <numeric>      <numeric> <numeric> <numeric>   <numeric>
## ENSG00000000003 1497.757504      -0.904490 0.0925777  -9.77006 1.51357e-22
## ENSG00000000419 1789.280526      -0.463760 0.0712266  -6.51105 7.46273e-11
## ENSG00000000457  304.687311      -0.132462 0.1169467  -1.13267 2.57352e-01
## ENSG00000000460  467.602306       1.197748 0.1136246  10.54128 5.57374e-26
## ENSG00000000938    0.618937       0.569104 2.8961996   0.19650 8.44219e-01
## ...                     ...            ...       ...       ...         ...
## ENSG00000288690    8.913454      -3.275897  0.903258 -3.626757 0.000287003
## ENSG00000288694    0.258397      -1.020007  4.407356 -0.231433 0.816978494
## ENSG00000288695   11.159112      -0.264531  0.698382 -0.378777 0.704853537
## ENSG00000288698   23.403296       1.095652  1.055949  1.037600 0.299456461
## ENSG00000288699    0.869679       0.823637  3.075925  0.267769 0.788877177
##                        padj
##                   <numeric>
## ENSG00000000003 8.25857e-22
## ENSG00000000419 2.39298e-10
## ENSG00000000457 3.36544e-01
## ENSG00000000460 3.44774e-25
## ENSG00000000938          NA
## ...                     ...
## ENSG00000288690 0.000587388
## ENSG00000288694          NA
## ENSG00000288695 0.771544281
## ENSG00000288698 0.383276923
## ENSG00000288699 0.839628917

Para modificar los threshold de expresión y estadístico:

##Emplea los argumentos lfcthreshold y alpha para especificar ambos puntos de corte de forma simultánea
res_LFC <- results(dds, lfcThreshold = 1, alpha = 0.05)

res_LFC
## log2 fold change (MLE): Treatment siRNA NRF2 vs siRNA control 
## Wald test p-value: Treatment siRNA NRF2 vs siRNA control 
## DataFrame with 21299 rows and 6 columns
##                    baseMean log2FoldChange     lfcSE      stat    pvalue
##                   <numeric>      <numeric> <numeric> <numeric> <numeric>
## ENSG00000000003 1497.757504      0.0982849 0.0902643  0.000000  1.000000
## ENSG00000000419 1789.280526     -0.0107935 0.0702472  0.000000  1.000000
## ENSG00000000457  304.687311     -0.0243559 0.1151762  0.000000  1.000000
## ENSG00000000460  467.602306      0.5142010 0.1159843  0.000000  1.000000
## ENSG00000000938    0.618937      1.3390032 2.8076874  0.120741  0.903896
## ...                     ...            ...       ...       ...       ...
## ENSG00000288690    8.913454      0.0326943  0.663455  0.000000  1.000000
## ENSG00000288694    0.258397      0.8179394  4.365460  0.000000  1.000000
## ENSG00000288695   11.159112     -0.3769746  0.701013  0.000000  1.000000
## ENSG00000288698   23.403296      0.5109625  1.060520  0.000000  1.000000
## ENSG00000288699    0.869679      1.5093784  3.007207  0.169386  0.865493
##                      padj
##                 <numeric>
## ENSG00000000003         1
## ENSG00000000419         1
## ENSG00000000457         1
## ENSG00000000460         1
## ENSG00000000938         1
## ...                   ...
## ENSG00000288690         1
## ENSG00000288694         1
## ENSG00000288695         1
## ENSG00000288698         1
## ENSG00000288699         1
summary(res_LFC)
## 
## out of 21299 with nonzero total read count
## adjusted p-value < 0.05
## LFC > 1.00 (up)    : 7, 0.033%
## LFC < -1.00 (down) : 8, 0.038%
## outliers [1]       : 26, 0.12%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Lo anterior nos indica que hay 15 genes que muestran expresión diferencial significativa (padj < 0.05) expresados 2 veces más o 0.5 veces menos respecto al control…

Un procedimiento muy común es realizar la prueba estadística con un lfcthreshold = 0 y posterioremente obtener (subset) los genes que muestran un log2FoldChange > |1| y reportarlos como diferencialmente expresados. Este procedimiento es incorrecto, ya que se el valor de padj pierde todo su significado y corremos el riesgo de reportar resultados falsos positivos

##Del objeto res (lfcthreshold = 0 y alpha = 0.1) filtra los genes que tienen un lfc > |1| y padj < 0.05
as.data.frame(results(dds, lfcThreshold = 0, alpha = 0.05)) %>% dplyr::filter(log2FoldChange > 1 & padj < 0.05 | log2FoldChange < -1 & padj < 0.05) %>%
  nrow()

##Compara los resultados con el objeto res_LFC
as.data.frame(res_LFC) %>% dplyr::filter(log2FoldChange > 1 & padj < 0.05 | log2FoldChange < -1 & padj < 0.05) %>%
  nrow()

Visualización de los resultados

MA-plot

Es una gráfica útil para visualizar la distribución de los genes que se encuentran diferencialmente expresados tomando en cuenta:

  • log2FoldChange en el eje “Y” (M de minus o contraste)
  • El promedio de las cuentas normalizadas de cada gen en todas las muestras analizadas en el eje “X” (A de average)
##Corre la función plotMA sobre el objeto de los resultados
plotMA(res)

Nota: Existe una sobre-estimación del log2FoldChange para genes poco abundantes (media de cuentas normalizadas baja). Para corregir la sobre-estimación del log2FoldChange para dichos genes sin afectar los valores de aquellos con mayor abundancia usamos la función lfcShrink() de la librería apeglm

res_shrink <- lfcShrink(dds = dds, res = res, coef = "Treatment_siRNA_NRF2_vs_siRNA_control"
, type = "apeglm")
## using 'apeglm' for LFC shrinkage. If used in published research, please cite:
##     Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
##     sequence count data: removing the noise and preserving large differences.
##     Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
summary(res_shrink)
## 
## out of 21299 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2271, 11%
## LFC < 0 (down)     : 2660, 12%
## outliers [1]       : 26, 0.12%
## low counts [2]     : 5778, 27%
## (mean count < 8)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
plotMA(res_shrink,  ylim = c(-3, 3))

Volcano plot

De manera similar al MA-plot, el volcano plot muestra la distribución de los genes diferencialmente expresados graficando el -log10(padj) en el eje “Y” y el log2FoldChange en el eje “X.” Usa la función personalizada volcanoplotR del script de funciones:

##En los argumentos logfc y padj utiliza los mismos valores empleados en la función results
res_shrink <- as.data.frame(res_shrink)
volcanoplotR(res_shrink, logfc = 0, p.adj = 0.1)

Utiliza la función glimmaVolcano() para generar un volcano plot interactivo

Heatmap

El mapa de calor o Heatmap permite visualizar los genes diferencialmente expresados, a nivel de cuentas normalizadas, y el agrupamiento de las muestras con respecto a la expresión de dichos genes.

Para generar un Heatmap, sigue estos pasos:

  1. Obten la matriz de cuentas normalizadas:
norm_counts <- counts(dds, normalized = T)
  1. Guarda un data.frame con los resultados de los genes mostrando expresión diferencial significativa:
deg <- res_shrink %>% dplyr::filter(log2FoldChange > 0 & padj < 0.1 |
                                                     log2FoldChange < 0 & padj < 0.1)

head(deg)
##                  baseMean log2FoldChange      lfcSE       pvalue         padj
## ENSG00000000460  467.6023      0.4721205 0.11836789 9.277210e-06 0.0001217192
## ENSG00000001084 1889.3028     -0.2991140 0.09464064 4.368589e-04 0.0031869719
## ENSG00000001167 1068.3113      0.2413441 0.08978195 2.877301e-03 0.0149359373
## ENSG00000001461  564.0667     -0.2490246 0.08409560 1.186278e-03 0.0073086983
## ENSG00000001561  612.2146     -0.2536187 0.10372375 5.197697e-03 0.0240905329
## ENSG00000001617  996.1836      0.2917830 0.09054617 3.973259e-04 0.0029527889
  1. Usando la librería de pheatmap() obten la gráfica
##Genera una paleta de colores 
RedBu <- brewer.pal(n = 10, name = "RdBu")
RedBu <- rev(RedBu)

pheatmap(norm_counts[rownames(deg), 1:6], ##Corta los genes con expresión diferencial de la matriz normalizada
         scale = "row", ##Calcula los valores z por fila (gen)
         color = RedBu, ##Nombre de la paleta de colores
         show_rownames = F, ##No mostrar el nombre de los genes en la gráfica
         show_colnames = F, ##No mostrar el nombre de las muestras en la gráfica
         clustering_distance_rows = "euclidean", ##Especifíca la medida de la distancia entre filas
         clustering_distance_cols = "euclidean", ##O columnas
         clustering_method = "single", ##Método de agrupamiento jerárquico
         border_color = NA, 
         annotation_col = coldata[1:6, c(2, 6)]) ##Añadir información relacionada a las columnas)