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

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

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

4- Correr fastqc

##Imprimir el manual del programa usando el flag help
fastqc --help

1- Checar los datos fastq

pwd
ls ../data/Archivos_fastq

2- Crear un directorio para guardar el output

mkdir ../output/QC

3- Correr FastQC

# la opción -o es para elegir el directorio de salida
fastqc ../data/Archivos_fastq/*.fastq.gz -o ../output/QC/

4- Checar los archivos generados

#Buscar los archivos generados por fastqc en el directorio de resultados
ls ../output/QC

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.

1. La instalación se hace con conda, así que primero activamos un ambiente conda

conda activate base
pip install multiqc

1. Correr MultiQC

cd ../output/.
multiQC .

1. Checar el output.html

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 Transcriptome/

• Código para genear el índice

salmon index -t --genecode path_to_transcriptome.fa.gz -i path_to_save_index

-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

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

salmon quant -i path_to_index 
             -l A -1 path_to_R1 -2 path_to_R2 
             -o path_to_store_results

-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

Importamos la información de cada muestra

##Leemos los datos
coldata <- read.delim(here::here("output/salmon_quants/metadata.txt"))
##Generamos la columna "names"
coldata <- coldata %>% dplyr::rename(names = unique_id)
coldata

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

##Ahora las cuentas generadas por `Salmon`, para ello 
##Ubicamos el dir de trabajo
dir <- file.path(here::here("output/salmon_quants"))
##Verificar que estén los directorios  
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?? Verificar que el path se generó correctamente
#file.path(dir,paste0(info$Sample,"_quant"),"quant.sf")
coldata$files <- file.path(dir,paste0(coldata$Sample,"_quant"),"quant.sf")
##Corroboramos 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

Importar los datos usando tximeta

#BiocManager::install("tximeta")
library("tximeta")
se <- tximeta(coldata)

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

• Para acceder a ColData usar el siguiente comando:

colData(se)

## DataFrame with 12 rows and 6 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
##               Condition
##                <factor>
## A_Control_1     Control
## A_Control_2     Control
## A_Control_3     Control
## M_Control_1     Control
## M_Control_2     Control
## ...                 ...
## A_Verafinib_2 Verafinib
## A_Verafinib_3 Verafinib
## M_Verafinib_1 Verafinib
## M_Verafinib_2 Verafinib
## M_Verafinib_3 Verafinib

• 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

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

assayNames(se)

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

• 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

¿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"

colnames(assay(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"

## 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

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

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

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

## # 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

1. Se calcula la fración muestra/pseudo-referencia

## # 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

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

## [1] 1.230235

## [1] 0.8222689

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

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

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

Se pueden ocupar otras transformaciones en la paquetería de DESeq2 pero 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")

# Definir los grupos que se desean comparar
se$Condition <- factor(se$Condition)
# Checar cantidad de muestras
table(se$Condition)
# Generar el objeto de DESeq 
dds <- DESeqDataSet(se, design = ~ Condition)
# Prefiltrado para quitar genes con cuentas bajas
keep <- rowSums(counts(dds)) >= 6
dds <- dds[keep, ]
# Funcion que normaliza los datos y realiza el analisis de expresión differencial
dds <- DESeq(dds)

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

# 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)

Primero instalen las siguientes librerías:

library(DESeq2)
library(PCAtools)

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

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

Si tuvieramos un número mayor de muestras (n > 30) entonces conviene:

vsd <- vst(dds, blind = T)

• 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")

DESeq2 evalúa los 500 genes con mayor varianza del set de datos

• 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)

Generemos un MDS plot interactivo con los datos. Para ello usemos la librería de Glimma y la función de GlimmaMDS

glimmaMDS(dds)

La librería de DESeq2 normaliza y realiza el análisis de expresión diferencial en una sola función. Así que ya realizamos el análisis, vamos a verlo:

library("DESeq2")
library("dplyr")
library("ggplot2")

Para ver los resultados podemos usar la función results

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

## log2 fold change (MLE): Condition Verafinib vs Control 
## Wald test p-value: Condition 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

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: Cond..
## stat                results Wald statistic: Cond..
## pvalue              results Wald test p-value: C..
## padj                results   BH adjusted p-values

Podemos hacer cortes:

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

## 
## FALSE  TRUE 
## 34309  1156

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

M eje y de minus:

A eje x de average -> Promedio de las cuentas normalizadas para cada gen en todas las muestras

Para generar el MAplot usemos el siguiente comando:

##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)

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

##Para los objetos creados con DESeq2 es importante que incluyas un nivel llamado group
group <- info$Treatment
dds$group <- group
glimmaMA(dds)