Objetivo: Trazer a camada de proteína para a análise. O fluxo de trabalho de RNA até o Bloco 5 ignorou completamente o ADT. O valor agregado do CITE-seq é a capacidade de resolver o dropout, validar a anotação e ponderar as modalidades por célula.
Ordem das operações: inspecionar a camada de ADT e contrastá-la com o RNA, normalizar, fazer o controle de qualidade (QC) do painel, contrastar RNA vs ADT para marcadores correspondentes, fazer o gating, verificar canais trocados, construir o PCA de ADT, e então integrar com WNN.
Carregar o arquivo no R:
library(here) # install.packages("here")here() starts at /private/var/folders/zz/cdcfnwts145c0xjqs_c0y6mm0000gn/T/RtmppqY2dt/file177bf72942b49sobj_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 Resumo por proteína: escala e distribuição:
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.8Contagem 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 saudável: mediana próxima de 0, média algumas unidades mais alta, máximo nas centenas. A proteína não está corada na maioria das células e é fortemente positiva em seu subconjunto alvo. Padrão de anticorpo morto: média abaixo de 1, máximo também baixo (abaixo de 10), pct_zero muito alto. Nada sobe acima do fundo porque o anticorpo não é funcional. No conjunto de dados injetado, CD56 mostra o padrão morto; CD3, CD4, CD8a, CD14, CD19 mostram o padrão bimodal. Examine a tabela ordenada por média para detectar a anomalia.
DEMONSTRAÇÃO AO 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 / SAÍDA:
Também poderia usar:
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
O Bloco 1 (Seção 4.4) mediu a esparsidade (sparsity) do RNA isoladamente. Agora que a camada de ADT também foi inspecionada, a comparação é o ponto principal: as mesmas células, dois assays, taxas de zero muito 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 %O ADT permaneceu intocado desde o Bloco 1 (Seção 4.4). Confirme seu estado explicitamente 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 Neste ponto esperamos: apenas counts. Sem data, sem scale.data. A normalização é a próxima etapa.
Normalização CLR (Centered Log-Ratio, log-razão centralizada). O argumento margin define a direção: + margin = 1: normaliza cada proteína em todas as células (médias de linha ~0) + margin = 2: normaliza cada célula em todas as proteínas (médias de coluna ~0)
margin = 2 é o correto para CITE-seq. Ele remove a variação por célula na captura total de anticorpos, análogo à normalização por tamanho de biblioteca em RNA. margin = 1 roda sem erro, mas destrói o sinal biológico. Você vai diagnosticar um objeto real normalizado da forma errada no Bloco 8, Cenário 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.Uma conjugação falha ou um anticorpo degradado produz um “canal morto”: a proteína lê perto de zero em todas as células, com quase nenhuma variância. Não gera erro. Se você fizer gating ou anotar sobre um canal morto, perde silenciosamente uma população inteira.
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"No conjunto de dados injetado: CD56 (codificado por NCAM1). CD56 marca as células NK. Com CD56 ADT morto, qualquer anotação de células NK baseada apenas em ADT falha. O bloco de diagnóstico compara CD56 ADT (quase zero) com NCAM1 RNA (sinal claro no cluster de NK) e expõe a inconsistência. Hábito: sempre verifique cruzadamente as proteínas de ADT com menor média contra suas contrapartes de RNA.
DEMONSTRAÇÃO AO 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 / SAÍDA: Painel FeaturePlot: CD56 ADT (plano) ao lado de NCAM1 RNA (positivo em um cluster)
Compare a proteína suspeita (ADT) com seu gene de RNA lado a lado.
Edite suspect_adt / suspect_rna para a proteína sinalizada acima.
suspect_adt <- "CD56" # protein flagged as near-zero in panel QC
suspect_rna <- "NCAM1" # its RNA counterpart (NK marker)Decida o título a partir dos números reais de panel_qc em vez de afirmar uma afirmação fixa. No objeto limpo, esta proteína está saudável; no objeto injetado, está morta. O título deve dizer o que for verdadeiro para o objeto que está realmente carregado.
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"Comparação direta. O mRNA de CD4 lê como zero em muitas células T CD4+ (dropout), enquanto a proteína CD4 do ADT é bimodal nas mesmas 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 com uma contraparte de RNA, calcule a correlação de Spearman entre todas as células. Correlação baixa = dropout alto em 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, encontre as células que são claramente positivas para a proteína (ADT > limiar) mas têm zero contagens de RNA. Essas são as células que seriam mal-anotadas por uma análise baseada apenas em 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%O gating digital replica os gráficos de dispersão biaxiais da citometria de fluxo usando dados de ADT. Ele permite a classificação computacional de células com a mesma lógica que os imunologistas usam na bancada.
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"Uma chamada de função que “funciona” no assay errado produz um erro silencioso.
FeaturePlot("CD4") com DefaultAssay = "RNA" plota o gene CD4 de RNA.FeaturePlot("CD4") com DefaultAssay = "ADT" plota a proteína CD4.O gráfico é RENDERIZADO em ambos os casos, e os dois gráficos podem parecer enganosamente semelhantes: min.cutoff = "q05" reescala o gradiente de cor de cada painel para seu próprio intervalo de valores, então um sinal de RNA esparso e um sinal de ADT denso acabam ambos esticados por uma escala de cor de aparência similar. A falha é silenciosa precisamente porque nada no gráfico renderizado avisa qual assay o produziu.
Quantifique a diferença que o gráfico por si só oculta: qual fração de células lê exatamente zero em cada modalidade 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.Demonstre desenhando a mesma chamada sob ambos os assays. Os títulos são mantidos curtos e o tamanho da fonte é restrito, já que os dois gráficos são renderizados lado a lado com metade da largura cada.
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"Você teria qualquer assay que estivesse ativo no momento da chamada do gráfico, e um leitor não conseguiria distinguir qual apenas pela figura: o título do gráfico diz apenas ‘CD4’, os eixos são coordenadas UMAP idênticas, e min.cutoff='q05' reescala o gradiente de cor próprio de cada painel de forma independente, então um sinal de RNA esparso e um sinal de ADT mais denso acabam esticados por escalas de aparência visualmente similar. É exatamente por isso que os números de taxa de zero importam: eles mostram a diferença real em qual fração de células lê como positiva em cada modalidade, uma diferença que o gradiente de cor reescalado oculta. A figura publicada poderia mostrar silenciosamente CD4 RNA quando CD4 proteína era pretendida, ou vice-versa, e nada na própria imagem sinalizaria a discrepância. Hábito: defina DefaultAssay explicitamente no topo de qualquer bloco de plotagem, E/OU passe o assay dentro da chamada; melhor ainda, edite o título do gráfico para incluir o nome do assay e considere relatar a taxa de zero junto com a figura para que a diferença de modalidade fique documentada, não apenas visível para quem já sabe procurar por isso.
DEMONSTRAÇÃO AO 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 / SAÍDA: Dois FeaturePlots lado a lado, um por assay, com títulos distintos; os números de taxa de zero explicam o que os gráficos sozinhos não mostram.
Leia o código abaixo. Ele parece correto. O que está errado?
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 sessionSua resposta:
REVELAR: a chamada pede marcadores de células T CD4, mas DefaultAssay é RNA. Ela retorna marcadores de RNA, não marcadores de ADT. A saída está correta para RNA, mas não responde à pergunta “quais proteínas distinguem as células T CD4”. Ou defina DefaultAssay = "ADT" primeiro, ou passe assay = "ADT" dentro da chamada FindMarkers. O min.pct = 0.25 também está errado para ADT (projetado para RNA esparso); para ADT use um valor como 0.5 ou 0.0 com logfc.threshold.
Um canal trocado ou mal rotulado roda perfeitamente. O CLR o normaliza, o FeaturePlot o desenha, o gating o agrupa. Nada dá erro. O único sinal é que a biologia é impossível. A defesa é verificar cruzadamente cada proteína contra um rótulo independente (seu marcador de RNA).
Cada proteína deveria se correlacionar com seu próprio gene de RNA entre as células.
Um marcador canônico com rho próximo de 0 é um canal mal rotulado.
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 Dois canais de ADT mostram rho próximo de 0 contra seu próprio gene de RNA, enquanto o par cruzado é alto. Confirme a troca com a visualização abaixo.
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 e sua contraparte de RNA entre as células. Marcadores saudáveis mostram rho de 0.3-0.7. Impressão digital da troca: ADT_X vs RNA_X próximo de zero E ADT_X vs RNA_Y alto; o mesmo com X,Y invertidos. No conjunto de dados injetado, os rownames de CD14 e CD19 em ADT estão trocados. O bloco de correção reconstrói o assay ADT com os rownames re-trocados e reexecuta o CLR. É condicional, então em dados limpos não faz nada.
DEMONSTRAÇÃO AO 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 / SAÍDA: Dois gráficos de dispersão lado a lado. Saudável: os dispersogramas de ADT e RNA coincidem. Com uma troca: eles divergem
Seguro para v5: reconstrói o assay ADT a partir de sua matriz de contagens com os nomes corrigidos.Condicional: só age se a troca for de fato detectada (seguro em dados limpos).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"A integração WNN (Weighted Nearest Neighbor, vizinho mais próximo ponderado) aprende a contribuição relativa de cada modalidade por célula. Células onde o RNA é informativo recebem mais peso de RNA; células onde a proteína é mais discriminante recebem mais peso de ADT. O WNN requer um PCA por modalidade. O PCA de RNA existe desde o Bloco 4; agora construímos o 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
)Com 24 proteínas, solicitar 20 PCs dispara um aviso do Seurat sobre calcular valores singulares demais via SVD truncado. Ambas as soluções funcionam:
Usamos (a) porque 15 PCs é sinal suficiente 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 O Seurat espera um nome de peso por modalidade. Passar apenas “RNA.weight” dispara um aviso de que ADT.weight foi auto-atribuído. Passe ambos explicitamente para silenciar o aviso e documentar os nomes de coluna pretendidos.
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.umapO WNN não necessariamente move todos os tipos celulares para uma região visualmente diferente; os layouts de UMAP podem rotacionar ou espelhar entre execuções mesmo quando as vizinhanças subjacentes mal mudam. Os dois painéis abaixo podem parecer semelhantes à primeira vista. O que realmente importa é se o conjunto de células agrupadas mudou, não se a imagem parece diferente.
A Etapa 6.11c (mais adiante) verifica isso diretamente por tipo celular através de uma tabela cruzada; aqui, uma comparação rápida de clustering dá uma primeira leitura antes disso.
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 82Uma pontuação de concordância simples: para cada cluster apenas-RNA, qual fração de suas células cai no único cluster de WNN que captura a maioria delas.
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 próximo de 1 = a célula é melhor caracterizada por RNARNA.weight próximo de 0 = a célula é melhor caracterizada por ADTRNA.weight e ADT.weight somam 1 para cada célula (FindMultiModalNeighbors os normaliza dessa forma), então 0.5 é o ponto de referência natural: abaixo dele, o ADT está contribuindo mais do que o RNA para a colocação daquela célula no grafo WNN, não apenas “alguma quantidade” de informação 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))
Células T CD4 e células T CD8, por uma ampla margem. Seus marcadores definidores (CD4, CD8A) têm o pior dropout de RNA. Células B (MS4A1, CD79A) e monócitos (CD14, LYZ) têm sinal de RNA mais forte e dependem menos do ADT. Células NK (NKG7, GNLY, NCAM1) são intermediárias. O histograma facetado de RNA.weight por tipo celular torna isso concreto, e a linha tracejada em RNA.weight = 0.5 marca o ponto em que o ADT supera o RNA para aquela célula, não apenas contribui com alguma quantidade de informação.
DEMONSTRAÇÃO AO 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 / SAÍDA: Histograma de RNA.weight facetado por tipo celular, com uma linha tracejada em 0.5. Células T atingem o pico à esquerda da linha, monócitos à direita.
Um cluster onde a maioria das células tem RNA.weight próximo de 0 é identificado quase inteiramente por proteína. Isso pode ser biologia (subtipos de células T, onde o dropout de RNA é severo) ou um artefato (um efeito de lote específico de proteína concentrado em um 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 99Marque os clusters onde a mediana de RNA.weight < 0.3 (a 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?Rótulos de SingleR semelhantes a célula T (qualquer caso)
t_cell_pattern <- "T[._ -]?cell|CD4|CD8|Tcell|T cell"Cluster mais impulsionado 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/isotypeExecute um clustering baseado em grafo sobre o grafo WNN e compare-o com a anotação de SingleR apenas-RNA, restrito às células onde o ADT dominou.
sobj_filt <- FindClusters(sobj_filt,
graph.name = "wsnn",
resolution = 0.5,
verbose = FALSE)
sobj_filt$wnn_clusters <- Idents(sobj_filt)Células onde o ADT dominou o 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. Mas com 3 doadores saudáveis e 4 com COVID-19, ele pode ser completamente impulsionado por um doador 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"))
Proporções por doador para UM tipo celular de exemplo, testadas entre condições. Execute isso para o rótulo que sua anotação produziu como o maior em COVID-19; o script tenta escolher automaticamente.
Proporções por doador para UM tipo celular de exemplo, testadas entre condições. Execute isso para o rótulo que sua anotação produziu como o maior em COVID-19; o script tenta escolher automaticamente. Restrito aos rótulos top-N: um rótulo sustentado por 1-2 células poderia mostrar um delta espúrio de 100% entre condições apenas por acaso, o que faria um “achado” sem sentido parecer o sinal mais forte do conjunto de dados.
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()Escolha o tipo celular com o maior delta entre grupos para fins de demonstração
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 poder: com n=3 vs n=4 e Wilcoxon em alpha=0.05, o menor tamanho de efeito padronizado detectável (d de Cohen) é aproximadamente 1.5 para 80% de poder.
Em termos simples: apenas diferenças de proporção maiores que ~1.5 desvios-padrão da distribuição em nível de doador são detectáveis. Qualquer valor-p aqui é exploratório, não confirmatório.
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.Visão final das proporções anotadas e o registro da sessão de R. A tabela de referência de armadilhas comuns vive no PDF do guia do instrutor.
Quatro resultados que juntos resumem o curso, reutilizando objetos já construídos nos Blocos 5 e 6 em vez de recalcular qualquer coisa:
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