Objetivo: Incorporar la capa de proteína al análisis. El flujo de trabajo de RNA hasta el Bloque 5 ignoró por completo el ADT. El valor añadido de CITE-seq es la capacidad de resolver el dropout, validar la anotación y ponderar las modalidades por célula.
Orden de las operaciones: inspeccionar la capa de ADT y contrastarla con RNA, normalizar, hacer control de calidad (QC) del panel, contrastar RNA vs ADT para marcadores coincidentes, realizar gating, verificar canales intercambiados, construir el PCA de ADT, y luego integrar con WNN.
Cargar el archivo en R:
library(here) # install.packages("here")here() starts at /private/var/folders/zz/cdcfnwts145c0xjqs_c0y6mm0000gn/T/RtmppqY2dt/file177bf1a6b81e1sobj_filt <- readRDS(
here("checkpoints", "outputs", "sobj_preprocessed.rds")
)
# top markers
top_markers <- readRDS(
here("checkpoints", "outputs", "top_markers.rds")
)DefaultAssay(sobj_filt) <- "ADT"
adt_counts <- LayerData(sobj_filt, assay = "ADT", layer = "counts")
cat("ADT count matrix:\n")ADT count matrix:cat(" Proteins (rows):", nrow(adt_counts), "\n") Proteins (rows): 24 cat(" Cells (cols) :", ncol(adt_counts), "\n") Cells (cols) : 2673 cat(" Class :", class(adt_counts), "\n") Class : dgCMatrix cat(" Layers in ADT :", paste(SafeLayers(sobj_filt[["ADT"]]), collapse = ", "), "\n\n") Layers in ADT : counts, data Resumen por proteína: escala y distribución:
adt_summary <- data.frame(
protein = rownames(adt_counts),
min = apply(adt_counts, 1, min),
median = apply(adt_counts, 1, median),
mean = round(Matrix::rowMeans(adt_counts), 2),
max = apply(adt_counts, 1, max),
pct_zero = round(Matrix::rowMeans(adt_counts == 0) * 100, 1)
)
cat("Per-protein distribution (first 8 proteins):\n")Per-protein distribution (first 8 proteins):print(head(adt_summary, 8), row.names = FALSE) protein min median mean max pct_zero
CD3 0 1 15.36 168 30.9
CD4 0 1 14.08 196 35.6
CD8a 0 1 6.00 160 40.9
CD14 0 1 21.64 294 32.7
CD16 0 1 9.03 186 37.3
CD19 0 1 10.30 209 38.4
CD20 0 1 8.16 233 39.9
CD25 0 1 0.80 5 45.8Conteo total de ADT por célula
total_adt_per_cell <- Matrix::colSums(adt_counts)
cat("\nPer-cell total ADT counts:\n")
Per-cell total ADT counts:print(round(quantile(total_adt_per_cell, c(0.05, 0.25, 0.5, 0.75, 0.95)), 0)) 5% 25% 50% 75% 95%
67 118 149 184 244 DefaultAssay(sobj_filt) <- "RNA"Bimodal saludable: mediana cercana a 0, media unas pocas unidades más alta, máximo en los cientos. La proteína no está teñida en la mayoría de las células y es fuertemente positiva en su subconjunto objetivo. Patrón de anticuerpo muerto: media por debajo de 1, máximo también bajo (por debajo de 10), pct_zero muy alto. Nada se eleva por encima del fondo porque el anticuerpo no es funcional. En el conjunto de datos inyectado, CD56 muestra el patrón muerto; CD3, CD4, CD8a, CD14, CD19 muestran el patrón bimodal. Examina la tabla ordenada por media para detectar la anomalía.
DEMOSTRACIÓN EN VIVO:
DefaultAssay(sobj_filt) <- "ADT"
adt <- LayerData(sobj_filt, assay="ADT", layer="counts")
qc <- data.frame(
protein = rownames(adt),
median = apply(adt, 1, median),
mean = round(Matrix::rowMeans(adt), 2),
max = apply(adt, 1, max),
pct_zero = round(Matrix::rowMeans(adt == 0)*100, 1)
)
qc[order(qc$mean), ] protein median mean max pct_zero
PD1 PD1 1 0.79 5 46.0
CD25 CD25 1 0.80 5 45.8
CD45RO CD45RO 1 0.80 5 44.0
LAG3 LAG3 1 0.80 6 43.7
CD69 CD69 1 0.81 6 44.1
CD197 CD197 1 0.81 5 44.1
TIGIT TIGIT 1 0.81 5 44.4
CD62L CD62L 1 0.82 6 43.0
CD86 CD86 1 2.66 117 44.7
CD45RA CD45RA 1 2.86 84 40.4
CD57 CD57 1 3.38 66 39.5
CD27 CD27 1 3.66 163 42.0
IgD IgD 1 4.17 79 39.5
CD8a CD8a 1 6.00 160 40.9
CD38 CD38 1 6.07 204 41.3
CD127 CD127 1 7.56 87 36.1
CD20 CD20 1 8.16 233 39.9
CD56 CD56 1 8.37 181 38.5
CD16 CD16 1 9.03 186 37.3
CD19 CD19 1 10.30 209 38.4
CD4 CD4 1 14.08 196 35.6
CD3 CD3 1 15.36 168 30.9
CD14 CD14 1 21.64 294 32.7
HLADR HLADR 1 21.82 219 29.3DefaultAssay(sobj_filt) <- "RNA"GRÁFICO / SALIDA:
También podría usarse:
RidgePlot(sobj_filt, features=rownames(sobj_filt[['ADT']])[1:6], assay='ADT')Picking joint bandwidth of 0.346Picking joint bandwidth of 0.174Picking joint bandwidth of 0.11Picking joint bandwidth of 0.449Picking joint bandwidth of 0.149Picking joint bandwidth of 0.141
El Bloque 1 (Sección 4.4) midió la dispersión (sparsity) de RNA por sí sola. Ahora que la capa de ADT también ha sido inspeccionada, la comparación es el punto clave: las mismas células, dos ensayos, tasas de ceros muy diferentes.
rna_counts_b6 <- LayerData(sobj_filt, assay = "RNA", layer = "counts")
rna_sparsity <- 1 - Matrix::nnzero(rna_counts_b6) / prod(dim(rna_counts_b6))
adt_sparsity <- 1 - Matrix::nnzero(adt_counts) / prod(dim(adt_counts))
cat("RNA sparsity:", round(rna_sparsity * 100, 1), "%\n")RNA sparsity: 90.5 %cat("ADT sparsity:", round(adt_sparsity * 100, 1), "%\n")ADT sparsity: 40.1 %ADT ha permanecido intacta desde el Bloque 1 (Sección 4.4). Confirma su estado explícitamente antes de normalizar.
DefaultAssay(sobj_filt) <- "ADT"
cat("ADT assay current state:\n")ADT assay current state:cat(" Active assay : ADT\n") Active assay : ADTcat(" Proteins :", nrow(sobj_filt[["ADT"]]), "\n") Proteins : 24 cat(" Cells (post-QC):", ncol(sobj_filt), "\n") Cells (post-QC): 2673 cat(" Layers :", paste(SafeLayers(sobj_filt[["ADT"]]), collapse = ", "), "\n") Layers : counts, data En este punto esperamos: solo counts. Sin data, sin scale.data. La normalización es el siguiente paso.
Normalización CLR (Centered Log-Ratio, log-ratio centrado). El argumento margin establece la dirección: + margin = 1: normaliza cada proteína a través de todas las células (medias de fila ~0) + margin = 2: normaliza cada célula a través de todas las proteínas (medias de columna ~0)
margin = 2 es correcto para CITE-seq. Elimina la variación por célula en la captura total de anticuerpos, análogo a la normalización por tamaño de biblioteca en RNA. margin = 1 se ejecuta sin error pero destruye la señal biológica. Diagnosticarás un objeto real normalizado de forma incorrecta en el Bloque 8, Escenario 3.
DefaultAssay(sobj_filt) <- "ADT"
sobj_filt <- NormalizeData(
sobj_filt,
normalization.method = "CLR",
margin = 2 # CORRECT: normalizes each cell across proteins
)Normalizing layer: countsNormalizing across cellsadt_m2 <- LayerData(sobj_filt, layer = "data")
cat("Column means after margin=2 (per cell), expected ~0:\n")Column means after margin=2 (per cell), expected ~0:print(round(colMeans(adt_m2)[1:8], 4))CELL000011 CELL000021 CELL000031 CELL000041 CELL000051 CELL000061 CELL000071
0.5889 0.6480 0.4619 0.6050 0.5824 0.6226 0.5801
CELL000081
0.5344 cat("\nRow means (per protein), should NOT be ~0 (biological variation preserved):\n")
Row means (per protein), should NOT be ~0 (biological variation preserved):print(round(rowMeans(adt_m2), 4)) CD3 CD4 CD8a CD14 CD16 CD19 CD20 CD25 CD27 CD38 CD45RA
1.0745 0.8865 0.5040 1.1483 0.7338 0.7100 0.5979 0.2592 0.4233 0.4623 0.4362
CD45RO CD56 CD57 CD62L CD69 CD86 CD127 CD197 PD1 TIGIT LAG3
0.2617 0.6312 0.5058 0.2677 0.2645 0.3738 0.7574 0.2646 0.2555 0.2631 0.2601
IgD HLADR
0.5157 1.2983 DefaultAssay(sobj_filt) <- "RNA"
cat("\nADT normalized correctly.\n")
ADT normalized correctly.Una conjugación fallida o un anticuerpo degradado producen un “canal muerto”: la proteína se lee cerca de cero en cada célula, con casi ninguna varianza. No genera error. Si haces gating o anotas sobre un canal muerto, pierdes silenciosamente toda una población.
DefaultAssay(sobj_filt) <- "ADT"
adt_counts <- LayerData(sobj_filt, layer = "counts")
panel_qc <- data.frame(
protein = rownames(adt_counts),
median = round(apply(adt_counts, 1, median), 2),
mean = round(Matrix::rowMeans(adt_counts), 2),
pct_zero = round(Matrix::rowMeans(adt_counts == 0) * 100, 1),
max = apply(adt_counts, 1, max)
)
panel_qc <- panel_qc[order(panel_qc$mean), ]
cat("ADT panel QC (sorted by mean; suspicious = bottom rows):\n")ADT panel QC (sorted by mean; suspicious = bottom rows):print(head(panel_qc, 6), row.names = FALSE) protein median mean pct_zero max
PD1 1 0.79 46.0 5
CD25 1 0.80 45.8 5
CD45RO 1 0.80 44.0 5
LAG3 1 0.80 43.7 6
CD69 1 0.81 44.1 6
CD197 1 0.81 44.1 5DefaultAssay(sobj_filt) <- "RNA"En el conjunto de datos inyectado: CD56 (codificada por NCAM1). CD56 marca las células NK. Con CD56 ADT muerta, cualquier anotación de células NK basada solo en ADT falla. El bloque de diagnóstico compara CD56 ADT (casi cero) con NCAM1 RNA (señal clara en el cluster NK) y expone la inconsistencia. Hábito: siempre verifica las proteínas de ADT con menor media contra sus contrapartes de RNA.
DEMOSTRACIÓN EN VIVO:
FeaturePlot(sobj_filt, c("CD56", "NCAM1"), reduction = "umap")Warning: Could not find CD56 in the default search locations, found in 'ADT'
assay instead
GRÁFICO / SALIDA: Panel FeaturePlot: CD56 ADT (plano) junto a NCAM1 RNA (positivo en un cluster)
Compara la proteína sospechosa (ADT) con su gen de RNA lado a lado.
Edita suspect_adt / suspect_rna a la proteína marcada anteriormente.
suspect_adt <- "CD56" # protein flagged as near-zero in panel QC
suspect_rna <- "NCAM1" # its RNA counterpart (NK marker)Decide el título a partir de los números reales de panel_qc en lugar de afirmar una afirmación fija. En el objeto limpio esta proteína está saludable; en el objeto inyectado está muerta. El título debe indicar lo que sea cierto para el objeto que realmente está cargado.
suspect_row <- panel_qc[panel_qc$protein == suspect_adt, ]
is_dead <- nrow(suspect_row) == 1 && suspect_row$pct_zero > 90 && suspect_row$mean < 1
adt_title <- if (is_dead) {
paste0(suspect_adt, " ADT\n(flat: likely failed antibody)")
} else {
paste0(suspect_adt, " ADT\n(signal present: not a dead channel here)")
}
DefaultAssay(sobj_filt) <- "ADT"
p_dead <- FeaturePlot(sobj_filt, suspect_adt, min.cutoff = "q05") +
ggtitle(adt_title)
DefaultAssay(sobj_filt) <- "RNA"
p_alive <- FeaturePlot(sobj_filt, suspect_rna, min.cutoff = "q05") +
ggtitle(paste0(suspect_rna, " RNA\n(clear NK signal: protein should exist)"))
p_dead | p_alive
DefaultAssay(sobj_filt) <- "RNA"Comparación directa. El mRNA de CD4 se lee como cero en muchas células T CD4+ (dropout), mientras que la proteína CD4 de ADT es bimodal a través de las mismas células.
DefaultAssay(sobj_filt) <- "RNA"
p_rna <- FeaturePlot(sobj_filt, "CD4", min.cutoff = "q05") +
ggtitle("CD4 via RNA\n(dropout: many CD4+ cells read as zero)") +
scale_color_gradient(low = "lightgrey", high = "#991b1b")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.DefaultAssay(sobj_filt) <- "ADT"
p_adt <- FeaturePlot(sobj_filt, "CD4", min.cutoff = "q05") +
ggtitle("CD4 via ADT protein\n(bimodal, reliable)") +
scale_color_gradient(low = "lightgrey", high = "#0d7377")Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.p_rna | p_adt
DefaultAssay(sobj_filt) <- "RNA"Para cada proteína de ADT con una contraparte de RNA, calcula la correlación de Spearman a través de todas las células. Correlación baja = dropout alto en RNA.
adt_rna_pairs <- list(
CD3 = "CD3E",
CD4 = "CD4",
CD8a = "CD8A",
CD14 = "CD14",
CD19 = "CD19",
CD56 = "NCAM1"
)
DefaultAssay(sobj_filt) <- "ADT"
adt_vals <- t(LayerData(sobj_filt, layer = "data"))
DefaultAssay(sobj_filt) <- "RNA"
rna_vals <- t(LayerData(sobj_filt, layer = "data"))
cat("Spearman correlation (RNA vs ADT) per marker:\n")Spearman correlation (RNA vs ADT) per marker:cat(sprintf(" %-8s %s\n", "Marker", "Spearman rho")) Marker Spearman rhocat(sprintf(" %-8s %s\n", "------", "------------")) ------ ------------for (adt_name in names(adt_rna_pairs)) {
rna_gene <- adt_rna_pairs[[adt_name]]
if (rna_gene %in% colnames(rna_vals) &&
adt_name %in% colnames(adt_vals)) {
rho <- cor(rna_vals[, rna_gene],
adt_vals[, adt_name],
method = "spearman")
cat(sprintf(" %-8s %.3f\n", adt_name, rho))
}
} CD3 0.753
CD4 0.585
CD8a 0.225
CD14 0.727
CD19 0.502
CD56 0.540DefaultAssay(sobj_filt) <- "RNA"Para cada marcador, encuentra las células que son claramente positivas para la proteína (ADT > umbral) pero tienen cero conteos de RNA. Estas son las células que se anotarían incorrectamente con un análisis basado solo en RNA.
DefaultAssay(sobj_filt) <- "ADT"
adt_data <- LayerData(sobj_filt, layer = "data")
DefaultAssay(sobj_filt) <- "RNA"
rna_counts <- LayerData(sobj_filt, layer = "counts")
cat("Dropout analysis: cells positive by ADT but zero by RNA\n")Dropout analysis: cells positive by ADT but zero by RNAcat(sprintf(" %-8s %-12s %-12s %s\n",
"Marker", "ADT+ cells", "RNA==0 among ADT+", "Dropout rate")) Marker ADT+ cells RNA==0 among ADT+ Dropout ratecat(sprintf(" %-8s %-12s %-12s %s\n",
"------", "----------", "-----------------", "------------")) ------ ---------- ----------------- ------------marker_pairs <- list(CD4 = "CD4", CD14 = "CD14", CD19 = "CD19")
for (adt_name in names(marker_pairs)) {
rna_gene <- marker_pairs[[adt_name]]
if (rna_gene %in% rownames(rna_counts) && adt_name %in% rownames(adt_data)) {
adt_pos <- which(adt_data[adt_name, ] > 1.0)
rna_zero <- sum(rna_counts[rna_gene, adt_pos] == 0)
pct <- round(rna_zero / length(adt_pos) * 100, 1)
cat(sprintf(" %-8s %-12d %-12d %s%%\n",
adt_name, length(adt_pos), rna_zero, pct))
}
} CD4 595 12 2%
CD14 811 126 15.5%
CD19 480 171 35.6%El gating digital replica los gráficos de dispersión biaxiales de la citometría de flujo usando datos de ADT. Permite la clasificación computacional de células con la misma lógica que usan los inmunólogos en el laboratorio.
DefaultAssay(sobj_filt) <- "ADT"
color_by <- if ("singler_label_top" %in% colnames(sobj_filt@meta.data)) {
"singler_label_top"
} else {
"cell_type"
}
gate_data <- FetchData(sobj_filt, vars = c("CD4", "CD8a", color_by))
names(gate_data)[3] <- "label"
ggplot(gate_data, aes(CD4, CD8a, color = label)) +
geom_point(alpha = 0.4, size = 0.8) +
geom_density_2d(color = "grey50", linewidth = 0.3) +
geom_vline(xintercept = 1.5, linetype = "dashed", color = "#991b1b") +
geom_hline(yintercept = 1.5, linetype = "dashed", color = "#991b1b") +
annotate("text", x = 3.5, y = 0.4, label = "CD4+ T cells", size = 3.5) +
annotate("text", x = 0.4, y = 3.5, label = "CD8+ T cells", size = 3.5) +
labs(title = "Digital gating: CD4 vs CD8a (ADT protein)",
subtitle = "Only possible with ADT; RNA dropout makes this scatter uninformative",
x = "CD4 (CLR normalized)", y = "CD8a (CLR normalized)",
color = NULL) +
theme_classic(base_size = 13) +
theme(legend.text = element_text(size = 7))
DefaultAssay(sobj_filt) <- "RNA"Una llamada a función que “funciona” en el ensayo equivocado produce un error silencioso.
FeaturePlot("CD4") con DefaultAssay = "RNA" grafica el gen CD4 de RNA.FeaturePlot("CD4") con DefaultAssay = "ADT" grafica la proteína CD4.El gráfico se RENDERIZA en ambos casos, y los dos gráficos pueden parecer engañosamente similares: min.cutoff = "q05" reescala el gradiente de color de cada panel a su propio rango de valores, de modo que una señal de RNA dispersa y una señal de ADT densa terminan ambas estiradas a través de una escala de color de apariencia similar. El fallo es silencioso precisamente porque nada en el gráfico renderizado te advierte de qué ensayo lo produjo.
Cuantifica la diferencia que el gráfico por sí solo oculta: qué fracción de células lee exactamente cero en cada modalidad para este marcador.
DefaultAssay(sobj_filt) <- "RNA"
cd4_rna_zero_pct <- mean(LayerData(sobj_filt, layer = "counts")["CD4", ] == 0) * 100
DefaultAssay(sobj_filt) <- "ADT"
cd4_adt_zero_pct <- mean(LayerData(sobj_filt, layer = "counts")["CD4", ] == 0) * 100
DefaultAssay(sobj_filt) <- "RNA"
cat(sprintf("CD4 zero rate: %.1f%% of cells in RNA, %.1f%% of cells in ADT\n",
cd4_rna_zero_pct, cd4_adt_zero_pct))CD4 zero rate: 67.3% of cells in RNA, 35.6% of cells in ADTcat("The plots below can look similar at a glance; the zero rates above are\n")The plots below can look similar at a glance; the zero rates above arecat("the actual difference the two assays are reporting for the same gene.\n")the actual difference the two assays are reporting for the same gene.Demuéstralo dibujando la misma llamada bajo ambos ensayos. Los títulos se mantienen cortos y el tamaño de fuente está restringido, ya que los dos gráficos se renderizan lado a lado a la mitad del ancho cada uno.
DefaultAssay(sobj_filt) <- "RNA"
p_as_rna <- FeaturePlot(sobj_filt, "CD4", min.cutoff = "q05") +
ggtitle("'CD4' | DefaultAssay = RNA") +
theme(plot.title = element_text(size = 11))
DefaultAssay(sobj_filt) <- "ADT"
p_as_adt <- FeaturePlot(sobj_filt, "CD4", min.cutoff = "q05") +
ggtitle("'CD4' | DefaultAssay = ADT") +
theme(plot.title = element_text(size = 11))
p_as_rna | p_as_adt
DefaultAssay(sobj_filt) <- "RNA"Tendrías el ensayo que estuviera activo en el momento de la llamada al gráfico, y un lector no podría distinguir cuál a partir de la figura sola: el título del gráfico dice solo ‘CD4’, los ejes son coordenadas UMAP idénticas, y min.cutoff='q05' reescala el gradiente de color propio de cada panel de forma independiente, de modo que una señal de RNA dispersa y una señal de ADT más densa terminan estiradas a través de escalas de apariencia visual similar. Esta es exactamente la razón por la que importan los números de tasa de cero: muestran la diferencia real en qué fracción de células se lee como positiva en cada modalidad, una diferencia que el gradiente de color reescalado oculta. La figura publicada podría mostrar silenciosamente CD4 RNA cuando se pretendía CD4 proteína, o viceversa, y nada en la imagen misma señalaría la discrepancia. Hábito: establece DefaultAssay explícitamente al inicio de cualquier bloque de graficado, Y/O pasa el ensayo dentro de la llamada; mejor aún, edita el título del gráfico para incluir el nombre del ensayo y considera reportar la tasa de cero junto con la figura para que la diferencia de modalidad quede documentada, no solo visible para quien ya sepa buscarla.
DEMOSTRACIÓN EN VIVO:
DefaultAssay(sobj_filt) <- "RNA"
mean(LayerData(sobj_filt, layer="counts")["CD4", ] == 0) * 100 # RNA zero rate[1] 67.34007DefaultAssay(sobj_filt) <- "ADT"
mean(LayerData(sobj_filt, layer="counts")["CD4", ] == 0) * 100 # ADT zero rate[1] 35.61541DefaultAssay(sobj_filt) <- "RNA"
FeaturePlot(sobj_filt, "CD4") + ggtitle("CD4 protein (ADT)")
FeaturePlot(sobj_filt, "CD4") + ggtitle("CD4 mRNA (RNA)")
GRÁFICO / SALIDA: Dos FeaturePlots lado a lado, uno por ensayo, con títulos distintos; los números de tasa de cero explican lo que los gráficos solos no muestran.
Lee el código de abajo. Parece correcto. ¿Qué está mal?
DefaultAssay(sobj_filt) <- "RNA"
adt_markers <- FindMarkers(sobj_filt, ident.1 = "CD4Tcell",
group.by = "cell_type",
min.pct = 0.25)For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
(default method for FindMarkers) please install the presto package
--------------------------------------------
install.packages('devtools')
devtools::install_github('immunogenomics/presto')
--------------------------------------------
After installation of presto, Seurat will automatically use the more
efficient implementation (no further action necessary).
This message will be shown once per sessionTu respuesta:
REVELAR: la llamada pide marcadores de células T CD4 pero DefaultAssay es RNA. Devuelve marcadores de RNA, no marcadores de ADT. La salida es correcta para RNA pero no responde la pregunta “qué proteínas distinguen a las células T CD4”. O bien establece DefaultAssay = "ADT" primero, o pasa assay = "ADT" dentro de la llamada a FindMarkers. El min.pct = 0.25 también es incorrecto para ADT (diseñado para RNA disperso); para ADT usa un valor como 0.5 o 0.0 con logfc.threshold.
Un canal intercambiado o mal etiquetado se ejecuta perfectamente. CLR lo normaliza, FeaturePlot lo dibuja, el gating lo agrupa. Nada da error. La única señal es que la biología es imposible. La defensa es verificar cada proteína contra una etiqueta independiente (su marcador de RNA).
Cada proteína debería correlacionarse con su propio gen de RNA a través de las células.
Un marcador canónico con rho cercano a 0 es un canal mal etiquetado.
DefaultAssay(sobj_filt) <- "ADT"
adt_d <- t(LayerData(sobj_filt, layer = "data"))
DefaultAssay(sobj_filt) <- "RNA"
rna_d <- t(LayerData(sobj_filt, layer = "data"))
check_pairs <- list(CD4 = "CD4", CD8a = "CD8A", CD19 = "CD19",
CD14 = "CD14", CD3 = "CD3E")
cat(sprintf(" %-6s %-12s %-10s\n", "ADT", "rho(self)", "verdict")) ADT rho(self) verdict cat(sprintf(" %-6s %-12s %-10s\n", "---", "---------", "-------")) --- --------- ------- for (a in names(check_pairs)) {
g <- check_pairs[[a]]
if (a %in% colnames(adt_d) && g %in% colnames(rna_d)) {
rho <- round(cor(adt_d[, a], rna_d[, g], method = "spearman"), 3)
verdict <- if (rho < 0.05) "SUSPECT" else "ok"
cat(sprintf(" %-6s %-12s %-10s\n", a, rho, verdict))
}
} CD4 0.585 ok
CD8a 0.225 ok
CD19 0.502 ok
CD14 0.727 ok
CD3 0.753 ok Dos canales de ADT muestran rho cercano a 0 contra su propio gen de RNA, mientras que el par cruzado es alto. Confirma el intercambio con la visualización de abajo.
DefaultAssay(sobj_filt) <- "ADT"
g <- FetchData(sobj_filt, vars = c("CD14", "CD19"))
DefaultAssay(sobj_filt) <- "RNA"
g$CD19_rna <- FetchData(sobj_filt, vars = "CD19")[, 1] # B cell RNA marker
ggplot(g, aes(CD14, CD19, color = CD19_rna)) +
geom_point(alpha = 0.5, size = 0.8) +
scale_color_gradient(low = "lightgrey", high = "#991b1b") +
labs(title = "CD14 vs CD19 (ADT), colored by CD19 RNA",
subtitle = "If CD19 RNA piles up on the CD14 axis, the labels are swapped",
x = "CD14 (ADT)", y = "CD19 (ADT)", color = "CD19 RNA") +
theme_classic(base_size = 12)
DefaultAssay(sobj_filt) <- "RNA"Rho de Spearman entre cada proteína de ADT y su contraparte de RNA a través de las células. Los marcadores saludables muestran rho de 0.3-0.7. Huella del intercambio: ADT_X vs RNA_X cerca de cero Y ADT_X vs RNA_Y alto; lo mismo con X,Y invertidos. En el conjunto de datos inyectado, los rownames de CD14 y CD19 en ADT están intercambiados. El bloque de corrección reconstruye el ensayo ADT con los rownames reintercambiados y vuelve a ejecutar CLR. Es condicional, así que en datos limpios no hace nada.
DEMOSTRACIÓN EN VIVO:
DefaultAssay(sobj_filt) <- "ADT"
FeatureScatter(sobj_filt, "CD14", "CD19") + ggtitle("ADT CD14 vs ADT CD19")
DefaultAssay(sobj_filt) <- "RNA"
FeatureScatter(sobj_filt, "CD14", "CD19") + ggtitle("RNA CD14 vs RNA CD19")
GRÁFICO / SALIDA: Dos gráficos de dispersión lado a lado. Saludable: los dispersogramas de ADT y RNA coinciden. Con un intercambio: divergen
Seguro para v5: reconstruye el ensayo ADT a partir de su matriz de conteos con los nombres corregidos.Condicional: solo actúa si el intercambio se detecta realmente (seguro en datos limpios).DefaultAssay(sobj_filt) <- "ADT"
adt_d <- t(LayerData(sobj_filt, layer = "data"))
rna_d <- t(LayerData(sobj_filt, assay = "RNA", layer = "data"))
swapped <- FALSE
if (all(c("CD14", "CD19") %in% colnames(adt_d))) {
rho14 <- cor(adt_d[, "CD14"], rna_d[, "CD14"], method = "spearman")
swapped <- is.na(rho14) || rho14 < 0.05
}
if (swapped) {
adt_counts <- LayerData(sobj_filt, assay = "ADT", layer = "counts")
rn <- rownames(adt_counts)
i14 <- which(rn == "CD14"); i19 <- which(rn == "CD19")
rn[c(i14, i19)] <- rn[c(i19, i14)]
rownames(adt_counts) <- rn
# CreateAssay5Object() builds a fresh assay with only the counts layer, so
# the [[<- replacement method warns that it differs in structure (no
# data/scale.data layers yet) from the ADT assay it replaces. Expected and
# harmless here; suppressWarnings() wraps the assignment itself, since the
# notice comes from the replacement method, not from CreateAssay5Object().
new_adt_assay <- CreateAssay5Object(counts = adt_counts)
suppressWarnings(sobj_filt[["ADT"]] <- new_adt_assay)
sobj_filt <- NormalizeData(sobj_filt, normalization.method = "CLR",
margin = 2, verbose = FALSE)
cat("ADT labels corrected (CD14 <-> CD19) and re-normalized.\n")
} else {
cat("No swap detected; ADT labels left unchanged.\n")
}No swap detected; ADT labels left unchanged.# Re-check
adt_d <- t(LayerData(sobj_filt, layer = "data"))
rna_d <- t(LayerData(sobj_filt, assay = "RNA", layer = "data"))
for (a in c("CD14", "CD19")) {
rho <- round(cor(adt_d[, a], rna_d[, a], method = "spearman"), 3)
cat(sprintf(" %-6s rho(self) = %s\n", a, rho))
} CD14 rho(self) = 0.727
CD19 rho(self) = 0.502DefaultAssay(sobj_filt) <- "RNA"La integración WNN (Weighted Nearest Neighbor, vecino más cercano ponderado) aprende la contribución relativa de cada modalidad por célula. Las células donde el RNA es informativo reciben más peso de RNA; las células donde la proteína es más discriminante reciben más peso de ADT. WNN requiere un PCA por modalidad. El PCA de RNA existe desde el Bloque 4; ahora construimos el PCA de ADT.
DefaultAssay(sobj_filt) <- "ADT"
adt_features <- rownames(sobj_filt[["ADT"]]) # all 24 proteins
sobj_filt <- ScaleData(
sobj_filt,
features = adt_features,
verbose = FALSE
)Con 24 proteínas, solicitar 20 PCs desencadena una advertencia de Seurat sobre calcular demasiados valores singulares vía SVD truncado. Ambas soluciones funcionan:
Usamos (a) porque 15 PCs es suficiente señal para 24 proteínas.
sobj_filt <- RunPCA(
sobj_filt,
features = adt_features,
npcs = 15,
reduction.name = "pca.adt",
reduction.key = "pcaADT_",
verbose = FALSE
)Warning in svd.function(A = t(x = object), nv = npcs, ...): You're computing
too large a percentage of total singular values, use a standard svd instead.cat("Reductions now available:", paste(SafeReductions(sobj_filt), collapse = ", "), "\n")Reductions now available: pca, umap, harmony, umap.harmony, pca.adt Seurat espera un nombre de peso por modalidad. Pasar solo “RNA.weight” desencadena una advertencia de que ADT.weight se ha asignado automáticamente. Pasa ambos explícitamente para silenciar el aviso y documentar los nombres de columna previstos.
sobj_filt <- FindMultiModalNeighbors(
sobj_filt,
reduction.list = list("pca", "pca.adt"),
dims.list = list(1:20, 1:15),
modality.weight.name = c("RNA.weight", "ADT.weight"),
verbose = FALSE
)
sobj_filt <- RunUMAP(
sobj_filt,
nn.name = "weighted.nn",
reduction.name = "wnn.umap",
reduction.key = "wnnUMAP_",
verbose = FALSE
)Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per sessionDefaultAssay(sobj_filt) <- "RNA"
cat("WNN complete. New reduction: wnn.umap\n")WNN complete. New reduction: wnn.umapWNN no necesariamente mueve cada tipo celular a una región visualmente distinta; los diseños de UMAP pueden rotar o reflejarse entre ejecuciones incluso cuando las vecindades subyacentes apenas cambian. Los dos paneles de abajo pueden parecer similares a primera vista. Lo que realmente importa es si el conjunto de células agrupadas cambió, no si la imagen se ve diferente.
El Paso 6.11c (más adelante) verifica esto directamente por tipo celular mediante una tabla cruzada; aquí, una comparación rápida de clustering da una primera lectura antes de eso.
sobj_filt <- FindClusters(sobj_filt, resolution = 0.5,
cluster.name = "rna_only_clusters_preview",
verbose = FALSE)
rna_only_clusters <- sobj_filt$rna_only_clusters_preview
sobj_filt <- FindClusters(sobj_filt, graph.name = "wsnn", resolution = 0.5,
cluster.name = "wnn_clusters_preview",
verbose = FALSE)
wnn_clusters_preview <- sobj_filt$wnn_clusters_preview
agreement_tab <- table(RNA_only = rna_only_clusters, WNN = wnn_clusters_preview)
cat("Cluster membership: RNA-only vs WNN (rows: RNA-only clusters, cols: WNN clusters)\n")Cluster membership: RNA-only vs WNN (rows: RNA-only clusters, cols: WNN clusters)print(agreement_tab) WNN
RNA_only 0 1 2 3 4 5 6 7 8
0 527 0 0 0 0 0 0 0 0
1 0 329 0 0 131 0 0 0 0
2 0 254 0 0 93 0 0 0 0
3 0 0 221 0 0 0 0 0 0
4 0 0 0 198 0 0 0 0 0
5 9 0 1 0 1 166 0 4 16
6 165 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 134 0 1
8 0 0 118 0 0 0 0 0 0
9 0 0 0 0 0 0 0 116 0
10 0 0 0 107 0 0 0 0 0
11 0 0 0 0 0 0 0 0 82Una puntuación de acuerdo simple: para cada cluster solo-RNA, qué fracción de sus células cae en el único cluster de WNN que captura la mayoría de ellas.
best_match_frac <- apply(agreement_tab, 1, function(row) max(row) / sum(row))
cat(sprintf("\nMedian per-cluster agreement: %.1f%% of cells stay grouped together\n",
100 * median(best_match_frac)))
Median per-cluster agreement: 100.0% of cells stay grouped togetherp_rna <- DimPlot(sobj_filt, reduction = "umap",
group.by = "singler_label_top",
label = TRUE, label.size = 2.5, repel = TRUE) +
ggtitle("RNA-only UMAP") + NoLegend()
p_wnn <- DimPlot(sobj_filt, reduction = "wnn.umap",
group.by = "singler_label_top",
label = TRUE, label.size = 2.5, repel = TRUE) +
ggtitle("WNN UMAP (RNA + ADT)") + NoLegend()
p_rna | p_wnn
RNA.weight cercano a 1 = la célula está mejor caracterizada por RNARNA.weight cercano a 0 = la célula está mejor caracterizada por ADTRNA.weight y ADT.weight suman 1 para cada célula (FindMultiModalNeighbors los normaliza de esa forma), así que 0.5 es el punto de referencia natural: por debajo de él, ADT está contribuyendo más que RNA a la ubicación de esa célula en el grafo de WNN, no solo “alguna cantidad” de información de proteína.pct_below_half <- round(mean(sobj_filt$RNA.weight < 0.5) * 100, 1)
cat(sprintf("Cells where ADT contributes more than RNA (RNA.weight < 0.5): %.1f%%\n",
pct_below_half))Cells where ADT contributes more than RNA (RNA.weight < 0.5): 74.8%ggplot(sobj_filt@meta.data, aes(x = RNA.weight, fill = singler_label_top)) +
geom_histogram(bins = 40, alpha = 0.8) +
geom_vline(xintercept = 0.5, linetype = "dashed", color = "#2c3142") +
facet_wrap(~singler_label_top, scales = "free_y") +
labs(title = "Per-cell RNA modality weight from WNN",
subtitle = "Dashed line at 0.5: left of it, ADT outweighs RNA for that cell",
x = "RNA weight (1 = fully RNA, 0 = fully ADT)") +
theme_classic(base_size = 10) +
theme(legend.position = "none",
strip.text = element_text(size = 7))
Las células T CD4 y las células T CD8, por un amplio margen. Sus marcadores definitorios (CD4, CD8A) tienen el peor dropout de RNA. Las células B (MS4A1, CD79A) y los monocitos (CD14, LYZ) tienen una señal de RNA más fuerte y dependen menos de ADT. Las células NK (NKG7, GNLY, NCAM1) son intermedias. El histograma facetado de RNA.weight por tipo celular hace esto concreto, y la línea discontinua en RNA.weight = 0.5 marca el punto en el que ADT supera a RNA para esa célula, no solo contribuye con alguna cantidad de información.
DEMOSTRACIÓN EN VIVO:
mean(sobj_filt$RNA.weight < 0.5) * 100 # cells where ADT outweighs RNA[1] 74.8223ggplot(sobj_filt@meta.data, aes(x = RNA.weight, fill = singler_label_top)) +
geom_histogram(bins = 40) + geom_vline(xintercept = 0.5, linetype = "dashed") +
facet_wrap(~singler_label_top)
GRÁFICO / SALIDA: Histograma de RNA.weight facetado por tipo celular, con una línea discontinua en 0.5. Las células T alcanzan su pico a la izquierda de la línea, los monocitos a la derecha.
Un cluster donde la mayoría de las células tienen RNA.weight cerca de 0 se identifica casi enteramente por proteína. Eso puede ser biología (subtipos de células T, donde el dropout de RNA es severo) o un artefacto (un efecto de lote específico de proteína concentrado en un cluster).
rna_w_by_cluster <- sobj_filt@meta.data %>%
group_by(seurat_clusters) %>%
summarise(median_rna_w = median(RNA.weight),
mean_rna_w = mean(RNA.weight),
n_cells = n(),
.groups = "drop") %>%
arrange(median_rna_w)
cat("RNA.weight per cluster (sorted; low = protein-driven):\n")RNA.weight per cluster (sorted; low = protein-driven):print(rna_w_by_cluster)# A tibble: 9 × 4
seurat_clusters median_rna_w mean_rna_w n_cells
<fct> <dbl> <dbl> <int>
1 4 0.142 0.158 225
2 1 0.398 0.402 583
3 0 0.435 0.445 701
4 5 0.439 0.431 166
5 3 0.467 0.468 305
6 6 0.476 0.476 134
7 2 0.484 0.490 340
8 7 0.491 0.487 120
9 8 0.775 0.712 99Marca los clusters donde la mediana de RNA.weight < 0.3 (la proteína domina)
protein_driven <- rna_w_by_cluster$seurat_clusters[rna_w_by_cluster$median_rna_w < 0.3]
if (length(protein_driven) > 0) {
cat("\nClusters where protein dominates (median RNA.weight < 0.3):\n")
cat(" ", paste(protein_driven, collapse = ", "), "\n")
cat("Inspect these clusters: are they T-cell subsets (expected) or",
"something else?\n")
}
Clusters where protein dominates (median RNA.weight < 0.3):
4
Inspect these clusters: are they T-cell subsets (expected) or something else?Etiquetas de SingleR similares a célula T (cualquier caso)
t_cell_pattern <- "T[._ -]?cell|CD4|CD8|Tcell|T cell"Cluster más impulsado por proteína = menor mediana de RNA.weight
md <- sobj_filt@meta.data
cluster_summary <- md %>%
group_by(seurat_clusters) %>%
summarise(median_rna_w = median(RNA.weight, na.rm = TRUE),
n_cells = n(),
t_cell_frac = mean(grepl(t_cell_pattern, singler_label,
ignore.case = TRUE), na.rm = TRUE),
.groups = "drop") %>%
arrange(median_rna_w)
cat("\nClusters sorted by median RNA.weight (lowest = most protein-driven):\n")
Clusters sorted by median RNA.weight (lowest = most protein-driven):print(cluster_summary)# A tibble: 9 × 4
seurat_clusters median_rna_w n_cells t_cell_frac
<fct> <dbl> <int> <dbl>
1 4 0.142 225 0.276
2 1 0.398 583 0.336
3 0 0.435 701 0.287
4 5 0.439 166 0.241
5 3 0.467 305 0.272
6 6 0.476 134 0.351
7 2 0.484 340 0.3
8 7 0.491 120 0.217
9 8 0.775 99 0.253target_cluster <- cluster_summary$seurat_clusters[1]
tfrac <- cluster_summary$t_cell_frac[1]
mrw <- cluster_summary$median_rna_w[1]
cat(sprintf("\nMost protein-driven cluster: %s (median RNA.weight = %.2f)\n",
target_cluster, mrw))
Most protein-driven cluster: 4 (median RNA.weight = 0.14)cat(sprintf("Fraction of its cells with T-cell label: %.1f%%\n", tfrac * 100))Fraction of its cells with T-cell label: 27.6%cat("\nDecision rule (concrete thresholds):\n")
Decision rule (concrete thresholds):cat(" t_cell_frac >= 70% -> protein dominance reflects RNA dropout (expected,\n") t_cell_frac >= 70% -> protein dominance reflects RNA dropout (expected,cat(" CD4/CD8 RNA drops out heavily). No action.\n") CD4/CD8 RNA drops out heavily). No action.cat(" t_cell_frac 30-70% -> mixed. Inspect ADT panel for this cluster: is one\n") t_cell_frac 30-70% -> mixed. Inspect ADT panel for this cluster: is onecat(" specific protein driving the dominance? Run:\n") specific protein driving the dominance? Run:cat(" FeaturePlot(sobj_filt, '<protein>',\n") FeaturePlot(sobj_filt, '<protein>',cat(" cells = WhichCells(sobj_filt, idents = '<cluster>'))\n") cells = WhichCells(sobj_filt, idents = '<cluster>'))cat(" t_cell_frac < 30% -> technical explanation. Check:\n") t_cell_frac < 30% -> technical explanation. Check:cat(" 1. ADT batch (Block 6.3, panel QC)\n") 1. ADT batch (Block 6.3, panel QC)cat(" 2. CLR margin (Block 6.2; verify row means != 0)\n") 2. CLR margin (Block 6.2; verify row means != 0)cat(" 3. Isotype background (Block 8 Scenario 4)\n") 3. Isotype background (Block 8 Scenario 4)verdict <- if (tfrac >= 0.70) {
"RNA dropout (expected)"
} else if (tfrac >= 0.30) {
"mixed; inspect ADT panel"
} else {
"technical artifact; check batch/CLR/isotype"
}
cat(sprintf("\nVerdict for cluster %s: %s\n", target_cluster, verdict))
Verdict for cluster 4: technical artifact; check batch/CLR/isotypeEjecuta un clustering basado en grafos sobre el grafo de WNN y compáralo con la anotación de SingleR solo-RNA, restringido a las células donde ADT dominó.
sobj_filt <- FindClusters(sobj_filt,
graph.name = "wsnn",
resolution = 0.5,
verbose = FALSE)
sobj_filt$wnn_clusters <- Idents(sobj_filt)Células donde ADT dominó en WNN (RNA.weight < 0.3)
low_rna_mask <- sobj_filt$RNA.weight < 0.3
n_low <- sum(low_rna_mask)
cat("Cells where ADT dominated WNN (RNA.weight < 0.3):", n_low,
sprintf("(%.1f%% of total)\n", 100 * n_low / ncol(sobj_filt)))Cells where ADT dominated WNN (RNA.weight < 0.3): 301 (11.3% of total)if (n_low >= 20) {
cat("\nCross-tab: RNA-based SingleR label vs WNN cluster (ADT-dominated cells)\n")
print(table(
SingleR_label = sobj_filt$singler_label[low_rna_mask],
WNN_cluster = sobj_filt$wnn_clusters[low_rna_mask]
))
}
Cross-tab: RNA-based SingleR label vs WNN cluster (ADT-dominated cells)
WNN_cluster
SingleR_label 0 1 2 3 4 5 6 7 8
B_cell 0 11 0 0 28 1 0 0 1
DC 1 4 0 0 8 1 0 0 0
Endothelial_cells 2 4 0 0 10 0 0 0 0
Epithelial_cells 0 2 0 0 6 0 0 0 0
Gametocytes 0 0 0 0 1 0 0 0 0
Hepatocytes 0 0 0 0 2 0 0 0 0
Macrophage 1 0 0 0 8 0 0 0 1
Monocyte 2 15 0 0 26 4 0 0 1
Myelocyte 0 0 0 0 3 0 0 0 0
Neurons 0 0 0 0 2 0 0 0 0
Neutrophils 0 1 0 0 12 0 0 0 3
NK_cell 3 13 0 0 24 1 0 0 1
Platelets 0 1 0 0 5 0 0 0 0
Pre-B_cell_CD34- 0 2 0 0 10 0 0 0 0
T_cells 5 18 0 0 53 4 0 0 0Este gráfico parece informativo. Pero con 3 donantes saludables y 4 con COVID-19, puede estar completamente impulsado por un donante atípico.
group_props <- sobj_filt@meta.data %>%
filter(!is.na(singler_label_top), !is.na(condition)) %>%
group_by(condition, singler_label_top) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(condition) %>%
mutate(prop = n / sum(n))
ggplot(group_props, aes(x = condition, y = prop, fill = singler_label_top)) +
geom_bar(stat = "identity", width = 0.6) +
scale_y_continuous(labels = percent_format()) +
labs(title = "Cell type proportions by condition",
subtitle = "Looks clear at the group level; check whether one donor drives it",
x = "Condition", y = "Proportion", fill = "Cell type") +
theme_classic(base_size = 12) +
theme(legend.text = element_text(size = 8))
donor_props <- sobj_filt@meta.data %>%
filter(!is.na(singler_label_top)) %>%
group_by(donor_id, condition, singler_label_top) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(donor_id) %>%
mutate(prop = n / sum(n))
ggplot(donor_props, aes(x = donor_id, y = prop, fill = singler_label_top)) +
geom_bar(stat = "identity", width = 0.7) +
scale_y_continuous(labels = percent_format()) +
facet_wrap(~condition, scales = "free_x") +
labs(title = "Cell type proportions per donor",
subtitle = "Consistent pattern across donors within each group; not driven by one outlier",
x = NULL, y = "Proportion (%)", fill = "Cell type") +
theme_classic(base_size = 11) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.text = element_text(size = 7),
strip.text = element_text(face = "bold"))
Proporciones por donante para UN tipo celular de ejemplo, probadas entre condiciones. Ejecuta esto para la etiqueta que tu anotación produjo como la más grande en COVID-19; el script intenta elegir automáticamente.
Proporciones por donante para UN tipo celular de ejemplo, probadas entre condiciones. Ejecuta esto para la etiqueta que tu anotación produjo como la más grande en COVID-19; el script intenta elegir automáticamente. Restringido a las etiquetas top-N: una etiqueta sostenida por 1-2 células podría mostrar un delta espurio del 100% entre condiciones solo por azar, lo que haría que un “hallazgo” sin sentido pareciera la señal más fuerte del conjunto de datos.
per_donor <- sobj_filt@meta.data %>%
filter(!is.na(condition), !is.na(singler_label_top), condition %in% c("Healthy", "COVID19")) %>%
group_by(donor_id, condition, singler_label_top) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(donor_id, condition) %>%
mutate(prop = n / sum(n)) %>%
ungroup()Elige el tipo celular con el delta más grande entre grupos para fines de demostración
delta_by_lbl <- per_donor %>%
group_by(singler_label_top, condition) %>%
summarise(mean_prop = mean(prop), .groups = "drop") %>%
tidyr::pivot_wider(names_from = condition, values_from = mean_prop) %>%
mutate(delta = abs(COVID19 - Healthy)) %>%
arrange(desc(delta))
target_lbl <- delta_by_lbl$singler_label_top[1]
cat("Largest mean-proportion delta is for label:", as.character(target_lbl), "\n")Largest mean-proportion delta is for label: T_cells cat(sprintf(" Healthy mean: %.1f%% COVID-19 mean: %.1f%% delta: %.1f%%\n",
100 * delta_by_lbl$Healthy[1],
100 * delta_by_lbl$COVID19[1],
100 * delta_by_lbl$delta[1])) Healthy mean: 39.1% COVID-19 mean: 24.4% delta: 14.7%subset_df <- per_donor %>% filter(singler_label_top == target_lbl)
n_h <- sum(subset_df$condition == "Healthy")
n_c <- sum(subset_df$condition == "COVID19")
if (n_h >= 2 && n_c >= 2) {
wt <- wilcox.test(prop ~ condition, data = subset_df, exact = FALSE)
cat(sprintf("\nWilcoxon rank-sum (n=%d Healthy vs n=%d COVID-19):\n", n_h, n_c))
cat(sprintf(" W = %g, p = %.3f\n", wt$statistic, wt$p.value))
}
Wilcoxon rank-sum (n=3 Healthy vs n=4 COVID-19):
W = 0, p = 0.052Piso de potencia: con n=3 vs n=4 y Wilcoxon a alpha=0.05, el tamaño de efecto estandarizado más pequeño detectable (d de Cohen) es aproximadamente 1.5 para un 80% de potencia.
En términos simples: solo diferencias de proporción mayores a ~1.5 desviaciones estándar de la distribución a nivel de donante son detectables. Cualquier valor p aquí es exploratorio, no confirmatorio.
cat("\n--- Power note ---\n")
--- Power note ---cat("n=3 vs n=4 donors. Minimum detectable effect size (Cohen's d) ~1.5 for\n")n=3 vs n=4 donors. Minimum detectable effect size (Cohen's d) ~1.5 forcat("80% power at alpha=0.05. Treat any p-value here as exploratory only.\n")80% power at alpha=0.05. Treat any p-value here as exploratory only.Vista final de las proporciones anotadas y el registro de la sesión de R. La tabla de referencia de errores comunes vive en el PDF de la guía del instructor.
Cuatro resultados que en conjunto resumen el curso, reutilizando objetos ya construidos en los Bloques 5 y 6 en lugar de recalcular nada:
p_final_umap <- DimPlot(sobj_filt, reduction = "wnn.umap",
group.by = "singler_label_clean",
label = TRUE, label.size = 2.5, repel = TRUE) +
ggtitle("1. Final WNN UMAP (confident annotation)") + NoLegend()
p_final_markers <- DotPlot(sobj_filt,
features = unique(top_markers$gene),
idents = setdiff(levels(Idents(sobj_filt)), "Ambiguous"),
group.by = "singler_label_clean") +
RotatedAxis() +
theme(axis.text.x = element_text(size = 7), legend.position = "none") +
ggtitle("2. Canonical markers behind the labels")
if ("singler_label_top" %in% colnames(sobj_filt@meta.data)) {
final_props <- sobj_filt@meta.data %>%
filter(!is.na(condition), !is.na(singler_label_top)) %>%
group_by(condition, singler_label_top) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(condition) %>%
mutate(prop = n / sum(n))
p_final_props <- ggplot(final_props, aes(condition, prop, fill = singler_label_top)) +
geom_bar(stat = "identity", width = 0.65) +
scale_y_continuous(labels = percent_format()) +
labs(x = "Condition", y = "Proportion", fill = NULL) +
theme_classic(base_size = 10) +
theme(legend.text = element_text(size = 6),
legend.position = "right",
axis.text.x = element_text(angle = 45, hjust = 1)) +
ggtitle("3. Cell type proportions by condition")
}
DefaultAssay(sobj_filt) <- "RNA"
p_cd4_rna <- FeaturePlot(sobj_filt, "CD4", min.cutoff = "q05") +
ggtitle("4a. CD4 (RNA)") +
theme(plot.title = element_text(size = 10), legend.position = "none")
DefaultAssay(sobj_filt) <- "ADT"
p_cd4_adt <- FeaturePlot(sobj_filt, "CD4", min.cutoff = "q05") +
ggtitle("4b. CD4 (ADT)") +
theme(plot.title = element_text(size = 10), legend.position = "none")
DefaultAssay(sobj_filt) <- "RNA"
p_final_cd4 <- p_cd4_rna | p_cd4_adt
(p_final_umap | p_final_markers) / (p_final_props | p_final_cd4)
sessionInfo()R version 4.6.0 (2026-04-24)
Platform: aarch64-apple-darwin23
Running under: macOS Tahoe 26.5.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Mexico_City
tzcode source: internal
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] future_1.70.0 here_1.0.2
[3] ggrepel_0.9.8 scales_1.4.0
[5] harmony_2.0.5 Rcpp_1.1.1-1.1
[7] BiocParallel_1.46.0 SingleR_2.14.0
[9] SummarizedExperiment_1.42.0 Biobase_2.72.0
[11] GenomicRanges_1.64.0 Seqinfo_1.2.0
[13] IRanges_2.46.0 S4Vectors_0.50.1
[15] BiocGenerics_0.58.1 generics_0.1.4
[17] MatrixGenerics_1.24.0 matrixStats_1.5.0
[19] Matrix_1.7-5 tidyr_1.3.2
[21] dplyr_1.2.1 patchwork_1.3.2
[23] ggplot2_4.0.3 Seurat_5.5.1
[25] SeuratObject_5.4.0 sp_2.2-1
loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 rstudioapi_0.19.0 jsonlite_2.0.0
[4] magrittr_2.0.5 ggbeeswarm_0.7.3 spatstat.utils_3.2-3
[7] farver_2.1.2 rmarkdown_2.31 vctrs_0.7.3
[10] ROCR_1.0-12 spatstat.explore_3.8-1 htmltools_0.5.9
[13] S4Arrays_1.12.0 SparseArray_1.12.2 sctransform_0.4.3
[16] parallelly_1.48.0 KernSmooth_2.23-26 htmlwidgets_1.6.4
[19] ica_1.0-3 plyr_1.8.9 plotly_4.12.0
[22] zoo_1.8-15 igraph_2.3.3 mime_0.13
[25] lifecycle_1.0.5 pkgconfig_2.0.3 R6_2.6.1
[28] fastmap_1.2.0 fitdistrplus_1.2-6 shiny_1.14.0
[31] digest_0.6.39 rprojroot_2.1.1 tensor_1.5.1
[34] RSpectra_0.16-2 irlba_2.3.7 beachmat_2.28.0
[37] labeling_0.4.3 progressr_0.19.0 spatstat.sparse_3.2-0
[40] httr_1.4.8 polyclip_1.10-7 abind_1.4-8
[43] compiler_4.6.0 withr_3.0.3 S7_0.2.2
[46] fastDummies_1.7.6 MASS_7.3-65 DelayedArray_0.38.2
[49] tools_4.6.0 vipor_0.4.7 lmtest_0.9-40
[52] otel_0.2.0 beeswarm_0.4.0 httpuv_1.6.17
[55] future.apply_1.20.2 goftest_1.2-3 glue_1.8.1
[58] nlme_3.1-169 promises_1.5.0 grid_4.6.0
[61] Rtsne_0.17 cluster_2.1.8.2 reshape2_1.4.5
[64] isoband_0.3.0 gtable_0.3.6 spatstat.data_3.1-9
[67] data.table_1.18.4 utf8_1.2.6 XVector_0.52.0
[70] spatstat.geom_3.8-1 RcppAnnoy_0.0.23 RANN_2.6.2
[73] pillar_1.11.1 stringr_1.6.0 limma_3.68.4
[76] spam_2.11-4 RcppHNSW_0.7.0 later_1.4.8
[79] splines_4.6.0 lattice_0.22-9 survival_3.8-6
[82] deldir_2.0-4 tidyselect_1.2.1 miniUI_0.1.2
[85] pbapply_1.7-4 knitr_1.51 gridExtra_2.3.1
[88] scattermore_1.2 xfun_0.59 statmod_1.5.2
[91] stringi_1.8.7 lazyeval_0.2.3 evaluate_1.0.5
[94] codetools_0.2-20 tibble_3.3.1 cli_3.6.6
[97] uwot_0.2.4 xtable_1.8-8 reticulate_1.46.0
[100] globals_0.19.1 spatstat.random_3.5-0 png_0.1-9
[103] ggrastr_1.0.2 spatstat.univar_3.2-0 parallel_4.6.0
[106] dotCall64_1.2 listenv_1.0.0 viridisLite_0.4.3
[109] ggridges_0.5.7 purrr_1.2.2 rlang_1.2.0
[112] cowplot_1.2.0