Esta práctica está basada en el flujo de trabajo o pipeline de Michael Love por lo que requerirémos de:
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:
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:
Corriendo el comando:
conda install -c bioconda <paquete>
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")
::install(version = 3.13) BiocManager
Instalado Bioconductor, descarga e instala los siguientes paquetes:
Corriendo el siguiente comando:
::install("paquete") BiocManager
Los paquetes que descargaremos desde CRAN son:
Corriendo el comando:
install.packages("paquete")
Para el análisis bioinformático vamos a seguir este pipeline:
Si decides realizar el ejercicio desde el inicio, empieza desde este punto. De lo contrario salta a la sección 4.2
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:
Pasos para correr FastQC
Abre una ventana de la terminal y ubicate en el directorio bin/
del presente repositorio
Activa tu ambiente de conda corriendo:
conda activate <ruta al directorio miniconda2 o minicoda3>
fastqc -h
fastqc ../data/*.fastq.gz -o ../Output/FastQC/
Para concatenar o juntar los resultados de todas las muestras analizadas por FastQC, usemos MultiQC
multiqc -h
.zip
generados por FastQCmultiqc ../Output/FastQC/*.zip -o ../Output/MultiQC/
Explora y revisa los archivos generados por FastQC y MultiQC
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.
Pasos para correr Cutadapt:
bin/
, tienes activo el ambiente de conda y cutadapt
se encuentra correctamente instalado:cutadapt -h
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.
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).
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
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
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:
salmon
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"
<- read.table(here("results/salmon_quants/metadata.txt"), he = T, stringsAsFactors = T, sep = "\t")
coldata
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
salmon
:<- dplyr::mutate(coldata, files = file.path(here("results/salmon_quants"), paste0(coldata$key, "_quant"), "quant.sf")) coldata
<- mutate(coldata, file_exist = file.exists(coldata$files))
coldata
2:8] coldata[,
## 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
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)<- mutate(coldata, names = paste0(coldata$key, "_quant")) coldata
tximeta
<- tximeta(coldata) se
## 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.
se
tiene clase de SummarizedExperiment y las cuentas están reportadas a nivel de transcrito (isoforma). Para obtener las cuentas a nivel de gen:<- summarizeToGene(se) gse
## 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))
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:
Previamente al análisis de expresión diferencial, es necesario explorar 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:
DESeq
(SummarizedExperiment) a partir del objeto gse
<- DESeqDataSet(se = gse, ##Objeto SummarizedExperiment que contiene la matriz de cuentas crudas
dds design = ~ Treatment) ##Variable con los grupos experimentales a contrastar. Forma parte de las columnas de coldata
## using counts and average transcript lengths from tximeta
##El primer nivel o elemento es el grupo control ("CT")
levels(dds$Treatment)
## [1] "CT" "mut_NRF2" "siRNA_control" "siRNA_NRF2"
$Treatment <- relevel(dds$Treatment, "siRNA_control")
dds
levels(dds$Treatment)
## [1] "siRNA_control" "CT" "mut_NRF2" "siRNA_NRF2"
##Recupera los genes que muestran más de 1 cuenta en al menos 3 muestras
<- rowSums(counts(dds) >= 1) >=3
keep
##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[keep, ] dds
<- rlog(dds, blind = F) rld
## using 'avgTxLength' from assays(dds), correcting for library size
PCAtools
:##Indicarle a la función pca que los datos deben de ser escalados para que sean más comparables entre sí
<- pca(assay(rld), metadata = colData(rld), scale = T)
pca
##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")
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.
El análisis de expresión diferencial en DESeq2
se encuentra completamente automatizado. Al correr la función DESeq()
, el programa:
El análisis de expresión diferencial en DESeq2
se lleva a cabo de la siguiente forma:
DESeq()
para normalizar, estimar la dispersión y ajustar los datos al modelo binomial negativo.<- DESeq(dds) 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
##Usemos la función results con los parámetros por default
<- results(dds)
res
##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:
…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
<- results(dds, lfcThreshold = 1, alpha = 0.05)
res_LFC
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()
Es una gráfica útil para visualizar la distribución de los genes que se encuentran diferencialmente expresados tomando en cuenta:
##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
<- lfcShrink(dds = dds, res = res, coef = "Treatment_siRNA_NRF2_vs_siRNA_control"
res_shrink 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))
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
<- as.data.frame(res_shrink)
res_shrink volcanoplotR(res_shrink, logfc = 0, p.adj = 0.1)
Utiliza la función glimmaVolcano()
para generar un volcano plot interactivo
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:
<- counts(dds, normalized = T) norm_counts
<- res_shrink %>% dplyr::filter(log2FoldChange > 0 & padj < 0.1 |
deg < 0 & padj < 0.1)
log2FoldChange
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
pheatmap()
obten la gráfica##Genera una paleta de colores
<- brewer.pal(n = 10, name = "RdBu")
RedBu <- rev(RedBu)
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)