library(tidyverse)
library(matrixStats)
library(patchwork)
library(RColorBrewer)
library(gridExtra)
library(tximport)
library(DESeq2)
library(apeglm)
library(sva)
library(umap)
library(Rtsne)
library(limma)
library(ConsensusClusterPlus)
library(pheatmap)
library(ggrepel)
library(data.table)
library(survival)
library(survminer)
library(clusterProfiler)
library(msigdbr)
library(fgsea)
library(ggridges)
library(enrichplot)
library(BiocParallel)
# Unsupervised Analysis Functions-----------------------------------------------
plot_pca <- function(matriz, variables, objeto_vsd, continua = FALSE) {
pca_res <- prcomp(t(matriz), scale. = FALSE)
df_plot <- as.data.frame(pca_res$x)
metadatos <- as.data.frame(colData(objeto_vsd)[colnames(matriz), ])
df_plot <- cbind(df_plot, metadatos)
percentVar <- round(100 * (pca_res$sdev^2 / sum(pca_res$sdev^2)))
p <- ggplot(df_plot, aes(PC1, PC2, color = .data[[variables]])) +
geom_point(size = 3, alpha = 0.8) +
labs(x = paste0("PC1: ", percentVar[1], "%"),
y = paste0("PC2: ", percentVar[2], "%"),
color = variables) +
theme_bw() +
theme(
axis.title = element_text(face = "bold", size = 12),
axis.text = element_text(size = 12, color = "black"),
legend.position = "top",
legend.title = element_text(face = "bold", size = 12),
legend.text = element_text(size = 12)
)
if (continua) {
p +
scale_color_viridis_c(name = variables) +
theme(legend.position = "right")
} else {
p + scale_color_brewer(palette = "Set1")
}
}
consensus_clustering <- function(matriz, objeto_vsd, Tejido,k)
{
genes_mad <- rowMads(matriz)
mat_cc <- matriz[genes_mad > quantile(genes_mad, 0.75), ]
mat_cc_scaled <- t(scale(t(mat_cc)))
cat( nrow(mat_cc), "genes retendios \n")
results <- ConsensusClusterPlus(mat_cc_scaled, maxK = 7, reps = 1000,
pItem = 0.8, pFeature = 1, clusterAlg = "hc",
distance = "pearson", plot = "png",title = paste0("CC_", Tejido), seed = 123)
vars <- intersect(c("Subtipo", "proliferation_score"), colnames(colData(objeto_vsd)))
ann_col <- as.data.frame(colData(objeto_vsd)[, vars, drop = FALSE])
colnames(ann_col)[colnames(ann_col) == "proliferation_score"] <- "Proliferación"
cm <- results[[k]]$consensusMatrix
colnames(cm) <- rownames(cm) <- colnames(mat_cc)
ann_colors <- list(
Subtipo = c(LCMc = "red", LCMng = "blue", Indeterminado = "lightgray"),
Proliferación = c("white", "darkcyan")
)
p <- pheatmap(cm,
annotation_col = ann_col, annotation_colors = ann_colors,
show_rownames = FALSE,fontsize_row = 10, show_colnames = FALSE, annotation_names_col = FALSE,silent = TRUE)
return(list(plot = p, res= results))
}
# DEA functions------------------------------------------------------------------
plot_volcano <- function(resLFC_df, top_n , lfc_cut = 0.5, padj_cut = 0.1)
{
deg_clean <- resLFC_df %>%
mutate(
diff_exp = case_when(
log2FoldChange > lfc_cut & padj < padj_cut ~ "Sobreexpresados",
log2FoldChange < -lfc_cut & padj < padj_cut ~ "Infraexpresados",
.default = "NS"
)
)
top_labels <- bind_rows(
filter(deg_clean, diff_exp == "Sobreexpresados") %>% slice_min(padj, n = top_n),
filter(deg_clean, diff_exp == "Infraexpresados") %>% slice_min(padj, n = top_n)
)
p <- ggplot(deg_clean, aes (x = log2FoldChange,y = -log10(padj),color = diff_exp)) +
geom_point(alpha = 0.5, size = 1) +
geom_vline(xintercept = c(-lfc_cut, lfc_cut),linetype = "dashed", color = "grey60") +
geom_hline(yintercept = -log10(padj_cut), linetype = "dashed", color = "grey60") +
scale_color_manual(values = c("Sobreexpresados" = "red", "Infraexpresados" = "blue","NS" = "grey80")) +
labs(
x = expression(log[2] ~ "Fold Change"),
y = expression(-log[10] ~ "(adj. p-value)"),
color = ""
) +
geom_text(data=top_labels,
aes(label = gene_label),
size = 3, color = "black",
vjust = -0.5,
check_overlap = TRUE
) +
theme_classic() + theme(legend.position = "top",legend.text = element_text(size=10))
return(p)
}
# Heatmap
plot_combined_heatmap <- function(contrast_results,
vsd,
dds,
top_n,
lfc_threshold = 0.5,
padj_threshold = 0.1,
cont_filter = FALSE) {
stopifnot(
is.list(contrast_results),
length(contrast_results) >= 1,
top_n >= 1
)
# Match the number of contrasts
cont_filter <- rep_len(cont_filter, length(contrast_results))
# Strip Ensembl version suffix
strip_version <- function(ids) sub("\\..*", "", ids)
# Select the top up- and down-regulated genes of a single contrast
select_top_genes <- function(results_df, n_genes, apply_cont_filter) {
if (is.null(results_df) || nrow(results_df) == 0) {
return(character(0))
}
classified <- results_df %>%
dplyr::mutate(
regulation = dplyr::case_when(
log2FoldChange > lfc_threshold & padj < padj_threshold ~ "up",
log2FoldChange < -lfc_threshold & padj < padj_threshold ~ "down",
.default = "ns"
)
)
if (apply_cont_filter) {
classified <- dplyr::filter(classified, regulation != "ns")
}
top_for_direction <- function(direction) {
classified %>%
dplyr::filter(regulation == direction) %>%
dplyr::slice_min(padj, n = n_genes, with_ties = FALSE) %>%
dplyr::pull(Gene.stable.ID)
}
strip_version(unique(c(top_for_direction("up"), top_for_direction("down"))))
}
# Gene selection across all contrasts
selected_genes <- purrr::imap(contrast_results, function(df, i) {
select_top_genes(df, top_n, cont_filter[i])
}) %>%
unlist() %>%
unique()
# Expression matrix, with version-stripped row names for clean matching
expression_matrix <- SummarizedExperiment::assay(vsd)
rownames(expression_matrix) <- strip_version(rownames(expression_matrix))
selected_genes <- selected_genes[selected_genes %in% rownames(expression_matrix)]
message(sprintf("Genes mapped to the heatmap: %d", length(selected_genes)))
if (length(selected_genes) == 0) {
stop("No matching genes found")
}
heatmap_matrix <- expression_matrix[selected_genes, ]
# Map Ensembl IDs to readable gene labels, pooling all contrasts
label_lookup <- contrast_results %>%
purrr::map(~ stats::setNames(.x$gene_label, strip_version(.x$Gene.stable.ID))) %>%
unlist()
label_lookup <- label_lookup[!duplicated(names(label_lookup))]
gene_labels <- label_lookup[selected_genes]
missing_labels <- is.na(gene_labels) | gene_labels == ""
gene_labels[missing_labels] <- selected_genes[missing_labels]
rownames(heatmap_matrix) <- gene_labels
column_annotation <- as.data.frame(
SummarizedExperiment::colData(dds)[, c("Cluster", "Subtipo"), drop = FALSE]
)
# Order samples by cluster instead of clustering the columns
cluster_order <- order(column_annotation$Cluster)
heatmap_matrix <- heatmap_matrix[, cluster_order]
column_annotation <- column_annotation[cluster_order, , drop = FALSE]
pheatmap::pheatmap(
heatmap_matrix,
cluster_rows = TRUE,
cluster_cols = FALSE,
show_rownames = TRUE,
show_colnames = FALSE,
fontsize_row = 8,
fontsize_col = 8,
annotation_col = column_annotation,
scale = "row",
color = colorRampPalette(c("blue", "white", "red"))(100),
border_color = NA
)
}
# GSEA functions---------------------------------------------------------------
run_gsea <- function(df_res, gene_sets_db) {
# Pre-ranked
gene_list <- df_res %>%
filter(!is.na(pvalue)) %>%
mutate(rank = -log10(pvalue) * sign(log2FoldChange)) %>%
arrange(desc(abs(rank))) %>%
distinct(Gene.stable.ID, .keep_all = TRUE) %>%
arrange(desc(rank)) %>%
select(Gene.stable.ID, rank) %>%
deframe()
gse_res <- GSEA(
geneList = gene_list,
TERM2GENE = gene_sets_db,
minGSSize = 10,
maxGSSize = 250,
pvalueCutoff = 0.05,
BPPARAM = SerialParam()
)
return(gse_res)
}
plot_gsea_dotplot <- function(gse_obj, show_n = 5, plot_title = "") {
if (nrow(as.data.frame(gse_obj)) == 0) {
return(NULL)
}
p <- dotplot(gse_obj, showCategory = show_n, split = ".sign", font.size = 10, title = plot_title) +
facet_grid(. ~ .sign) +
theme(plot.title = element_text(face = "bold", hjust = 0.5)) +
scale_y_discrete(labels = function(x) {
x %>%
stringr::str_replace_all("^(REACTOME_|WP_|HALLMARK_)", "") %>%
stringr::str_replace_all("_", " ") %>%
stringr::str_wrap(width = 44)
})
return(p)
}
################################################################################
# Import metadata
################################################################################
load("info.Paula.v2.RData")
sample.info <- sample.info %>%
dplyr::rename(Subtipo=RNAseq_subtype ,Pureza=Purity,Tejido=Tissue, Centro_Secuenciacion=SequencingCenter ) %>%
mutate(Date.of.sample = iconv(Date.of.sample, from = "", to = "UTF-8", sub = "")) %>%
mutate(Centro_Secuenciacion = fct_other(Centro_Secuenciacion, keep = "IDIBAPS", other_level = "Otros")) %>%
mutate(Pureza = replace_na(Pureza, "Desconocida")) %>%
mutate(across(c(-RIN, -Date.of.sample), as.factor)) %>%
drop_na(Tejido,Subtipo)
summary(sample.info)
## Sample Mcode Date.of.sample Centro_Secuenciacion
## M001 : 1 M389 : 6 Length:191 IDIBAPS:176
## M001-T2: 1 M401 : 4 Class :character Otros : 15
## M003 : 1 M032 : 3 Mode :character
## M004 : 1 M432 : 3
## M006 : 1 M001 : 2
## M007 : 1 M014 : 2
## (Other):185 (Other):171
## SequencingCenterProjecteID Sequencer RNAseqDoneBy
## Projecte_22_0496_NGS_EC:42 CNAG_Illumina: 12 BG: 18
## Projecte_25_0162 :42 NextSeq2000 :176 BP: 3
## Projecte_22_0363_NGS_EC:24 Other : 3 CL:158
## Projecte_v376 :15 FN: 12
## CLL_55 :12
## (Other) :53
## NA's : 3
## RIN Material Pureza Tejido Subytpe
## Min. : 2.900 Cryotube: 21 Desconocida: 19 colon : 1 cMCL :84
## 1st Qu.: 6.150 FFPE? : 1 sorted : 68 LN : 62 nnMCL :32
## Median : 7.600 Fresh :123 unsorted :104 Normal: 0 unclass: 4
## Mean : 7.335 OCT : 46 PB :127 NA's :71
## 3rd Qu.: 8.900 Spleen: 1
## Max. :10.000
##
## Subtipo
## cMCL :130
## nnMCL : 53
## Undetermined: 8
##
##
##
##
sample.info$RIN_cat <- cut(
sample.info$RIN,
breaks = c(1, 5, 8, 10),
labels = c("Degradado", "Parcialmente degradado", "Integro")
)
# Remove technical replicates and longitudinal samples + RNA highly degraded samples | RIN_cat == "RNA highly degraded")
samples_excluded <- sample.info %>%
filter(str_detect(Sample, "[_-]")| RIN_cat == "Degradado")%>%
select(Sample,RIN_cat) %>%
distinct()
# Remove technical replicates and longitudinal samples
df <- sample.info %>%
filter(!Sample %in% samples_excluded$Sample) %>%
droplevels()
df$Subtipo<-factor(df$Subtipo,levels = c("cMCL","nnMCL","Undetermined"),labels=c("LCMc","LCMng","Indeterminado"))
df$Material<-factor(df$Material,levels = c("Cryotube","FFPE?","Fresh","OCT"),labels=c("Criotubo","FFPE","Fresco","OCT"))
df$Tejido<-factor(df$Tejido,levels = c("LN","PB","Spleen"),labels=c("TL","SP","Bazo"))
df$Pureza<-factor(df$Pureza,levels = c("Desconocida","sorted","unsorted"),labels=c("Desconocida","Purificada","No purificada"))
summary(df)
## Sample Mcode Date.of.sample Centro_Secuenciacion
## M001 : 1 M001 : 1 Length:135 IDIBAPS:124
## M003 : 1 M003 : 1 Class :character Otros : 11
## M004 : 1 M004 : 1 Mode :character
## M006 : 1 M006 : 1
## M007 : 1 M007 : 1
## M008 : 1 M008 : 1
## (Other):129 (Other):129
## SequencingCenterProjecteID Sequencer RNAseqDoneBy
## Projecte_25_0162 :34 CNAG_Illumina: 8 BG: 16
## Projecte_22_0363_NGS_EC:23 NextSeq2000 :124 BP: 3
## Projecte_22_0496_NGS_EC:20 Other : 3 CL:108
## Projecte_v346 :10 FN: 8
## Projecte_v376 :10
## (Other) :35
## NA's : 3
## RIN Material Pureza Tejido Subytpe
## Min. :5.10 Criotubo:19 Desconocida :12 TL :41 cMCL :55
## 1st Qu.:6.85 FFPE : 1 Purificada :41 SP :93 nnMCL :19
## Median :8.00 Fresco :83 No purificada:82 Bazo: 1 unclass: 3
## Mean :7.75 OCT :32 NA's :58
## 3rd Qu.:8.90
## Max. :9.70
##
## Subtipo RIN_cat
## LCMc :94 Parcialmente degradado:69
## LCMng :36 Integro :66
## Indeterminado: 5
##
##
##
##
################################################################################
# Import counts and gene annotation
################################################################################
# Build paths to kallisto abundance files
files <- list.files("kallistos_tsv", full.names = TRUE)
names(files) <- gsub("_abundance\\.tsv","", basename(files))
# Filter active samples
files <- files[rownames(df)]
# Load Ensembl gene annotation
gene_annotation <- read.delim("Homo_sapiens.GRCh38_genes_transcripts_release-105.txt")
gene_names <- gene_annotation %>%
select(Gene.stable.ID, HGNC.symbol) %>%
distinct() %>%
mutate(gene_label = ifelse(HGNC.symbol == "" | is.na(HGNC.symbol),
Gene.stable.ID, HGNC.symbol))
tx2gene <- gene_annotation %>%
select(Transcript.stable.ID.version, Gene.stable.ID) %>%
distinct()
# Import transcript counts and aggregate to gene level
txi.kallisto <- tximport(files, type = "kallisto", tx2gene = tx2gene)
################################################################################
# Pre-processing of counts
################################################################################
# Filter metadata to match available counts
all(rownames(df) == colnames(txi.kallisto$counts))
## [1] TRUE
# Build general DESeqDataSet
dds <- DESeqDataSetFromTximport(txi.kallisto, colData = df, design = ~ 1)
dds
## class: DESeqDataSet
## dim: 66892 135
## metadata(1): version
## assays(2): counts avgTxLength
## rownames(66892): ENSG00000000003 ENSG00000000005 ... ENSG00000289643
## ENSG00000289644
## rowData names(0):
## colnames(135): M001 M003 ... M545 M546
## colData names(14): Sample Mcode ... Subtipo RIN_cat
# Pre-filtering
# keep only rows that have a count of at least 10 for a minimal number of samples
smallestGroupSize <- min(table(df$Subtipo))
keep <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep,]
# Variance stabilizing transformation for unsupervised analysis
vsd <- vst(dds, blind = TRUE)
mat_vsd <-assay(vsd)
# Split dds---------------------------------------------------------------------
dds_SP <- dds[, dds$Tejido == "SP"]
smallestGroupSize_SP <- min(table(dds_SP$Subtipo))
keep_SP <- rowSums(counts(dds_SP) >= 10) >= smallestGroupSize_SP
dds_SP <- dds_SP[keep_SP,]
dds_SP <- estimateSizeFactors(dds_SP)
vsd_SP <- vst(dds_SP, blind =TRUE)
mat_vsd_SP <- assay(vsd_SP)
dds_TL <- dds[, dds$Tejido %in% c("TL","Bazo")]
smallestGroupSize_TL <- min(table(dds_TL$Subtipo))
keep_TL <- rowSums(counts(dds_TL) >= 10) >= smallestGroupSize_TL
dds_TL <- dds_TL [keep_TL,]
dds_TL <- estimateSizeFactors(dds_TL)
vsd_TL<- vst(dds_TL, blind = TRUE)
mat_vsd_TL <- assay(vsd_TL )
################################################################################
# Proliferation score calculate
################################################################################
str(proliferation.signature)
## Named num [1:17] -1 -1 -1 -1 1 1 1 1 1 1 ...
## - attr(*, "names")= chr [1:17] "GLIPR1" "ZDHHC21" "FMNL3" "SPG3A" ...
# gene weights
prolif.weights <- c(GLIPR1 = -29.91, ZDHHC21 = -23.47, FMNL3 = -21.46, ATL1 = -19.64,
ZWINT = 5.41, FAM83D = 5.92, CCNB2 = 6.01, E2F2 = 6.02, H2AX = 6.08,
KIF2C = 6.19, CDC20 = 6.35, CDKN3 = 6.4, NCAPG = 6.44, TOP2A = 6.46,
ESPL1 = 6.5, FOXM1 = 6.55, MKI67 = 6.65)
# Ensembl IDs mapping for signature
sig_ensembl <- gene_names %>%
filter(gene_label %in% names(prolif.weights)) %>%
select(Gene.stable.ID, gene_label)%>%
distinct(gene_label, .keep_all = TRUE) # E2EF replicate:ENSG00000282899 ENSG00000007968
v_weights <- prolif.weights[sig_ensembl$gene_label]
names(v_weights) <- sig_ensembl$Gene.stable.ID
# Multiply each gene by its corresponding weight
prolif_score_SP <- sweep(x = mat_vsd_SP[names(v_weights),],
MARGIN = 1,
STATS = v_weights,
FUN = "*")
# Add proliferation score to sample metadata
vsd_SP$Proliferacion <- colSums(prolif_score_SP)
dds_SP$Proliferacion <- vsd_SP$Proliferacion
prolif_score_TL <- sweep(x = mat_vsd_TL [names(v_weights),],
MARGIN = 1,
STATS = v_weights,
FUN = "*")
# Add proliferation score to sample metadata SP
vsd_TL$Proliferacion<- colSums(prolif_score_TL)
dds_TL$Proliferacion<- vsd_TL$Proliferacion
data_mki67_SP <- plotCounts(dds_SP,
gene = "ENSG00000148773",
intgroup = "Subtipo",
returnData = TRUE)
a<-ggplot(data_mki67_SP,
aes(x = Subtipo, y = count, fill = Subtipo)) +
geom_boxplot(outlier.shape = NA, alpha = 0.4) +
geom_jitter(width = 0.2, size = 3, aes(color = Subtipo)) +
scale_y_log10() +
theme_minimal() +
labs(
title = "Correlación subtipo y gen MKI67 en sangre périferica",
x = "Subtype", y = "Normalized Counts (log10)")
data_mki67_TL <- plotCounts(dds_TL ,
gene = "ENSG00000148773",
intgroup = "Subtipo",
returnData = TRUE)
b<-ggplot(data_mki67_TL ,
aes(x = Subtipo, y = count, fill = Subtipo)) +
geom_boxplot(outlier.shape = NA, alpha = 0.4) +
geom_jitter(width = 0.2, size = 3, aes(color = Subtipo)) +
scale_y_log10() +
theme_minimal() +
labs(
title = "Correlación subtipo y gen MKI67 en tejido linfoide secundario",
x = "Subtype", y = "Normalized Counts (log10)")
a+b
################################################################################
# Ánalisis exploratorio inicial
################################################################################
vars <- c("Tejido","Subtipo","Centro_Secuenciacion","Material","Pureza","RIN_cat")
# PCA global
p_before <- lapply(vars, function(v) plot_pca(mat_vsd, v, vsd))
wrap_plots(p_before, ncol = 2) + plot_annotation(tag_levels = "A")
# PCA según Tejido
vars <- (vars [-1])
TL_before <- lapply(vars, function(v) plot_pca(mat_vsd_TL, v, vsd_TL))
pca_prolif <- plot_pca(mat_vsd_TL, "Proliferacion", vsd_TL, continua=TRUE)
TL_before[[6]] <- pca_prolif
wrap_plots(TL_before, ncol = 2) + plot_annotation(tag_levels = "A")
SP_before <- lapply(vars, function(v) plot_pca(mat_vsd_SP, v, vsd_SP))
pca_prolif <- plot_pca(mat_vsd_SP, "Proliferacion", vsd_SP, continua=TRUE)
SP_before[[6]] <- pca_prolif
wrap_plots(SP_before, ncol = 2) + plot_annotation(tag_levels = "A")
# UMAP SP ---------------------------------------------------------------------------
umap_result <- umap(t(mat_vsd_SP))
umap_df <- as.data.frame(umap_result$layout) %>%
mutate(Subtipo = vsd_SP$Subtipo,
Pureza = vsd_SP$Pureza,
Centro = vsd_SP$Centro_Secuenciacion)
a <- ggplot(umap_df, aes(x =V1, y = V2, color = Subtipo)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_brewer(palette = "Set1") +
labs(x = "UMAP1", y = "UMAP2") +
theme_classic()
# t-SNE SP ---------------------------------------------------------------------
tsne_result <- Rtsne(t(mat_vsd_SP), perplexity = 10)
tsne_df <- as.data.frame(tsne_result$Y) %>%
mutate(Subtipo = vsd_SP$Subtipo,
Pureza = vsd_SP$Pureza)
b <- ggplot(tsne_df, aes(x = V1, y = V2, color = Subtipo)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_brewer(palette = "Set1") +
labs(x = "tSNE1", y = "tSNE2") +
theme_classic()
# MDS SP -------------------------------------------------------------------------
mds <- plotMDS(mat_vsd_SP, gene.selection = "pairwise", plot = FALSE)
mds_df <- data.frame(X = mds$x, Y = mds$y) %>%
mutate(Subtipo = vsd_SP$Subtipo,
Pureza = vsd_SP$Pureza)
per <- mds$var.explained * 100
c <- ggplot(mds_df, aes(X, Y, color = Subtipo)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_brewer(palette = "Set1") +
labs(x = paste0("MDS1 (", round(per[1], 1), "%)"),
y = paste0("MDS2 (", round(per[2], 1), "%)")) +
theme_classic()
# Representación conjunta de los plots
(a + b + c) +
plot_layout(ncol = 3, guides = "collect") +
plot_annotation(tag_levels = "A") &
theme_bw() &
theme(legend.position = "bottom",
axis.title = element_text(face = "bold", size = 14),
axis.text = element_text(size = 12, color="black"),
legend.title = element_text(size=14, face = "bold"),
legend.text = element_text(size=12))
################################################################################
# Análisis de Variables Subrrogadas y correción de facores latentes y efecto de lote
################################################################################
# Definición del modelo para TL
mod_TL <- model.matrix(~ Subtipo + Proliferacion, data = colData(vsd_TL))
mod0_TL <- model.matrix(~ 1, data = colData(vsd_TL))
n_sv_TL <- num.sv(mat_vsd_TL, mod_TL, method = "leek", seed = 123)
svobj_TL <- svaseq(mat_vsd_TL, mod_TL, mod0_TL ,n.sv=n_sv_TL )
## Number of significant surrogate variables is: 1
## Iteration (out of 5 ):1 2 3 4 5
# Correlación SVs con los metadatos
vars2 <- setdiff(c(vars, "RIN", "Proliferacion"), "RIN_cat")
for(i in 1:ncol(svobj_TL$sv)) {
sv_vector <- svobj_TL$sv[,i]
cat(sprintf("\n---------------------- SV%d ------------------------\n", i))
cat(sprintf("%-20s | %-10s | %-6s\n", "Variable", "p-val", "rho"))
cat("----------------------------------------------------\n")
for(var_name in vars2) {
x <- colData(vsd_TL )[[var_name]]
tryCatch({
if (is.numeric(x)) {
res <- cor.test(sv_vector, x, method = "spearman")
p_val <- if(res$p.value < 0.001) "<0.001" else sprintf("%.3f", res$p.value)
cat(sprintf("%-20s | %-10s | %5.3f\n",
var_name, p_val, res$estimate))
} else {
res <- kruskal.test(sv_vector, as.factor(x))
p_val <- if(res$p.value < 0.001) "<0.001" else sprintf("%.3f", res$p.value)
cat(sprintf("%-20s | %-10s | %-6s\n",
var_name, p_val, ""))
}
}, error = function(e) NULL)
}
cat ("\n")
}
##
## ---------------------- SV1 ------------------------
## Variable | p-val | rho
## ----------------------------------------------------
## Subtipo | 0.049 |
## Centro_Secuenciacion | 0.516 |
## Material | 0.035 |
## Pureza | 0.288 |
## RIN | 0.047 | -0.309
## Proliferacion | 0.002 | -0.461
# Correción SV en TL
mat_uns_TL <- removeBatchEffect(
mat_vsd_TL,
covariates = svobj_TL$sv,
design = mod_TL
)
# Definición del modelo para SP
mod_SP <- model.matrix(~ Subtipo + Proliferacion + Centro_Secuenciacion + Pureza, data = colData(vsd_SP))
mod0_SP <- model.matrix(~ Centro_Secuenciacion + Pureza, data = colData(vsd_SP))
n_sv_SP <- num.sv(mat_vsd_SP, mod_SP, method = "leek", seed = 123)
svobj_SP <- svaseq(mat_vsd_SP, mod_SP, mod0_SP,n.sv=n_sv_SP)
## Number of significant surrogate variables is: 1
## Iteration (out of 5 ):1 2 3 4 5
# Correlación SVs con los metadatos
for(i in 1:ncol(svobj_SP$sv)) {
sv_vector <- svobj_SP$sv[, i]
cat(sprintf("\n---------------------- SV%d ------------------------\n", i))
cat(sprintf("%-20s | %-10s | %-6s\n", "Variable", "p-val", "rho"))
cat("----------------------------------------------------\n")
for(var_name in vars2) {
x <- colData(vsd_SP)[[var_name]]
tryCatch({
if (is.numeric(x)) {
res <- cor.test(sv_vector, x, method = "spearman")
p_val <- if(res$p.value < 0.001) "<0.001" else sprintf("%.3f", res$p.value)
cat(sprintf("%-20s | %-10s | %5.3f \n",
var_name, p_val, res$estimate))
} else {
res <- kruskal.test(sv_vector, as.factor(x))
p_val <- if(res$p.value < 0.001) "<0.001" else sprintf("%.3f", res$p.value)
cat(sprintf("%-20s | %-10s | %-6s \n",
var_name, p_val, ""))
}
}, error = function(e) NULL)
}
cat ("\n")
}
##
## ---------------------- SV1 ------------------------
## Variable | p-val | rho
## ----------------------------------------------------
## Subtipo | 0.774 |
## Centro_Secuenciacion | 0.003 |
## Material | 0.016 |
## Pureza | <0.001 |
## RIN | <0.001 | 0.383
## Proliferacion | 0.883 | 0.015
# Correción SV1 y efecto de lote en SP
mat_uns_SP <- removeBatchEffect(
mat_vsd_SP,
covariates = cbind(svobj_SP$sv, vsd_SP$Pureza,vsd_SP$Centro_Secuenciacion),
design = model.matrix(~ Subtipo + Proliferacion, data = colData(vsd_SP))
)
################################################################################
# Métodos no supervisados
################################################################################
# PCA plots---------------------------------------------------------------------
TL_after <- lapply(vars, function(v) plot_pca(mat_uns_TL, v, vsd_TL))
pca_prolif <- plot_pca(mat_uns_TL, "Proliferacion", vsd_TL, continua=TRUE)
TL_after[[6]] <- pca_prolif
wrap_plots(TL_after, ncol = 2) + plot_annotation(tag_levels = "A")
SP_after <- lapply(vars, function(v) plot_pca(mat_uns_SP, v, vsd_SP))
pca_prolif <- plot_pca(mat_uns_SP, "Proliferacion", vsd_SP, continua=TRUE)
SP_after[[6]] <- pca_prolif
wrap_plots(SP_after, ncol = 2) + plot_annotation(tag_levels = "A")
# UMAP SP ---------------------------------------------------------------------------
umap_result <- umap(t(mat_uns_SP))
umap_df <- as.data.frame(umap_result$layout) %>%
mutate(Subtipo = vsd_SP$Subtipo,
Pureza = vsd_SP$Pureza,
Centro = vsd_SP$Centro_Secuenciacion)
d <- ggplot(umap_df, aes(x =V1, y = V2, color = Subtipo)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_brewer(palette = "Set1") +
labs(x = "UMAP1", y = "UMAP2") +
theme_classic()
# t-SNE SP ---------------------------------------------------------------------
tsne_result <- Rtsne(t(mat_uns_SP), perplexity = 10)
tsne_df <- as.data.frame(tsne_result$Y) %>%
mutate(Subtipo = vsd_SP$Subtipo,
Pureza= vsd_SP$Pureza)
e <- ggplot(tsne_df, aes(x = V1, y = V2, color = Subtipo)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_brewer(palette = "Set1") +
labs(x = "tSNE1", y = "tSNE2") +
theme_classic()
# MDS SP -------------------------------------------------------------------------
mds <- plotMDS(mat_uns_SP, gene.selection = "pairwise", plot = FALSE)
mds_df <- data.frame(X = mds$x, Y = mds$y) %>%
mutate(Subtipo = vsd_SP$Subtipo,
Pureza= vsd_SP$Pureza)
per <- mds$var.explained * 100
f <- ggplot(mds_df, aes(X, Y, color = Subtipo)) +
geom_point(size = 3, alpha = 0.8) +
scale_color_brewer(palette = "Set1") +
labs(x = paste0("MDS1 (", round(per[1], 1), "%)"),
y = paste0("MDS2 (", round(per[2], 1), "%)")) +
theme_classic()
# Representación conjunta de los plots
(a + b + c + d + e + f ) +
plot_layout(ncol = 3, guides = "collect") +
plot_annotation(tag_levels = "A") &
theme_bw() &
theme(legend.position = "bottom",
axis.title = element_text(face = "bold", size = 14),
axis.text = element_text(size = 12, color="black"),
legend.title = element_text(size=14, face = "bold"),
legend.text = element_text(size=12))
# Hierarchical clustering -------------------------------------------------------
# Matriz que corrige SV1 y efectos de lote, preservando la señal biológica
p<-consensus_clustering (mat_uns_SP,vsd_SP,"SP_mat_uns25",k=3)
## 8613 genes retendios
grid.arrange(p$plot$gtable)
# PCA de los clústeres que preservan las variables biológicas
vsd_SP0 <- vsd_SP
assay(vsd_SP0) <- mat_uns_SP
clusteres0 <- p$res[[3]]$consensusClass
clusteres0 <- clusteres0[ rownames(colData(vsd_SP0))]
colData(vsd_SP0)$Cluster <- as.factor(clusteres0)
plotPCA(vsd_SP0, intgroup = "Cluster",ntop=8613)
t_cont0<-table(vsd_SP0$Cluster,vsd_SP0$Subtipo)
fisher.test(t_cont0)
##
## Fisher's Exact Test for Count Data
##
## data: t_cont0
## p-value = 1.296e-10
## alternative hypothesis: two.sided
vsd_SP$SV1 <- svobj_SP$sv
# Matriz corregida por efectos de lote, preservando las variables biológicas conocidas
mat_adjusted_SP0.1 <- removeBatchEffect(assay(vsd_SP),
covariates = cbind(vsd_SP$Centro_Secuenciacion, vsd_SP$Pureza),
design = model.matrix(~ vsd_SP$Subtipo + vsd_SP$Proliferacion)
)
# Identificación y observación de clusteres
p0<-consensus_clustering (mat_adjusted_SP0.1,vsd_SP,"SP.025",k=3)
## 8613 genes retendios
grid.arrange(p0$plot$gtable)
# Matriz corregida por efectos de lote y variables biológicas conocidos para identificar nuevos clústeres
mat_adjusted_SP1 <- removeBatchEffect(
assay(vsd_SP),
covariates = cbind(
vsd_SP$Proliferacion,vsd_SP$Subtipo,
vsd_SP$Centro_Secuenciacion,vsd_SP$Pureza
)
)
p1<-consensus_clustering (mat_adjusted_SP1,vsd_SP,"SP1.25",k=3)
## 8613 genes retendios
grid.arrange(p1$plot$gtable)
# Matriz corregida por centro de secuenciación,pureza, efectos biológicos conocidos y SVs identificadas para identificar nuevos clusteres
mat_adjusted_SP2 <- removeBatchEffect(
assay(vsd_SP),
covariates = cbind(
vsd_SP$Proliferacion,vsd_SP$Subtipo,
vsd_SP$Centro_Secuenciacion,vsd_SP$Pureza,
vsd_SP$SV1
)
)
p2<-consensus_clustering (mat_adjusted_SP2,vsd_SP,"SP2.25",k=3)
## 8613 genes retendios
grid.arrange(p2$plot$gtable)
# PCA de los clústeres ajustados por las variables biológicas
vsd_SP1 <- vsd_SP
assay(vsd_SP1) <- mat_adjusted_SP1
clusteres1 <- p1$res[[3]]$consensusClass
clusteres1 <- clusteres1[ rownames(colData(vsd_SP1)) ]
colData(vsd_SP1)$Cluster <- as.factor(clusteres1)
plotPCA(vsd_SP1, intgroup = "Cluster",ntop=8613)
t_cont1<-table(vsd_SP1$Cluster,vsd_SP1$Subtipo)
fisher.test(t_cont1)
##
## Fisher's Exact Test for Count Data
##
## data: t_cont1
## p-value = 0.8342
## alternative hypothesis: two.sided
################################################################################
# Análisis de supervivencia
################################################################################
surv_df <- read.delim("survival.patients.txt", header = TRUE, stringsAsFactors = FALSE)
# Clúesteres que preservan las variables biológicas
cluster_data0 <- data.frame(
id = rownames(colData(vsd_SP0)),
Cluster = colData(vsd_SP0)$Cluster,
Subtipo = colData(vsd_SP0)$Subtipo,
Proliferacion = colData(vsd_SP0)$Proliferacion
)
head(cluster_data0)
## id Cluster Subtipo Proliferacion
## M001 M001 1 LCMc -446.1330
## M003 M003 1 LCMng -259.5152
## M004 M004 2 LCMng -546.5540
## M006 M006 3 LCMc -517.2294
## M008 M008 3 LCMc -538.1799
## M009 M009 3 LCMng -540.2379
surv_data0 <- merge(surv_df, cluster_data0, by = "id")
head(surv_data0)
## id os_yrs event Cluster Subtipo Proliferacion
## 1 M001 0.6926762 1 1 LCMc -446.1330
## 2 M003 7.4168378 1 1 LCMng -259.5152
## 3 M004 4.2710472 1 2 LCMng -546.5540
## 4 M008 1.4839151 0 3 LCMc -538.1799
## 5 M009 6.3299110 1 3 LCMng -540.2379
## 6 M019 5.1033539 1 1 LCMc -474.7440
fit0 <- survfit(Surv(surv_data0$os_yrs, surv_data0$event) ~ Cluster, data = surv_data0)
summary(fit0)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## Cluster=1 38 38 38 18 7.218764 1.311128 5.103354 3.433265
## Cluster=2 16 16 16 6 7.737965 2.009521 8.043806 4.271047
## Cluster=3 23 23 23 8 9.909080 1.509100 8.156057 6.329911
## 0.95UCL
## Cluster=1 NA
## Cluster=2 NA
## Cluster=3 NA
ggsurvplot(fit0,
pval = TRUE,
risk.table = TRUE,
xlab = "Tiempo (Años)",
surv.median.line = "hv",
ylab = "Probabilidad de Supervivencia",
ggtheme = theme_minimal(),
palette = "Set4")
cox0 <- coxph(Surv(surv_data0$os_yrs, surv_data0$event) ~ Cluster+ Subtipo + Proliferacion, data = surv_data0)
summary(cox0)
## Call:
## coxph(formula = Surv(surv_data0$os_yrs, surv_data0$event) ~ Cluster +
## Subtipo + Proliferacion, data = surv_data0)
##
## n= 77, number of events= 32
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Cluster2 0.576676 1.780111 0.653983 0.882 0.3779
## Cluster3 -0.684824 0.504179 0.482738 -1.419 0.1560
## SubtipoLCMng -0.751088 0.471853 0.504321 -1.489 0.1364
## SubtipoIndeterminado -0.225571 0.798061 1.109374 -0.203 0.8389
## Proliferacion 0.004169 1.004178 0.001899 2.195 0.0281 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Cluster2 1.7801 0.5618 0.49406 6.414
## Cluster3 0.5042 1.9834 0.19574 1.299
## SubtipoLCMng 0.4719 2.1193 0.17560 1.268
## SubtipoIndeterminado 0.7981 1.2530 0.09073 7.020
## Proliferacion 1.0042 0.9958 1.00045 1.008
##
## Concordance= 0.69 (se = 0.047 )
## Likelihood ratio test= 9.65 on 5 df, p=0.09
## Wald test = 8.94 on 5 df, p=0.1
## Score (logrank) test = 9.47 on 5 df, p=0.09
# Clusteres ajustados por las variables biológicas
cluster_data1 <- data.frame(
id = rownames(colData(vsd_SP1)),
Cluster = colData(vsd_SP1)$Cluster,
Subtipo = colData(vsd_SP1)$Subtipo,
Proliferacion = colData(vsd_SP1)$Proliferacion
)
head(cluster_data1)
## id Cluster Subtipo Proliferacion
## M001 M001 1 LCMc -446.1330
## M003 M003 2 LCMng -259.5152
## M004 M004 3 LCMng -546.5540
## M006 M006 1 LCMc -517.2294
## M008 M008 1 LCMc -538.1799
## M009 M009 3 LCMng -540.2379
surv_data1 <- merge(surv_df, cluster_data1, by = "id")
head(surv_data1)
## id os_yrs event Cluster Subtipo Proliferacion
## 1 M001 0.6926762 1 1 LCMc -446.1330
## 2 M003 7.4168378 1 2 LCMng -259.5152
## 3 M004 4.2710472 1 3 LCMng -546.5540
## 4 M008 1.4839151 0 1 LCMc -538.1799
## 5 M009 6.3299110 1 3 LCMng -540.2379
## 6 M019 5.1033539 1 3 LCMc -474.7440
fit1 <- survfit(Surv(surv_data1$os_yrs, surv_data1$event) ~ Cluster, data = surv_data1)
summary(fit1)$table
## records n.max n.start events rmean se(rmean) median 0.95LCL
## Cluster=1 28 28 28 13 8.186282 1.423389 6.809035 3.433265
## Cluster=2 25 25 25 11 8.029012 1.623715 7.416838 3.474333
## Cluster=3 24 24 24 8 8.347471 1.688781 5.837098 4.271047
## 0.95UCL
## Cluster=1 NA
## Cluster=2 NA
## Cluster=3 NA
ggsurvplot(fit1,
pval = TRUE,
risk.table = TRUE,
xlab = "Tiempo (Años)",
surv.median.line = "hv",
ylab = "Probabilidad de Supervivencia",
ggtheme = theme_minimal(),
palette = "Set1")
cox1 <- coxph(Surv(surv_data1$os_yrs, surv_data1$event) ~ Cluster + Subtipo + Proliferacion, data = surv_data1)
summary(cox1)
## Call:
## coxph(formula = Surv(surv_data1$os_yrs, surv_data1$event) ~ Cluster +
## Subtipo + Proliferacion, data = surv_data1)
##
## n= 77, number of events= 32
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Cluster2 -0.113326 0.892860 0.421942 -0.269 0.7883
## Cluster3 -0.202695 0.816527 0.464248 -0.437 0.6624
## SubtipoLCMng -0.657962 0.517906 0.391013 -1.683 0.0924 .
## SubtipoIndeterminado -0.916847 0.399778 1.061494 -0.864 0.3877
## Proliferacion 0.002541 1.002544 0.001712 1.484 0.1378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Cluster2 0.8929 1.1200 0.39050 2.041
## Cluster3 0.8165 1.2247 0.32870 2.028
## SubtipoLCMng 0.5179 1.9309 0.24067 1.115
## SubtipoIndeterminado 0.3998 2.5014 0.04992 3.202
## Proliferacion 1.0025 0.9975 0.99919 1.006
##
## Concordance= 0.678 (se = 0.055 )
## Likelihood ratio test= 5.7 on 5 df, p=0.3
## Wald test = 5.64 on 5 df, p=0.3
## Score (logrank) test = 5.77 on 5 df, p=0.3
################################################################################
# DEA
################################################################################
# Clíesres que preservan las variables biológicas
clusteres0 <- clusteres0[ rownames(colData(dds_SP)) ]
colData(dds_SP)$Cluster <- as.factor(clusteres0)
dds_SP$SV1 <- svobj_SP$sv[, 1]
design(dds_SP) <- ~ Centro_Secuenciacion + Pureza + SV1 + Proliferacion + Subtipo + Cluster
dds_SP0<- DESeq(dds_SP)
resultsNames(dds_SP0)
## [1] "Intercept"
## [2] "Centro_Secuenciacion_Otros_vs_IDIBAPS"
## [3] "Pureza_Purificada_vs_Desconocida"
## [4] "Pureza_No.purificada_vs_Desconocida"
## [5] "SV1"
## [6] "Proliferacion"
## [7] "Subtipo_LCMng_vs_LCMc"
## [8] "Subtipo_Indeterminado_vs_LCMc"
## [9] "Cluster_2_vs_1"
## [10] "Cluster_3_vs_1"
# Clíesres ajustados por las variables biológicas
clusteres1 <- clusteres1[ rownames(colData(dds_SP))]
colData(dds_SP)$Cluster <- as.factor(clusteres1)
design(dds_SP) <- ~ Centro_Secuenciacion + Pureza + SV1 + Proliferacion + Subtipo + Cluster
dds_SP1 <- DESeq(dds_SP)
resultsNames(dds_SP1)
## [1] "Intercept"
## [2] "Centro_Secuenciacion_Otros_vs_IDIBAPS"
## [3] "Pureza_Purificada_vs_Desconocida"
## [4] "Pureza_No.purificada_vs_Desconocida"
## [5] "SV1"
## [6] "Proliferacion"
## [7] "Subtipo_LCMng_vs_LCMc"
## [8] "Subtipo_Indeterminado_vs_LCMc"
## [9] "Cluster_2_vs_1"
## [10] "Cluster_3_vs_1"
# Volcano plot de los clústeres que preservan las variables biológicas
# Clúster 1 VS Clúster 2
res_C2_vs_C1 <- results(dds_SP0, contrast = c("Cluster", "2", "1"), alpha = 0.1, lfcThreshold = 0.5)
summary(res_C2_vs_C1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 4, 0.012%
## LFC < -0.50 (down) : 9, 0.026%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
resLFC_SP0_2_1 <- lfcShrink(dds_SP0, contrast = c("Cluster", "2", "1"), res = res_C2_vs_C1, type = "ashr")
summary(resLFC_SP0_2_1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 4, 0.012%
## LFC < -0.50 (down) : 9, 0.026%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
df_SP0.1 <- as.data.frame(resLFC_SP0_2_1) %>%
rownames_to_column("Gene.stable.ID") %>%
filter(!is.na(pvalue) & !is.na(log2FoldChange)) %>%
left_join(gene_names, by = "Gene.stable.ID")
plot_volcano(df_SP0.1, top_n = 20)
# Clúster 2 VS Clúster 3
res_C2_vs_C3 <- results(dds_SP0, contrast = c("Cluster", "2", "3"), alpha = 0.1, lfcThreshold = 0.5)
summary(res_C2_vs_C3)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 214, 0.62%
## LFC < -0.50 (down) : 98, 0.28%
## outliers [1] : 0, 0%
## low counts [2] : 668, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
resLFC_SP0_2_3 <- lfcShrink(dds_SP0, contrast = c("Cluster", "2", "3"), res = res_C2_vs_C3, type = "ashr")
summary(resLFC_SP0_2_3)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 214, 0.62%
## LFC < -0.50 (down) : 98, 0.28%
## outliers [1] : 0, 0%
## low counts [2] : 668, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
df_SP0.2 <- as.data.frame(resLFC_SP0_2_3) %>%
rownames_to_column("Gene.stable.ID") %>%
filter(!is.na(pvalue) & !is.na(log2FoldChange)) %>%
left_join(gene_names, by = "Gene.stable.ID")
plot_volcano(df_SP0.2, top_n = 20)
# Clúster 3 VS Clúster 1
res_C3_vs_1 <- results(dds_SP0, contrast = c("Cluster", "3", "1"), alpha = 0.1, lfcThreshold = 0.5)
summary(res_C3_vs_1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 503, 1.5%
## LFC < -0.50 (down) : 589, 1.7%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
resLFC_SP0_3_1 <- lfcShrink(dds_SP0, contrast = c("Cluster", "3", "1"), res = res_C3_vs_1, type = "ashr")
summary(resLFC_SP0_3_1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 503, 1.5%
## LFC < -0.50 (down) : 589, 1.7%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
df_SP0.3 <- as.data.frame(resLFC_SP0_3_1) %>%
rownames_to_column("Gene.stable.ID") %>%
filter(!is.na(pvalue) & !is.na(log2FoldChange)) %>%
left_join(gene_names, by = "Gene.stable.ID")
plot_volcano(df_SP0.3, top_n = 20)
# Heatmap de clústeres que preservan las variables biológicas
vsd_SP0 <- vst(dds_SP0, blind = FALSE)
plot_combined_heatmap(list(df_SP0.1,df_SP0.2,df_SP0.3),
vsd = vsd_SP0,
dds = dds_SP0,
top_n = 20
)
# Volcano plot de clústeres ajustados por variables biológicas
# Clúster 2 VS Clúster 1
res_C2_vs_C1 <- results(dds_SP1, contrast = c("Cluster", "2", "1"), alpha = 0.1, lfcThreshold = 0.5)
summary(res_C2_vs_C1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 829, 2.4%
## LFC < -0.50 (down) : 459, 1.3%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
resLFC_SP1_2_1 <- lfcShrink(dds_SP1, contrast = c("Cluster", "2", "1"), res = res_C2_vs_C1, type = "ashr")
summary(resLFC_SP1_2_1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 829, 2.4%
## LFC < -0.50 (down) : 459, 1.3%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
df_SP1.1 <- as.data.frame(resLFC_SP1_2_1) %>%
rownames_to_column("Gene.stable.ID") %>%
filter(!is.na(pvalue) & !is.na(log2FoldChange)) %>%
left_join(gene_names, by = "Gene.stable.ID")
plot_volcano(df_SP1.1, top_n = 20)
# Clúster 2 VS Clúster 3
res_C2_vs_C3 <- results(dds_SP1, contrast = c("Cluster", "2", "3"), alpha = 0.1, lfcThreshold = 0.5)
summary(res_C2_vs_C3)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 74, 0.21%
## LFC < -0.50 (down) : 67, 0.19%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
resLFC_SP1_2_3 <- lfcShrink(dds_SP1, contrast = c("Cluster", "2", "3"), res = res_C2_vs_C3, type = "ashr")
summary(resLFC_SP1_2_3)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 74, 0.21%
## LFC < -0.50 (down) : 67, 0.19%
## outliers [1] : 0, 0%
## low counts [2] : 0, 0%
## (mean count < 0)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
df_SP1.2 <- as.data.frame(resLFC_SP1_2_3) %>%
rownames_to_column("Gene.stable.ID") %>%
filter(!is.na(pvalue) & !is.na(log2FoldChange)) %>%
left_join(gene_names, by = "Gene.stable.ID")
plot_volcano(df_SP1.2, top_n = 20)
# Clúster 3 VS Clúster 1
res_C3_vs_1 <- results(dds_SP1, contrast = c("Cluster", "3", "1"), alpha = 0.1, lfcThreshold = 0.5)
summary(res_C3_vs_1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 5, 0.015%
## LFC < -0.50 (down) : 32, 0.093%
## outliers [1] : 0, 0%
## low counts [2] : 668, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
resLFC_SP1_3_1 <- lfcShrink(dds_SP1, contrast = c("Cluster", "3", "1"), res = res_C3_vs_1, type = "ashr")
summary(resLFC_SP1_3_1)
##
## out of 34452 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0.50 (up) : 5, 0.015%
## LFC < -0.50 (down) : 32, 0.093%
## outliers [1] : 0, 0%
## low counts [2] : 668, 1.9%
## (mean count < 2)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
df_SP1.3 <- as.data.frame(resLFC_SP1_3_1) %>%
rownames_to_column("Gene.stable.ID") %>%
filter(!is.na(pvalue) & !is.na(log2FoldChange)) %>%
left_join(gene_names, by = "Gene.stable.ID")
plot_volcano(df_SP1.3, top_n = 20)
# Heatmap de clústeres ajustado por las variables biológicas
vsd_SP1 <- vst(dds_SP1, blind = FALSE)
plot_combined_heatmap(list(df_SP1.1,df_SP1.2,df_SP1.3),
vsd = vsd_SP1,
dds = dds_SP1,
top_n = 20
)
################################################################################
# GSEA
################################################################################
set.seed(123)
# Cargar bases de datos
gene_sets_c2 <- msigdbr(species = "Homo sapiens", category = "C2", subcollection = "CP:REACTOME") %>%
select(gs_name, ensembl_gene)
gene_sets_H <- msigdbr(species = "Homo sapiens", category = "H" ) %>%
select(gs_name, ensembl_gene)
gene_sets <- rbind(gene_sets_H, gene_sets_c2)
# GSEA para los contrastes enre los clústeres que preservan las variables biológicas
gse0.1 <- run_gsea(df_SP0.1, gene_sets)
gse0.2 <- run_gsea(df_SP0.2, gene_sets)
gse0.3 <- run_gsea(df_SP0.3, gene_sets)
plot_gsea_dotplot(gse0.1, plot_title = "SP0 Contraste C2vsC1")
plot_gsea_dotplot(gse0.2, plot_title = "SP0 Contraste C2vsC3")
plot_gsea_dotplot(gse0.3, plot_title = "SP0 Contraste C3vsC1")
# GSEA para los contrastes enre los clústeres ajustados por las variables biológicas
gse1.1 <- run_gsea(df_SP1.1, gene_sets)
gse1.2 <- run_gsea(df_SP1.2, gene_sets)
gse1.3 <- run_gsea(df_SP1.3, gene_sets)
plot_gsea_dotplot(gse1.1, plot_title = "SP1 Contraste C2vsC1")
plot_gsea_dotplot(gse1.2, plot_title = "SP1 Contraste C2vsC3")
plot_gsea_dotplot(gse1.3, plot_title = "SP1 Contraste C3vsC1")
sessionInfo()
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8 LC_CTYPE=Spanish_Spain.utf8
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Spain.utf8
##
## time zone: Europe/Madrid
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] enrichplot_1.30.5 ggridges_0.5.7
## [3] fgsea_1.36.2 msigdbr_25.1.1
## [5] clusterProfiler_4.18.4 survminer_0.5.2
## [7] ggpubr_0.6.3 survival_3.8-6
## [9] data.table_1.18.2.1 ggrepel_0.9.7
## [11] pheatmap_1.0.13 ConsensusClusterPlus_1.74.0
## [13] limma_3.66.0 Rtsne_0.17
## [15] umap_0.2.10.0 sva_3.58.0
## [17] BiocParallel_1.44.0 genefilter_1.92.0
## [19] mgcv_1.9-4 nlme_3.1-168
## [21] apeglm_1.32.0 DESeq2_1.50.2
## [23] SummarizedExperiment_1.40.0 Biobase_2.70.0
## [25] MatrixGenerics_1.22.0 GenomicRanges_1.62.1
## [27] Seqinfo_1.0.0 IRanges_2.44.0
## [29] S4Vectors_0.48.0 BiocGenerics_0.56.0
## [31] generics_0.1.4 tximport_1.38.2
## [33] gridExtra_2.3 RColorBrewer_1.1-3
## [35] patchwork_1.3.2 matrixStats_1.5.0
## [37] lubridate_1.9.5 forcats_1.0.1
## [39] stringr_1.6.0 dplyr_1.2.0
## [41] purrr_1.2.1 readr_2.2.0
## [43] tidyr_1.3.2 tibble_3.3.1
## [45] ggplot2_4.0.2 tidyverse_2.0.0
##
## loaded via a namespace (and not attached):
## [1] splines_4.5.2 ggplotify_0.1.3 R.oo_1.27.1
## [4] polyclip_1.10-7 XML_3.99-0.23 lifecycle_1.0.5
## [7] mixsqp_0.3-54 rstatix_0.7.3 edgeR_4.8.2
## [10] vroom_1.7.0 lattice_0.22-9 MASS_7.3-65
## [13] backports_1.5.1 magrittr_2.0.4 sass_0.4.10
## [16] rmarkdown_2.31 jquerylib_0.1.4 yaml_2.3.12
## [19] otel_0.2.0 ggtangle_0.1.2 askpass_1.2.1
## [22] reticulate_1.45.0 cowplot_1.2.0 DBI_1.3.0
## [25] abind_1.4-8 R.utils_2.13.0 yulab.utils_0.2.4
## [28] tweenr_2.0.3 rappdirs_0.3.4 gdtools_0.5.1
## [31] irlba_2.3.7 tidytree_0.4.7 RSpectra_0.16-2
## [34] annotate_1.88.0 codetools_0.2-20 DelayedArray_0.36.0
## [37] xml2_1.5.2 ggtext_0.1.2 ggforce_0.5.0
## [40] DOSE_4.4.0 tidyselect_1.2.1 aplot_0.2.9
## [43] farver_2.1.2 jsonlite_2.0.0 Formula_1.2-5
## [46] systemfonts_1.3.2 bbmle_1.0.25.1 tools_4.5.2
## [49] ggnewscale_0.5.2 treeio_1.34.0 Rcpp_1.1.1
## [52] glue_1.8.0 SparseArray_1.10.9 xfun_0.57
## [55] qvalue_2.42.0 withr_3.0.2 numDeriv_2016.8-1.1
## [58] fastmap_1.2.0 openssl_2.4.2 truncnorm_1.0-9
## [61] digest_0.6.39 timechange_0.4.0 R6_2.6.1
## [64] gridGraphics_0.5-1 GO.db_3.22.0 RSQLite_2.4.6
## [67] R.methodsS3_1.8.2 fontLiberation_0.1.0 htmlwidgets_1.6.4
## [70] httr_1.4.8 S4Arrays_1.10.1 scatterpie_0.2.6
## [73] pkgconfig_2.0.3 gtable_0.3.6 blob_1.3.0
## [76] S7_0.2.1 XVector_0.50.0 htmltools_0.5.9
## [79] fontBitstreamVera_0.1.1 carData_3.0-6 scales_1.4.0
## [82] png_0.1-8 ashr_2.2-63 ggfun_0.2.0
## [85] knitr_1.51 rstudioapi_0.18.0 tzdb_0.5.0
## [88] reshape2_1.4.5 curl_7.0.0 coda_0.19-4.1
## [91] bdsmatrix_1.3-7 cachem_1.1.0 parallel_4.5.2
## [94] AnnotationDbi_1.72.0 pillar_1.11.1 grid_4.5.2
## [97] vctrs_0.7.1 tidydr_0.0.6 car_3.1-5
## [100] xtable_1.8-8 cluster_2.1.8.2 evaluate_1.0.5
## [103] invgamma_1.2 mvtnorm_1.3-3 cli_3.6.5
## [106] locfit_1.5-9.12 compiler_4.5.2 rlang_1.1.7
## [109] crayon_1.5.3 SQUAREM_2026.1 ggsignif_0.6.4
## [112] labeling_0.4.3 emdbook_1.3.14 plyr_1.8.9
## [115] fs_2.1.0 ggiraph_0.9.6 stringi_1.8.7
## [118] viridisLite_0.4.3 babelgene_22.9 assertthat_0.2.1
## [121] Biostrings_2.78.0 lazyeval_0.2.3 GOSemSim_2.36.0
## [124] fontquiver_0.2.1 Matrix_1.7-4 hms_1.1.4
## [127] bit64_4.6.0-1 KEGGREST_1.50.0 statmod_1.5.1
## [130] gridtext_0.1.6 igraph_2.3.1 broom_1.0.13
## [133] memoise_2.0.1 bslib_0.11.0 ggtree_4.0.5
## [136] fastmatch_1.1-8 bit_4.6.0 gson_0.1.0
## [139] ape_5.8-1