library(readr)
  library(DESeq2)
  library(GEOquery)
  library(canvasXpress)
  library(ggplot2)
  library(clinfun)
  library(GGally)
  library(factoextra)
  library(tidyr)
  library(tidyverse)
  library(dplyr)
  library(tibble)
  library(apeglm)
  library(EnhancedVolcano)
  library(pheatmap)
  library(biomartr)
  library(clusterProfiler)
  library(enrichplot)
  library(org.At.tair.db)
  library(biomaRt)
  library(gplots)

INPUT DATA

# download.file("https://cursoronline.usal.es/GSE141547_count_matrix.txt","GSE141547_count_matrix.txt")
raw_counts <- read.delim("GSE141547_count_matrix.txt", 
                         header = T, 
                         stringsAsFactors = F) %>% 
  column_to_rownames("Gene") 

# download.file("https://cursoronline.usal.es/coldata.csv","coldata.csv")
xp_design <- read.csv("coldata.csv", sep = ";", header = T, 
                      row.names = 1,
                      stringsAsFactors = F)

Compruebo que los nombres de las muestras son iguales en ambos archivos.

all(colnames(raw_counts) %in% xp_design$group)
## [1] TRUE
all(colnames(raw_counts) == xp_design$group)
## [1] TRUE

Análisis de componentes principales (PCA)

Creo el objeto de clase DESeqDataSet (objeto S4 DESeq2).

dds <- DESeqDataSetFromMatrix(countData = raw_counts, 
                              colData = xp_design, 
                              design = ~ 1)          # sin especificar modelo

Gráfico de la distribucion de conteo por cada muestra.

as_tibble(assay(dds), rownames = "gene") %>%
  pivot_longer(names_to = "group", values_to = "counts", cols = 2:13) %>%
  ggplot(aes(x = log2(counts + 1), fill = group)) +
  geom_histogram(bins = 20) +
  facet_wrap(~ group)

Transformación estabilizadora de varianza (vst).

Primero hay que hacer la estimación de factores de tamaño y dispersiones.

dds = estimateSizeFactors(dds)
dds = estimateDispersions(object = dds, 
                          fitType = "parametric", 
                          quiet = TRUE)

vsd = varianceStabilizingTransformation(object = dds)

Genero una matriz de los conteos transformados.

variance_stabilised_counts <- assay(vsd)

Visualización de la relación media - desviación estándar.

p_mean_sd_vst <- 
  variance_stabilised_counts %>% 
  as.data.frame() %>% 
  rownames_to_column("gene") %>% 
  pivot_longer(cols = - gene, names_to = "sample", values_to = "counts") %>% 
  group_by(gene) %>% 
  summarise(gene_average = mean(counts), gene_stdev = sd(counts)) %>% 
  ungroup() %>% 
  ggplot(., aes(x = gene_average, y = gene_stdev)) +
  geom_point(alpha = 0.5, fill = "grey", colour = "black") +
  labs(x = "Promedio del recuento de genes (varianza estabilizada)", 
       y = "Desviación estándar del recuento de genes (varianza estabilizada)") +
  ggtitle("Relación media - desviación estándar\n(tras estabilización de varianza)")
p_mean_sd_vst

Con los datos transformados se realiza el PCA. Cargamos la función “mypca”() para generar gráficos PCA.

mypca <- function(x, center = TRUE, scale = TRUE){
  # Muestras deben ir en filas
  # Variables deben ir en columnas
  
  # elimina variables constantes
  constant_val = apply(x,2,'sd')
  x_reduced = x[,constant_val>0]
  
  # realiza SVD
  SVD <- svd(scale(x_reduced,center = center, scale = scale))
  
  # crea data frame de scores
  scores <- as.data.frame(SVD$u %*% diag(SVD$d))
  rownames(scores) <- rownames(x)
  colnames(scores) <- paste0("PC", c(1:dim(scores)[2]))
  
  # crea data frame de 'loadings'
  loadings <- data.frame(SVD$v)
  colnames(loadings) <- paste0("PC", c(1:dim(loadings)[2]))
  rownames(loadings) <- colnames(x_reduced)
  
  # crea data frame de 'explained variances'
  explained_var <- as.data.frame(round((SVD$d^2) / sum(SVD$d^2)*100, digits = 1))
  rownames(explained_var) <- paste0("PC", c(1:dim(loadings)[2]))
  colnames(explained_var) <- "exp_var"
  
  # return resultados
  return (list("scores" = scores, "loadings" = loadings, "explained_var" = explained_var))
}

Hacemos el PCA usando mypca(), que devuelve una lista con 3 objetos: scores, loadings y explained_var dataframes. En esta función, las muestras deben ir en filas, y las variables en columnas. Transponemos los datos, cambiando filas por columnas:

t_variance_stabilised_counts <- t(variance_stabilised_counts)
dim(t_variance_stabilised_counts)
## [1]    12 21543
pca_results <- mypca(t_variance_stabilised_counts, 
                     center = TRUE, 
                     scale = TRUE)

Hacemos un screen plot con el dataframe ‘explained variance’.

ggplot(pca_results$explained_var, 
       aes(x = seq(from = 1, to = nrow(pca_results$explained_var)), 
           y = exp_var)) +
  ylab('explained variance (%)') + 
  ggtitle('Explained variance per component') + 
  geom_bar(stat = "identity") +
  labs(x = "Principal Component number") +
  scale_x_continuous(breaks = seq(
    from = 1, 
    to = nrow(pca_results$explained_var)))

cumsum(pca_results$explained_var)[2,1] # PC1 y PC2 explican el 53.5% de la varianza
## [1] 53.5
cumsum(pca_results$explained_var)[3,1] # PC1, PC2 y PC3 explican el 70.5% de la varianza
## [1] 70.5

Extraemos los scores de los resultados del PCA.

scores <- pca_results$scores
scores[1:5,1:5]
##               PC1       PC2       PC3       PC4       PC5
## Col0_C1 -83.04716  29.37289 -29.88341 -70.48800 -10.73317
## Col0_C2 -64.49673 -63.67578  31.84949 -51.93314 -19.40872
## Col0_C3 -78.81453 -58.02116  31.14457 -45.66848  19.97232
## Col0_P1  51.16797 -30.27315 -94.82714  22.02300 -36.29980
## Col0_P2  25.43227 -88.40931 -49.74875  43.45586 -10.63132

La matriz de scores contiene las muestras en filas con sus coordenadas en diferentes PCs. Añadimos los factores del diseño experimental al dataframe de scores.

scores_with_conditions <- 
  scores %>% 
  rownames_to_column("group") %>% 
  left_join(x = .,                 # cargamos el dataframe 'scores'
            y = xp_design,         # cargamos el dataframe del diseño experimental
            by = "group")
scores_with_conditions[1:5, 11:15]
##         PC10      PC11         PC12  sample Genotype
## 1  -9.976663  1.498587 4.934371e-13 sample1      Col
## 2   9.236434  7.181170 4.934371e-13 sample2      Col
## 3   5.896607 -8.330307 4.934371e-13 sample3      Col
## 4  45.724442 -3.714991 4.934371e-13 sample4      Col
## 5 -23.077050 -1.606868 4.934371e-13 sample5      Col
explained_variance <- 
  pca_results$explained_var %>% 
  pull("exp_var")

PCA con las 2 variables (Genotype, Condition).

ggplot(scores_with_conditions, 
          aes(PC1, PC2, color = Condition, shape = Genotype)) +
  geom_point(size = 4) +
  xlab(paste0("PC1: ",explained_variance[1],"% variance")) +
  ylab(paste0("PC2: ",explained_variance[2],"% variance")) + 
  coord_fixed(ratio = 1) +
  ggtitle("PCA")

Las muestras de la misma condición han sido agrupadas juntas en el PCA. Muestras de diferentes condiciones se han agrupado por separado. PC1 y PC2 explican gran porcentaje de la variación total presente en el dataset. La variable ‘Condition’ explica gran parte de la variabilidad en los datos.

Heatmap de todas las muestras

dds <- DESeqDataSetFromMatrix(countData = raw_counts, 
                              colData = xp_design, 
                              design = ~Genotype * Condition)   # diseño con interacción
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
dds$Genotype
##  [1] Col  Col  Col  Col  Col  Col  med8 med8 med8 med8 med8 med8
## Levels: Col med8
dds$Condition
##  [1] CT CT CT MV MV MV CT CT CT MV MV MV
## Levels: CT MV
resultsNames(dds)
## [1] "Intercept"                "Genotype_med8_vs_Col"    
## [3] "Condition_MV_vs_CT"       "Genotypemed8.ConditionMV"
res <- results(dds)
all_genes_results <- results(dds, contrast = c("Condition", "MV", "CT"))

Filtrado DEGs por p-value < 0.01 y log2foldchange > 1 y < -1

diff_genes_lfc1 = all_genes_results %>% 
  as.data.frame() %>% 
  rownames_to_column("genes") %>% 
  filter(padj < 0.01 & abs(log2FoldChange) > 1) %>% 
  arrange(desc(log2FoldChange), 
          desc(padj))
dim(diff_genes_lfc1) 
## [1] 1314    7

Obtener la matriz counts_scaled, con los conteos normalizados

Cargamos la funcion para normalizar.

mor_normalization = function(data){
  
  # take the log
  log_data = log(data) 
  
  # find the psuedo-references per sample by taking the geometric mean
  log_data = log_data %>% 
    rownames_to_column('gene') %>% 
    mutate (gene_averages = rowMeans(log_data)) %>% 
    filter(gene_averages != "-Inf")
  
  # the last columns is the pseudo-reference column 
  pseudo_column = ncol(log_data)
  
  # where to stop before the pseudo column 
  before_pseduo = pseudo_column - 1
  
  # find the ratio of the log data to the pseudo-reference
  ratios = sweep(log_data[,2:before_pseduo], 1, log_data[,pseudo_column], "-")
  
  # find the median of the ratios
  sample_medians = apply(ratios, 2, median)
  
  # convert the median to a scaling factor
  scaling_factors = exp(sample_medians)
  
  # use scaling factors to scale the original data
  manually_normalized = sweep(data, 2, scaling_factors, "/")
  return(manually_normalized)
}

counts_normalised_only_diff_genes = 
  mor_normalization(raw_counts) %>%             # normalización con la función creada
  rownames_to_column("genes") %>%               
  pivot_longer(- genes,                         
               names_to = "sample", 
               values_to = "counts") %>% 
  filter(genes %in% diff_genes_lfc1$genes) %>%     # filtrado, para quedarnos solo con diff_genes        
  pivot_wider(names_from = "sample",            
              values_from = "counts")  %>%      
  column_to_rownames("genes") 
dim(counts_normalised_only_diff_genes) 
## [1] 1314   12
counts_scaled = 
  counts_normalised_only_diff_genes %>% 
  t(.) %>%                               
  scale() %>%                            
  t(.)

Heatmap con clustering.

pheatmap(counts_scaled, 
         cluster_rows = TRUE,                      
         cluster_cols = TRUE, 
         show_rownames = FALSE, 
         show_colnames = TRUE,
         main = "Heatmap")

Añadir barras de colores con las condiciones.

xp_design_filtered = xp_design[,2:4]
anno_col_info = xp_design_filtered %>% column_to_rownames("group")
anno_info_colors = list(
  Condition = c(CT = "lightgrey", 
                MV = "black"))
pheatmap(counts_scaled, 
         cluster_rows = TRUE,                       
         cluster_cols = TRUE, 
         show_rownames = FALSE, 
         show_colnames = TRUE,
         annotation_col = anno_col_info,
         annotation_colors = anno_info_colors,
         main = "Clustering de genes diferencialmente expresados")

Análisis de expresión diferencial: comparaciones a pares

[Col_MV vs Col_CT] -> DEGs en Col_MV

head(raw_counts) 
##           Col0_C1 Col0_C2 Col0_C3 Col0_P1 Col0_P2 Col0_P3 med8_C1 med8_C2
## AT1G01010     209     298     292     443     457     417     204     217
## AT1G01020     148     161     148     205     163     155     176     159
## AT1G01030      81      83      74      78      71      65      64      75
## AT1G01040    1809    1711    1506    2030    1817    1948    1672    1676
## AT1G01046       6       6       7      12      15       9       6      13
## AT1G01050    1675    1810    1842    2480    2228    2305    1789    1780
##           med8_C3 med8_P1 med8_P2 med8_P3
## AT1G01010     118     394     458     449
## AT1G01020     142     138     151     217
## AT1G01030     102      44      73      75
## AT1G01040    1932    1943    2047    2155
## AT1G01046       7       9      11      15
## AT1G01050    1807    2155    2173    2295
head(xp_design)
##    sample   group Genotype Condition
## 1 sample1 Col0_C1      Col        CT
## 2 sample2 Col0_C2      Col        CT
## 3 sample3 Col0_C3      Col        CT
## 4 sample4 Col0_P1      Col        MV
## 5 sample5 Col0_P2      Col        MV
## 6 sample6 Col0_P3      Col        MV

Selección de las muestras Col (CT y MV).

xp_design_CT_vs_MV = xp_design %>% 
  filter(Genotype == "Col")

Selección en la matriz de conteo de acuerdo a las muestras seleccionadas en el diseño experimental.

raw_counts_filtered = raw_counts[, colnames(raw_counts) %in% xp_design_CT_vs_MV$group]

all(colnames(raw_counts_filtered) %in% xp_design_CT_vs_MV$group)
## [1] TRUE
all(colnames(raw_counts_filtered) == xp_design_CT_vs_MV$group)
## [1] TRUE

Objeto DESeq2 para realizar la comparacion del efecto de la condición sobre plantas Col.

dds <- DESeqDataSetFromMatrix(countData = raw_counts_filtered, 
                              colData = xp_design_CT_vs_MV, 
                              design = ~ Condition)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
dds
## class: DESeqDataSet 
## dim: 21543 6 
## metadata(1): version
## assays(1): counts
## rownames(21543): AT1G01010 AT1G01020 ... AT5G67630 AT5G67640
## rowData names(0):
## colnames(6): Col0_C1 Col0_C2 ... Col0_P2 Col0_P3
## colData names(4): sample group Genotype Condition
dds$Condition
## [1] CT CT CT MV MV MV
## Levels: CT MV
dds$Condition <- relevel(dds$Condition, ref = "CT")
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
resultsNames(dds)
## [1] "Intercept"          "Condition_MV_vs_CT"
res <- results(dds)
res
## log2 fold change (MLE): Condition MV vs CT 
## Wald test p-value: Condition MV vs CT 
## DataFrame with 21543 rows and 6 columns
##             baseMean log2FoldChange     lfcSE       stat      pvalue
##            <numeric>      <numeric> <numeric>  <numeric>   <numeric>
## AT1G01010  347.10841     0.52698620  0.160349  3.2864982  0.00101441
## AT1G01020  162.57772    -0.00564394  0.176636 -0.0319524  0.97451002
## AT1G01030   75.70065    -0.35213995  0.232689 -1.5133514  0.13019042
## AT1G01040 1797.56339     0.00557764  0.106464  0.0523900  0.95821794
## AT1G01046    8.99624     0.72812682  0.656379  1.1093078  0.26729741
## ...              ...            ...       ...        ...         ...
## AT5G67600   900.2656      0.2043385  0.129693    1.57556 1.15128e-01
## AT5G67610   489.4217     -0.0989258  0.110191   -0.89777 3.69308e-01
## AT5G67620    53.4717      0.4373731  0.382165    1.14446 2.52432e-01
## AT5G67630   744.7855      0.6748510  0.122901    5.49102 3.99608e-08
## AT5G67640   173.8583      0.1232257  0.205233    0.60042 5.48226e-01
##                  padj
##             <numeric>
## AT1G01010  0.00576689
## AT1G01020  0.98720468
## AT1G01030  0.28428510
## AT1G01040  0.97930588
## AT1G01046  0.46641178
## ...               ...
## AT5G67600 2.60283e-01
## AT5G67610 5.74769e-01
## AT5G67620 4.48503e-01
## AT5G67630 6.29697e-07
## AT5G67640 7.28307e-01
res0.01 <- results(dds, alpha = 0.01)
summary(res0.01)
## 
## out of 21542 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 0 (up)       : 2160, 10%
## LFC < 0 (down)     : 1995, 9.3%
## outliers [1]       : 1, 0.0046%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
dim(res)          # 21543 genes 
## [1] 21543     6

False Discovery Rates

Filtro por threshold de padj = 0.01

res %>% 
  as.data.frame() %>% 
  filter(padj < 0.01) %>% 
  dim()                         # 4155 genes 
## [1] 4155    6

Filtro por threshold de padj = 0.001

res %>% 
  as.data.frame() %>% 
  filter(padj < 0.001) %>% 
  dim()                         # 2890 genes
## [1] 2890    6

Distribution of adjusted p-values.

hist(res$padj, col="lightblue", main = "Distribución del p-value ajustado")

Distribution of non-adjusted p-values.

hist(res$pvalue, col="grey", main = "Distribución del p-value no ajustado")

La distribucion de los p-values es bastante similar, indicando que una buena proporcion de los tests tienen un p-value significativo. Esto sugiere que la mayor parte serán ‘true positives’ (DEGs reales).

Obtención tablas DEGs todas las comparaciones

[Col_MV vs Col_CT] -> DEGs en Col_MV Filtrados por p-value < 0.01 y log2foldchange > 1 y < -1 (up y down-regulated).

DEG_ColMV_vs_ColCT = res %>% 
  as.data.frame() %>% 
  rownames_to_column("genes") %>% 
  filter(padj < 0.01 & abs(log2FoldChange) > 1) %>% 
  arrange(desc(log2FoldChange), 
          desc(padj))
dim(DEG_ColMV_vs_ColCT) 
## [1] 1263    7

Volcano plot # reducción del log2foldchange para eliminar ruido.

resLFC <- lfcShrink(dds = dds, 
                    res = res,
                    type = "apeglm",
                    coef = "Condition_MV_vs_CT")
## 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
EnhancedVolcano(toptable = resLFC,
                x = "log2FoldChange",
                y = "padj",
                lab = rownames(resLFC),
                xlim = c(-8, +8),
                ylim = c(0,120),
                pCutoff = 1e-06,
                FCcutoff = 2,
                title = "Volcano Col_MV vs Col_CT"
)

MA plot

plotMA(resLFC, alpha = 0.01)

[med8_CT vs Col_CT] -> DEGs en med8_CT

xp_design_med8_vs_Col = xp_design %>% 
  filter(Condition == "CT")
raw_counts_filtered = raw_counts[, colnames(raw_counts) %in% xp_design_med8_vs_Col$group]
all(colnames(raw_counts_filtered) %in% xp_design_med8_vs_Col$group)
## [1] TRUE
all(colnames(raw_counts_filtered) == xp_design_med8_vs_Col$group)
## [1] TRUE
dds <- DESeqDataSetFromMatrix(countData = raw_counts_filtered, 
                              colData = xp_design_med8_vs_Col, 
                              design = ~ Genotype)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
dds$Genotype
## [1] Col  Col  Col  med8 med8 med8
## Levels: Col med8
dds$Genotype <- relevel(dds$Genotype, ref = "Col")
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
resultsNames(dds)
## [1] "Intercept"            "Genotype_med8_vs_Col"
res <- results(dds)

res0.01 <- results(dds, alpha = 0.01)
summary(res0.01)
## 
## out of 21539 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 0 (up)       : 264, 1.2%
## LFC < 0 (down)     : 304, 1.4%
## outliers [1]       : 1, 0.0046%
## low counts [2]     : 418, 1.9%
## (mean count < 3)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
DEG_med8CT_vs_ColCT = res %>% 
  as.data.frame() %>% 
  rownames_to_column("genes") %>% 
  filter(padj < 0.01 & abs(log2FoldChange) > 1) %>% 
  arrange(desc(log2FoldChange), 
          desc(padj))
dim(DEG_med8CT_vs_ColCT) 
## [1] 164   7

Volcano plot

resLFC <- lfcShrink(dds = dds, 
                    res = res,
                    type = "apeglm",
                    coef = "Genotype_med8_vs_Col")
## 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
EnhancedVolcano(toptable = resLFC,
                x = "log2FoldChange",
                y = "padj",
                lab = rownames(resLFC),
                xlim = c(-8, +8),
                ylim = c(0,120),
                pCutoff = 1e-06,
                FCcutoff = 2,
                title = "Volcano med8_CT vs Col_CT"
)

MA plot

plotMA(resLFC, alpha = 0.01)

[med8_MV vs med8_CT] -> DEGs en med8_MV

xp_design_med8MV_vs_med8CT = xp_design %>% 
  filter(Genotype == "med8")

raw_counts_filtered = raw_counts[, colnames(raw_counts) %in% xp_design_med8MV_vs_med8CT$group]
all(colnames(raw_counts_filtered) %in% xp_design_med8MV_vs_med8CT$group)
## [1] TRUE
all(colnames(raw_counts_filtered) == xp_design_med8MV_vs_med8CT$group)
## [1] TRUE
dds <- DESeqDataSetFromMatrix(countData = raw_counts_filtered, 
                              colData = xp_design_med8MV_vs_med8CT, 
                              design = ~ Condition)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
dds$Condition <- relevel(dds$Condition, ref = "CT")
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
resultsNames(dds)
## [1] "Intercept"          "Condition_MV_vs_CT"
res <- results(dds)

res0.01 <- results(dds, alpha = 0.01)
summary(res0.01)
## 
## out of 21541 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 0 (up)       : 2917, 14%
## LFC < 0 (down)     : 3356, 16%
## outliers [1]       : 4, 0.019%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
DEG_med8MV_vs_med8CT = res %>% 
  as.data.frame() %>% 
  rownames_to_column("genes") %>% 
  filter(padj < 0.01 & abs(log2FoldChange) > 1) %>% 
  arrange(desc(log2FoldChange), 
          desc(padj))
dim(DEG_med8MV_vs_med8CT)
## [1] 1889    7

Volcano plot

resLFC <- lfcShrink(dds = dds, 
                    res = res,
                    type = "apeglm",
                    coef = "Condition_MV_vs_CT")
## 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
EnhancedVolcano(toptable = resLFC,
                x = "log2FoldChange",
                y = "padj",
                lab = rownames(resLFC),
                xlim = c(-8, +8),
                ylim = c(0,120),
                pCutoff = 1e-06,
                FCcutoff = 2,
                title = "Volcano med8_MV vs med8_CT"
)

MA plot

plotMA(resLFC, alpha = 0.01)

[med8_MV vs Col_MV] -> DEGs en med8_MV vs Col_MV

xp_design_med8MV_vs_ColMV = xp_design %>% 
  filter(Condition == "MV")

raw_counts_filtered = raw_counts[, colnames(raw_counts) %in% xp_design_med8MV_vs_ColMV$group]
all(colnames(raw_counts_filtered) %in% xp_design_med8MV_vs_ColMV$group)
## [1] TRUE
all(colnames(raw_counts_filtered) == xp_design_med8MV_vs_ColMV$group)
## [1] TRUE
dds <- DESeqDataSetFromMatrix(countData = raw_counts_filtered, 
                              colData = xp_design_med8MV_vs_ColMV, 
                              design = ~ Genotype)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
dds$Genotype <- relevel(dds$Genotype, ref = "Col")
dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
resultsNames(dds)
## [1] "Intercept"            "Genotype_med8_vs_Col"
res <- results(dds)

res0.01 <- results(dds, alpha = 0.01)
summary(res0.01)
## 
## out of 21542 with nonzero total read count
## adjusted p-value < 0.01
## LFC > 0 (up)       : 2008, 9.3%
## LFC < 0 (down)     : 2519, 12%
## outliers [1]       : 0, 0%
## low counts [2]     : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
DEG_med8MV_vs_ColMV = res %>% 
  as.data.frame() %>% 
  rownames_to_column("genes") %>% 
  filter(padj < 0.01 & abs(log2FoldChange) > 1) %>% 
  arrange(desc(log2FoldChange), 
          desc(padj))
dim(DEG_med8MV_vs_ColMV)
## [1] 1270    7

Volcano plot

resLFC <- lfcShrink(dds = dds, 
                    res = res,
                    type = "apeglm",
                    coef = "Genotype_med8_vs_Col")
## 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
EnhancedVolcano(toptable = resLFC,
                x = "log2FoldChange",
                y = "padj",
                lab = rownames(resLFC),
                xlim = c(-8, +8),
                ylim = c(0,120),
                pCutoff = 1e-06,
                FCcutoff = 2,
                title = "Volcano med8_MV vs Col_MV"
)

MA plot

plotMA(resLFC, alpha = 0.01)

Visualización de los DEGs obtenidos

Diagrama de Venn

dim(DEG_ColMV_vs_ColCT)
## [1] 1263    7
dim(DEG_med8CT_vs_ColCT)
## [1] 164   7
dim(DEG_med8MV_vs_med8CT)
## [1] 1889    7
dim(DEG_med8MV_vs_ColMV)
## [1] 1270    7

Hay que introducir listas de genes.

DEG_1 <- DEG_ColMV_vs_ColCT [1:1263, 1]
DEG_2 <- DEG_med8CT_vs_ColCT [1:164, 1]
DEG_3 <- DEG_med8MV_vs_med8CT [1:1889, 1]
DEG_4 <- DEG_med8MV_vs_ColMV [1:1270, 1]

x <- list(
  ColMV = DEG_1, 
  med8CT = DEG_2, 
  med8MV = DEG_3,
  med8MVvsCT = DEG_4
)
library(ggvenn)
ggvenn(
  x, 
  fill_color = c("#0073C2FF", "#EFC000FF", "#868686FF", "#CD534CFF"),
  stroke_size = 0.5, set_name_size = 4
)

Se pueden mostrar sólo las listas de interés, para encontrar genes comunes

ggvenn(x,c("ColMV","med8MV"),show_percentage=T, 
       fill_color=c("red","orange"))

ggvenn(x,c("ColMV","med8MVvsCT"),show_percentage=T, 
       fill_color=c("red","orange"))

Heatmaps

Se empieza con los raw counts, y se normalizan (método “median of ratios”). Usaremos la función customizada para normalización (mor_normalization). Generamos tabla de datos normalizados y filtrados (por los genes de la tabla DEG_ColMV_vs_ColCT).

counts_normalised_only_diff_genes = 
  mor_normalization(raw_counts) %>%             
  rownames_to_column("genes") %>%               
  pivot_longer(- genes,                         
               names_to = "sample", 
               values_to = "counts") %>% 
  filter(genes %in% DEG_ColMV_vs_ColCT$genes) %>%             
  pivot_wider(names_from = "sample",            
              values_from = "counts")  %>%      
  column_to_rownames("genes")                   
dim(counts_normalised_only_diff_genes)  
## [1] 1263   12
counts_scaled = 
  counts_normalised_only_diff_genes %>% 
  t(.) %>%                               
  scale() %>%                           
  t(.)

Filtrado para quedarnos con las 6 columnas que se corresponden con xp_design_CT_vs_MV.

counts_scaled_filtered = 
  counts_scaled %>% 
  as.data.frame() %>%
  dplyr::select(xp_design_CT_vs_MV$group) 
head(counts_scaled_filtered)
##              Col0_C1    Col0_C2    Col0_C3    Col0_P1    Col0_P2    Col0_P3
## AT1G01060 -0.9997198 -1.0647231 -0.9824216  1.6055872  0.9838184  0.7512321
## AT1G01190 -0.1124159 -0.5307067 -0.2462124 -1.1679173 -1.4619436 -1.2630718
## AT1G01520 -0.9865496 -0.9731675 -0.9093698  2.0663596  1.4244424  1.1433754
## AT1G01570 -0.6524403 -0.4173666 -0.6250506  2.0938080  1.7416185  0.9382720
## AT1G01720 -0.8721725 -0.5413346 -0.6267859  1.4278793  1.6597364  1.2348087
## AT1G02065  0.5530190 -0.2562367  1.7379532 -0.6898748 -0.4803629 -0.6889708
xp_design_CT_vs_MV_filtered = xp_design_CT_vs_MV[,2:4]

anno_col_info = xp_design_CT_vs_MV_filtered %>% column_to_rownames("group")
anno_info_colors = list(
  Genotype = c(Col = "#d8b365"),
  Condition = c(CT = "lightgrey", 
                MV = "black"))
pheatmap(counts_scaled_filtered, 
         cluster_rows = TRUE,                       
         cluster_cols = TRUE, 
         show_rownames = FALSE, 
         show_colnames = TRUE,
         annotation_col = anno_col_info,
         annotation_colors = anno_info_colors,
         main = "Clustering DEGs Col_MV vs Col_CT")

Otra comparación: DEG_med8CT_vs_ColCT.

counts_normalised_only_diff_genes = 
  mor_normalization(raw_counts) %>%             
  rownames_to_column("genes") %>%               
  pivot_longer(- genes,                         
               names_to = "sample", 
               values_to = "counts") %>% 
  filter(genes %in% DEG_med8CT_vs_ColCT$genes) %>%             
  pivot_wider(names_from = "sample",            
              values_from = "counts")  %>%      
  column_to_rownames("genes")                   
dim(counts_normalised_only_diff_genes)  
## [1] 164  12
counts_scaled = 
  counts_normalised_only_diff_genes %>% 
  t(.) %>%                               
  scale() %>%                           
  t(.)

Filtrado para quedarnos con las 6 columnas que se corresponden con xp_design_CT_vs_MV.

counts_scaled_filtered = 
  counts_scaled %>% 
  as.data.frame() %>%
  dplyr::select(xp_design_med8_vs_Col$group) 
head(counts_scaled_filtered)
##               Col0_C1     Col0_C2     Col0_C3    med8_C1    med8_C2     med8_C3
## AT1G01520 -0.98654955 -0.97316748 -0.90936979 -0.4395027 -0.5112365 -0.07829492
## AT1G03820 -0.28358506  1.78973736  0.20113660 -1.4210811 -1.6619227 -1.04997813
## AT1G05300 -0.22400173  0.25051240  0.73385540 -0.8533846 -0.9227341 -1.04633476
## AT1G07610  1.27873414  0.90770052  1.30277772 -1.0651969 -0.9198594 -2.06935875
## AT1G10070  0.32570143  2.32261914  1.59495368 -0.1763710 -0.2015724 -0.65431891
## AT1G11070  0.03119132 -0.08371861 -0.01166175  0.9213960  1.8381867  1.78720424
xp_design_med8_vs_Col_filtered = xp_design_med8_vs_Col[,2:4]

anno_col_info = xp_design_med8_vs_Col_filtered %>% column_to_rownames("group")
anno_info_colors = list(
  Condition = c(CT = "#d8b365"),
  Genotype = c(Col = "lightgrey", 
                med8 = "black"))
pheatmap(counts_scaled_filtered, 
         cluster_rows = TRUE,                       
         cluster_cols = TRUE, 
         show_rownames = FALSE, 
         show_colnames = TRUE,
         annotation_col = anno_col_info,
         annotation_colors = anno_info_colors,
         main = "Clustering DEGs med8_CT vs Col_CT")

Análisis de clustering de los DEGs

Cargamos una función que devuelve una lista de todos los DEGs, con sus valores de conteo normalizados. Esta función requiere: dataframe de raw_counts, dataframe con diseño experimental, especificar umbral de log2foldchange (f) y p-value máximo (p)

normalize<-function(countdata, xp_design, f, p){
  col.names = colnames(xp_design)                                                
  columns.to.discard = c("sample","group")                                    
  colsForConditions = col.names[! col.names %in% columns.to.discard]              
  # una condición
  if (length(colsForConditions) == 1){
    condition <- factor(xp_design[,colsForConditions])
    # dos condiciones
  } else if (length(colsForConditions) == 2) {
    # tres condiciones --> tercera columna que es el producto de las dos
    xp_design$conditions = paste(xp_design[,colsForConditions[1]],xp_design[,colsForConditions[2]],sep = ".")
    condition <- factor(x = xp_design[,"conditions"],levels = unique(xp_design[,"conditions"]))
  }
  # Analysis with DESeq2
  coldata <- data.frame(row.names=colnames(countdata), condition)
  dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)   
  dds <- DESeq(dds)
  # crea dataframe con los conteos normalizados
  resdata <- as.data.frame(counts(dds, normalized=TRUE))
  # crea valores de expresión diferencial 
  total_selec <- list()
  x <- 1
  for(i in levels(condition)){
    x <- x + 1
    if (x <= length(levels(condition))){
      for(j in x:length(levels(condition))){
        res <- results(dds, contrast=c("condition", i,  levels(condition)[j]))                      # i = first in pair, levels(condition)[j] is the second in pair.
        d <- paste(i, levels(condition)[j], sep="&")                                                # paste the two conditions in one character, to be used as the pair name
        res$genenames <- rownames(res)
        resul <- as.data.frame(res)
        significantDE <- resul %>% filter(padj<p & (log2FoldChange>f | log2FoldChange<(-1*f)) )
        selec <- as.list(significantDE$genenames)
        total_selec <- append(total_selec, selec)
      }
    }
  }
  total_selec <- c(unique(total_selec))
  total_selec <- t(as.data.frame(total_selec))
  selection <- resdata[total_selec[,1],]
  
  return(selection)
}

Filtramos por el factor de interés.

xp_design_CT_vs_MV = xp_design %>% 
  filter(Genotype == "Col")
raw_counts_filtered = raw_counts[, colnames(raw_counts) %in% xp_design_CT_vs_MV$group]

Ejecutamos la función, con los mismos valores de log2foldchange y p-value que antes.

DEgenes_ColMV_ColCT <- normalize(raw_counts_filtered,xp_design_CT_vs_MV,1,0.01) 
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
dim(DEgenes_ColMV_ColCT)
## [1] 1263    6
dim(DEG_ColMV_vs_ColCT)
## [1] 1263    7

El resultado tras filtrar es el mismo número de DEGs, pero la matriz de valores es diferente, aquí están con los valores normalizados y con el nombre de las muestras.

Se puede repetir este proceso con otras comparaciones (ej. med8_CT vs Col_CT).

xp_design_med8_vs_Col = xp_design %>% 
  filter(Condition == "CT")
raw_counts_filtered = raw_counts[, colnames(raw_counts) %in% xp_design_med8_vs_Col$group]
DEgenes_med8_Col <- normalize(raw_counts_filtered,xp_design_med8_vs_Col,1,0.01) 
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

Para realizar el clustering de DEGs, es necesario escalar los datos, para evitar que genes altamente expresados se agrupen juntos.

DEgenes_ColMV_ColCT <- as.matrix(DEgenes_ColMV_ColCT)
scaledata <- t(scale(t(DEgenes_ColMV_ColCT)))

Visualización de los DEGs

Construcción de los dendrogramas: Requiere calcular la distancia entre las diferentes muestras y genes.

Árbol jerárquico de las muestras.

hc <- hclust(as.dist(1-cor(scaledata, method="spearman")), method="complete") # Clusters columns by Spearman correlation.

sampleTree = as.dendrogram(hc, method="average")
plot(sampleTree,
     main = "Clustering de muestras",
     ylab = "Height")

Árbol jerárquico de los genes.

hr <- hclust(as.dist(1-cor(t(scaledata), method="pearson")), method="complete") # Cluster rows by Pearson correlation.

geneTree = as.dendrogram(hr, method="average")
plot(geneTree,
     leaflab = "none",             
     main = "Clustering de genes",
     ylab = "Height")

Clustering jerárquico de los datos.

library(dendextend)

Árbol con número de clusters de interés (ej. k=5).

hclustk5 = cutree(hr, k=5)
plot(geneTree,
     leaflab = "none",
     main = "Clustering de genes",
     ylab = "Height")
colored_bars(hclustk5, geneTree, sort_by_labels_order = T, y_shift=-0.1, rowLabels = c("k=5"),cex.rowLabels=0.7)

clustColBar <- rainbow(length(unique(hclustk5)), start=0.1, end=0.9)
clustColBar <- clustColBar[as.vector(hclustk5)]

heatmap.2(DEgenes_ColMV_ColCT,
          Rowv=as.dendrogram(hr), 
          Colv=as.dendrogram(hc),
          col=redgreen(100),
          scale="row",
          margins = c(7, 7),
          cexCol = 0.7,
          labRow = F,
          main = "Clustering DEGs Col_MV vs Col_CT",
          trace = "none",
          RowSideColors=clustColBar,
          key = FALSE)

Visualizar los ‘cores’ de los clusters.

clust.core = function(i, dat, clusters) {
  ind = (clusters == i)
  colMeans(dat[ind,])
}
clusters <- hclustk5
cores <- sapply(unique(clusters), clust.core, scaledata, clusters)

library(reshape2)
moltenCores <- melt(cores) 
colnames(moltenCores) <- c('sample','cluster','value')
ggplot(moltenCores, aes(x=sample,y=value, group=cluster, colour=as.factor(cluster))) + 
  geom_point() + 
  geom_line() +
  xlab("Sample") +
  ylab("Expression") +
  labs(title= "Cluster Expression of the samples",color = "Cluster") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Aislar un cluster determinado.

clust5 <- t(scaledata[hclustk5==5,])  # cluster 5
clust1 <- t(scaledata[hclustk5==1,])  # cluster 1

Cluster 5

clust5Molten <- melt(clust5, id.vars = "Sample")
colnames(clust5Molten) <- c('sample','gene','value')
core <- moltenCores[moltenCores$cluster==5,]

ggplot(clust5Molten, aes(x=sample,y=value)) + 
  geom_line(color="grey", aes(color="grey", group=gene)) +
  geom_line(data=core, aes(sample,value, group=cluster), color="blue",inherit.aes=FALSE) +
  geom_point(data=core, aes(sample,value, group=cluster), color="blue",inherit.aes=FALSE) +
  xlab("Muestras") +
  ylab("Expresion") +
  labs(title= paste0("Cluster 5 compuesto de ", ncol(clust5), " genes"),color = "Score") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Cluster 1

clust1Molten <- melt(clust1, id.vars = "Sample")
colnames(clust1Molten) <- c('sample','gene','value')
core <- moltenCores[moltenCores$cluster==1,]

ggplot(clust1Molten, aes(x=sample,y=value)) + 
  geom_line(color="grey", aes(color="grey", group=gene)) +
  geom_line(data=core, aes(sample,value, group=cluster), color="blue",inherit.aes=FALSE) +
  geom_point(data=core, aes(sample,value, group=cluster), color="blue",inherit.aes=FALSE) +
  xlab("Muestras") +
  ylab("Expresion") +
  labs(title= paste0("Cluster 1 compuesto de ", ncol(clust1), " genes"),color = "Score") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Análisis de enriquecimiento funcional con ClusterProfiler

Cogemos la lista de DEGs [Col0_MV vs Col0_CT] filtrados.

dim(DEG_ColMV_vs_ColCT) 
## [1] 1263    7

Anotación de DEGs con Ensembl y biomartr.

biomartr::organismBM(organism = "Arabidopsis thaliana")
## Starting retrieval of all available BioMart datasets for Arabidopsis thaliana ...
## # A tibble: 2 × 5
##   organism_name description                                mart  dataset version
##   <chr>         <chr>                                      <chr> <chr>   <chr>  
## 1 athaliana     Arabidopsis thaliana genes (TAIR10)        plan… athali… TAIR10 
## 2 athaliana     Arabidopsis thaliana Short Variants (SNPs… plan… athali… TAIR10
arabido_attributes = 
  biomartr::organismAttributes("Arabidopsis thaliana") %>% 
  filter(dataset == "athaliana_eg_gene")
## Starting retrieval of all available BioMart datasets for Arabidopsis thaliana ...
## 
## Starting retrieval of all available BioMart attributes for Arabidopsis thaliana ...
arabido_attributes
## # A tibble: 3,408 × 4
##    name                  description              dataset           mart       
##    <chr>                 <chr>                    <chr>             <chr>      
##  1 ensembl_gene_id       Gene stable ID           athaliana_eg_gene plants_mart
##  2 ensembl_transcript_id Transcript stable ID     athaliana_eg_gene plants_mart
##  3 ensembl_peptide_id    Protein stable ID        athaliana_eg_gene plants_mart
##  4 ensembl_exon_id       Exon stable ID           athaliana_eg_gene plants_mart
##  5 description           Gene description         athaliana_eg_gene plants_mart
##  6 chromosome_name       Chromosome/scaffold name athaliana_eg_gene plants_mart
##  7 start_position        Gene start (bp)          athaliana_eg_gene plants_mart
##  8 end_position          Gene end (bp)            athaliana_eg_gene plants_mart
##  9 strand                Strand                   athaliana_eg_gene plants_mart
## 10 band                  Karyotype band           athaliana_eg_gene plants_mart
## # ℹ 3,398 more rows
attributes_to_retrieve = c("tair_symbol", "uniprotswissprot","entrezgene_id")

DEG_ColMV_vs_ColCT_annotated <- biomartr::biomart(genes = DEG_ColMV_vs_ColCT$genes,         
                                mart       = "plants_mart",                     
                                dataset    = "athaliana_eg_gene",               
                                attributes = attributes_to_retrieve,            
                                filters =   "ensembl_gene_id" )
## Starting BioMart query ...
## 
## 
## Please cite: Drost HG, Paszkowski J. Biomartr: genomic data retrieval with R. Bioinformatics (2017) 33(8): 1216-1217. doi:10.1093/bioinformatics/btw821.
head(DEG_ColMV_vs_ColCT_annotated)
##   ensembl_gene_id tair_symbol uniprotswissprot entrezgene_id
## 1       AT1G01060         LHY                         839341
## 2       AT1G01060         LHY           Q6R0H1        839341
## 3       AT1G01190     CYP78A8                         839233
## 4       AT1G01520                                     839337
## 5       AT1G01520                       Q6R0H0        839337
## 6       AT1G01570                                     837183

Obtenemos los DEGs anotados con símbolos de bases de datos.

Análisis ORA-GO con ClusterProfiler

Construcción del universo de genes (selecciona directamente los nombres de las filas).

all_genes <- rownames(raw_counts)

Con los correspondientes símbolos TAIR/Ensembl, EntrezGeneids.

attributes_to_retrieve = c("tair_symbol", "uniprotswissprot","entrezgene_id")

Query Ensembl API

all_genes_annotated <- biomartr::biomart(genes = all_genes,
                                                     mart       = "plants_mart",                 
                                                     dataset    = "athaliana_eg_gene",           
                                                     attributes = attributes_to_retrieve,        
                                                     filters =  "ensembl_gene_id" ) 
## Starting BioMart query ...
## 
## 
## Please cite: Drost HG, Paszkowski J. Biomartr: genomic data retrieval with R. Bioinformatics (2017) 33(8): 1216-1217. doi:10.1093/bioinformatics/btw821.

Genes convertidos a caracteres, para la compatibilidad con enrichGO.

all_genes_annotated$entrezgene_id = as.character(
  all_genes_annotated$entrezgene_id) 

head(all_genes_annotated)
##   ensembl_gene_id tair_symbol uniprotswissprot entrezgene_id
## 1       AT1G01010     ANAC001           Q0WV96        839580
## 2       AT1G01020        ARV1                         839569
## 3       AT1G01020        ARV1           Q5MK24        839569
## 4       AT1G01030        NGA3           Q9MAN1        839321
## 5       AT1G01030        NGA3                         839321
## 6       AT1G01040        DCL1           Q9SP32        839574
head(DEG_ColMV_vs_ColCT_annotated)
##   ensembl_gene_id tair_symbol uniprotswissprot entrezgene_id
## 1       AT1G01060         LHY                         839341
## 2       AT1G01060         LHY           Q6R0H1        839341
## 3       AT1G01190     CYP78A8                         839233
## 4       AT1G01520                                     839337
## 5       AT1G01520                       Q6R0H0        839337
## 6       AT1G01570                                     837183

ORA para la clase GO biological process.

ora_analysis_bp <- enrichGO(gene = DEG_ColMV_vs_ColCT_annotated$entrezgene_id, 
                            universe = all_genes_annotated$entrezgene_id, 
                            OrgDb = org.At.tair.db, 
                            keyType = "ENTREZID",
                            ont = "BP",              
                            pAdjustMethod = "BH",
                            qvalueCutoff = 0.05,
                            readable = TRUE, 
                            pool = FALSE)

ORA para la clase GO cellular component.

ora_analysis_cc <- enrichGO(gene = DEG_ColMV_vs_ColCT_annotated$entrezgene_id, 
                            universe = all_genes_annotated$entrezgene_id, 
                            OrgDb = org.At.tair.db, 
                            keyType = "ENTREZID",
                            ont = "CC",              
                            pAdjustMethod = "BH",
                            qvalueCutoff = 0.05,
                            readable = TRUE, 
                            pool = FALSE)

ORA para la clase GO molecular function.

ora_analysis_mf <- enrichGO(gene = DEG_ColMV_vs_ColCT_annotated$entrezgene_id, 
                            universe = all_genes_annotated$entrezgene_id, 
                            OrgDb = org.At.tair.db,  
                            keyType = "ENTREZID",
                            ont = "MF",              
                            pAdjustMethod = "BH",
                            qvalueCutoff = 0.05,
                            readable = TRUE, 
                            pool = FALSE)

Simplificamos la clasificacion GO.

ora_analysis_bp_simplified <- clusterProfiler::simplify(ora_analysis_bp) 
ora_analysis_cc_simplified <- clusterProfiler::simplify(ora_analysis_cc) 
ora_analysis_mf_simplified <- clusterProfiler::simplify(ora_analysis_mf)

View(ora_analysis_bp_simplified@result)

Obtengo como resultado una lista de los procesos biológicos enriquecidos en los DEGs identificados

Extraer tabla de resultados:

write_delim(x = as.data.frame(ora_analysis_bp_simplified@result),
            "go_results.tsv", 
            delim = "\t")

Gráficos a partir del Gene Ontology ORA analysis

dotplot(ora_analysis_bp_simplified)

dotplot(ora_analysis_mf_simplified)

Mapa de enriquecimiento que conecta términos GO con puentes entre conjuntos de genes.

ora_analysis_bp <- pairwise_termsim(ora_analysis_bp, method = "JC")
emapplot(ora_analysis_bp, color = "qvalue")