## Checar el tiempo de cada chunk
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      # record the current time before each chunk
      now <<- Sys.time()
    } else {
      # calculate the time difference after a chunk
      res <- difftime(Sys.time(), now)
      # return a character string to show the time
      paste("Tiempo: ", round(res, digits = 2))
    }
  }
}))

knitr::opts_chunk$set(time_it = TRUE)

Descripción general

Este documento encontrarás el código que se ocupará en el mini curso de RNA-seq, el código usado para generar las slides las puedes encontrar en Github.

	Tema	Práctica	Tiempo (min)	Encargado
1	Introducción a RNA-seq	NA	5	AnaB
2	Organización del proyecto (folder de datos, fig, scripts	NA	5	Rodolfo
3	Flujo de trabajo para el análisis de RNA-seq	NA	5	AnaB
4	Tipos de archivos	NA	5	Rodolfo
5	Control de calidad de los datos	Si	10	AnaB
6	Alineamiento (Salmon)	Si	10	Rodolfo
7	tximeta, tximport	Si	10	AnaB
8	Objeto SummarizedExperiment en R	Si	10	Rodolfo
9	Normalización de los datos (RPKM, TPM)	Si	15	AnaB
10	Exploración de los datos	Si	15	Rodolfo
11	DE (analysis)	Si	15	AnaB
12	DE (visualization/resultados)	Si	15	Rodolfo

Requisitos

Está práctica está basada en Mike love vignette y ocuparemos los siguientes programas para analizar los datos del GSE152699.

Experimento

Número de acceso GEO: GSE152699

Tiempo: 0.02

Descargar los datos que ocuparemos en está práctica
FastQC
MultiQC (opcional)
Salmon (opcional)
R > 4.0
Rstudio
Bioconductor >= 3.12
Paquetería de R:

Analizar archivos crudos de RNAseq requiere poder de computo, si quieres hacer el análisis completo desde cero, necesitarás una computadora con al menos RAM 8 GB y espacio para almacenar los archivos crudos (.fastq), aún así todo el proceso puede tomar más tiempo del que tendremos en la práctica. Normalmente estos análisis se realizan en un cluster o servidor que suelen tener mayor capacidad que una laptop. Por lo que se verán los comandos para hacer el alineamiento en Salmon pero NO correremos este paso.

Archivos

Si decides realizar toda la práctica necesitarás estos archivos. Descargar el archivo en formato zip llamado Files.zip desde este vínculo y colocalo en el folder principal de este repositiorio. Descomprime el archivo y sigue las indicaciones de los instructores. En el archivo Files.zip encontrarás:

Archivos_fastq/: Folder que contiene los archivos en formato .fastq
Salmon_quants/: Folder que contiene las matrices de cuentas generadas por Salmon
genecode.v37.transcripts.fa: Transcriptoma de referencia para realizar el alineamiento y conteo con Salmon
ObjetosR.RData: Archivo de R que contiene los objetos de SummarizedExperiment y DESeq para realizar el análisis de expresión diferencial. Descarga este objeto si no tienes suficiente espacio en tu computadora

Instalación

Instalación de FastQC

1- Descargar el programa de su página web fastQC

2- Hacer el archivo ejecutable

chmod 755 /Applications/FastQC.app/Contents/MacOS/fastqc
ls -l /Applications/FastQC.app/Contents/MacOS/fastqc

Tiempo: 0

3- Crear una liga simbolica para ejecutar fastqc desde cualquier ubicación

ln -s /Applications/FastQC.app/Contents/MacOS/fastqc /usr/local/bin/fastqc

Tiempo: 0.14

4- Correr fastqc

fastqc --help

Tiempo: 0

Instalación de MultiQC

Esta parte es opcional, pero si trabajas con RNA-seq te puede intersar. MultiQC es un programa que te perimte juntar los archivos de fastQCen uno solo y facilita su visualización.

Conda: La instalación se puede hacer con conda, así que primero activamos un ambiente conda

conda activate base
conda install -c bioconda -c conda-forge multiqc

Tiempo: 0

Pip: o con pip

pip install multiqc

Tiempo: 0

Instalación de Salmon

Descargar el archivo salmon-1.4.0_linux_x86_64.tar.gz desde la página de Salmon y seguir las instrucciones.

Si lo prefieres, la instalación la puedes realizar empleando conda:

#Descarga e instala el gestor de paquetes de python miniconda desde https://conda.io/projects/conda/en/latest/user-guide/install/index.html
conda activate miniconda
conda install -c bioconda salmon

Tiempo: 0

Instalación de R y Rstudio

Puedes instalar R a través de Rstudio, asegurate de tener una versión de R > 4.0. Puedes checar la versión de R que tienes si ejecutas version en tu consola de R:

version

##                _                           
## platform       x86_64-apple-darwin17.0     
## arch           x86_64                      
## os             darwin17.0                  
## system         x86_64, darwin17.0          
## status                                     
## major          4                           
## minor          0.2                         
## year           2020                        
## month          06                          
## day            22                          
## svn rev        78730                       
## language       R                           
## version.string R version 4.0.2 (2020-06-22)
## nickname       Taking Off Again

Tiempo: 0.01

Instalación de Bioconductor

Durante esta práctica se usarán paqueterías que se enuentran en el repositorio de Bioconductor, usaremos BiocManager para instalar los diferentes paquetes pero asegurate de tener la versión más actual de

# Instalar BiocManager
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Checar la versión
BiocManager::version()

# Si se requiere, se puede instalar la versión 3.12 manualmente
# BiocManager::install(version = "3.12")

Tiempo: 0

Instalación de la paquetería de R

Instalar la paquetería utilizando BiocManager

paquetes = c("DESeq2", "tximeta", "PCAtools", "SummarizedExperiment")
BiocManager::install()

Tiempo: 0

Adicionalmente, ocuparemos algunas funciones de tidyverse como la paquetería de ggplot2 y pheatmap. La instalaremos utilizando el repositorio de CRAN

install.packages(c('tidyverse','pheatmap'))

Tiempo: 0

# 2: Organización del proyecto

Slides

Tiempo: 0.01

## Experimento

Número de acceso GEO: GSE152699

Tiempo: 0.01

# 3: Flujo de trabajo para el análisis de RNA-seq

Slides

En este mini curso veremos brevemente como analizar datos de RNA-seq, tendremos una breve introdución a RNA-seq, así como los pasos a seguir en el análisis bioinformático (como se observa en el diagrama).

r knitr::include_graphics("https://biocorecrg.github.io/RNAseq_course_2019/images/RNAseq_workflow.png")

Tiempo: 0.01

Imagen tomada de BiCore RNAseq course 2019

5: Control de calidad de los datos

Slides

Muchas veces a esté paso se le conoce como QC (Quality Check).

Pasos para correr FastQC:

1- Checar los datos fastq

pwd
ls data/Archivos_fastq

Tiempo: 0

2- Crear un directorio para guardar el output

mkdir Output/FastQC

Tiempo: 0

3- Correr FastQC

# la opción -o es para elejir el directorio de salida
fastqc data/Archivos_fastq/*.fastq.gz -o Output/FastQC/

Tiempo: 0

4- Checar los archivos generados

ls Output/FastQC

Tiempo: 0

MultiQC

MultiQC te permite juntar los output en un solo archivo y poder visualizarlos mejor.

*Esta parte es opcional, pero si trabajas con RNA-seq te puede intersar.

Instalar MultiQC
Correr MultiQC

cd Output/.
multiQC FastQC/ -o MultiQC/

Tiempo: 0

Checar el output.html

6: Alineamiento (Salmon)

Slides

Si deseas realizar el pseudo-alineamiento con Salmon te recomendamos descargar los archivos .fastq de las líneas celulares así como el archivo .fasta del transcriptoma humano.

Generar índice del transcriptoma

Para ello vamos a requerir:

Archivo Fasta del transcriptoma de referencia del humano. Descargado del sitio de GeneCode. Lo pueden encontrar en el folder de Reference/

Código para genear el índice

-t: Ubicación al archivo Fasta del transcriptoma

-i: Ubicación para salvar el índice

--genecode: El Fasta del transcriptoma de referencia está en formato de GENECODE

salmon index -t ../transcriptome/gencode.v37.transcripts.fa.gz -i ../transcriptome/genecode.v37.salmon141 -p 6

Tiempo: 0

Cuantificar la abundancia de los transcritos

Requerimientos:

Archivos Fastq de las lecturas -> Paired end R1 y R2

Ubicados en el folder de data

Índice generado en el paso anterior

Código para producir cuantificar los transcritos

-i: Ubicación del índice

-l: Tipo de librería

-1 y -2: Ubicación a las lecturas R1 y R2

-o: Ubicación para almacenar los resultados

##Crear un directorio para almacenar los datos de los conteos
mkdir -p ../salmon_quant
##Llamar salmon para realizar el conteo
salmon quant -i ../transcriptome/genecode.v37.salmon141 \
             -l A \
             -1 ../Data/s1_R1.fastq.gz -2 ../Data/s1_R2.fastq.gz \
             -p 6 --validateMappings \
             -o ../salmon_quant/s1_quant

Tiempo: 0

7: Tximeta

Slides

knitr::include_graphics("https://journals.plos.org/ploscompbiol/article/figure/image?size=large&id=10.1371/journal.pcbi.1007664.g001")

Tiempo: 0.01

Love, Michael I., et al. “Tximeta: Reference sequence checksums for provenance identification in RNA-seq.” PLoS computational biology 16.2 (2020): e1007664.

Instalación de tximeta

Instalar tximeta desde R usando BiocManager

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Puedes instalar la paqueteria de tximeta usando BiocManager::install() 
BiocManager::install("tximeta")

Tiempo: 0

Importar datos

Utilizamos tximeta para importar los datos.

Nota: tximeta espera dos columnas:

files: con la rut a quant.sf
names: con los nombres de las muestras

Importamos la información de cada muestra

# Leemos los datos
coldata <- read.delim(here::here("output/salmon_quants/metadata.txt"))

# Generamos la columna "names"
# info$names <- info$Unique_id # Puedes agregar una columna nueva
coldata <- dplyr::rename(coldata, names = unique_id) # Renombrar la columna 
rownames(coldata) <- coldata$names
coldata

##               Sample Cell Treatment Replicate         names
## A_Control_1       s1    A   Control         1   A_Control_1
## A_Control_2       s2    A   Control         2   A_Control_2
## A_Control_3       s3    A   Control         3   A_Control_3
## M_Control_1       s4    M   Control         1   M_Control_1
## M_Control_2       s5    M   Control         2   M_Control_2
## M_Control_3       s6    M   Control         3   M_Control_3
## A_Verafinib_1     s7    A Verafinib         1 A_Verafinib_1
## A_Verafinib_2     s8    A Verafinib         2 A_Verafinib_2
## A_Verafinib_3     s9    A Verafinib         3 A_Verafinib_3
## M_Verafinib_1    s10    M Verafinib         1 M_Verafinib_1
## M_Verafinib_2    s11    M Verafinib         2 M_Verafinib_2
## M_Verafinib_3    s12    M Verafinib         3 M_Verafinib_3

Tiempo: 0.63

Ahora las cuentas generadas por Salmon

# Ubicamos el dir de trabajo
dir <- file.path(here::here("output/salmon_quants"))
#Checamos que esten los folders  
list.files(dir)

##  [1] "Icon\r"       "metadata.txt" "s1_quant"     "s10_quant"    "s11_quant"   
##  [6] "s12_quant"    "s2_quant"     "s3_quant"     "s4_quant"     "s5_quant"    
## [11] "s6_quant"     "s7_quant"     "s8_quant"     "s9_quant"

# Que hacemos aqui??  Checar
# file.path(dir,paste0(info$Sample,"_quant"),"quant.sf")
coldata$files <- file.path(dir,paste0(coldata$Sample,"_quant"),"quant.sf")

#Checamos que los archivos existan
data.frame(coldata$Sample, file.exists(coldata$files))

##    coldata.Sample file.exists.coldata.files.
## 1              s1                       TRUE
## 2              s2                       TRUE
## 3              s3                       TRUE
## 4              s4                       TRUE
## 5              s5                       TRUE
## 6              s6                       TRUE
## 7              s7                       TRUE
## 8              s8                       TRUE
## 9              s9                       TRUE
## 10            s10                       TRUE
## 11            s11                       TRUE
## 12            s12                       TRUE

Tiempo: 0.02

Importar los datos usando tximeta

# Importamos la paqueteria 
library(tximeta)

## Warning: package 'tximeta' was built under R version 4.0.5

# Checamos que el df tenga la columna "names" y "files"
coldata

##               Sample Cell Treatment Replicate         names
## A_Control_1       s1    A   Control         1   A_Control_1
## A_Control_2       s2    A   Control         2   A_Control_2
## A_Control_3       s3    A   Control         3   A_Control_3
## M_Control_1       s4    M   Control         1   M_Control_1
## M_Control_2       s5    M   Control         2   M_Control_2
## M_Control_3       s6    M   Control         3   M_Control_3
## A_Verafinib_1     s7    A Verafinib         1 A_Verafinib_1
## A_Verafinib_2     s8    A Verafinib         2 A_Verafinib_2
## A_Verafinib_3     s9    A Verafinib         3 A_Verafinib_3
## M_Verafinib_1    s10    M Verafinib         1 M_Verafinib_1
## M_Verafinib_2    s11    M Verafinib         2 M_Verafinib_2
## M_Verafinib_3    s12    M Verafinib         3 M_Verafinib_3
##                                                                                                                   files
## A_Control_1    /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s1_quant/quant.sf
## A_Control_2    /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s2_quant/quant.sf
## A_Control_3    /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s3_quant/quant.sf
## M_Control_1    /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s4_quant/quant.sf
## M_Control_2    /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s5_quant/quant.sf
## M_Control_3    /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s6_quant/quant.sf
## A_Verafinib_1  /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s7_quant/quant.sf
## A_Verafinib_2  /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s8_quant/quant.sf
## A_Verafinib_3  /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s9_quant/quant.sf
## M_Verafinib_1 /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s10_quant/quant.sf
## M_Verafinib_2 /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s11_quant/quant.sf
## M_Verafinib_3 /Users/user/Documents/Doctorado/Courses/Taught/minicurso_abr_2021/output/salmon_quants/s12_quant/quant.sf

# Generamos el objeto 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 
## found matching transcriptome:
## [ Ensembl - Homo sapiens - release 103 ]
## loading existing EnsDb created: 2021-04-14 19:47:23
## loading existing transcript ranges created: 2021-04-14 19:49:06

## Warning in checkAssays2Txps(assays, txps): 
## 
## Warning: the annotation is missing some transcripts that were quantified.
## 12741 out of 184844 txps were missing from GTF/GFF but were in the indexed FASTA.
## (This occurs sometimes with Ensembl txps on haplotype chromosomes.)
## In order to build a ranged SummarizedExperiment, these txps were removed.
## To keep these txps, and to skip adding ranges, use skipMeta=TRUE
## 
## Example missing txps: [ENST00000631435, ENST00000632524, ENST00000633009, ...]

Tiempo: 20.94

# Opción si no funciona tximeta
#load("../data/ObjetosR.RData")

Tiempo: 0

8: Objeto SummarizedExperiment en R

Slides

Tiempo: 0

library(SummarizedExperiment)

Tiempo: 0.25

se

## class: RangedSummarizedExperiment 
## dim: 172103 12 
## metadata(6): tximetaInfo quantInfo ... txomeInfo txdbInfo
## assays(3): counts abundance length
## rownames(172103): ENST00000415118 ENST00000434970 ... ENST00000606902
##   ENST00000427827
## rowData names(6): tx_id tx_biotype ... gene_id tx_name
## colnames(12): A_Control_1 A_Control_2 ... M_Verafinib_2 M_Verafinib_3
## colData names(5): Sample Cell Treatment Replicate names

Tiempo: 0.06

ColData

El cajón o slot correspondiente a ColData contiene la tabla de metadatos creada para importar las cuentas con tximeta.

Para acceder a ColData usar el siguiente comando:

colData(se)

## DataFrame with 12 rows and 5 columns
##                    Sample        Cell   Treatment Replicate         names
##               <character> <character> <character> <integer>   <character>
## A_Control_1            s1           A     Control         1   A_Control_1
## A_Control_2            s2           A     Control         2   A_Control_2
## A_Control_3            s3           A     Control         3   A_Control_3
## M_Control_1            s4           M     Control         1   M_Control_1
## M_Control_2            s5           M     Control         2   M_Control_2
## ...                   ...         ...         ...       ...           ...
## A_Verafinib_2          s8           A   Verafinib         2 A_Verafinib_2
## A_Verafinib_3          s9           A   Verafinib         3 A_Verafinib_3
## M_Verafinib_1         s10           M   Verafinib         1 M_Verafinib_1
## M_Verafinib_2         s11           M   Verafinib         2 M_Verafinib_2
## M_Verafinib_3         s12           M   Verafinib         3 M_Verafinib_3

Tiempo: 0.03

El slot ColData es un objeto con clase de DataFrame

class(colData(se))

## [1] "DFrame"
## attr(,"package")
## [1] "S4Vectors"

Tiempo: 0.01

Rownames de ColData corresponden a los Colnames en el slot Assay- Importante para análisis con DESeq2

rownames(colData(se))

##  [1] "A_Control_1"   "A_Control_2"   "A_Control_3"   "M_Control_1"  
##  [5] "M_Control_2"   "M_Control_3"   "A_Verafinib_1" "A_Verafinib_2"
##  [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3"

Tiempo: 0.01

rowRanges

El cajón de rowRanges hace referencia a las coordenadas de cada transcrito

Para acceder al rowRanges usar:

rowRanges(se)

## GRanges object with 172103 ranges and 6 metadata columns:
##                   seqnames              ranges strand |           tx_id
##                      <Rle>           <IRanges>  <Rle> |     <character>
##   ENST00000415118       14   22438547-22438554      + | ENST00000415118
##   ENST00000434970       14   22439007-22439015      + | ENST00000434970
##   ENST00000448914       14   22449113-22449125      + | ENST00000448914
##   ENST00000604642       15   20003840-20003862      - | ENST00000604642
##   ENST00000603326       15   20004797-20004815      - | ENST00000603326
##               ...      ...                 ...    ... .             ...
##   ENST00000444456        1 247974741-247975653      + | ENST00000444456
##   ENST00000416377        1 248003129-248004068      + | ENST00000416377
##   ENST00000318104        1 248121932-248122882      + | ENST00000318104
##   ENST00000606902        1 248498308-248498649      + | ENST00000606902
##   ENST00000427827        1 248549444-248549762      - | ENST00000427827
##                               tx_biotype tx_cds_seq_start tx_cds_seq_end
##                              <character>        <integer>      <integer>
##   ENST00000415118              TR_D_gene         22438547       22438554
##   ENST00000434970              TR_D_gene         22439007       22439015
##   ENST00000448914              TR_D_gene         22449113       22449125
##   ENST00000604642              IG_D_gene         20003840       20003862
##   ENST00000603326              IG_D_gene         20004797       20004815
##               ...                    ...              ...            ...
##   ENST00000444456 unprocessed_pseudogene             <NA>           <NA>
##   ENST00000416377 unprocessed_pseudogene             <NA>           <NA>
##   ENST00000318104 unprocessed_pseudogene             <NA>           <NA>
##   ENST00000606902 unprocessed_pseudogene             <NA>           <NA>
##   ENST00000427827 unprocessed_pseudogene             <NA>           <NA>
##                           gene_id         tx_name
##                       <character>     <character>
##   ENST00000415118 ENSG00000223997 ENST00000415118
##   ENST00000434970 ENSG00000237235 ENST00000434970
##   ENST00000448914 ENSG00000228985 ENST00000448914
##   ENST00000604642 ENSG00000270961 ENST00000604642
##   ENST00000603326 ENSG00000271317 ENST00000603326
##               ...             ...             ...
##   ENST00000444456 ENSG00000237492 ENST00000444456
##   ENST00000416377 ENSG00000232215 ENST00000416377
##   ENST00000318104 ENSG00000177233 ENST00000318104
##   ENST00000606902 ENSG00000271934 ENST00000606902
##   ENST00000427827 ENSG00000227102 ENST00000427827
##   -------
##   seqinfo: 47 sequences from GRCh38 genome

Tiempo: 0.04

El cajón rowRanges almacena metadatos de cada secuencia (en este caso cromosomas)

seqinfo(rowRanges(se))

## Seqinfo object with 47 sequences from GRCh38 genome:
##   seqnames   seqlengths isCircular genome
##   1           248956422       <NA> GRCh38
##   10          133797422       <NA> GRCh38
##   11          135086622       <NA> GRCh38
##   12          133275309       <NA> GRCh38
##   13          114364328       <NA> GRCh38
##   ...               ...        ...    ...
##   KI270744.1     168472       <NA> GRCh38
##   KI270750.1     148850       <NA> GRCh38
##   MT              16569       <NA> GRCh38
##   X           156040895       <NA> GRCh38
##   Y            57227415       <NA> GRCh38

Tiempo: 0.02

Assay

El slot assay almacena la información de las cuentas para cada transcrito dividida en tres niveles:

assayNames(se)

## [1] "counts"    "abundance" "length"

Tiempo: 0

Para acceder a la matriz de cuentas estimadas por Salmon, correr:

head(assay(se), 5)

##                 A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2
## ENST00000415118           0           0           0           0           0
## ENST00000434970           0           0           0           0           0
## ENST00000448914           0           0           0           0           0
## ENST00000604642           0           0           0           0           0
## ENST00000603326           0           0           0           0           0
##                 M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3
## ENST00000415118           0             0             0             0
## ENST00000434970           0             0             0             0
## ENST00000448914           0             0             0             0
## ENST00000604642           0             0             0             0
## ENST00000603326           0             0             0             0
##                 M_Verafinib_1 M_Verafinib_2 M_Verafinib_3
## ENST00000415118             0             0             0
## ENST00000434970             0             0             0
## ENST00000448914             0             0             0
## ENST00000604642             0             0             0
## ENST00000603326             0             0             0

Tiempo: 0.02

Las matriz de abundancia (TPM) puede obtenerse:

## Obtener matriz de TPM
head(se@assays@data$abundance, 5)

##                 A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2
## ENST00000415118           0           0           0           0           0
## ENST00000434970           0           0           0           0           0
## ENST00000448914           0           0           0           0           0
## ENST00000604642           0           0           0           0           0
## ENST00000603326           0           0           0           0           0
##                 M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3
## ENST00000415118           0             0             0             0
## ENST00000434970           0             0             0             0
## ENST00000448914           0             0             0             0
## ENST00000604642           0             0             0             0
## ENST00000603326           0             0             0             0
##                 M_Verafinib_1 M_Verafinib_2 M_Verafinib_3
## ENST00000415118             0             0             0
## ENST00000434970             0             0             0
## ENST00000448914             0             0             0
## ENST00000604642             0             0             0
## ENST00000603326             0             0             0

Tiempo: 0.01

¿Recuerdan a qué tiene que ser igual Rownames del slot colData?

rownames(colData(se))

##  [1] "A_Control_1"   "A_Control_2"   "A_Control_3"   "M_Control_1"  
##  [5] "M_Control_2"   "M_Control_3"   "A_Verafinib_1" "A_Verafinib_2"
##  [9] "A_Verafinib_3" "M_Verafinib_1" "M_Verafinib_2" "M_Verafinib_3"

Tiempo: 0.01

## Comprobar que rownames de colData es igual a colnames de assay
row.names(colData(se)) == colnames(assay(se))

##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Tiempo: 0.02

9: Normalización de los datos (RPKM, TPM)

Slides

Es el primer paso del análisis de expresión diferencial y es necesario para hacer comparaciones acertadas entre muestras.
Las cuentas crudas están conformadas por un componente “interesante” (la expresión de RNA) y componentes “no interesantes” (como los batch effects, ruido de la plataforma, etc.).
La normalzación escala las cuentas para tratar de reducir los componentes “no interesantes” y poder comparar las muestras entre si.

Métodos comunes

Método de normalización	Descripción	Factores de evaluación	Recomendaciones de uso
CPM (cuentas por millón): cuentas escalados por el número total de lecturas	Profundidad de secuenciación	Comparaciones de cuentas de genes entre réplicas del mismo grupo de muestras	NO para comparaciones dentro de la muestra o análisis de DE.
TPM (transcritos por millón de lecturas): cuentas escalados por el número total de lecturas	Profundidad de secuenciación	Comparaciones de cuentas de genes entre réplicas del mismo grupo de muestras	NO para análisis de DE.
RPKM/FPKM (lecturas/fragmentos por kilobase de exón por millón de lecturas/fragmentos mapeados	Similar a TPM, profundidad de secuenciación y longitud del gen	Comparaciones de cuentas entre genes dentro de una muestra	NO para comparaciones entre muestras o análisis de DE.
Factor de normalización de DESeq2	Cuentas divididas por factores de tamaño específicos de la muestra determinados por la mediana del ratio de cuentas de genes en relación con la media geométrica por gen	Profundidad de secuenciación y composición del RNA	Comparaciones de cuentas de genes entre muestras y para el análisis de DE; NO para comparaciones dentro de la muestra
La media cortada de los valores M de EdgeR (TMM)	Utiliza una media recortada ponderada de los ratios de expresión logarítmica entre las muestras	Profundidad de secuenciación	Composición de RNA y longitud de los genes.

Tabla tomada de Intro to DGE

DESeq2

Normalización para análisis de expresión diferencial:
- Factor de normalización de DESeq2
Normalización para visualización u otras aplicaciones:
- Variance stabilizing transformation (VST)
- Regularized-logarithm transformation (rlog)

DESeq2ajusta a un modelo lineal generalizado (GLM) de la familia binomial negativa (NB).

Factor de normalización en DESeq2

Ejemplo de como se calcula el factor de normalización en DESeq2.

Primero generamos una tabla con cuentas simuladas

# tabla con cuentas 
df <- tibble(gene = c("gene1", "gene2"),
                 muestraA = c(1749, 35),
                 muestraB = c(943, 29)
                 )
df

## # A tibble: 2 x 3
##   gene  muestraA muestraB
##   <chr>    <dbl>    <dbl>
## 1 gene1     1749      943
## 2 gene2       35       29

Tiempo: 0.04

Crea una pseudo-referencia por muestra (promedio geometrico por fila) sqrt(muestraA * muestra B)

# Calcular el promedio geometrico
df <- df %>%
  rowwise() %>% 
  mutate(prom_geom = sqrt(muestraA * muestraB))
df

## # A tibble: 2 x 4
## # Rowwise: 
##   gene  muestraA muestraB prom_geom
##   <chr>    <dbl>    <dbl>     <dbl>
## 1 gene1     1749      943    1284. 
## 2 gene2       35       29      31.9

Tiempo: 0.05

Se calcula la fración muestra/pseudo-referencia

# Dividir las cuentas entre el promedio geometrico
df <- df %>% 
  rowwise() %>% 
  mutate(muestraA_pseudo_ref = muestraA / prom_geom) %>% 
  mutate(muestraB_pseudo_ref = muestraB / prom_geom)

df

## # A tibble: 2 x 6
## # Rowwise: 
##   gene  muestraA muestraB prom_geom muestraA_pseudo_ref muestraB_pseudo_ref
##   <chr>    <dbl>    <dbl>     <dbl>               <dbl>               <dbl>
## 1 gene1     1749      943    1284.                 1.36               0.734
## 2 gene2       35       29      31.9                1.10               0.910

Tiempo: 0.07

Se calcula un factor de normalización (size factor) utilizando la median por columnas.

# Se calcula el factor de normalizacion usando la mediana para cada muestra
norm_factor_muestraA <- median(df$muestraA_pseudo_ref)
norm_factor_muestraA

## [1] 1.230235

# Repetimos el proceso para la muestra B
norm_factor_muestraB <- median(df$muestraB_pseudo_ref)
norm_factor_muestraB

## [1] 0.8222689

Tiempo: 0.04

Se dividen las cuentas crudas/size factor para calcular las cuentas normalizadas.

# Columnas con las cuentas crudas
df %>% 
  select(gene, muestraA, muestraB)

## # A tibble: 2 x 3
## # Rowwise: 
##   gene  muestraA muestraB
##   <chr>    <dbl>    <dbl>
## 1 gene1     1749      943
## 2 gene2       35       29

# Dividir las cuentras entre el factor de normalización
df$norm_muestraA <- df$muestraA/norm_factor_muestraA
df$norm_muestraB <- df$muestraB/norm_factor_muestraB

# Columnas con las cuentas normalizadas
df %>%  
  select(norm_muestraA, norm_muestraB)

## # A tibble: 2 x 2
## # Rowwise: 
##   norm_muestraA norm_muestraB
##           <dbl>         <dbl>
## 1        1422.         1147. 
## 2          28.4          35.3

Tiempo: 0.2

Ejercicio

Se pueden ocupar otras transformaciones en la paquetería de DESeq2pero no son las mas recomendadas para DE.

Para normalizar, primero construimos el objeto DESeqDataSet, para ello necesitamos definir el modelo de DE (Quienes son Controles y quienes tratamiento).

# Importamos las librerías que necesitemos
library("DESeq2")
library("dplyr")
library("ggplot2")
library("vsn")

Tiempo: 1.25

# Definir los grupos que se desean comparar
se$Treatment <- factor(se$Treatment)

# Checar cantidad de muestras
table(se$Treatment)

## 
##   Control Verafinib 
##         6         6

# Generar el objeto de DESeq 
dds <- DESeqDataSet(se, design = ~ Treatment)

## using counts and average transcript lengths from tximeta

# Prefiltrado para quitar genes con cuentas bajas
keep <- rowSums(counts(dds) >= 3) >= 6
dds <- dds[keep, ]

# Funcion que normaliza los datos y realiza el analisis de expresión differencial
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

Tiempo: 22.44

Otras transformaciones

Puedes realizar otras transformaciones en DESeq2 para estabilizar la varianza a través de los differentes valores promedio de expresión.

# Grafica de cuentas crudas
meanSdPlot(counts(dds), ranks = F)

# Gráfica de cuentas en escala log2
meanSdPlot(log2(counts(dds) +1), ranks = F)

Tiempo: 1.77

# variance stabilizing transformation (VST), (Anders and Huber 2010)
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)

##                 A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2
## ENST00000641506    5.243738    5.022432    5.236241    5.121340    5.055195
## ENST00000611116    5.058149    5.063223    4.758866    4.737813    4.871766
## ENST00000390289    4.900795    5.154803    5.002771    4.909267    4.423150
##                 M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3
## ENST00000641506    5.176386      5.508017      5.293662      5.474252
## ENST00000611116    5.289621      4.874653      4.423150      4.423150
## ENST00000390289    5.040044      5.058810      4.423150      5.039269
##                 M_Verafinib_1 M_Verafinib_2 M_Verafinib_3
## ENST00000641506      5.716058      5.895262      5.482647
## ENST00000611116      5.185938      5.045729      5.180529
## ENST00000390289      5.571450      5.182835      5.294481

# regularized-logarithm transformation (rlog), (Love, Huber, and Anders 2014)
rld <- rlog(dds, blind = FALSE)
head(assay(rld), 3)

##                 A_Control_1 A_Control_2 A_Control_3 M_Control_1 M_Control_2
## ENST00000641506    2.905283    2.239445   2.8848975   2.5548691   2.4040501
## ENST00000611116    1.829812    1.846142   0.6284644   0.5186003   1.1642110
## ENST00000390289    1.416904    2.289820   1.8015912   1.4496085   0.1791245
##                 M_Control_3 A_Verafinib_1 A_Verafinib_2 A_Verafinib_3
## ENST00000641506    2.739348      3.526028    3.03947595    3.45477910
## ENST00000611116    2.444677      1.174766   -0.05979834   -0.07010645
## ENST00000390289    1.918379      1.975739    0.15929486    1.91598057
##                 M_Verafinib_1 M_Verafinib_2 M_Verafinib_3
## ENST00000641506      3.938412      4.256634      3.472666
## ENST00000611116      2.166800      1.752071      2.152501
## ENST00000390289      3.267274      2.339360      2.635467

# Normalización con el factor de normalizacion 
dds <- estimateSizeFactors(dds)

## using 'avgTxLength' from assays(dds), correcting for library size

# Juntar los datos de las tres normalizaciones
df <- bind_rows(
  as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>%
         mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"),
  as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"))

## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.

# Renombrar columnas
colnames(df)[1:2] <- c("x", "y")  

# Nombre de las graficas
lvls <- c("log2(x + 1)", "vst", "rlog")

# Agrupar los tres tipos de normalizacion en grupos como factores
df$transformation <- factor(df$transformation, levels=lvls)

# Plotear los datos
ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
  coord_fixed() + facet_grid( . ~ transformation)

Tiempo: 12.62

10: Exploración de los datos

Slides

PCA

Análisis de componentes principales
Método algebraico para reducir la dimensionalidad de sets de datos complejos (múltiples variables)
Reducción de dimensionalidad o variables permite un análisis exploratorio más intuitivo
Reducción de variables implica preservar y captar la mayor información sobre los datos

¿Cómo crear nuestro propio PCA?

Primero instalen las siguientes librerías:

library(DESeq2)
library(PCAtools)

## Warning: package 'PCAtools' was built under R version 4.0.3

## Loading required package: ggrepel

## 
## Attaching package: 'PCAtools'

## The following objects are masked from 'package:stats':
## 
##     biplot, screeplot

Tiempo: 0.9

Usemos la siguiente normalización logarítmica de las cuentas (toma en cuenta el tamaño de

## Transformación recomendada para sets de con menos de 30 muestras
rld <- rlog(dds, blind = F)

Tiempo: 9.7

Usemos la función interna de DESeq2 para graficar nuestro PCA:

## El argumento intgroup permite especificar mediante cuál variable colorear los datos
plotPCA(rld, intgroup = "Treatment")

Tiempo: 0.48

Con la librería de PCAtools:

## Crear un objeto que contenga los datos del PCA
PCA <- pca(assay(rld), scale = T, metadata = colData(rld)) 
## Graficar en 2D los resultados
biplot(PCA, colby = "Treatment")

## Generar un screeplot para visualizar la varianza asociada a cada componente
screeplot(PCA)

Tiempo: 1.65

¿Cómo ven los resultados?

plotPCA(rld, intgroup = "Treatment")+
  geom_label_repel(aes(label = colnames(assay(rld))), 
                   segment.color = "grey50", 
                   box.padding = 0.35, 
                   point.padding = 0.5)

Tiempo: 1.46

biplot(PCA, lab = colnames(assay(rld)), colby = "Treatment", legendPosition = "right")

Tiempo: 1.03

11: Análisis de expresión diferencial

Slides

Paso para obtener cuales son los genes que varían entre las condiciones

knitr::include_graphics("https://raw.githubusercontent.com/hbctraining/DGE_workshop/master/img/de_theory.png")

Tiempo: 0

Imagen tomada de Intro to DGE

Recordemos el objeto que creeamos

# Checar objeto generado 
dds

## class: DESeqDataSet 
## dim: 35483 12 
## metadata(7): tximetaInfo quantInfo ... txdbInfo version
## assays(7): counts abundance ... H cooks
## rownames(35483): ENST00000641506 ENST00000611116 ... ENST00000500626
##   ENST00000314835
## rowData names(28): tx_id tx_biotype ... deviance maxCooks
## colnames(12): A_Control_1 A_Control_2 ... M_Verafinib_2 M_Verafinib_3
## colData names(5): Sample Cell Treatment Replicate names

Tiempo: 0.02

Para ver los resultados podemos usar la función results

# Obtener los resultados del DE
res <- results(dds)
res

## log2 fold change (MLE): Treatment Verafinib vs Control 
## Wald test p-value: Treatment Verafinib vs Control 
## DataFrame with 35483 rows and 6 columns
##                   baseMean log2FoldChange     lfcSE       stat      pvalue
##                  <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## ENST00000641506   10.05792      1.3610503  0.360096  3.7796900 1.57024e-04
## ENST00000611116    3.26143     -0.0421221  0.808389 -0.0521062 9.58444e-01
## ENST00000390289    4.55373      0.9128357  0.601744  1.5169831 1.29271e-01
## ENST00000361390 3379.18027      2.1135181  0.379169  5.5740792 2.48843e-08
## ENST00000361453 6522.68449      1.5906491  0.385459  4.1266327 3.68114e-05
## ...                    ...            ...       ...        ...         ...
## ENST00000622095    7.43638       0.769911  0.620186   1.241418 2.14451e-01
## ENST00000605551   40.04273      -0.846855  0.397376  -2.131115 3.30797e-02
## ENST00000579209    9.61591       1.597767  0.386827   4.130447 3.62059e-05
## ENST00000500626    3.46026       0.798818  0.733797   1.088610 2.76326e-01
## ENST00000314835   17.48904       0.166453  0.674825   0.246661 8.05171e-01
##                        padj
##                   <numeric>
## ENST00000641506 1.62737e-03
## ENST00000611116 9.82689e-01
## ENST00000390289 3.24642e-01
## ENST00000361390 7.31775e-07
## ENST00000361453 4.68935e-04
## ...                     ...
## ENST00000622095 0.451633725
## ENST00000605551 0.124474321
## ENST00000579209 0.000463051
## ENST00000500626 0.525737614
## ENST00000314835 0.915847174

summary(res)

## 
## out of 35483 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 4420, 12%
## LFC < 0 (down)     : 4315, 12%
## outliers [1]       : 18, 0.051%
## low counts [2]     : 0, 0%
## (mean count < 1)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

Tiempo: 1.65

Para ver que hay en cada columna:

mcols(res, use.names = TRUE)

## DataFrame with 6 rows and 2 columns
##                        type            description
##                 <character>            <character>
## baseMean       intermediate mean of normalized c..
## log2FoldChange      results log2 fold change (ML..
## lfcSE               results standard error: Trea..
## stat                results Wald statistic: Trea..
## pvalue              results Wald test p-value: T..
## padj                results   BH adjusted p-values

Tiempo: 0.01

Podemos hacer cortes:

resLFC1 <- results(dds, lfcThreshold=1)
table(resLFC1$padj < 0.1)

## 
## FALSE  TRUE 
## 34309  1156

Tiempo: 1.55

12: Visualización de los resultados

Slides

MAplot

Gráfico que representa la distribución de los genes o transcritos en las comparaciones de interés

res_lfc <- lfcShrink(dds, coef = "Treatment_Verafinib_vs_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

group <- coldata$Treatment
dds$group <- group

Tiempo: 20.36

##Es importante que recuerden que la hipótesis nula que se probó fue
##"El lfc del gen n es igual a 0" por lo tanto los genes coloreados son...
plotMA(res)

Tiempo: 2.77

También pueden generar un MAplot interactivo con la librería de Glimma

library(Glimma)
glimmaMA(dds)

Tiempo: 0

Volcano plot

De manera similar al MAplot con el volcano plot visualizamos los genes que muestran expresión diferencial en nuestra condición de interés

En el eje y se grafica el -log10 de padj
En el eje x se grafica el lfc o log2foldchange

##Para crear un volcano plot necesitas convertir los resultados de DESeq a un data frame
DEG <- as.data.frame(res)
##En el script de funciones encontrarás la función de volcanoplotR para generar tu gráfico
#los valores de los argumentos logfc y p.adj deben ser iguales al threshold utilizado para generar los resultados
#en type debes indicar que los resultados provienen de DESeq
source("./functions.R")
volcanoplotR(DEG, logfc = 0, p.adj = 0.1, type = "DESeq")

Tiempo: 4.56

Heatmap

El heatmap nos permite visualizar la expresión de los genes diferencialmente expresados en terminos de las cuentas normalizadas

Consideraciones:

Usar los valores de las cuentas normalizadas para una mejor comparación entre muestras
Escalar los valores de las cuentas (renglones) para visualizar las diferencias en la expresión

Usaremos la librería de pheatmap

library(pheatmap)

##Guardar la lista de transcritos que mostraron expresión diferencial significativa
significant <- DEG %>% filter(log2FoldChange > 0 & padj < 0.1 |
                                   log2FoldChange < 0 & padj < 0.1)
##Generar la matriz de cuentas normalizadas
norm_counts <- counts(dds, normalized = T)
##Cortar la matriz de cuentas normalizadas con los id de los transcritos diferencialmente expresados
norm_counts <- norm_counts[rownames(significant), ]
##Generar una tabla de anotaciones preservando el tratamiento y el tipo de células
annotation_col <- coldata[, c("names","Cell", "Treatment")]
##Generar el heatmap empleando clustering jerarquico
pheatmap(norm_counts, 
         border_color = NA, 
         scale = "row",
         clustering_distance_rows = "euclidean", 
         clustering_distance_cols = "euclidean", 
         clustering_method = "average", 
         show_colnames = F, 
         show_rownames = F, 
         annotation_col = annotation_col)

Tiempo: 12.11

library(pheatmap)
pheatmap(norm_counts, 
         border_color = NA, 
         scale = "row",
         clustering_distance_rows = "euclidean", 
         clustering_distance_cols = "euclidean", 
         clustering_method = "average", 
         show_colnames = F, 
         show_rownames = F, 
         annotation_col = annotation_col)

Tiempo: 9.81

Paquetería

sessionInfo()

## R version 4.0.2 (2020-06-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS  10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pheatmap_1.0.12             PCAtools_2.2.0             
##  [3] ggrepel_0.9.1               vsn_3.58.0                 
##  [5] DESeq2_1.30.1               forcats_0.5.1              
##  [7] stringr_1.4.0               dplyr_1.0.5                
##  [9] purrr_0.3.4                 readr_1.4.0                
## [11] tidyr_1.1.3                 tibble_3.1.0               
## [13] ggplot2_3.3.3               tidyverse_1.3.0            
## [15] SummarizedExperiment_1.20.0 Biobase_2.50.0             
## [17] GenomicRanges_1.42.0        GenomeInfoDb_1.26.7        
## [19] IRanges_2.24.1              S4Vectors_0.28.1           
## [21] BiocGenerics_0.36.0         MatrixGenerics_1.2.1       
## [23] matrixStats_0.58.0          tximeta_1.8.5              
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1                  backports_1.2.1              
##   [3] AnnotationHub_2.22.0          BiocFileCache_1.14.0         
##   [5] plyr_1.8.6                    lazyeval_0.2.2               
##   [7] splines_4.0.2                 BiocParallel_1.24.1          
##   [9] digest_0.6.27                 ensembldb_2.14.0             
##  [11] htmltools_0.5.1.1             fansi_0.4.2                  
##  [13] magrittr_2.0.1                memoise_2.0.0                
##  [15] limma_3.46.0                  Biostrings_2.58.0            
##  [17] annotate_1.68.0               modelr_0.1.8                 
##  [19] bdsmatrix_1.3-4               askpass_1.1                  
##  [21] prettyunits_1.1.1             colorspace_2.0-0             
##  [23] apeglm_1.12.0                 blob_1.2.1                   
##  [25] rvest_1.0.0                   rappdirs_0.3.3               
##  [27] haven_2.3.1                   xfun_0.22                    
##  [29] hexbin_1.28.2                 crayon_1.4.1                 
##  [31] RCurl_1.98-1.3                jsonlite_1.7.2               
##  [33] tximport_1.18.0               genefilter_1.72.1            
##  [35] survival_3.2-10               glue_1.4.2                   
##  [37] gtable_0.3.0                  zlibbioc_1.36.0              
##  [39] XVector_0.30.0                DelayedArray_0.16.3          
##  [41] BiocSingular_1.6.0            scales_1.1.1                 
##  [43] mvtnorm_1.1-1                 DBI_1.1.1                    
##  [45] Rcpp_1.0.6                    emdbook_1.3.12               
##  [47] xtable_1.8-4                  progress_1.2.2               
##  [49] dqrng_0.2.1                   rsvd_1.0.3                   
##  [51] bit_4.0.4                     preprocessCore_1.52.1        
##  [53] httr_1.4.2                    RColorBrewer_1.1-2           
##  [55] ellipsis_0.3.1                pkgconfig_2.0.3              
##  [57] XML_3.99-0.6                  farver_2.1.0                 
##  [59] sass_0.3.1                    dbplyr_2.1.1                 
##  [61] locfit_1.5-9.4                utf8_1.2.1                   
##  [63] here_1.0.1                    reshape2_1.4.4               
##  [65] labeling_0.4.2                tidyselect_1.1.0             
##  [67] rlang_0.4.10                  later_1.1.0.1                
##  [69] AnnotationDbi_1.52.0          munsell_0.5.0                
##  [71] BiocVersion_3.12.0            cellranger_1.1.0             
##  [73] tools_4.0.2                   cachem_1.0.4                 
##  [75] cli_2.4.0                     generics_0.1.0               
##  [77] RSQLite_2.2.6                 broom_0.7.6                  
##  [79] evaluate_0.14                 fastmap_1.1.0                
##  [81] yaml_2.2.1                    knitr_1.32                   
##  [83] bit64_4.0.5                   fs_1.5.0                     
##  [85] AnnotationFilter_1.14.0       sparseMatrixStats_1.2.1      
##  [87] mime_0.10                     xml2_1.3.2                   
##  [89] biomaRt_2.46.3                BiocStyle_2.18.1             
##  [91] compiler_4.0.2                rstudioapi_0.13              
##  [93] curl_4.3                      interactiveDisplayBase_1.28.0
##  [95] affyio_1.60.0                 reprex_2.0.0                 
##  [97] geneplotter_1.68.0            bslib_0.2.4                  
##  [99] stringi_1.5.3                 highr_0.8                    
## [101] GenomicFeatures_1.42.3        lattice_0.20-41              
## [103] ProtGenerics_1.22.0           Matrix_1.3-2                 
## [105] vctrs_0.3.7                   pillar_1.6.0                 
## [107] lifecycle_1.0.0               BiocManager_1.30.12          
## [109] jquerylib_0.1.3               irlba_2.3.3                  
## [111] cowplot_1.1.1                 bitops_1.0-6                 
## [113] httpuv_1.5.5                  rtracklayer_1.50.0           
## [115] R6_2.5.0                      affy_1.68.0                  
## [117] promises_1.2.0.1              MASS_7.3-53.1                
## [119] assertthat_0.2.1              openssl_1.4.3                
## [121] rprojroot_2.0.2               withr_2.4.1                  
## [123] GenomicAlignments_1.26.0      Rsamtools_2.6.0              
## [125] GenomeInfoDbData_1.2.4        hms_1.0.0                    
## [127] beachmat_2.6.4                grid_4.0.2                   
## [129] coda_0.19-4                   DelayedMatrixStats_1.12.3    
## [131] rmarkdown_2.7                 bbmle_1.0.23.1               
## [133] numDeriv_2016.8-1.1           shiny_1.6.0                  
## [135] lubridate_1.7.10

Tiempo: 0.28

Introdución a RNA-seq

Ana BVA, Rodolfo

15 April 2021

Descripción general

Requisitos

Experimento

Archivos

Instalación

Instalación de FastQC

Instalación de MultiQC

Instalación de Salmon

Instalación de R y Rstudio

Instalación de Bioconductor

Instalación de la paquetería de R

5: Control de calidad de los datos

Pasos para correr FastQC:

MultiQC

6: Alineamiento (Salmon)

Generar índice del transcriptoma

Cuantificar la abundancia de los transcritos

7: Tximeta

Instalación de tximeta

Importar datos

8: Objeto SummarizedExperiment en R

ColData

rowRanges

Assay

¿Recuerdan a qué tiene que ser igual Rownames del slot colData?

9: Normalización de los datos (RPKM, TPM)

Métodos comunes

DESeq2

Factor de normalización en DESeq2

Ejercicio

Otras transformaciones

10: Exploración de los datos

PCA

¿Cómo crear nuestro propio PCA?

¿Cómo ven los resultados?

11: Análisis de expresión diferencial

12: Visualización de los resultados

MAplot

Volcano plot

Heatmap

Paquetería