¡OJO! Es importante actualizar y descargar el paquete DESeq2, ya que lo utilizaremos en esta práctica.
Datos
Utilizaremos las medidas de expresión génica del archivo salmon.merged.gene_counts.ts (descargado desde el aula virtual). Los datos proceden de 12 genotecas de RNA-seq de lecturas emparejadas, de muestras de Drosophila melanogaster de dos orígenes geográficos (Panamá y Maine) y dos tratamientos de temperatura (low y high) (Zhao 2015). Hay 3 réplicas biológicas para cada combinación de región y temperatura. Además, necesitáis el archivo dme_elev_samples.tsv, una tabla que indica a qué región y tratamiento pertenece cada muestra (descargado desde el aula virtual). Los dos archivos deben estar guardados en la carpeta de datos.
sample population temp
1 SRR1576457 maine low
2 SRR1576458 maine low
3 SRR1576459 maine low
4 SRR1576460 maine high
5 SRR1576461 maine high
6 SRR1576462 maine high
7 SRR1576463 panama low
8 SRR1576464 panama low
9 SRR1576465 panama low
10 SRR1576514 panama high
11 SRR1576515 panama high
12 SRR1576516 panama high
Ejercicio 1
Recuerda dejar constancia en un bloque de código de cómo obtienes las respuestas.
A) Usa la función paste() para añadir una columna a muestras que contenga la población y la temperatura unidas por _.
B) ¿Cuál es el número total de lecturas por muestra?
El dataframe genex tiene el dataframe del gen y una columna por cada muestra, los números que vemos son “counts”. Es un recuento del número de lecturas que mapean en cada uno de los genes/muestras proveniendo del RNASeq. Si pone un punto en el número es porque el programa no está “seguro” acerca del recuento y pone un decimal. No lo tenemos que aplicar a toda la sumas de columnas, ya que las dos primeras no son números, entonces no deben ser sumadas. Vemos que hay entre 20-30 millones de lecturas por muestra. Si hay 20.000 genes y entre 20-30 millones de lecturas por muestra nos hacemos la idea de lecturas por gen.
C) ¿Los recuentos tienen distribuciones parecidas entre muestras?
Para separar las muestras metemos las columnas una bajo de otra, para ello usamos pivot_longer. Enviamos los nombres a la variable muestra y los valores a la de expresion. Aparecen algunos errores porque son valores (de expresión) con cero, es decir, no se expresan los genes (y nosotros estamos usando logaritmos). A continuación vemos un violin plot.
Nota: he añadido en el siguiente código warning:false y message:false, ya que salía un error a la hora de representar.
Otra opción sería la función pairs acepta una matriz/dataframe de datos numéricos.
pairs(genex[, -c(1,2)])
La función pairs muestra la nube de puntos de cada pareja. La mayor parte de muestras tienen una expresión baja, se concentran todos en la esquina inferior izquierda. Tenemos muchos paneles, podríamos haberlo hecho condición por condición.
D) ¿Cuántos genes no se expresan en absoluto en ninguna de las muestras?
Igualamos a cero para crear un vector de verdadero y falso, así nos dicen en cuántos genes la expresión es cero. En este caso, hay 4722 genes sin expresión.
Para crear un objeto de clase DESeqDataSet necesitamos que las filas de la tabla muestras estén en el mismo orden que las columnas en genex. La función DESeqDataSetFromMatrix() necesita que los recuentos sean una matriz, no un marco de datos. Y en muestras debe haber factores.
# Pasamos la columna "sample" a nombres de fila:row.names(muestras) <- muestras$sample# Eliminamos la columna "sample":muestras$sample <-NULL# Y convertimos "population" y "temp" en factores:muestras$population <-factor(muestras$population,levels =c('maine', 'panama'))muestras$temp <-factor(muestras$temp,levels =c('low', 'high'))muestras
population temp nueva_columna
SRR1576457 maine low maine_low
SRR1576458 maine low maine_low
SRR1576459 maine low maine_low
SRR1576460 maine high maine_high
SRR1576461 maine high maine_high
SRR1576462 maine high maine_high
SRR1576463 panama low panama_low
SRR1576464 panama low panama_low
SRR1576465 panama low panama_low
SRR1576514 panama high panama_high
SRR1576515 panama high panama_high
SRR1576516 panama high panama_high
Aunque no hacía falta definir los niveles de los factores, este es un buen momento para definir manualmente su orden (por defecto se ordenan alfabéticamente). El orden de los niveles determina qué nivel actúa como referencia en los contrastes estadísticos (el primero). Esto resultaría importante, por ejemplo, al comparar tratamientos respecto de un control.
Los datos de recuento deben ser números enteros y en forma de matriz.
dds1 <-DESeqDataSetFromMatrix(countData = genex,colData = muestras,design =~ population + temp)
converting counts to integer mode
dds1
class: DESeqDataSet
dim: 24254 12
metadata(1): version
assays(1): counts
rownames(24254): 7SLRNA:CR32864 a ... RR51477 Bari1{}RR51478
rowData names(0):
colnames(12): SRR1576457 SRR1576458 ... SRR1576515 SRR1576516
colData names(3): population temp nueva_columna
Nota: los objetos DESeqDataSet pueden contener información sobre las coordenadas cromosómicas de los genes o tránscritos, si se construyen a partir de otro tipo de datos.
Pre-filtrado
Una fracción importante de genes no se están expresando en ningúna muestra. Aunque el filtrado no es necesario, ayuda a reducir la memoria ocupada por dds1:
Por defecto, la función results() produce la tabla de resultados para el contraste entre el nivel último y el de referencia del último factor en la tabla que describe las muestras. En este caso, el factor es temp, el nivel de referencia es low y el otro nivel es high. Por tanto, los valores de log2FoldChange representan el logaritmo en base 2 de la relación entre niveles de expresión en muestras high respecto de muestras low.
En referencia a estos resultados: log2foldchange - logaritmo en base 2 de la relación de cambio. Valores positivos son genes que a temperatura alta se expresan más que en baja. El nivel de expresión es mayor a temperatura baja para los valores negativos. Tenemos un dataframe con 19.508 filas (esto es porque hemos prefiltrado y hemos eliminado los que tienen una expresión muy baja).
Ejercicio 3
A) Mira la ayuda de results() y genera también la tabla de resultados para el contraste entre las dos regiones de procedencia de las moscas.
log2 fold change (MLE): population panama vs maine
Wald test p-value: population panama vs maine
DataFrame with 19508 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue
<numeric> <numeric> <numeric> <numeric> <numeric>
7SLRNA:CR32864 346.333 -0.3405225 1.2828671 -0.265439 0.7906716
a 555.596 -0.1061901 0.1079132 -0.984033 0.3250994
abd-A 531.996 0.0468148 0.0780163 0.600065 0.5484628
Abd-B 634.760 -0.0386942 0.0827328 -0.467701 0.6399985
Abl 1654.842 -0.1435019 0.0846880 -1.694478 0.0901746
... ... ... ... ... ...
RR51007 4.90906 1.0441599 2.913048 0.3584424 0.72001228
RR51048 16.60710 -1.2841124 0.397258 -3.2324379 0.00122739
RR51093 36.93279 -0.8907347 1.591194 -0.5597901 0.57562259
RR51475 159.62678 -0.0858284 1.937892 -0.0442896 0.96467360
RR51477 70.83394 -0.2603940 0.461897 -0.5637486 0.57292520
padj
<numeric>
7SLRNA:CR32864 0.966622
a 0.790762
abd-A 0.907772
Abd-B 0.928931
Abl 0.523390
... ...
RR51007 0.9498358
RR51048 0.0472064
RR51093 0.9181683
RR51475 0.9967946
RR51477 0.9169334
También lo podemos hacer así (no lo pongo en código de bloque de R para evitar sobreescribir los datos de antes): res2 <- results(dds1, contrast = c(‘population’, ‘panama’, ‘maine’))
Otro: p-valor ajustado para resolver el problema de los tests múltiples.
B) Usa summary() sobre los resultados para ver el número de genes diferencialmente expresados entre cada par de condiciones.
summary(res)
out of 19508 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 2418, 12%
LFC < 0 (down) : 2336, 12%
outliers [1] : 90, 0.46%
low counts [2] : 3018, 15%
(mean count < 2)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
Tenemos un 12% de genes up-regulated o sobreexpresados ene temperaturas altas sobre bajas, también un 12% de subexpresados. Outliers tenemos 90 y low counts 3018.
summary(res2)
out of 19508 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 350, 1.8%
LFC < 0 (down) : 351, 1.8%
outliers [1] : 46, 0.24%
low counts [2] : 2270, 12%
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
Tenemos porcentaje bastante equilibrado entre main y panama. La temperatura causa la sobreexpresión. Nota: En clase hemos probado este bloque de código pero me da distinto. Debería dar LFC > 0 de 623 y LFC < 0 de 634. Los outliers deberían ser de 90 y low count de 3018.
C) ¿Cuál es la tasa de falsos positivos (FDR)?
En summary podemos apreciar que el FDR es menor de 0,1 en ambos casos. Esto podemos conocerlo de la siguiente forma: en la función results tenemos un parametro alfa que es por defecto 0,1. Por su definición alfa es el umbral de significancia utilziado para filtrar el filtrado independiente. Si somos más exigentes y no lo toleramos y lo queremos hacer con 0,5 lo cambiamos. Este no ha sido el caso de esta práctica.
Gráficas
plotMA(res)
Este plot es simétrico y centrado en cero. Vemos el logFC según el nivel de expresión. Vemos que hay una dispersión grande.
Los genes con niveles de expresión más bajos producen resultados más ruidosos (la varianza es mayor). La función lfcShrink() contrae (shrink) las estimaciones de los efectos (log fold change) para reducir el ruido y hacer el resultado más realista: Aplicamos empirical bayes que es una forma de descomprimir. Lo comprime asumiendo que el logFC lo calcula con error. Vemos que ha disminuido la escala del logFC. Es una forma de ser más conservadores, estimamos la variación pero con muy pocas réplicas. Si asumimos que la tasas de cambio está entre 0,5 y -0,5, vemos que hay dispersión pero la media está muy concentrada. Suponemos que los valores elevados o bajos con pocas muestra se parecen más a la muestra, lo que hace es comprimirlo.
using 'normal' for LFC shrinkage, the Normal prior from Love et al (2014).
Note that type='apeglm' and type='ashr' have shown to have less bias than type='normal'.
See ?lfcShrink for more details on shrinkage type, and the DESeq2 vignette.
Reference: https://doi.org/10.1093/bioinformatics/bty895
plotMA(resShrink)
Ejercicio 4
¿Sabrías hacer una gráfica volcán?
res_df <-as.data.frame(res)volcano <-ggplot(res_df, aes(x = log2FoldChange, y =-log10(padj))) +geom_point(alpha =0.5) +labs(title ="Volcano Plot", x ="resShrink", y ="-log10(P-value)") +theme_minimal()print(volcano)
Extraer datos transformados
Aunque los datos de entrada no deben estar normalizados y la función DESeq() ya tiene en cuenta las diferencias de tamaño de genoteca entre muestras y la relación entre media y varianza, puede ser útil disponer de los datos normalizados para representarlos gráficamente. Para ello podemos usar las funciones vst() y assay():
¿Podrías explicar qué hace el bloque de código anterior?
vst tiene la función de estimar la tendencia de la dispersió (Quickly estimate dispersion trend and apply a variance stabilizing transformation).
assay contiene una matriz donde las filas representan objetos de interés y las columnas representan las muestras.
Nota: He consultado help(vst) y help(assay) respectivamente.
Ejercicio 6
A) En este enlace analizan exactamente los mismos datos con otro paquete de R: limma. Una de sus conclusiones es se sobreexpresan más genes en temperatura baja que en temperatura alta. ¿Confirma tu análisis esta conclusión?
Recapitulando con la primera parte de la práctica, yo he visualizado estos valores mediante summary(res). De esta forma, llegué a la conclusión de se están expresando más genes a temperatura alta: 2418 se sobreexpresan en temperatura alta 2336 se sobreexpresan en temperatura baja No obstante, el enlace no me redirige, así que no puedo comparar los resultados. Mi análisis se confirmaría en el caso de que en el otro análisis también se llegue a la conclusión de que hay más genes sobreexpresados en temperatura alta.
B) Estos son los 10 genes más diferencialmente expresados según aquel análisis entre las condiciones de baja y alta temperatura:
He elaborado la tabla a partir de los resultados que se ven en el guion del Aula Virtual
logFC
AveExpr
t
P.Value
adj.P.Val
B
Cpr92F
2.473527
6.201774
24.28881
1.306910e-13
1.739105e-09
21.29258
Vajk3
1.675184
4.029281
21.90298
6.069753e-13
3.524528e-09
19.66813
lncRNA:CR46017
-3.974938
1.044005
-21.50807
7.945880e-13
3.524528e-09
16.80878
CG7214
2.198970
7.945664
20.75054
1.349821e-12
4.490516e-09
19.14847
CG3739
-2.869172
6.514652
-20.49734
4.414988e-12
1.175005e-08
18.01975
Acp1
2.231834
7.638125
17.79441
2.089853e-11
3.997808e-08
16.54909
¿Cómo se comparan estos resultados con los tuyos?
Mediante el siguiente código voy a comprobar mis valores para estos genes:
log2 fold change (MLE): temp high vs low
Wald test p-value: temp high vs low
DataFrame with 5 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
Cpr92F 2718.979 -2.39126 0.1172256 -20.3988 1.71374e-92 1.40527e-88
Vajk3 513.411 -1.65423 0.0886464 -18.6610 1.02701e-77 3.36859e-74
Acp1 7042.192 -2.17394 0.1217046 -17.8624 2.31289e-71 4.21460e-68
CG3739 3912.044 2.93354 0.1497441 19.5903 1.86963e-85 7.66550e-82
CG7214 8615.934 -2.14594 0.1178994 -18.2015 5.02408e-74 1.17707e-70
A través de esto, y retomando el apartado A, veo que los resultados están “invertidos”. Esto lo observo mediante el FC, ya que cuando me da positivo a mi, en el análisis por limma da negativo. Esto quiere decir que yo para calcular el FC lo he hecho temp alta/temp bajo, mientras que en el analisis por limma se ha hecho temp baja/temp alta. Para poder comprobar los resultados probaríamos a calcular el FC de la misma forma que en el otro estudio.
