4  scRNA‑seq — Errores y Supervivencia

ImportanteSobre los conjuntos de datos

Este conjunto de datos proviene de GSE149689. Los datos han sido reducidos y se introdujeron errores deliberados con fines educativos.

Errores deliberados

Se diseñaron cinco pasos para que fallaran (por ejemplo, acceso al slot de v4, uso incorrecto de GetAssayData slot=, umbrales de filtrado incorrectos, dimensiones mayores al número de componentes principales, graficar una etiqueta antes de agregarla).

Estos pasos están envueltos en try(), de modo que el error se imprime en la consola sin detener la ejecución completa. El objetivo es que el usuario lea el mensaje de error, lo interprete y continúe con el análisis.

4.1 Opciones de directorio de trabajo

Al iniciar un análisis, es necesario asegurarse de que R sepa dónde encontrar tus scripts y datos. Hay dos enfoques comunes:

• Establecer el directorio de trabajo Puedes establecer manualmente el directorio de trabajo en RStudio mediante: Session > Set Working Directory > To Source File Location. Esta opción hace que R se ejecute desde la carpeta que contiene tu script y la subcarpeta data/. Es rápido y útil para scripts pequeños o puntuales, pero requiere reconfigurarse cada vez que abres el proyecto.

• Crear un Proyecto de R El enfoque recomendado para la reproducibilidad y la colaboración. Un archivo .Rproj establece automáticamente la raíz del proyecto como directorio de trabajo. Esto permite usar rutas relativas (por ejemplo, data/file.csv) sin ajustes manuales. También se integra perfectamente con Git/GitHub, Quarto y RMarkdown, haciendo que tu flujo de trabajo sea más organizado y consistente.

TipBuena práctica

Para proyectos de largo plazo, especialmente aquellos compartidos en GitHub o usados en la enseñanza, crear un Proyecto de R es la opción más confiable.

4.2 Cómo usar este documento

Ejecuta el código de forma interactiva, sección por sección (Ctrl+Enter / Cmd+Enter), de principio a fin. NO uses source() para todo el archivo de una sola vez: varios bloques son errores deliberados pensados para ser leídos, y el Bloque 8 es trabajo opcional para casa que depende de un archivo que se publica por separado después de la sesión.

4.3 Bloque 0 - Instalación de paquetes y configuración

4.3.1 Paso 0.1 - Instalar paquetes

Verifica cada paquete antes de instalarlo. Si install.packages() falla, reintenta a través de BiocManager. Es seguro ejecutarlo varias veces.

# Install BiocManager
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

# All packages
all_pkgs <- c(
  "Seurat", "SeuratObject", "ggplot2", "patchwork",  "dplyr", "tidyr", "Matrix", "harmony", "scales", "ggrepel",
  "SingleR", "celldex", "BiocParallel", "scDblFinder",
  "SingleCellExperiment"
)

for (pkg in all_pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    message("Installing: ", pkg)
    tryCatch(install.packages(pkg, quiet = TRUE),
             error   = function(e) NULL,
             warning = function(w) NULL)
    if (!requireNamespace(pkg, quietly = TRUE)) {
      message("  CRAN failed. Retrying via BiocManager...")
      tryCatch(
        BiocManager::install(pkg, ask = FALSE, update = FALSE),
        error   = function(e) message("  FAILED: ", e$message),
        warning = function(w) NULL
      )
    }
    if (requireNamespace(pkg, quietly = TRUE))
      message("  OK: ", pkg, " (", as.character(packageVersion(pkg)), ")")
    else
      message("  FAILED: ", pkg, ". Run the fallback chunk below.")
  } else {
    message("Already installed: ", pkg,
            " (", as.character(packageVersion(pkg)), ")")
  }
}

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Already installed: Seurat (5.5.0)
Already installed: SeuratObject (5.4.0)
Already installed: ggplot2 (4.0.3)
Already installed: patchwork (1.3.2)
Already installed: dplyr (1.2.1)
Already installed: tidyr (1.3.2)
Already installed: Matrix (1.7.5)
Already installed: harmony (2.0.4)
Already installed: scales (1.4.0)
Already installed: ggrepel (0.9.8)
Installing: SingleR
  CRAN failed. Retrying via BiocManager...
'getOption("repos")' replaces Bioconductor standard repositories, see
'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
Bioconductor version 3.23 (BiocManager 1.30.27), R 4.6.0 (2026-04-24)
Installing package(s) 'BiocVersion', 'SingleR'
also installing the dependencies ‘XVector’, ‘MatrixGenerics’, ‘GenomicRanges’, ‘Biobase’, ‘IRanges’, ‘Seqinfo’, ‘S4Arrays’, ‘SparseArray’, ‘SummarizedExperiment’, ‘BiocGenerics’, ‘S4Vectors’, ‘DelayedArray’, ‘beachmat’, ‘assorthead’

  FAILED: SingleR. Run the fallback chunk below.
Installing: celldex
  CRAN failed. Retrying via BiocManager...
'getOption("repos")' replaces Bioconductor standard repositories, see
'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
Bioconductor version 3.23 (BiocManager 1.30.27), R 4.6.0 (2026-04-24)
Installing package(s) 'BiocVersion', 'celldex'
also installing the dependencies ‘dir.expiry’, ‘dbplyr’, ‘Biostrings’, ‘XVector’, ‘rhdf5filters’, ‘V8’, ‘biocmake’, ‘h5mread’, ‘MatrixGenerics’, ‘GenomicRanges’, ‘Biobase’, ‘BiocGenerics’, ‘IRanges’, ‘Seqinfo’, ‘S4Arrays’, ‘BiocFileCache’, ‘BiocBaseUtils’, ‘KEGGREST’, ‘SparseArray’, ‘sparseMatrixStats’, ‘alabaster.schemas’, ‘rhdf5’, ‘jsonvalidate’, ‘assorthead’, ‘Rhdf5lib’, ‘HDF5Array’, ‘alabaster.ranges’, ‘blob’, ‘SummarizedExperiment’, ‘ExperimentHub’, ‘AnnotationHub’, ‘AnnotationDbi’, ‘S4Vectors’, ‘DelayedArray’, ‘DelayedMatrixStats’, ‘gypsum’, ‘alabaster.base’, ‘alabaster.matrix’, ‘alabaster.se’, ‘DBI’, ‘RSQLite’

  FAILED: celldex. Run the fallback chunk below.
Installing: BiocParallel
  CRAN failed. Retrying via BiocManager...
'getOption("repos")' replaces Bioconductor standard repositories, see
'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
Bioconductor version 3.23 (BiocManager 1.30.27), R 4.6.0 (2026-04-24)
Installing package(s) 'BiocVersion', 'BiocParallel'
also installing the dependencies ‘formatR’, ‘lambda.r’, ‘futile.options’, ‘futile.logger’, ‘snow’

  FAILED: BiocParallel. Run the fallback chunk below.
Installing: scDblFinder
  CRAN failed. Retrying via BiocManager...
'getOption("repos")' replaces Bioconductor standard repositories, see
'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
Bioconductor version 3.23 (BiocManager 1.30.27), R 4.6.0 (2026-04-24)
Installing package(s) 'BiocVersion', 'scDblFinder'
also installing the dependencies ‘formatR’, ‘lambda.r’, ‘futile.options’, ‘locfit’, ‘beeswarm’, ‘vipor’, ‘Cairo’, ‘cigarillo’, ‘RCurl’, ‘rjson’, ‘futile.logger’, ‘snow’, ‘assorthead’, ‘beachmat’, ‘ScaledMatrix’, ‘rsvd’, ‘MatrixGenerics’, ‘Biobase’, ‘Seqinfo’, ‘S4Arrays’, ‘edgeR’, ‘limma’, ‘statmod’, ‘metapod’, ‘SparseArray’, ‘ggbeeswarm’, ‘viridis’, ‘RcppML’, ‘pheatmap’, ‘ggrastr’, ‘UCSC.utils’, ‘Biostrings’, ‘XVector’, ‘Rhtslib’, ‘XML’, ‘GenomicAlignments’, ‘BiocIO’, ‘restfulr’, ‘SingleCellExperiment’, ‘BiocGenerics’, ‘BiocParallel’, ‘BiocNeighbors’, ‘BiocSingular’, ‘S4Vectors’, ‘SummarizedExperiment’, ‘scran’, ‘scater’, ‘scuttle’, ‘bluster’, ‘DelayedArray’, ‘xgboost’, ‘IRanges’, ‘GenomicRanges’, ‘GenomeInfoDb’, ‘Rsamtools’, ‘rtracklayer’

  FAILED: scDblFinder. Run the fallback chunk below.
Installing: SingleCellExperiment
  CRAN failed. Retrying via BiocManager...
'getOption("repos")' replaces Bioconductor standard repositories, see
'help("repositories", package = "BiocManager")' for details.
Replacement repositories:
    CRAN: https://cran.rstudio.com/
Bioconductor version 3.23 (BiocManager 1.30.27), R 4.6.0 (2026-04-24)
Installing package(s) 'BiocVersion', 'SingleCellExperiment'
also installing the dependencies ‘XVector’, ‘MatrixGenerics’, ‘Biobase’, ‘IRanges’, ‘Seqinfo’, ‘S4Arrays’, ‘SparseArray’, ‘SummarizedExperiment’, ‘S4Vectors’, ‘BiocGenerics’, ‘GenomicRanges’, ‘DelayedArray’

  FAILED: SingleCellExperiment. Run the fallback chunk below.

4.3.2 Paso 0.2 - Verificar y cargar todas las librerías

Todos los paquetes se cargan aquí. No aparecen llamadas a library() más adelante en el script. Si ves “namespace ggplot2 is imported by Seurat…” eso no es un error; el paquete ya está activo. Siempre realiza Session > Restart R antes de abrir el script.

pkgs <- c("Seurat", "SeuratObject", "ggplot2", "patchwork",
          "dplyr", "tidyr", "Matrix", "SingleR", "celldex",
          "harmony", "scales", "ggrepel", "BiocParallel",
          "scDblFinder", "SingleCellExperiment")

for (pkg in pkgs)
  cat(sprintf("  %-18s %s\n", pkg,
              ifelse(requireNamespace(pkg, quietly = TRUE), "OK", "MISSING")))

  Seurat             OK
  SeuratObject       OK
  ggplot2            OK
  patchwork          OK
  dplyr              OK
  tidyr              OK
  Matrix             OK
  SingleR            OK
  celldex            MISSING
  harmony            OK
  scales             OK
  ggrepel            OK
  BiocParallel       OK
  scDblFinder        OK
  SingleCellExperiment OK

cat("\nSeurat version:", as.character(packageVersion("Seurat")), "\n")

Seurat version: 5.5.1 

suppressPackageStartupMessages({
  library(Seurat);   library(SeuratObject); library(ggplot2)
  library(patchwork); library(dplyr);     library(tidyr);       library(Matrix)
  library(SingleR);  library(celldex);      library(BiocParallel)
  library(harmony);  library(scales);       library(ggrepel)
})

Warning: package 'Seurat' was built under R version 4.6.1

Error in `library()`:
! there is no package called 'celldex'

cat("\nAll libraries loaded.\n")

All libraries loaded.

Caso 1: scDblFinder

El comando usado fue:

BiocManager::install("scDblFinder", type="binary", dependencies=TRUE)

👉 En este caso, la instalación fue exitosa porque Bioconductor proporciona compilaciones binarias para macOS ARM64. Al forzar type="binary", R evitó compilar código C++, y todas las dependencias se instalaron sin errores.

Caso 2: celldex

El comando fue:

BiocManager::install("celldex", type="binary", dependencies=TRUE)

👉 Aquí el error persistió porque una dependencia crítica (alabaster.base) aún no tiene una compilación binaria disponible. R intentó compilarla desde el código fuente, pero falló con el siguiente error:

ld: library 'ssl' not found
clang++: error: linker command failed with exit code 1
ERROR: compilation failed for package ‘alabaster.base’
ERROR: dependency ‘alabaster.base’ is not available for package ‘celldex’

Esto significa que el compilador no pudo encontrar la librería OpenSSL requerida para el enlazado. Aunque solicitamos la instalación binaria, Bioconductor aún no ha publicado un binario para esta dependencia, por lo que el error continúa.

4.3.3 Paso 0.3 - Estructura de carpetas y descarga de archivos

Crear y eliminar archivos

# Helper function to recreate a folder with messages
recreate_dir <- function(dir_name) {
  if (dir.exists(dir_name)) {
    cat(sprintf("The folder '%s' already exists. It will be deleted and recreated.\n", dir_name))
    unlink(dir_name, recursive = TRUE)  # delete folder and contents
  } else {
    cat(sprintf("The folder '%s' does not exist. It will be created.\n", dir_name))
  }
  dir.create(dir_name, showWarnings = FALSE)
}

# Apply the function to your folders
recreate_dir("data")

The folder 'data' already exists. It will be deleted and recreated.

recreate_dir("outputs")

The folder 'outputs' already exists. It will be deleted and recreated.

cat("Process completed. Folders are ready.\n")

Process completed. Folders are ready.

Descargar archivos

Aquí encontrarás los archivos reducidos disponibles para descargar desde Google Drive.

FILE_ID <- "1mCg4VD91qK7pg-zL1iYUUBFWBr7fmmCM"
options(timeout = 600)
download.file(
  url      = paste0("https://drive.google.com/uc?export=download&confirm=t&id=",
                    FILE_ID),
  destfile = "data/gse149689_subset_3k.rds",
  mode     = "wb"
)
cat("File exists:", file.exists("data/gse149689_subset_3k.rds"), "\n")

File exists: TRUE 

cat("Size (MB)  :", round(file.size("data/gse149689_subset_3k.rds") / 1e6, 1), "\n")

Size (MB)  : 14.1 

4.3.4 Paso 0.4 - Ayudante reveal() (solo para uso del instructor)

Imprime la respuesta + código de demostración en vivo para cualquier PREGUNTA o ACERTIJO de este script. Las respuestas NO están en este archivo. Se encuentran en un archivo separado instructor_answers.R que solo el instructor carga con source() antes de la clase. Los estudiantes que ejecuten reveal() en una sesión nueva verán un breve aviso y nada más.

👉 Descargar: instructor_answers.R y reveal_function.R

#Instructor only
source("scripts/instructor_answers.R")

Instructor answers loaded. 47 entries. Try: reveal("1.2b")

source("scripts/reveal_function.R")

Uso: reveal("number")

reveal("5.2")

4.3.5 Paso 0.5 - Envoltorios de acceso seguros a espacios de nombres

En R, los nombres que comienzan con un punto (.) generalmente se usan para funciones auxiliares internas. La función .safe_accessor no está pensada para ser llamada directamente por el usuario. En cambio, construye las versiones “seguras” de las funciones de acceso de Seurat: SafeAssays, SafeLayers, y SafeReductions.

Estos envoltorios seguros aseguran que se use la función correcta incluso si diferentes versiones de Seurat u otros paquetes de Bioconductor introducen conflictos. Sin ellos, podrías obtener errores confusos al llamar directamente a Assays(), Layers(), o Reductions().

👉 Descargar: safe_accessor_function.R

Cárgala en tu sesión de R con:

#Instructor only
source("scripts/safe_accessor_function.R")

TipExplicación de cada función

Usa las funciones seguras en tu análisis:

• SafeAssays(seurat_object)

  • Devuelve una lista de todos los ensayos (assays) presentes en el objeto Seurat.
  • Más seguro que llamar directamente a Assays() porque evita problemas si el objeto tiene slots de ensayo inusuales o corruptos.
  • Útil para verificar qué tipos de datos (RNA, ATAC, proteína, etc.) están disponibles antes de ejecutar el análisis posterior.

SafeAssays(seurat_object)

• SafeLayers(seurat_object)

  • Enumera todas las capas (layers) dentro de un ensayo dado (por ejemplo, conteos crudos, datos normalizados, datos escalados).
  • Te ayuda a confirmar qué representaciones de los datos están almacenadas y previene errores al cambiar entre capas.
  • Importante en flujos de trabajo multimodales donde coexisten múltiples capas.

SafeLayers(seurat_object)

• SafeReductions(seurat_object)

  • Muestra todos los resultados de reducción de dimensionalidad (PCA, UMAP, t-SNE, etc.) almacenados en el objeto.
  • Asegura que sepas qué reducciones están disponibles antes de graficar o agrupar (clustering).
  • Evita errores como llamar a DimPlot() sobre una reducción que no existe.

SafeReductions(seurat_object)

4.4 Bloque 1 - Inspección del objeto (RNA + ADT)

Objetivo: Inspeccionar el ensayo RNA, contrastarlo con el ensayo ADT, y recorrer el sistema de capas, los metadatos y el slot de reducciones.

4.4.1 Paso 1.1 - Cargar el objeto

Carga el objeto inyectado (con errores incorporados) cuando está presente; de lo contrario, carga el objeto limpio.

main_candidates <- c("data/gse149689_subset_3k_errors.rds",
                     "data/gse149689_subset_3k.rds")
main_path <- main_candidates[file.exists(main_candidates)][1]
if (is.na(main_path)) stop("No dataset found in data/. Expected one of: ",
                           paste(main_candidates, collapse = ", "))
sobj <- readRDS(main_path)
cat("Loaded:", main_path, "\n")

Loaded: data/gse149689_subset_3k.rds 

print(sobj)

An object of class Seurat 
524 features across 3000 samples within 2 assays 
Active assay: RNA (500 features, 500 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: ADT
 2 dimensional reductions calculated: pca, umap

TipPREGUNTA 1.1

• ¿Cuántas células, cuántas características (features) y cuántos ensayos (assays) hay?
• ¿Cuál es el ensayo activo (default)? Ejecuta la línea de abajo para obtener los cuatro números.

cat("\nAnswer to QUESTION 1.1:\n")
cat("  Cells         :", ncol(sobj), "\n")
cat("  Features (RNA):", nrow(sobj[["RNA"]]), "\n")
cat("  Features (ADT):", nrow(sobj[["ADT"]]), "\n")
cat("  Assays        :", length(SafeAssays(sobj)), "(", paste(SafeAssays(sobj), collapse = ", "), ")\n")
cat("  Default assay :", DefaultAssay(sobj), "\n")

TipRESPUESTA

3,000 células a través de 2 ensayos. El ensayo RNA tiene 500 características, un subconjunto con fines didácticos (no el transcriptoma completo). El ensayo ADT tiene 24 proteínas. El ensayo por defecto es RNA. El subconjunto de 500 genes es relevante: cualquier umbral de control de calidad (QC) copiado de un tutorial de 33,000 genes será incorrecto. Esta es la primera oportunidad para anclar la discusión de que ‘los valores por defecto de los tutoriales no se transfieren directamente’.

DEMOSTRACIÓN EN VIVO:

ncol(sobj)                    # cells

[1] 3000

nrow(sobj[["RNA"]])           # RNA features (500)

[1] 500

nrow(sobj[["ADT"]])           # ADT features (24)

[1] 24

Assays(sobj)                  # c("RNA", "ADT")

An object of class "SimpleAssays"
Slot "data":
List of length 1

DefaultAssay(sobj)            # "RNA"

[1] "RNA"

4.4.2 Paso 1.1b - Agregar genes mitocondriales (este panel no incluye ninguno)

El panel de 500 genes usado para este curso fue curado en torno a marcadores de linaje y activación; no incluye genes MT-. Cada paso de QC que depende de percent.mt ( Sección 4.5, Bloque 2) necesita una señal real para ser útil, así que agregamos aquí 9 genes MT- sintéticos, con conteos correlacionados con el tamaño total de la biblioteca de cada célula y un componente de estrés para un subconjunto de células.

set.seed(7)
mt_gene_names <- c("MT-ND1", "MT-ND2", "MT-CO1", "MT-CO2", "MT-ATP8",
                    "MT-ATP6", "MT-CO3", "MT-ND3", "MT-CYB")

if (!any(grepl("^MT-", rownames(sobj[["RNA"]])))) {
  rna_counts  <- LayerData(sobj, assay = "RNA", layer = "counts")
  total_count <- Matrix::colSums(rna_counts)

  # A log-normal draw gives a continuous, right-skewed distribution: most
  # cells sit at a low, biologically typical fraction, with a gradual tail
  # toward the small subset of stressed/dying cells. This mirrors real PBMC
  # data, where percent.mt has no gap between "healthy" and "stressed"; it
  # is one continuous distribution that QC thresholds cut into at some
  # chosen point.
  n_cells     <- ncol(rna_counts)
  target_frac <- rlnorm(n_cells, meanlog = log(0.035), sdlog = 0.6)
  target_frac <- pmin(target_frac, 0.35)  # cap at a biologically plausible ceiling

  mt_total_per_cell <- round(total_count * target_frac)
  mt_total_per_cell[mt_total_per_cell < length(mt_gene_names)] <- length(mt_gene_names)

  # Spread each cell's MT total across the 9 genes with unequal weights
  # (mirrors real data, where MT-CO1/MT-CO3/MT-CYB are usually the highest).
  gene_weights <- c(0.10, 0.08, 0.22, 0.14, 0.03, 0.07, 0.18, 0.06, 0.12)
  mt_mat <- sapply(seq_len(n_cells), function(i) {
    as.integer(round(mt_total_per_cell[i] * gene_weights))
  })
  rownames(mt_mat) <- mt_gene_names
  colnames(mt_mat) <- colnames(rna_counts)
  mt_mat_sparse <- as(mt_mat, "CsparseMatrix")

  rna_counts_with_mt <- rbind(rna_counts, mt_mat_sparse)
  # The warning here ("Different cells and/or features from existing assay
  # RNA") comes from the [[<- replacement method noticing the new assay has
  # 509 features instead of 500, not from CreateAssay5Object itself. It is
  # the expected, harmless side effect of adding genes to an assay; wrap the
  # whole assignment so it does not look like an error in the console.
  new_rna_assay <- CreateAssay5Object(counts = rna_counts_with_mt)
  suppressWarnings(sobj[["RNA"]] <- new_rna_assay)
  DefaultAssay(sobj) <- "RNA"

  cat("Added", length(mt_gene_names), "synthetic MT- genes to the RNA assay.\n")
  cat("RNA assay now has", nrow(sobj[["RNA"]]), "features (was 500).\n")
} else {
  cat("MT- genes already present; skipping synthesis.\n")
}

MT- genes already present; skipping synthesis.

4.4.3 Paso 1.2 - Seurat v4 vs v5: cómo se almacena la matriz de conteos

Seurat v5 cambió la forma en que las matrices de conteos viven dentro de un ensayo. El código de v4 que accede directamente a @counts o usa slot= producirá un error. Los dos patrones a continuación son fallos deliberados: lee cada mensaje de error.

# >>> DELIBERATE ERROR - read the message, then continue. <<<
try({
  # In Seurat v4, this was the standard way to get the count matrix.
  counts_old <- sobj@assays$RNA@counts
})

Error in try({ : 
  no slot of name "counts" for this object of class "Assay5"

# >>> DELIBERATE ERROR - read the message, then continue. <<<
try({
  counts_old2 <- GetAssayData(sobj, slot = "counts")
})

Error : The `slot` argument of `GetAssayData()` was deprecated in SeuratObject
5.0.0 and is now defunct.
ℹ Please use the `layer` argument instead.

TipExplicación

En Seurat v5, el slot @counts ya no existe en los objetos Assay5. El argumento correcto ahora es layer = "counts". Estos errores son ejemplos intencionales: muestran lo que sucede cuando ejecutas scripts de pipelines antiguos o recibes objetos creados con versiones anteriores de Seurat. Lee el mensaje de error, entiende por qué ocurre, y luego continúa.

Patrón 1: LayerData(). Recomendado para Seurat v5.

counts_rna <- LayerData(sobj, assay = "RNA", layer = "counts")
cat("RNA counts matrix dimensions (genes x cells):\n")

RNA counts matrix dimensions (genes x cells):

cat(" ", dim(counts_rna), "\n\n")

  500 3000 

Patrón 1b: tres formas equivalentes de obtener los conteos de UN gen sin convertir toda la matriz dispersa (sparse) a densa.

fetched <- FetchData(sobj, vars = "CD4", layer = "counts")[, 1]
indexed <- LayerData(sobj, assay = "RNA", layer = "counts")["CD4", ]
get_ad  <- GetAssayData(sobj, assay = "RNA", layer = "counts")["CD4", ]
cat("Three equivalent accessors for one gene's counts:\n")

Three equivalent accessors for one gene's counts:

cat("  FetchData            length:", length(fetched), " class:", class(fetched), "\n")

  FetchData            length: 3000  class: numeric 

cat("  LayerData[gene, ]    length:", length(indexed), " class:", class(indexed), "\n")

  LayerData[gene, ]    length: 3000  class: numeric 

cat("  GetAssayData[gene, ] length:", length(get_ad),  " class:", class(get_ad),  "\n")

  GetAssayData[gene, ] length: 3000  class: numeric 

cat("  All identical        :", all(fetched == indexed) && all(indexed == get_ad), "\n\n")

  All identical        : TRUE 

Patrón 2: clase del assay

cat("RNA assay class :", class(sobj[["RNA"]]), "\n")

RNA assay class : Assay5 

cat("ADT assay class :", class(sobj[["ADT"]]), "\n\n")

ADT assay class : Assay5 

Patrón 3: capas disponibles

cat("RNA layers      :", paste(SafeLayers(sobj[["RNA"]]), collapse = ", "), "\n")

RNA layers      : counts, data, scale.data 

cat("ADT layers      :", paste(SafeLayers(sobj[["ADT"]]), collapse = ", "), "\n\n")

ADT layers      : counts, data 

Patrón 4: versión del objeto Seurat

cat("Object version  :", as.character(sobj@version), "\n")

Object version  : 5.1.0 

TipPREGUNTA 1.2

• ¿Cuál es la clase del ensayo RNA? ¿Y del ensayo ADT?
• ¿Por qué podrían ser de la misma clase aunque los datos se comporten de forma muy diferente?

TipRESPUESTA

Ambos son de clase Assay5. Assay5 es un contenedor de almacenamiento, no una especificación de normalización. El contenedor mantiene un sistema de capas (counts, data, scale.data) de la misma manera para cualquier modalidad. Las diferencias entre RNA y ADT (disperso vs denso, dropout vs bimodalidad, LogNormalize vs CLR, escala logarítmica vs escala CLR) son propiedades de los valores almacenados en las capas y de la normalización elegida, no de la clase del ensayo. Una función que opera sobre sobj[['ANY']] funciona mecánicamente en ambos. Esto también significa que ejecutar una normalización específica de RNA sobre ADT no produce un error de tipo.

DEMOSTRACIÓN EN VIVO:

class(sobj[["RNA"]])    # "Assay5"

[1] "Assay5"
attr(,"package")
[1] "SeuratObject"

class(sobj[["ADT"]])    # "Assay5"

[1] "Assay5"
attr(,"package")
[1] "SeuratObject"

# Same class, different content. Confirm by inspecting layers:
Layers(sobj[["RNA"]])   # counts (data and scale.data appear after preprocessing)

[1] "counts"     "data"       "scale.data"

Layers(sobj[["ADT"]])   # counts

[1] "counts" "data"  

TipPREGUNTA 1.2c

FetchData(), LayerData()[gene, ], y GetAssayData()[gene, ] son todas equivalentes aquí. ¿Por qué podrían NO ser equivalentes para datos normalizados en lugar de conteos crudos? (Piensa en cuál es el valor predeterminado de cada función para el argumento de capa/layer).

TipRESPUESTA

En la capa de conteos, las tres leen valores enteros crudos de la misma matriz dispersa subyacente, por lo que el resultado es idéntico. Para datos normalizados difieren en los valores predeterminados de sus argumentos: FetchData() usa por defecto layer='data' (normalizado); LayerData() requiere layer= explícito; GetAssayData() en versiones más nuevas de Seurat también espera layer= (slot= está obsoleto, ver E2). Aparecen discrepancias cuando una llamada usa por defecto data y otra counts. Siempre pasa layer= de forma explícita.

DEMOSTRACIÓN EN VIVO:

# After Block 4:
fd <- FetchData(sobj_filt, vars = "CD3E")[, 1] # data layer
lc <- LayerData(sobj_filt, assay="RNA", layer="counts")["CD3E", ] # counts
ld <- LayerData(sobj_filt, assay="RNA", layer="data")["CD3E", ] # data
all(fd == ld)   # TRUE
all(fd == lc)   # FALSE (different layers)

4.4.4 Paso 1.2b - Mapa completo de slots de un objeto Seurat

Un objeto Seurat tiene más slots de los que muestran la mayoría de los tutoriales. Algunos contienen datos crudos, otros contienen resultados derivados, otros son cachés, y otros son metadatos sobre el historial del análisis. Por diseño, el mismo valor a menudo vive en dos o tres lugares. Conocer el mapa evita tres clases de errores:

1. Leer de la copia equivocada después de que otra fue actualizada.
2. No encontrar datos que existen bajo un slot que no conocías.
3. Confiar en resultados posteriores cuando un slot anterior fue sobrescrito.

A.1 - Slots de nivel superior: imprime cada nombre de slot y una descripción de 1 línea de cada uno

cat("\n== TOP-LEVEL SLOTS ==\n")

== TOP-LEVEL SLOTS ==

top_slots <- slotNames(sobj)
for (s in top_slots) cat(sprintf("  @%-15s class: %s\n", s, class(slot(sobj, s))[1]))

  @assays          class: list
  @meta.data       class: data.frame
  @active.assay    class: character
  @active.ident    class: factor
  @graphs          class: list
  @neighbors       class: list
  @reductions      class: list
  @images          class: list
  @project.name    class: character
  @misc            class: list
  @version         class: package_version
  @commands        class: list
  @tools           class: list

A.2 - @assays: lista de ensayos. Cada ensayo es un Assay5 (v5) o Assay (v4)

cat("\n== @assays ==\n")

== @assays ==

cat("  Assay names              :", paste(names(sobj@assays), collapse = ", "), "\n")

  Assay names              : RNA, ADT 

cat("  Active assay             :", DefaultAssay(sobj), "\n")

  Active assay             : RNA 

cat("  Number of assays         :", length(sobj@assays), "\n")

  Number of assays         : 2 

Dentro de un ensayo (v5 Assay5), el mapa de slots es diferente al de v4:

cat("\n  [inside sobj[['RNA']] (v5 Assay5)]\n")

  [inside sobj[['RNA']] (v5 Assay5)]

rna_assay_slots <- slotNames(sobj[["RNA"]])
for (s in rna_assay_slots) {
  val <- slot(sobj[["RNA"]], s)
  cat(sprintf("    @%-15s class: %-15s length/dim: %s\n",
              s, class(val)[1],
              if (is.list(val))    paste0("list of ", length(val))
              else if (is.null(dim(val))) paste0("len ", length(val))
              else paste(dim(val), collapse = " x ")))
}

    @layers          class: list            length/dim: list of 3
    @cells           class: LogMap          length/dim: 3000 x 3
    @features        class: LogMap          length/dim: 500 x 3
    @default         class: integer         length/dim: len 1
    @assay.orig      class: character       length/dim: len 0
    @meta.data       class: data.frame      length/dim: list of 8
    @misc            class: list            length/dim: list of 0
    @key             class: character       length/dim: len 1

A.3 - @meta.data: data.frame de metadatos a nivel de célula

Filas = células. Columnas = lo que se haya agregado durante el preprocesamiento.

cat("\n== @meta.data (cell-level) ==\n")

== @meta.data (cell-level) ==

cat("  Dimensions       :", dim(sobj@meta.data)[1], "cells x",
                            dim(sobj@meta.data)[2], "columns\n")

  Dimensions       : 3000 cells x 14 columns

cat("  Column names     :\n")

  Column names     :

print(colnames(sobj@meta.data))

 [1] "orig.ident"      "nCount_RNA"      "nFeature_RNA"    "percent.mt"     
 [5] "donor_id"        "condition"       "severity"        "age"            
 [9] "cell_type"       "nCount_ADT"      "nFeature_ADT"    "percent.ribo"   
[13] "RNA_snn_res.0.5" "seurat_clusters"

A.4 - @active.ident: un FACTOR sobre las células. La "identidad" actual usada por las funciones posteriores de Seurat (FindMarkers, DimPlot group.by = NULL default)

Es independiente de las columnas de metadatos y se establece con SetIdent / Idents().

cat("\n== @active.ident ==\n")

== @active.ident ==

cat("  Class            :", class(sobj@active.ident), "\n")

  Class            : factor 

cat("  Levels currently :", paste(head(levels(sobj@active.ident), 5), collapse = ", "),
                             ifelse(length(levels(sobj@active.ident)) > 5, " ...", ""), "\n")

  Levels currently : 0, 1, 2, 3, 4  ... 

cat("  First 5 cells    :\n")

  First 5 cells    :

print(head(sobj@active.ident, 5))

CELL000001 CELL000011 CELL000021 CELL000031 CELL000041 
         8          9         10          5          2 
Levels: 0 1 2 3 4 5 6 7 8 9 10 11

A.5 - @reductions: lista de objetos DimReduc (pca, umap, harmony, etc.)

Cada uno tiene sus PROPIOS slots: cell.embeddings (células x dimensiones), feature.loadings (características x dimensiones), stdev (varianza por dimensión), key (prefijo de columna), assay.used.

cat("\n== @reductions ==\n")

== @reductions ==

cat("  Reductions present:", paste(names(sobj@reductions), collapse = ", "), "\n")

  Reductions present: pca, umap 

if (length(sobj@reductions) > 0) {
  for (red_name in names(sobj@reductions)) {
    red <- sobj@reductions[[red_name]]
    cat(sprintf("  [%s]\n", red_name))
    cat(sprintf("    cell.embeddings : %s\n", paste(dim(red@cell.embeddings), collapse = " x ")))
    cat(sprintf("    feature.loadings: %s\n", paste(dim(red@feature.loadings), collapse = " x ")))
    cat(sprintf("    stdev length    : %d\n", length(red@stdev)))
    cat(sprintf("    key             : %s\n", red@key))
    cat(sprintf("    assay.used      : %s\n", red@assay.used))
  }
} else {
  cat("  (none yet; PCA/UMAP happen in Block 4)\n")
}

  [pca]
    cell.embeddings : 3000 x 30
    feature.loadings: 500 x 30
    stdev length    : 30
    key             : PC_
    assay.used      : RNA
  [umap]
    cell.embeddings : 3000 x 2
    feature.loadings: 0 x 0
    stdev length    : 0
    key             : umap_
    assay.used      : RNA

A.6 - @graphs y @neighbors: cachés llenados por FindNeighbors. Vacíos aquí

cat("\n== @graphs and @neighbors ==\n")

== @graphs and @neighbors ==

cat("  Graphs    :", ifelse(length(sobj@graphs)    == 0, "empty (filled by FindNeighbors)",
                            paste(names(sobj@graphs), collapse = ", ")), "\n")

  Graphs    : RNA_nn, RNA_snn 

cat("  Neighbors :", ifelse(length(sobj@neighbors) == 0, "empty (filled by FindNeighbors)",
                            paste(names(sobj@neighbors), collapse = ", ")), "\n")

  Neighbors : empty (filled by FindNeighbors) 

A.7 - @commands: cada llamada a una función de Seurat realizada sobre este objeto, con todos sus argumentos, proporcionando un historial completo de análisis. Si alguna vez te preguntas “¿qué argumentos pasó el usuario anterior a NormalizeData?”, revisa aquí.

cat("\n== @commands (analysis history) ==\n")

== @commands (analysis history) ==

cat("  Commands recorded:", length(sobj@commands), "\n")

  Commands recorded: 8 

if (length(sobj@commands) > 0) {
  cat("  First 5 command names:\n")
  print(head(names(sobj@commands), 5))
  # Inspect one command in detail
  first_cmd <- sobj@commands[[1]]
  cat(sprintf("\n  [detail of first command: %s]\n", names(sobj@commands)[1]))
  cat("    call.string   :", first_cmd@call.string[1], "\n")
  cat("    time.stamp    :", as.character(first_cmd@time.stamp), "\n")
  cat("    assay.used    :", first_cmd@assay.used, "\n")
  cat("    params (names):", paste(names(first_cmd@params), collapse = ", "), "\n")
}

  First 5 command names:
[1] "NormalizeData.RNA"        "FindVariableFeatures.RNA"
[3] "ScaleData.RNA"            "RunPCA.RNA"              
[5] "FindNeighbors.RNA.pca"   

  [detail of first command: NormalizeData.RNA]
    call.string   : NormalizeData(sobj, verbose = FALSE) 
    time.stamp    : 2026-05-19 18:41:26.947809 
    assay.used    : RNA 
    params (names): assay, normalization.method, scale.factor, margin, verbose 

A.8 - @misc, @tools, @project.name, @version: pequeños slots administrativos

cat("\n== Administrative slots ==\n")

== Administrative slots ==

cat("  @project.name :", sobj@project.name, "\n")

  @project.name : GSE149689 

cat("  @version      :", as.character(sobj@version), "\n")

  @version      : 5.1.0 

cat("  @misc names   :", ifelse(length(sobj@misc) == 0, "empty",
                                paste(names(sobj@misc), collapse = ", ")), "\n")

  @misc names   : empty 

cat("  @tools names  :", ifelse(length(sobj@tools) == 0, "empty",
                                paste(names(sobj@tools), collapse = ", ")), "\n")

  @tools names  : empty 

A.9 - Mismos datos, diferentes rutas: una fuente común de confusión.

Los valores de abajo son IDÉNTICOS, pero se acceden mediante diferentes slots/accesores.

cat("\n== Same data, different paths (sanity proofs) ==\n")

== Same data, different paths (sanity proofs) ==

cat("Test 1: assay access\n")

Test 1: assay access

cat("  identical(sobj@assays$RNA, sobj[['RNA']])           :",
    identical(sobj@assays$RNA, sobj[["RNA"]]), "\n")

  identical(sobj@assays$RNA, sobj[['RNA']])           : TRUE 

cat("Test 2: metadata column access (4 equivalent paths)\n")

Test 2: metadata column access (4 equivalent paths)

v1 <- sobj$nFeature_RNA
v2 <- sobj@meta.data$nFeature_RNA
v3 <- sobj[[]]$nFeature_RNA
v4 <- FetchData(sobj, vars = "nFeature_RNA")[, 1]

identical() es sensible a los atributos (names, integer vs numeric storage) que difieren entre estos cuatro accesores incluso cuando cada valor es el mismo.

Elimina los nombres y homogeneiza el tipo antes de comparar, ya que el punto aquí es la igualdad de valores, no la igualdad de atributos.

all_match <- all(unname(as.numeric(v1)) == unname(as.numeric(v2))) &&
             all(unname(as.numeric(v2)) == unname(as.numeric(v3))) &&
             all(unname(as.numeric(v3)) == unname(as.numeric(v4)))
cat("  all four return the same values                     :", all_match, "\n")

  all four return the same values                     : TRUE 

cat("  (identical() can still report FALSE here; that compares attributes\n")

  (identical() can still report FALSE here; that compares attributes

cat("   like names or integer-vs-numeric storage, not the values themselves)\n")

   like names or integer-vs-numeric storage, not the values themselves)

cat("Test 3: cell name access (3 equivalent paths)\n")

Test 3: cell name access (3 equivalent paths)

cat("  identical(Cells(sobj), colnames(sobj))              :",
    identical(Cells(sobj), colnames(sobj)), "\n")

  identical(Cells(sobj), colnames(sobj))              : TRUE 

cat("  identical(Cells(sobj), rownames(sobj@meta.data))    :",
    identical(Cells(sobj), rownames(sobj@meta.data)), "\n")

  identical(Cells(sobj), rownames(sobj@meta.data))    : TRUE 

cat("Test 4: feature name access (varies with default assay)\n")

Test 4: feature name access (varies with default assay)

cat("  rownames(sobj) returns features of the ACTIVE assay only\n")

  rownames(sobj) returns features of the ACTIVE assay only

cat("  default assay is", DefaultAssay(sobj), "\n")

  default assay is RNA 

cat("  identical(rownames(sobj), Features(sobj))           :",
    identical(rownames(sobj), Features(sobj)), "\n")

  identical(rownames(sobj), Features(sobj))           : TRUE 

cat("  identical(rownames(sobj), rownames(sobj[['RNA']])) :",
    identical(rownames(sobj), rownames(sobj[["RNA"]])), "\n")

  identical(rownames(sobj), rownames(sobj[['RNA']])) : TRUE 

cat("  identical(rownames(sobj), rownames(sobj[['ADT']])) :",
    identical(rownames(sobj), rownames(sobj[["ADT"]])), "\n")

  identical(rownames(sobj), rownames(sobj[['ADT']])) : FALSE 

cat("  >> rownames depends on DefaultAssay. Always pass assay= explicitly.\n")

  >> rownames depends on DefaultAssay. Always pass assay= explicitly.

A.10 - ACERTIJOS DE SLOTS: ¿dónde buscarías para encontrar cada elemento a continuación?

Intenta escribir la respuesta mentalmente antes de ejecutar.

cat("\n== SLOT PUZZLES (try before running) ==\n")

== SLOT PUZZLES (try before running) ==

TipACERTIJO 1.2b/1

¿Dónde está la versión original de Seurat que creó este objeto?

cat("\nPUZZLE 1: original Seurat version that created this object?\n")

PUZZLE 1: original Seurat version that created this object?

cat("  Answer: sobj@version (top-level slot)\n")

  Answer: sobj@version (top-level slot)

cat("  Value :", as.character(sobj@version), "\n")

  Value : 5.1.0 

TipRESPUESTA

sobj@version. Slot de nivel superior que registra la versión de Seurat que originalmente construyó el objeto, no la versión actualmente cargada. Útil al depurar problemas de migración de v4 a v5.

DEMOSTRACIÓN EN VIVO:

sobj@version

[1] '5.1.0'

packageVersion("Seurat")   # version currently loaded

[1] '5.5.1'

TipPREGUNTA 1.2b

Para cada una de las 4 pruebas de verificación anteriores (A.9), ¿cuál está GARANTIZADA a ser idéntica sin importar el estado del análisis, y cuál DEPENDE de DefaultAssay? ¿Por qué es importante esta distinción al compartir código?

TipRESPUESTA

Las pruebas 1 y 3 están garantizadas idénticas vía identical(): el acceso al ensayo (sobj[['RNA']] vs sobj@assays$RNA) siempre devuelve el mismo objeto Assay5, y los nombres de células siempre viven en colnames(sobj). La prueba 2 (acceso a columnas de metadatos) es más sutil de lo que parece: sobj$nFeature_RNA, sobj@meta.data$nFeature_RNA, y sobj[[]]$nFeature_RNA devuelven el mismo vector con nombres, pero FetchData(sobj, vars='nFeature_RNA')[, 1] elimina el atributo de nombres de célula cuando se extrae la columna del data.frame con [, 1]. Los VALORES son idénticos; identical() no lo es, porque también compara el atributo de nombres. Este es un buen ejemplo de que identical() es demasiado estricto para la pregunta que realmente se está haciendo: al verificar la igualdad de valores entre accesores, elimina los nombres (por ejemplo unname()) o compara numéricamente (all(x == y)) en lugar de usar identical() directamente. La prueba 4 es la trampa con consecuencias reales: rownames(sobj) devuelve solo las características del ensayo ACTIVO. Si un colaborador escribe marker %in% rownames(sobj) asumiendo RNA, pero el ensayo activo se ha establecido en ADT, la verificación falla silenciosamente para cualquier gen que no sea también un nombre de proteína. Regla general al compartir o recibir código: pasa assay= explícitamente siempre que una llamada a función pueda resolverse de forma diferente dependiendo del ensayo activo, y no asumas que un fallo de identical() significa que los valores difieren; también puede significar que solo un atributo difiere.

DEMOSTRACIÓN EN VIVO:

# Run the 4 tests live:
identical(sobj@assays$RNA, sobj[["RNA"]])                  # TRUE always

[1] TRUE

# Test 2: identical() can be FALSE here even though values match,
# because FetchData()[, 1] drops cell names that $ and [[]] keep:
identical(sobj$nFeature_RNA, sobj@meta.data$nFeature_RNA)  # TRUE

[1] FALSE

v4 <- FetchData(sobj, vars = "nFeature_RNA")[, 1]
identical(sobj$nFeature_RNA, v4)                           # often FALSE (names)

[1] FALSE

all(unname(sobj$nFeature_RNA) == unname(v4))               # TRUE (values match)

[1] TRUE

identical(Cells(sobj), colnames(sobj))                     # TRUE always

[1] TRUE

# Now the trap:
DefaultAssay(sobj) <- "RNA"; head(rownames(sobj))          # gene symbols

[1] "CD3E" "CD3D" "CD3G" "CD4"  "CD8A" "CD8B"

DefaultAssay(sobj) <- "ADT"; head(rownames(sobj))          # protein names (different)

[1] "CD3"  "CD4"  "CD8a" "CD14" "CD16" "CD19"

DefaultAssay(sobj) <- "RNA"                                # restore

4.4.5 Paso 1.3 - Inspección del ensayo RNA (intermedio)

Tres propiedades de una matriz de conteos de RNA de célula única que determinan cada decisión posterior: dispersión (sparsity), huella de memoria y estado de las capas. Inspecciónalas ahora, antes de cualquier normalización o escalado. La contraparte de ADT de este paso abre el Bloque 6 (Sección 5.1), una vez que el flujo de trabajo de CITE-seq realmente lo necesite.

Dispersión (sparsity) = fracción de entradas cero. Determina las decisiones de normalización.

counts_mat <- LayerData(sobj, assay = "RNA", layer = "counts")
total_entries  <- prod(dim(counts_mat))
nonzero        <- Matrix::nnzero(counts_mat)
sparsity       <- 1 - nonzero / total_entries

cat("RNA count matrix:\n")

RNA count matrix:

cat("  Genes   :", nrow(counts_mat), "\n")

  Genes   : 500 

cat("  Cells   :", ncol(counts_mat), "\n")

  Cells   : 3000 

cat("  Total entries  :", format(total_entries, big.mark = ","), "\n")

  Total entries  : 1,500,000 

cat("  Non-zero entries:", format(nonzero, big.mark = ","), "\n")

  Non-zero entries: 142,375 

cat("  Sparsity        :", round(sparsity * 100, 1), "%\n\n")

  Sparsity        : 90.5 %

TipPREGUNTA 1.3a

Un conjunto de datos típico de PBMC con transcriptoma completo muestra una dispersión de RNA superior al 90%. Este panel de 500+9 genes es menor. ¿Qué te dice el valor de dispersión que acabas de imprimir sobre el tamaño del panel de genes versus la tasa de dropout?

TipRESPUESTA

La dispersión en este panel es menor que en un transcriptoma completo principalmente porque los 500 genes fueron curados como marcadores de linaje y activación, que tienden a expresarse de forma más alta y consistente que el gen promedio en un transcriptoma completo (la mayoría de los genes en un panel completo tienen baja expresión o están restringidos a tipos celulares específicos, lo cual impulsa la dispersión por encima del 90%). Un panel más pequeño y curado no es inmune al dropout: las mismas ineficiencias de captura y transcripción reversa aplican por molécula sin importar el tamaño del panel. Lo que cambia es el nivel promedio de expresión de los genes incluidos, no la biología subyacente del dropout. Esta distinción importa al leer un valor de dispersión de cualquier conjunto de datos: un valor de dispersión bajo puede significar un panel curado de genes bien expresados, no necesariamente un experimento técnicamente superior.

DEMOSTRACIÓN EN VIVO:

rna <- LayerData(sobj, assay="RNA", layer="counts")
mean(rna == 0) * 100        # RNA sparsity % for this panel

[1] 90.50833

Compara con la expectativa de un transcriptoma completo (comúnmente 90%+ para 10x PBMC)

PrecauciónPrecaución

El almacenamiento denso vs disperso importa a gran escala. Un experimento completo de 10x con 50,000 células y 33,000 genes como matriz densa supera los 50 GB.

4.4.6 Comparar memoria dispersa (counts) vs densa (después del scaling)

counts_rna_sz <- object.size(LayerData(sobj, assay = "RNA", layer = "counts"))
cat("RNA counts layer (sparse):", format(counts_rna_sz, units = "Mb"), "\n")

RNA counts layer (sparse): 1.9 Mb 

El objeto completo:

obj_size <- object.size(sobj)
cat("Full Seurat object       :", format(obj_size, units = "Mb"), "\n\n")

Full Seurat object       : 23.9 Mb 

cat("\nscale.data is a DENSE genes-by-cells matrix.\n")

scale.data is a DENSE genes-by-cells matrix.

cat("On full datasets, scale only the genes used in PCA to keep memory bounded.\n")

On full datasets, scale only the genes used in PCA to keep memory bounded.

TipPREGUNTA 1.3b

La capa de conteos es dispersa, la capa scale.data es densa. Proyecta esto a un experimento de 50,000 células y 33,000 genes: ¿cuál es la implicación en memoria, y qué te dice sobre cómo usar ScaleData?

TipRESPUESTA

Matriz densa de doble precisión: 33,000 características x 50,000 células x 8 bytes = 13.2 GB. Una representación dispersa con 90% de dispersión almacena aproximadamente 33,000 x 50,000 x 0.10 valores no cero x 16 bytes por valor no cero = 2.6 GB. ScaleData centra y escala cada característica, produciendo una matriz densa incluso cuando la entrada era dispersa. Ejecutar ScaleData sobre todas las características a esta escala agotará la RAM de cualquier laptop. Práctica estándar: escalar solo los genes altamente variables usados en PCA, típicamente 2,000-3,000 características. La matriz escalada densa se convierte en 3,000 x 50,000 x 8 = 1.2 GB, manejable. Pasa features = VariableFeatures(sobj) a ScaleData.

DEMOSTRACIÓN EN VIVO:

object.size(LayerData(sobj, assay="RNA", layer="counts"))   # sparse

1973416 bytes

object.size(as.matrix(LayerData(sobj, assay="RNA", layer="counts")))  # dense

12251800 bytes

Buena práctica (ya usada en Paso 4.4, Sección 4.7.4):

ScaleData(sobj, features = VariableFeatures(sobj))

4.4.7 Paso 1.6 - Metadatos

Las reducciones (PCA, UMAP, Harmony, etc.) se calculan a partir del Bloque 4 (Sección 4.7) en adelante. En este punto del script, aún no existe ninguna. Confirma esto explícitamente en lugar de asumirlo: es el mismo hábito que verificar DefaultAssay() antes de confiar en cualquier accesor. Los patrones de acceso para cell.embeddings y feature.loadings se cubren en el Paso 4.4 (Sección 4.7.4), una vez que el PCA realmente existe.

cat("Reductions in the object:\n")

Reductions in the object:

print(SafeReductions(sobj))

[1] "pca"  "umap"

cat("(empty; PCA happens in Block 4.)\n\n")

(empty; PCA happens in Block 4.)

cat("Metadata columns and types:\n")

Metadata columns and types:

str(sobj@meta.data)

'data.frame':   3000 obs. of  14 variables:
 $ orig.ident     : chr  "Donor01" "Donor01" "Donor01" "Donor01" ...
 $ nCount_RNA     : num  336 344 426 330 473 488 413 294 455 336 ...
 $ nFeature_RNA   : int  36 40 39 35 51 47 55 42 36 47 ...
 $ percent.mt     : num  18.15 15.12 0 3.33 5.07 ...
 $ donor_id       : chr  "Donor01" "Donor01" "Donor01" "Donor01" ...
 $ condition      : chr  "Healthy" "Healthy" "Healthy" "Healthy" ...
 $ severity       : chr  "" "" "" "" ...
 $ age            : int  35 35 35 35 35 35 35 35 35 35 ...
 $ cell_type      : chr  "NKcell" "DC" "Bcell" "Monocyte" ...
 $ nCount_ADT     : num  182 202 265 135 199 133 209 141 150 280 ...
 $ nFeature_ADT   : int  16 18 14 10 16 16 15 18 16 9 ...
 $ percent.ribo   : num  46.1 64.5 66.7 60.9 39.7 ...
 $ RNA_snn_res.0.5: Factor w/ 12 levels "0","1","2","3",..: 9 10 11 6 3 3 3 8 6 11 ...
 $ seurat_clusters: Factor w/ 12 levels "0","1","2","3",..: 9 10 11 6 3 3 3 8 6 11 ...

El análisis real requiere constantemente subconjuntar o consultar células por combinaciones de metadatos. Practica los patrones aquí.

• ¿Cuántas células por donante?

cat("Cells per donor:\n")

Cells per donor:

print(table(sobj$donor_id))

Donor01 Donor02 Donor03 Donor04 Donor05 Donor06 Donor07 
    350     300     350     450     450     550     550 

• ¿Cuántas células por condición?

cat("\nCells per condition:\n")

Cells per condition:

print(table(sobj$condition))

COVID19 Healthy 
   2000    1000 

Células de donantes con COVID-19 con más de 40 genes detectados. El umbral está calibrado para este panel: nFeature_RNA se ubica en el rango de 27-96 aquí, no en el rango de 200+ típico de un transcriptoma completo, así que un número tomado prestado de un tutorial de transcriptoma completo coincidiría silenciosamente con cero células.

covid_highgene <- sobj@meta.data %>%
  filter(condition == "COVID19", nFeature_RNA > 40)
cat("\nCOVID-19 cells with nFeature_RNA > 40:", nrow(covid_highgene), "\n")

COVID-19 cells with nFeature_RNA > 40: 1862 

• Tabulación cruzada: condición x severidad

cat("\nCondition x Severity:\n")

Condition x Severity:

print(table(sobj$condition, sobj$severity))

         
               Mild Moderate Severe
  COVID19    0  450      450   1100
  Healthy 1000    0        0      0

Misma idea, variable diferente: severidad por donante. La tabulación revela el diseño donante-condición (qué donantes son Sanos frente a qué nivel de severidad se asignó a cada donante con COVID-19).

cat("\nSeverity by donor:\n")

Severity by donor:

print(table(sobj$severity, sobj$donor_id, useNA = "ifany"))

          
           Donor01 Donor02 Donor03 Donor04 Donor05 Donor06 Donor07
               350     300     350       0       0       0       0
  Mild           0       0       0     450       0       0       0
  Moderate       0       0       0       0     450       0       0
  Severe         0       0       0       0       0     550     550

4.4.8 Paso 1.7 - Nomenclatura de las características de ADT

• Los nombres de las características de ADT a menudo difieren de los nombres de los genes de RNA.
• La proteína CD3 != el gen CD3E. La proteína CD8a != el gen CD8A.
• Esta discrepancia es una fuente frecuente de confusión.

cat("ADT protein names:\n")

ADT protein names:

print(rownames(sobj[["ADT"]]))

 [1] "CD3"    "CD4"    "CD8a"   "CD14"   "CD16"   "CD19"   "CD20"   "CD25"  
 [9] "CD27"   "CD38"   "CD45RA" "CD45RO" "CD56"   "CD57"   "CD62L"  "CD69"  
[17] "CD86"   "CD127"  "CD197"  "PD1"    "TIGIT"  "LAG3"   "IgD"    "HLADR" 

cat("\nRNA marker names for comparison:\n")

RNA marker names for comparison:

rna_markers <- c("CD3E", "CD4", "CD8A", "CD14", "CD19", "NCAM1")
cat(rna_markers, "\n")

CD3E CD4 CD8A CD14 CD19 NCAM1 

cat("\nImportant: 'CD8a' in ADT vs 'CD8A' in RNA.\n")

Important: 'CD8a' in ADT vs 'CD8A' in RNA.

cat("FetchData() handles this transparently, but FeaturePlot()\n")

FetchData() handles this transparently, but FeaturePlot()

cat("requires you to be on the correct DefaultAssay first.\n")

requires you to be on the correct DefaultAssay first.

Demostración: qué sucede cuando intentas acceder a una característica de ADT usando el nombre del gen de RNA:

cat("\nDoes 'CD8A' exist in ADT?\n")

Does 'CD8A' exist in ADT?

cat("'CD8A' in ADT rownames:", "CD8A" %in% rownames(sobj[["ADT"]]), "\n")

'CD8A' in ADT rownames: FALSE 

cat("'CD8a' in ADT rownames:", "CD8a" %in% rownames(sobj[["ADT"]]), "\n")

'CD8a' in ADT rownames: TRUE 

4.4.9 Paso 1.8 - Marcador almacenado bajo un alias

Los símbolos de genes tienen sinónimos: NCAM1 = CD56, FCGR3A = CD16, MS4A1 = CD20. Si un objeto almacena un gen bajo su alias, FeaturePlot("FCGR3A") devuelve “feature not found” y el marcador se lee como ausente cuando en realidad los datos están intactos.

Confirma que los marcadores canónicos existen bajo su símbolo esperado antes de graficar o anotar.

canonical <- c("CD3E", "CD4", "CD8A", "CD14", "CD19", "MS4A1", "NCAM1", "FCGR3A")
present   <- canonical %in% rownames(sobj[["RNA"]])

cat("Canonical RNA markers present under expected symbol:\n")

Canonical RNA markers present under expected symbol:

for (i in seq_along(canonical))
  cat(sprintf("  %-8s %s\n", canonical[i], ifelse(present[i], "OK", "MISSING")))

  CD3E     OK
  CD4      OK
  CD8A     OK
  CD14     OK
  CD19     OK
  MS4A1    OK
  NCAM1    OK
  FCGR3A   OK

Si algo está MISSING, busca alias conocidos en los rownames

aliases <- list(FCGR3A = "CD16", NCAM1 = "CD56", MS4A1 = "CD20")
missing <- canonical[!present]
if (length(missing) > 0) {
  cat("\nMISSING markers detected. Searching for aliases:\n")
  rna_counts <- LayerData(sobj, assay = "RNA", layer = "counts")
  rn <- rownames(rna_counts)
  changed <- FALSE
  for (m in missing) {
    al <- aliases[[m]]
    if (!is.null(al) && al %in% rn) {
      cat(sprintf("  %s is stored as alias '%s'. Renaming back.\n", m, al))
      rn[rn == al] <- m
      changed <- TRUE
    }
  }
  if (changed) {
    # v5-safe rename: rebuild the RNA assay from counts with corrected names.
    # Run before normalization (Block 4), so only the counts layer exists,
    # matching the structure of the assay it replaces. suppressWarnings()
    # wraps the assignment itself, since any "different features" notice
    # would come from the [[<- replacement method, not from
    # CreateAssay5Object().
    rownames(rna_counts) <- rn
    new_rna_assay <- CreateAssay5Object(counts = rna_counts)
    suppressWarnings(sobj[["RNA"]] <- new_rna_assay)
    DefaultAssay(sobj) <- "RNA"
  }
  cat("\nAfter fix, all canonical markers present:",
      all(canonical %in% rownames(sobj[["RNA"]])), "\n")
}

TipPREGUNTA 1.8

¿Qué otros alias de RNA verificarías rutinariamente en un conjunto de datos de PBMC?

TipRESPUESTA

Pares de alias comunes: MS4A1/CD20 (células B), NCAM1/CD56 (NK), FCGR3A/CD16 (NK y monocitos no clásicos), ITGAX/CD11c (DC, mono), ITGAM/CD11b (mieloide), PTPRC/CD45 (pan-inmune), IL3RA/CD123 (pDC, basófilo), CD3E/CD3, CD8A/CD8a (las mayúsculas importan: RNA en mayúsculas, ADT en minúsculas), FOXP3 (Treg), FCER1A (DC). Hábito defensivo: mantén una lista de marcadores canónicos de PBMC y ejecuta %in% rownames(sobj[['RNA']]) al inicio de cada paso de anotación.

DEMOSTRACIÓN EN VIVO:

canonical <- c("CD3E","CD3D","CD8A","CD4","MS4A1","CD79A","NCAM1",
               "FCGR3A","CD14","LYZ","FOXP3","FCER1A","KLRD1")
data.frame(marker = canonical,
           present = canonical %in% rownames(sobj[["RNA"]]))

   marker present
1    CD3E    TRUE
2    CD3D    TRUE
3    CD8A    TRUE
4     CD4    TRUE
5   MS4A1    TRUE
6   CD79A    TRUE
7   NCAM1    TRUE
8  FCGR3A    TRUE
9    CD14    TRUE
10    LYZ    TRUE
11  FOXP3    TRUE
12 FCER1A    TRUE
13  KLRD1    TRUE

# Check aliases for missing ones:
aliases <- c(MS4A1="CD20", NCAM1="CD56", FCGR3A="CD16")

4.5 Bloque 2 - Métricas de control de calidad (RNA)

Objetivo: Calcular métricas de control de calidad, observar sus distribuciones a través de los donantes, y filtrar células usando umbrales elegidos a partir de los datos, no de un tutorial.

4.5.1 Paso 2.1 - Calcular métricas de control de calidad

• La fracción mitocondrial (percent.mt) es un indicador de muerte/estrés celular.
• La fracción ribosomal (percent.ribo) señala células dominadas por transcritos de mantenimiento (housekeeping).

sobj[["percent.mt"]]   <- PercentageFeatureSet(sobj, pattern = "^MT-")
sobj[["percent.ribo"]] <- PercentageFeatureSet(sobj, pattern = "^RP[SL]")

summary(sobj@meta.data[, c("nFeature_RNA", "nCount_RNA", "percent.mt")])

  nFeature_RNA     nCount_RNA      percent.mt    
 Min.   :27.00   Min.   :172.0   Min.   : 0.000  
 1st Qu.:43.00   1st Qu.:381.0   1st Qu.: 2.788  
 Median :47.00   Median :451.0   Median : 5.290  
 Mean   :47.46   Mean   :451.8   Mean   : 6.463  
 3rd Qu.:52.00   3rd Qu.:517.0   3rd Qu.: 8.853  
 Max.   :69.00   Max.   :844.0   Max.   :28.605  

head(sobj@meta.data, 10)

           orig.ident nCount_RNA nFeature_RNA percent.mt donor_id condition
CELL000001    Donor01        336           36  18.154762  Donor01   Healthy
CELL000011    Donor01        344           40  15.116279  Donor01   Healthy
CELL000021    Donor01        426           39   0.000000  Donor01   Healthy
CELL000031    Donor01        330           35   3.333333  Donor01   Healthy
CELL000041    Donor01        473           51   5.073996  Donor01   Healthy
CELL000051    Donor01        488           47   9.631148  Donor01   Healthy
CELL000061    Donor01        413           55  12.106538  Donor01   Healthy
CELL000071    Donor01        294           42  15.986395  Donor01   Healthy
CELL000081    Donor01        455           36   3.956044  Donor01   Healthy
CELL000091    Donor01        336           47   5.357143  Donor01   Healthy
           severity age  cell_type nCount_ADT nFeature_ADT percent.ribo
CELL000001           35     NKcell        182           16     46.13095
CELL000011           35         DC        202           18     64.53488
CELL000021           35      Bcell        265           14     66.66667
CELL000031           35   Monocyte        135           10     60.90909
CELL000041           35   CD4Tcell        199           16     39.74630
CELL000051           35   CD4Tcell        133           16     48.36066
CELL000061           35   CD4Tcell        209           15     39.46731
CELL000071           35 ncMonocyte        141           18     59.86395
CELL000081           35   Monocyte        150           16     58.24176
CELL000091           35      Bcell        280            9     52.08333
           RNA_snn_res.0.5 seurat_clusters
CELL000001               8               8
CELL000011               9               9
CELL000021              10              10
CELL000031               5               5
CELL000041               2               2
CELL000051               2               2
CELL000061               2               2
CELL000071               7               7
CELL000081               5               5
CELL000091              10              10

4.5.2 Paso 2.2 | Verificación de sanidad: ¿se calculó realmente percent.mt?

• percent.mt depende completamente de que el patrón ‘^MT-’ coincida con nombres de genes reales.

Vale la pena confirmar esto explícitamente en lugar de asumir que funcionó: un prefijo MT- renombrado o ausente (que es exactamente lo que simula el objeto inyectado) devuelve percent.mt = 0 para cada célula, y cualquier filtro basado en él entonces no hace nada o rechaza todas las células.

mt_genes <- grep("^MT-", rownames(sobj), value = TRUE)
cat("Genes matching '^MT-':", length(mt_genes), "\n")

Genes matching '^MT-': 9 

if (length(mt_genes) > 0) print(mt_genes)

[1] "MT-CO1"  "MT-CO2"  "MT-CO3"  "MT-ND1"  "MT-ND2"  "MT-ATP6" "MT-CYB" 
[8] "MT-RNR1" "MT-RNR2"

cat("\nAre rownames symbols or Ensembl IDs? First 5 rownames:\n")

Are rownames symbols or Ensembl IDs? First 5 rownames:

print(head(rownames(sobj), 5))

[1] "CD3E" "CD3D" "CD3G" "CD4"  "CD8A"

cat("\nmax(percent.mt):", round(max(sobj$percent.mt), 4), "\n")

max(percent.mt): 28.6052 

if (max(sobj$percent.mt) == 0) {
  cat("\nWARNING: percent.mt is 0 for all cells. The '^MT-' pattern matched\n")
  cat("no genes. Do NOT filter on percent.mt here: it is uninformative in\n")
  cat("this panel. Report this limitation; do not claim an mt-based QC.\n")
} else {
  cat("MT- pattern matched", length(mt_genes), "genes. percent.mt is usable.\n")
}

MT- pattern matched 9 genes. percent.mt is usable.

TipPREGUNTA 2.2

Si max(percent.mt) es 0 en todas las células, ¿cuáles son las dos posibles explicaciones, y cómo las distingues?

TipRESPUESTA

Explicación 1 (casi siempre): el patrón MT- no coincidió con ningún gen porque los símbolos mitocondriales usan un prefijo diferente (mt-, Mt-, MTX, o una referencia no humana). Diagnóstico: grep('^MT-', rownames(sobj)) devuelve character(0). Inspecciona head(rownames(sobj)) y busca el prefijo real. Explicación 2 (poco plausible): las células se filtraron previamente de forma tan agresiva que no queda ningún transcrito MT. La discrepancia de patrón es el caso realista. El conjunto de datos inyectado renombra los genes MT- a MTX- para provocar exactamente este fallo.

DEMOSTRACIÓN EN VIVO:

grep("^MT-", rownames(sobj), value = TRUE)

[1] "MT-CO1"  "MT-CO2"  "MT-CO3"  "MT-ND1"  "MT-ND2"  "MT-ATP6" "MT-CYB" 
[8] "MT-RNR1" "MT-RNR2"

grep("^mt-", rownames(sobj), value = TRUE)

character(0)

grep("^MTX", rownames(sobj), value = TRUE)

character(0)

head(rownames(sobj), 20)

 [1] "CD3E"   "CD3D"   "CD3G"   "CD4"    "CD8A"   "CD8B"   "IL7R"   "CCR7"  
 [9] "TCF7"   "SELL"   "FOXP3"  "IL2RA"  "PDCD1"  "LAG3"   "HAVCR2" "CD19"  
[17] "MS4A1"  "CD79A"  "CD79B"  "IGHM"  

summary(sobj$percent.mt)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.788   5.290   6.463   8.853  28.605 

4.5.3 Paso 2.3 - Visualizar las distribuciones de control de calidad

VlnPlot(sobj,
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        ncol = 3, pt.size = 0.05, alpha = 0.3)

p1 <- FeatureScatter(sobj, "nCount_RNA", "nFeature_RNA") +
  ggtitle("Counts vs genes detected") +
  theme(plot.title = element_text(size = 10))
p2 <- FeatureScatter(sobj, "nCount_RNA", "percent.mt") +
  ggtitle("Counts vs mitochondrial %") +
  theme(plot.title = element_text(size = 10))
p1 | p2

4.5.4 Paso 2.4 - Detección de valores atípicos (outliers)

Las células que tienen un nFeature alto en relación con nCount podrían ser doublets.

• nCount bajo con nFeature bajo = gota vacía (empty droplet).
• percent.mt alto por sí solo = célula estresada/moribunda.

Células que son potenciales valores atípicos: nFeature > 2 DE por encima de la media

mean_feat <- mean(sobj$nFeature_RNA)
sd_feat   <- sd(sobj$nFeature_RNA)
high_feat <- sum(sobj$nFeature_RNA > mean_feat + 2 * sd_feat)
cat("Cells with nFeature > mean + 2 SD (potential doublets):",
    high_feat, "\n")

Cells with nFeature > mean + 2 SD (potential doublets): 58 

Células con mt alto: > percentil 95

mt_95 <- quantile(sobj$percent.mt, 0.95)
high_mt <- sum(sobj$percent.mt > mt_95)
cat("Cells with percent.mt > 95th percentile:", high_mt, "\n")

Cells with percent.mt > 95th percentile: 150 

Los dos conteos anteriores son fáciles de calcular y fáciles de malinterpretar: un conteo por sí solo no muestra si las células marcadas forman un grupo claro y separable, o si el umbral cortó a través del medio de una distribución continua. Marca ambas categorías en los metadatos y grafícalas directamente contra los mismos ejes usados en el Paso 2.3 (Sección 4.5.3), para que los valores atípicos sean visibles como puntos, no solo como un número.

sobj$outlier_flag <- "normal"
sobj$outlier_flag[sobj$nFeature_RNA > mean_feat + 2 * sd_feat] <- "high nFeature (possible doublet)"
sobj$outlier_flag[sobj$percent.mt > mt_95]                     <- "high percent.mt (stressed/dying)"
sobj$outlier_flag <- factor(sobj$outlier_flag,
                            levels = c("normal", "high nFeature (possible doublet)",
                                       "high percent.mt (stressed/dying)"))

p_out1 <- FeatureScatter(sobj, "nCount_RNA", "nFeature_RNA", group.by = "outlier_flag") +
  ggtitle("Outlier flags on counts vs genes detected") +
  theme(legend.position = "bottom", legend.text = element_text(size = 7))

p_out2 <- FeatureScatter(sobj, "nCount_RNA", "percent.mt", group.by = "outlier_flag") +
  ggtitle("Outlier flags on counts vs mitochondrial %") +
  theme(legend.position = "bottom", legend.text = element_text(size = 7))

p_out1 | p_out2

TipPREGUNTA 2.4

Observando los dos gráficos, ¿las células con nFeature alto y las células con percent.mt alto ocupan regiones distintas, o algunas células califican como ambas? ¿Qué sería probablemente una célula marcada en ambos ejes?

TipRESPUESTA

En la mayoría de las ejecuciones, los dos grupos marcados son en gran parte distintos: los valores atípicos de nFeature alto se agrupan hacia el lado derecho del gráfico nCount-vs-nFeature (más genes detectados que la mayoría de las células con un conteo de UMI similar), mientras que los valores atípicos de percent.mt alto se agrupan en la región superior del gráfico nCount-vs-percent.mt sin importar el nFeature. Se espera algo de superposición y es el caso más informativo: una célula marcada en ambos ejes (nFeature alto Y percent.mt alto) es la más difícil de interpretar a partir de una sola métrica. Podría ser un doublet que también resulta estar estresado, o podrían ser dos problemas técnicos no relacionados que coinciden en la misma célula. La medida práctica no es tratar de asignar una sola causa; marca la célula como de baja confianza y deja que los pasos posteriores (puntuación de doublets en el Bloque 3, confianza de anotación en el Bloque 5, Sección 4.8) tomen la decisión final con más evidencia.

DEMOSTRACIÓN EN VIVO:

table(sobj$outlier_flag)

                          normal high nFeature (possible doublet) 
                            2797                               53 
high percent.mt (stressed/dying) 
                             150 

# Cells flagged on both definitions (recompute the masks to check overlap):
both <- (sobj$nFeature_RNA > mean_feat + 2*sd_feat) & (sobj$percent.mt > mt_95)
sum(both)

[1] 5

Distribución por condición: ¿está elevado el %mt en COVID-19?

cat("\nMedian percent.mt by condition:\n")

Median percent.mt by condition:

sobj@meta.data %>%
  group_by(condition) %>%
  summarise(median_mt = round(median(percent.mt), 2),
            mean_mt   = round(mean(percent.mt),   2),
            .groups   = "drop") %>%
  print()

# A tibble: 2 × 3
  condition median_mt mean_mt
  <chr>         <dbl>   <dbl>
1 COVID19        5.05    6.08
2 Healthy        5.74    7.23

4.5.5 Paso 2.5 - Filtro incorrecto: umbrales de tutorial aplicados a ciegas

Los umbrales estándar de tutorial (nFeature > 200 & < 6000, percent.mt < 5%) se heredan de experimentos 10x completos con ~33,000 genes. Este conjunto de datos es un subconjunto de 500 genes con fines didácticos. Aplicar los umbrales del tutorial elimina todas las células.

# >>> DELIBERATE ERROR - read the message, then continue. <<<
try({
  sobj_filt <- subset(
    sobj,
    subset = nFeature_RNA > 200 &
             nFeature_RNA < 5000 &
             percent.mt   < 20
  )
})

Error in subset(sobj, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 &  : 
  No cells found

Importante

“No cells found” es el resultado canónico de copiar umbrales de un tutorial sin verificar la distribución de tus propios datos. El umbral nFeature_RNA > 200 se calibró para un transcriptoma de ~33,000 genes donde las células típicamente detectan 2,000-3,000 genes. En un panel de 500 genes, la mayoría de las células detectan entre 50 y 150 genes. El filtro elimina cada una de las células.

TipPREGUNTA 2.5

¿Cuál de los tres parámetros de filtro anteriores es el más incorrecto para este conjunto de datos de 500 genes, y qué valor probarías primero? Ejecuta quantile(sobj$nFeature_RNA, c(0.05, 0.95)) para ver el rango real antes de proponer un número.

TipRESPUESTA

nFeature_RNA > 200 es el más incorrecto. El valor del tutorial 200 se calibró para ~33,000 características donde las células detectan 2,000-3,000 genes; aquí el panel de 500 genes produce entre 50 y 150 detectados por célula, así que un piso de 200 elimina esencialmente todas las células. Un umbral inicial razonable aquí es el percentil 5 de nFeature_RNA, típicamente cerca de 30-50. nFeature_RNA < 6000 es técnicamente correcto porque ninguna célula tiene cerca de 6000 características en un panel de 500 genes; simplemente no hace nada. percent.mt < 5 es límite; el efecto depende de si el patrón MT coincidió.

DEMOSTRACIÓN EN VIVO:

quantile(sobj$nFeature_RNA, c(0.05, 0.50, 0.95))

 5% 50% 95% 
 37  47  58 

summary(sobj$nFeature_RNA)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  27.00   43.00   47.00   47.46   52.00   69.00 

lo <- quantile(sobj$nFeature_RNA, 0.05)
hi <- quantile(sobj$nFeature_RNA, 0.95)
subset(sobj, subset = nFeature_RNA >= lo & nFeature_RNA <= hi)

An object of class Seurat 
524 features across 2712 samples within 2 assays 
Active assay: RNA (500 features, 500 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: ADT
 2 dimensional reductions calculated: pca, umap

Una primera estimación razonable para este panel: nFeature_RNA > 50. Pruébalo y observa cuántas células sobreviven, antes de pasar a los umbrales totalmente basados en datos en los Pasos 2.6 y 2.7 (Sección 4.5.6, Sección 4.5.7).

n_above_50 <- sum(sobj$nFeature_RNA > 50)
cat("Cells with nFeature_RNA > 50:", n_above_50,
    sprintf("(%.1f%% of total)\n", 100 * n_above_50 / ncol(sobj)))

Cells with nFeature_RNA > 50: 982 (32.7% of total)

4.5.6 Paso 2.6 - Comparar umbrales de tutorial vs basados en datos

cat("--- Published thresholds (full 10x, 33k genes) ---\n")

--- Published thresholds (full 10x, 33k genes) ---

cat("  nFeature_RNA > 200   removes empty droplets\n")

  nFeature_RNA > 200   removes empty droplets

cat("  nFeature_RNA < 5000  removes likely doublets\n")

  nFeature_RNA < 5000  removes likely doublets

cat("  percent.mt   < 20%   removes dead/stressed cells\n\n")

  percent.mt   < 20%   removes dead/stressed cells

cat("--- This dataset (500 genes) ---\n")

--- This dataset (500 genes) ---

print(round(quantile(sobj$nFeature_RNA,
      probs = c(0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99))))

 1%  5% 25% 50% 75% 95% 99% 
 33  37  43  47  52  58  62 

cat("\n--- Data-driven thresholds ---\n")

--- Data-driven thresholds ---

min_features <- max(10, round(quantile(sobj$nFeature_RNA, 0.02)))
max_features <- round(quantile(sobj$nFeature_RNA, 0.98))
max_mt       <- round(quantile(sobj$percent.mt,   0.95), 1)

cat("  nFeature_RNA > ", min_features, "  (same intent: empty droplets)\n")

  nFeature_RNA >  34   (same intent: empty droplets)

cat("  nFeature_RNA < ", max_features, " (same intent: doublets)\n")

  nFeature_RNA <  60  (same intent: doublets)

cat("  percent.mt   < ", max_mt, "%  (same intent: dead cells)\n")

  percent.mt   <  16.5 %  (same intent: dead cells)

4.5.7 Paso 2.7 - Aplicar umbrales basados en datos

cat("Cells before filtering:", ncol(sobj), "\n")

Cells before filtering: 3000 

sobj_filt <- tryCatch(
  subset(
    sobj,
    subset = nFeature_RNA > min_features &
             nFeature_RNA < max_features &
             percent.mt   < max_mt
  ),
  error = function(e) {
    cat("\nsubset() failed:", conditionMessage(e), "\n")
    cat("This means the three thresholds together match zero cells.\n")
    cat("Check each threshold individually before combining them:\n")
    cat("  > min_features:", sum(sobj$nFeature_RNA > min_features), "cells\n")
    cat("  < max_features:", sum(sobj$nFeature_RNA < max_features), "cells\n")
    cat("  < max_mt      :", sum(sobj$percent.mt   < max_mt),       "cells\n")
    stop(e)
  }
)

cat("Cells after filtering :", ncol(sobj_filt), "\n")

Cells after filtering : 2683 

cat("Cells removed         :", ncol(sobj) - ncol(sobj_filt), "\n")

Cells removed         : 317 

cat("Retention rate        :",
    round(ncol(sobj_filt) / ncol(sobj) * 100, 1), "%\n")

Retention rate        : 89.4 %

if (ncol(sobj_filt) < 500) {
  cat("\nWARNING: fewer than 500 cells remaining.\n")
  cat("Options: relax max_mt or lower min_features,\n")
  cat("or request pre-filtered matrices from GEO.\n")
} else {
  cat("\nCell count adequate for downstream analysis.\n")
}

Cell count adequate for downstream analysis.

TipPREGUNTA 2.7

¿Cómo cambiarían tus umbrales si el conjunto de datos tuviera ~33,000 genes en lugar de 500?

TipRESPUESTA

Los umbrales sobre los conteos de características detectadas escalan de forma aproximadamente lineal con el espacio de características, pero los umbrales sobre métricas de calidad (percent.mt, percent.ribo) no lo hacen. Para 33,000 genes: nFeature_RNA inferior 200-500, superior 5,000-8,000; nCount_RNA superior en decenas de miles. percent.mt es independiente del número de características (es una fracción de UMIs), así que un límite superior de 10-20% está determinado por el tejido, no por el tamaño. Siempre inspecciona la distribución antes de fijar cualquier número.

DEMOSTRACIÓN EN VIVO:

p1 <- VlnPlot(sobj, "nFeature_RNA", group.by="donor_id") + ggtitle("500-gene panel")
print(p1)

4.5.8 Paso 2.8 - Fallo silencioso: percent.mt como fracción vs como porcentaje

Un fallo silencioso común: alguien copió un umbral de un tutorial que usaba percent.mt expresado como FRACCIÓN (0 a 1), pero Seurat devuelve percent.mt como PORCENTAJE (0 a 100). El filtro parece correcto y se ejecuta sin error, pero elimina casi todas las células.

# >>> DELIBERATE ERROR - read the message, then continue. <<<
try({
  # This filter expects 5 percent (i.e. percent.mt < 5). Written as 0.05 it
  # means "less than 0.05 percent", which removes almost everything.
  sobj_fraction <- subset(sobj, subset = percent.mt < 0.05)
  cat("Cells passing 'percent.mt < 0.05' filter :", ncol(sobj_fraction), "\n")
})

Cells passing 'percent.mt < 0.05' filter : 22 

cat("Compare to the correct threshold (percent.mt < 5):\n")

Compare to the correct threshold (percent.mt < 5):

sobj_pct <- subset(sobj, subset = percent.mt < 5)
cat("Cells passing 'percent.mt < 5' filter    :", ncol(sobj_pct), "\n")

Cells passing 'percent.mt < 5' filter    : 1423 

TipPREGUNTA 2.8

¿La distribución de percent.mt de tu conjunto de datos se parece más a una fracción o a un porcentaje? Ejecuta summary(sobj$percent.mt) para confirmar y recuerda esta trampa al leer código de otros grupos.

TipRESPUESTA

Es un porcentaje. PercentageFeatureSet devuelve valores en una escala de 0-100. Un conjunto de datos PBMC real típicamente se sitúa entre 0 y 15 por ciento mitocondrial. Si ves valores entre 0 y 1, estás viendo una codificación de fracción (alguien dividió entre 100), y cualquier umbral expresado como porcentaje será incorrecto. Diagnóstico: summary(sobj$percent.mt) - si el máximo es menor a 1, es una fracción; de lo contrario, es un porcentaje.

DEMOSTRACIÓN EN VIVO:

summary(sobj$percent.mt)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   2.788   5.290   6.463   8.853  28.605 

4.5.9 Paso 2.9 - Encuentra el error

Lee el código a continuación CUIDADOSAMENTE antes de ejecutarlo. ¿Qué tiene de incorrecto?

Escribe tu respuesta como un comentario en la siguiente línea.

covid_cells <- subset(sobj_filt, subset = condition == "COVID")

Error in `subset()`:
! No cells found

Tu respuesta:

REVELAR: el valor de condición en este conjunto de datos es “COVID19”, no “COVID”. - subset() devuelve cero células sin ningún error. Siempre inspecciona los valores únicos de una columna categórica antes de subconjuntar sobre ella.

cat("Unique values of condition in this dataset:\n")

Unique values of condition in this dataset:

print(unique(sobj_filt$condition))

[1] "Healthy" "COVID19"

TipConsejo

Hábito: imprime unique() de una columna factor o character antes de referenciar cualquier valor específico en subset() o filter().

4.6 Bloque 3 - Detección de doublets

Objetivo: Identificar y eliminar doublets técnicos que sobreviven a los umbrales de control de calidad. + scDblFinder simula doublets artificiales y puntúa cada célula real contra ellos. El resultado es una etiqueta de clase de doublet y una puntuación numérica por célula.

4.6.1 Paso 3.1 - Ejecutar scDblFinder

library(scDblFinder)

Ambos paquetes ya están cargados desde el Bloque 0; no es necesario volver a llamar library().

• scDblFinder llama a xgboost internamente; las versiones recientes de xgboost emiten advertencias de obsolescencia no relacionadas con nuestro análisis. Envuelve la llamada para mantener la consola enfocada en el resultado real.

set.seed(42)
sce <- SingleCellExperiment(
  assays = list(counts = LayerData(sobj_filt, assay = "RNA", layer = "counts"))
)
sce <- suppressWarnings(scDblFinder(sce))

Creating ~2147 artificial doublets...

Dimensional reduction

Evaluating kNN...

Training model...

iter=0, 78 cells excluded from training.

iter=1, 83 cells excluded from training.

iter=2, 92 cells excluded from training.

Threshold found:0.377

10 (0.4%) doublets called

sobj_filt$scDblFinder.class <- sce$scDblFinder.class
sobj_filt$scDblFinder.score <- sce$scDblFinder.score

cat("Doublet classification:\n")

Doublet classification:

print(table(sobj_filt$scDblFinder.class))

singlet doublet 
   2673      10 

cat("\nDoublet rate:",
    round(mean(sobj_filt$scDblFinder.class == "doublet") * 100, 1), "%\n")

Doublet rate: 0.4 %

4.6.2 Paso 3.2 - ¿Dónde se ubican los doublets marcados en el gráfico de dispersión de QC?

df <- sobj_filt@meta.data
p1 <- ggplot(df, aes(nCount_RNA, nFeature_RNA, color = scDblFinder.class)) +
  geom_point(alpha = 0.5, size = 0.7) +
  scale_color_manual(values = c("singlet" = "grey70", "doublet" = "#991b1b")) +
  labs(title = "Doublets sit on the high diagonal", color = NULL) +
  theme_classic(base_size = 11)

p2 <- ggplot(df, aes(scDblFinder.class, nFeature_RNA, fill = scDblFinder.class)) +
  geom_violin(alpha = 0.7) +
  scale_fill_manual(values = c("singlet" = "grey70", "doublet" = "#991b1b")) +
  labs(title = "Genes detected per class", x = NULL) +
  theme_classic(base_size = 11) + theme(legend.position = "none")

p1 | p2

TipPREGUNTA 3.2

¿Por qué se espera que los doublets NO se ubiquen todos en la parte más alta de nCount_RNA? ¿Qué implica esto sobre usar únicamente umbrales de nCount para eliminar doublets?

TipRESPUESTA

Un doublet de dos células similares (dos células T CD4, dos monocitos) tiene aproximadamente el mismo contenido transcripcional que una célula, solo con mayor captura; dependiendo de la eficiencia de captura y la preparación de la biblioteca, su nCount puede ubicarse en cualquier lugar dentro de la distribución de singletes. Los doublets fáciles de detectar son heterotípicos (T + monocito, B + DC) porque tienen transcriptomas híbridos. Los doublets difíciles de detectar son homotípicos. Los umbrales de nCount solo capturan la cola muy alta. scDblFinder captura los híbridos transcripcionales que los umbrales pasan por alto. Implicación: el filtrado por nCount es necesario pero no suficiente.

DEMOSTRACIÓN EN VIVO:

sobj_filt$is_dbl <- sce$scDblFinder.class == "doublet"
FeatureScatter(sobj_filt, "nCount_RNA", "nFeature_RNA", group.by = "is_dbl")

GRÁFICO / SALIDA: FeatureScatter() coloreado por clase de doublet; los doublets se distribuyen por toda la nube, no solo en la parte superior

4.6.3 Paso 3.3 - Eliminar doublets

cat("Cells before doublet removal:", ncol(sobj_filt), "\n")

Cells before doublet removal: 2683 

sobj_filt <- subset(sobj_filt, subset = scDblFinder.class == "singlet")
cat("Cells after doublet removal :", ncol(sobj_filt), "\n")

Cells after doublet removal : 2673 

4.6.4 Paso 3.4 - Reinspeccionar el objeto después del control de calidad y la eliminación de doublets

Después de cada paso que modifica el objeto, confirma qué tienes ahora. Los errores silenciosos solo salen a la luz cuando el objeto se reinspecciona, no cuando se vuelve a usar.

cat("Object after QC + doublet removal:\n")

Object after QC + doublet removal:

print(sobj_filt)

An object of class Seurat 
524 features across 2673 samples within 2 assays 
Active assay: RNA (500 features, 500 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: ADT
 2 dimensional reductions calculated: pca, umap

cat("\nCells lost   :", ncol(sobj) - ncol(sobj_filt), "\n")

Cells lost   : 327 

cat("Cells kept   :", ncol(sobj_filt), "\n")

Cells kept   : 2673 

cat("Median genes :", median(sobj_filt$nFeature_RNA), "\n")

Median genes : 47 

cat("Median UMIs  :", median(sobj_filt$nCount_RNA), "\n")

Median UMIs  : 452 

Dos vistas de la misma eliminación: el cambio general en la distribución, y el desglose por donante. Un filtro que parece razonable en conjunto todavía puede eliminar casi por completo a un donante; el gráfico de barras por donante es lo que detecta esto antes de que se convierta en una sorpresa en pasos posteriores.

before_after <- bind_rows(
  data.frame(nFeature_RNA = sobj$nFeature_RNA, stage = "Before filtering"),
  data.frame(nFeature_RNA = sobj_filt$nFeature_RNA, stage = "After filtering")
)
before_after$stage <- factor(before_after$stage,
                             levels = c("Before filtering", "After filtering"))

p_before_after <- ggplot(before_after, aes(x = nFeature_RNA, fill = stage)) +
  geom_histogram(bins = 40, alpha = 0.7, position = "identity") +
  labs(title = "nFeature_RNA distribution before vs after filtering",
       x = "nFeature_RNA", y = "Cell count", fill = NULL) +
  theme_classic(base_size = 11)

cell_counts <- data.frame(
  donor_id  = c(sobj$donor_id, sobj_filt$donor_id),
  condition = c(sobj$condition, sobj_filt$condition),
  stage     = rep(c("Before", "After"), c(ncol(sobj), ncol(sobj_filt)))
) %>%
  dplyr::count(donor_id, condition, stage) %>%
  dplyr::mutate(stage = factor(stage, levels = c("Before", "After")))

p_donor_counts <- ggplot(cell_counts, aes(x = donor_id, y = n, fill = stage)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~condition, scales = "free_x") +
  labs(title = "Cell count per donor, before vs after QC and doublet removal",
       x = NULL, y = "Cells", fill = NULL) +
  theme_classic(base_size = 11) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p_before_after

p_donor_counts

TipPREGUNTA 3.4

Observa el gráfico de barras por donante. ¿Algún donante en particular perdió una fracción de células mucho mayor que los demás? Si es así, ¿es eso una variación de control de calidad a nivel de donante, o una señal de que los umbrales del filtro se ajustaron en torno a la mayoría de los donantes a costa de un valor atípico?

TipRESPUESTA

Compara la altura antes/después para cada donante en lugar de la tasa de retención general. Un donante que pierde una fracción notablemente mayor que el resto merece una segunda mirada antes de continuar: podría ser variación genuina biológica o técnica a nivel de donante (un donante con calidad de RNA sistemáticamente menor, más células moribundas, o un lote de procesamiento diferente), o podría significar que los umbrales basados en datos del Paso 2.7 (Sección 4.5.7), calculados con todos los donantes agrupados, caen en un rango que penaliza desproporcionadamente la distribución de un donante. La solución es la misma en cualquier caso: si la pérdida de un donante parece extrema, calcula los umbrales por donante y compara, en lugar de asumir que un único umbral global sirve por igual a todos los donantes. Perder silenciosamente la mayoría de las células de un donante cambia lo que realmente mide cada comparación posterior (por condición, por severidad).

DEMOSTRACIÓN EN VIVO:

cell_counts %>% group_by(donor_id) %>%
  summarise(before = sum(n[stage=="Before"]),
            after  = sum(n[stage=="After"]),
            pct_kept = round(100*after/before, 1))

# A tibble: 7 × 4
  donor_id before after pct_kept
  <chr>     <int> <int>    <dbl>
1 Donor01     350   300     85.7
2 Donor02     300   254     84.7
3 Donor03     350   300     85.7
4 Donor04     450   404     89.8
5 Donor05     450   403     89.6
6 Donor06     550   504     91.6
7 Donor07     550   508     92.4

4.7 Bloque 4 - Normalización, reducción de dimensionalidad, clustering (RNA)

Objetivo: Llevar el ensayo RNA desde conteos crudos hasta un UMAP agrupado (clustered). La normalización de ADT deliberadamente NO se realiza aquí; pertenece al Bloque 6 (Sección 5.1) donde comienza el flujo de trabajo de CITE-seq. Mantén las modalidades separadas hasta WNN.

4.7.1 Paso 4.1 - Confirmar que la capa de conteos contiene conteos enteros crudos

Un fallo silencioso común: un objeto llega con la capa de datos copiada en la capa de conteos (ya normalizada logarítmicamente). NormalizeData entonces se ejecuta sobre valores ya normalizados. La corrección debe hacerse aquí, antes de cualquier llamada de normalización.

raw <- LayerData(sobj_filt, assay = "RNA", layer = "counts")
vals <- raw@x  # non-zero values
cat("Min non-zero count :", round(min(vals), 4), "\n")

Min non-zero count : 1 

cat("Max count          :", round(max(vals), 2), "\n")

Max count          : 145 

cat("All integers?      :", all(vals == round(vals)), "\n")

All integers?      : TRUE 

if (!all(vals == round(vals))) {
  cat("\nWARNING: counts layer contains non-integer values.\n")
  cat("This is not raw counts. Do NOT run NormalizeData on it.\n")
  cat("Obtain the original count matrix before continuing.\n")
} else {
  cat("\nCounts layer looks like raw counts. Safe to normalize.\n")
}

Counts layer looks like raw counts. Safe to normalize.

# >>> DELIBERATE ERROR - read the message, then continue. <<<
# What does this check look like on a CORRUPTED counts layer?
# Simulate: a well-meaning collaborator copied the data layer back into counts
# "for consistency". Run the integrity check on that simulated layer and
# compare with the real one above.
try({
  # Build a fake assay where counts holds log-normalized values
  fake_data   <- LayerData(sobj_filt, assay = "RNA", layer = "counts")
  fake_data@x <- log1p(fake_data@x / 1e4)  # roughly what LogNormalize produces
  fake_assay  <- CreateAssay5Object(counts = fake_data)
  fake_vals   <- LayerData(fake_assay, layer = "counts")@x

  cat("\n--- Integrity check on a CORRUPTED counts layer ---\n")
  cat("Min non-zero :", round(min(fake_vals), 4), "\n")
  cat("Max          :", round(max(fake_vals), 4), "\n")
  cat("All integers?:", all(fake_vals == round(fake_vals)), "\n")
  cat("This is what a contaminated counts layer looks like.\n")
  cat("If you see this on a real object, STOP and ask for the raw matrix.\n")
})

--- Integrity check on a CORRUPTED counts layer ---
Min non-zero : 1e-04 
Max          : 0.0144 
All integers?: FALSE 
This is what a contaminated counts layer looks like.
If you see this on a real object, STOP and ask for the raw matrix.

TipPREGUNTA 4.1

Si la verificación de integridad falla en un objeto real heredado, ¿qué le pides al colaborador? ¿Cuál es el mínimo absoluto que necesitas para reiniciar el análisis desde un estado limpio?

TipRESPUESTA

Pide la salida cruda de CreateSeuratObject(), O la salida original de 10x Genomics (filtered_feature_bc_matrix), O el .rds original guardado antes de que se llamara por primera vez a NormalizeData. Necesidad mínima: la matriz de conteos crudos con nombres de células y características que coincidan con el resto de los metadatos. Todo lo posterior (valores normalizados, PCA, UMAP, clustering, anotación) puede regenerarse. Sin los conteos crudos no puedes verificar ninguna afirmación cuantitativa. Deja explícito: todo proyecto orientado a publicación debe preservar un checkpoint en CreateSeuratObject, antes de cualquier normalización.

DEMOSTRACIÓN EN VIVO:

is_integer_counts <- function(layer) {
  v <- layer@x
  is.numeric(v) && all(v >= 0) && all(v == round(v))
}
is_integer_counts(LayerData(sobj_filt, assay="RNA", layer="counts"))

[1] TRUE

4.7.2 Paso 4.2 - Normalización de RNA: LogNormalize

sobj_filt <- NormalizeData(
  sobj_filt,
  normalization.method = "LogNormalize",
  scale.factor         = 10000
)

Normalizing layer: counts

cat("RNA normalized. Layer 'data' now populated.\n")

RNA normalized. Layer 'data' now populated.

cat("Layers present:", paste(SafeLayers(sobj_filt[["RNA"]]), collapse = ", "), "\n")

Layers present: counts, data, scale.data 

Verificar: los valores normalizados logarítmicamente deben ser no negativos

data_layer <- LayerData(sobj_filt, assay = "RNA", layer = "data")
cat("Any negative values in RNA data layer:", any(data_layer < 0), "\n")

Any negative values in RNA data layer: FALSE 

cat("Max value in RNA data layer          :",
    round(max(data_layer@x), 2), "\n")

Max value in RNA data layer          : 8.03 

TipPREGUNTA 4.2

¿Por qué se aplica LogNormalize a RNA pero se usa CLR (margin = 2) para ADT? ¿Qué propiedad de cada modalidad impulsa la diferencia?

TipRESPUESTA

RNA: tamaño de biblioteca variable por célula (decenas a decenas de miles de UMIs), conteos de genes fuertemente sesgados a la derecha, muy disperso. LogNormalize divide los conteos de cada célula entre su total, multiplica por 10,000, y luego aplica log1p. ADT: la carga de anticuerpos por célula es mucho más uniforme entre células (cada célula recibió la misma mezcla de tinción), las proteínas no son dispersas, y la señal significativa es la abundancia relativa de una proteína frente a otras dentro de la misma célula. CLR (log-ratio centrado) con margin=2 trata los conteos de proteínas de cada célula como una composición y normaliza dentro de la célula, centrando en cero. margin=1 normalizaría entre células por proteína, lo que elimina la señal biológica de qué células tiñen positivo.

DEMOSTRACIÓN EN VIVO:

summary(LayerData(sobj_filt, assay="RNA", layer="counts")@x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  1.000   1.000   4.000   9.541  14.000 145.000 

summary(LayerData(sobj_filt, assay="ADT", layer="counts")@x)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1.00    1.00    1.00   10.54    2.00  294.00 

Nota: conteos de RNA típicamente 0-200s; conteos de ADT típicamente 0-1000s.

4.7.3 Paso 4.3 - Genes altamente variables

sobj_filt <- FindVariableFeatures(
  sobj_filt,
  selection.method = "vst",
  nfeatures        = 2000
)

Finding variable features for layer counts

top10 <- head(VariableFeatures(sobj_filt), 10)
cat("Top 10 most variable genes:\n")

Top 10 most variable genes:

print(top10)

 [1] "ITGA2B" "PF4"    "GP1BB"  "TUBB1"  "PPBP"   "FCGR3A" "IRF8"   "MS4A7" 
 [9] "LYZ"    "CCR7"  

cat("\nTotal variable features selected:", length(VariableFeatures(sobj_filt)), "\n")

Total variable features selected: 500 

• Desafío: ¿qué fracción de todos los genes se selecciona como variable?

all_genes <- nrow(sobj_filt[["RNA"]])
n_var     <- length(VariableFeatures(sobj_filt))
cat("Variable fraction:", round(n_var / all_genes * 100, 1), "% of all genes\n")

Variable fraction: 100 % of all genes

4.7.4 Paso 4.4 - Escalar y ejecutar PCA

vars.to.regress elimina el efecto lineal de percent.mt de cada gen antes de escalar. Esto reduce la influencia del estrés/calidad celular sobre los componentes principales.

sobj_filt <- ScaleData(
  sobj_filt,
  vars.to.regress = "percent.mt",
  verbose         = FALSE
)

sobj_filt <- RunPCA(sobj_filt, npcs = 30, verbose = FALSE)

¿Dónde viven las coordenadas de los PC, ahora que el PCA realmente existe? El Paso 1.6 (Sección 4.4.7) confirmó que las reducciones estaban vacías antes de este punto; aquí están los tres patrones de acceso en la práctica.

cat("Reductions now available:", paste(SafeReductions(sobj_filt), collapse = ", "), "\n\n")

Reductions now available: pca, umap 

cat("Cell embeddings (cells x dims matrix), first 6 cells:\n")

Cell embeddings (cells x dims matrix), first 6 cells:

print(head(sobj_filt@reductions$pca@cell.embeddings))

                 PC_1       PC_2       PC_3     PC_4       PC_5         PC_6
CELL000011 -0.3873038  0.5668010 -0.8098380 5.706257 -8.1384301  3.928232203
CELL000021 -3.1533253 -7.3538408 -0.9866441 3.724133  2.2049994 -0.217801983
CELL000031 -2.5582802  2.4435804  3.4516301 5.465073  1.2167644 -0.002586255
CELL000041  6.4807853 -0.3325541 -0.3057807 1.854186  1.4233016  0.228833361
CELL000051  6.8346946 -0.3618520 -0.3202430 1.728192 -0.1796808 -2.190339236
CELL000061  7.3818265 -0.8512410  0.1201711 1.197970  0.7898951 -0.202019481
                 PC_7       PC_8       PC_9      PC_10      PC_11       PC_12
CELL000011  1.5609123 -1.8382579 -1.4743282 -0.2543423 -1.6961021  1.45573562
CELL000021  0.8133849 -1.1628040 -1.0338714 -0.7173800 -1.6887688 -0.97375034
CELL000031  0.6921562 -1.2557761  0.1039847 -0.2911183  0.2403640  2.13562491
CELL000041  0.6663791 -1.6453516  1.1500996  0.7569188  1.5069090  0.05940694
CELL000051  1.7836942  0.9288688  0.7521531  1.1371763 -0.6041043 -0.31607613
CELL000061 -0.8747301  1.8619551  0.1415833 -2.5307814 -0.3209138  2.32895553
                PC_13         PC_14      PC_15      PC_16      PC_17
CELL000011 -1.1607951 -2.1497776263 -1.6291037  0.2034311  0.3943008
CELL000021 -1.6112321  1.2585421285 -1.8412905 -0.2469064  0.8337801
CELL000031  0.7086251 -0.0002333597 -2.6518755  0.1517192  0.4822051
CELL000041  0.8131166 -0.4732009502 -1.0927415  0.2928089 -0.9084758
CELL000051 -1.0054588 -0.7089813023 -0.8314306  1.2667540 -1.7013513
CELL000061 -0.5625249 -3.3401948895  2.6490380  1.9061582 -0.6421049
                 PC_18      PC_19      PC_20      PC_21      PC_22      PC_23
CELL000011  2.47724967  0.6225494 -1.9762847  1.2580835 -0.7622494  1.3256096
CELL000021  0.43637139 -2.5473309 -1.6926438 -0.8692788 -1.2675537 -2.5943491
CELL000031  3.24255953  0.2145089  0.3356989 -2.5228945 -0.1019087 -2.3193680
CELL000041 -0.22792707  1.1687828 -1.7467221  0.8999425 -1.2552311  0.2056092
CELL000051  1.50945823 -0.3530648 -0.9523098  1.2861083  0.7915990 -2.0324108
CELL000061  0.03839866  1.2027816  2.5029708 -0.5543171 -2.1359935 -5.4095704
                PC_24      PC_25      PC_26      PC_27       PC_28      PC_29
CELL000011 -2.2012503 -1.8859956  2.4210827  2.4203230  2.33220401 -0.9544978
CELL000021 -0.8276845 -2.6342547  1.4857436  2.4008630 -2.89468619  1.6235719
CELL000031  1.8525981 -1.1689952 -0.9711090  0.6461746  0.09605737  1.5547858
CELL000041  0.7222750  2.3428183  0.3006477 -0.9846822 -0.21528628  0.4340726
CELL000051  0.2612944  0.3925716 -0.2098121 -0.5995705  1.35395647  0.3964083
CELL000061  3.4441244 -3.0830608  2.3359469 -0.2317923 -1.88698247 -1.0545951
                PC_30
CELL000011 -0.8733779
CELL000021 -0.8554295
CELL000031 -1.3149205
CELL000041 -0.1685520
CELL000051  2.4027365
CELL000061  2.2413082

cat("\nSame thing via the recommended accessor:\n")

Same thing via the recommended accessor:

print(head(Embeddings(sobj_filt, reduction = "pca")))

                 PC_1       PC_2       PC_3     PC_4       PC_5         PC_6
CELL000011 -0.3873038  0.5668010 -0.8098380 5.706257 -8.1384301  3.928232203
CELL000021 -3.1533253 -7.3538408 -0.9866441 3.724133  2.2049994 -0.217801983
CELL000031 -2.5582802  2.4435804  3.4516301 5.465073  1.2167644 -0.002586255
CELL000041  6.4807853 -0.3325541 -0.3057807 1.854186  1.4233016  0.228833361
CELL000051  6.8346946 -0.3618520 -0.3202430 1.728192 -0.1796808 -2.190339236
CELL000061  7.3818265 -0.8512410  0.1201711 1.197970  0.7898951 -0.202019481
                 PC_7       PC_8       PC_9      PC_10      PC_11       PC_12
CELL000011  1.5609123 -1.8382579 -1.4743282 -0.2543423 -1.6961021  1.45573562
CELL000021  0.8133849 -1.1628040 -1.0338714 -0.7173800 -1.6887688 -0.97375034
CELL000031  0.6921562 -1.2557761  0.1039847 -0.2911183  0.2403640  2.13562491
CELL000041  0.6663791 -1.6453516  1.1500996  0.7569188  1.5069090  0.05940694
CELL000051  1.7836942  0.9288688  0.7521531  1.1371763 -0.6041043 -0.31607613
CELL000061 -0.8747301  1.8619551  0.1415833 -2.5307814 -0.3209138  2.32895553
                PC_13         PC_14      PC_15      PC_16      PC_17
CELL000011 -1.1607951 -2.1497776263 -1.6291037  0.2034311  0.3943008
CELL000021 -1.6112321  1.2585421285 -1.8412905 -0.2469064  0.8337801
CELL000031  0.7086251 -0.0002333597 -2.6518755  0.1517192  0.4822051
CELL000041  0.8131166 -0.4732009502 -1.0927415  0.2928089 -0.9084758
CELL000051 -1.0054588 -0.7089813023 -0.8314306  1.2667540 -1.7013513
CELL000061 -0.5625249 -3.3401948895  2.6490380  1.9061582 -0.6421049
                 PC_18      PC_19      PC_20      PC_21      PC_22      PC_23
CELL000011  2.47724967  0.6225494 -1.9762847  1.2580835 -0.7622494  1.3256096
CELL000021  0.43637139 -2.5473309 -1.6926438 -0.8692788 -1.2675537 -2.5943491
CELL000031  3.24255953  0.2145089  0.3356989 -2.5228945 -0.1019087 -2.3193680
CELL000041 -0.22792707  1.1687828 -1.7467221  0.8999425 -1.2552311  0.2056092
CELL000051  1.50945823 -0.3530648 -0.9523098  1.2861083  0.7915990 -2.0324108
CELL000061  0.03839866  1.2027816  2.5029708 -0.5543171 -2.1359935 -5.4095704
                PC_24      PC_25      PC_26      PC_27       PC_28      PC_29
CELL000011 -2.2012503 -1.8859956  2.4210827  2.4203230  2.33220401 -0.9544978
CELL000021 -0.8276845 -2.6342547  1.4857436  2.4008630 -2.89468619  1.6235719
CELL000031  1.8525981 -1.1689952 -0.9711090  0.6461746  0.09605737  1.5547858
CELL000041  0.7222750  2.3428183  0.3006477 -0.9846822 -0.21528628  0.4340726
CELL000051  0.2612944  0.3925716 -0.2098121 -0.5995705  1.35395647  0.3964083
CELL000061  3.4441244 -3.0830608  2.3359469 -0.2317923 -1.88698247 -1.0545951
                PC_30
CELL000011 -0.8733779
CELL000021 -0.8554295
CELL000031 -1.3149205
CELL000041 -0.1685520
CELL000051  2.4027365
CELL000061  2.2413082

cat("\nFeature loadings (genes x dims matrix), first 6 genes:\n")

Feature loadings (genes x dims matrix), first 6 genes:

print(head(sobj_filt@reductions$pca@feature.loadings))

              PC_1         PC_2          PC_3        PC_4        PC_5
ITGA2B -0.01687587  0.005211474  6.426538e-03 -0.03257683 -0.08838266
PF4    -0.01338218  0.002973087 -5.743934e-05 -0.03797926 -0.08850547
GP1BB  -0.01004316 -0.002913585  3.037135e-03 -0.03195500 -0.08878522
TUBB1  -0.01120083  0.007213298 -8.998449e-04 -0.04374412 -0.09544160
PPBP   -0.01075531  0.008453901  1.316460e-03 -0.03349028 -0.08304103
FCGR3A -0.02012235  0.005710962  3.067386e-03  0.02628131 -0.09162497
             PC_6       PC_7       PC_8         PC_9        PC_10         PC_11
ITGA2B -0.3160597  0.2327805 0.02334507 -0.005077348  0.003825966 -0.0020273926
PF4    -0.3246330  0.2459271 0.02262608  0.013201597  0.008296717  0.0165600177
GP1BB  -0.3265169  0.2375194 0.01481990 -0.013251636 -0.007482809  0.0002795763
TUBB1  -0.3230936  0.2269098 0.02363318  0.010104419  0.007482969  0.0088062067
PPBP   -0.3303714  0.2418020 0.01789425  0.014695685 -0.005422940 -0.0252453357
FCGR3A -0.2025325 -0.3795360 0.02596463  0.013722655  0.046430924 -0.0217529256
               PC_12       PC_13        PC_14        PC_15        PC_16
ITGA2B -0.0008544128 0.013983617 -0.001266814 -0.001163291 -0.009135415
PF4     0.0307799520 0.009394718 -0.007484448  0.006896253  0.007304150
GP1BB   0.0001825089 0.019984007 -0.003350069  0.020869346  0.013690731
TUBB1   0.0146127690 0.010005869 -0.007599485  0.017866731  0.028596820
PPBP    0.0178133160 0.010643787 -0.009738071  0.024940378 -0.024243129
FCGR3A -0.0002877477 0.009727228 -0.007321154  0.010032782 -0.031800122
             PC_17        PC_18         PC_19        PC_20         PC_21
ITGA2B 0.007393687  0.002171734 -0.0044473244  0.011770104  0.0144165173
PF4    0.006250811 -0.003313275 -0.0031773166 -0.025356773  0.0009774063
GP1BB  0.021720083  0.027710983  0.0085306345 -0.014444921 -0.0005432358
TUBB1  0.002864765 -0.002214925 -0.0167476185 -0.005458695  0.0074148475
PPBP   0.014494628  0.006629371 -0.0151095545 -0.008093622  0.0285454781
FCGR3A 0.028889905 -0.015324592 -0.0007436423 -0.008751208 -0.0044508591
              PC_22        PC_23        PC_24        PC_25        PC_26
ITGA2B  0.005491894  0.004131900 -0.007293246 -0.005144066 -0.021643538
PF4    -0.006284602  0.011715859  0.002935507  0.008369986  0.003068438
GP1BB   0.011093066  0.002879251 -0.012739674  0.006344725 -0.012250609
TUBB1  -0.002009008 -0.002331423 -0.016289368  0.019886108 -0.007190101
PPBP    0.010497415  0.023477789 -0.011954778  0.006518588 -0.006737003
FCGR3A -0.021800046 -0.007231774 -0.011825965 -0.016529404  0.005007931
              PC_27        PC_28        PC_29        PC_30
ITGA2B -0.010288036  0.020678739  0.011486635 -0.007136209
PF4     0.001360859 -0.005729753  0.014267929 -0.010136190
GP1BB   0.007057034 -0.004169101  0.005576511 -0.008762779
TUBB1  -0.001575694  0.019432641 -0.013219720 -0.011897489
PPBP    0.001144814  0.004305241 -0.019385258 -0.001162891
FCGR3A -0.016231220  0.002229713 -0.016115785 -0.008175363

Inspecciona las cargas (loadings) principales de PC1 y PC2

cat("\nTop genes loading on PC1:\n")

Top genes loading on PC1:

print(head(sobj_filt@reductions$pca@feature.loadings[
  order(abs(sobj_filt@reductions$pca@feature.loadings[,1]),
        decreasing = TRUE), 1], 10))

      CD4      CCR7      IL7R    HAVCR2      CD8B     PDCD1     IL2RA      CD3D 
0.2342448 0.2338181 0.2331904 0.2314558 0.2309892 0.2305871 0.2297123 0.2296829 
     SELL      CD8A 
0.2295633 0.2295542 

TipPREGUNTA 4.4c

Heredas un objeto donde Reductions(sobj) lista “pca” y “umap” pero Layers(sobj[["RNA"]]) muestra solo “counts”. La capa data no está. ¿Puedes confiar en el UMAP? ¿Cuál es tu próximo paso?

TipRESPUESTA

No, no puedes confiar en él. El PCA se calculó a partir de la capa data (normalizada logarítmicamente), y el UMAP a partir del PCA. Si la capa data ha sido eliminada, la entrada previa al PCA desapareció, así que no puedes verificar que el PCA usó los valores, la normalización o las características correctas. Las reducciones sin su capa de origen no son verificables. Dos próximos pasos: (1) pedirle al colaborador el objeto antes de que se eliminara la capa data, O (2) volver a ejecutar NormalizeData, FindVariableFeatures, ScaleData, RunPCA a partir de la capa de conteos y comparar tu nuevo PCA con el heredado. Un desacuerdo sustancial significa que el UMAP heredado es sospechoso.

DEMOSTRACIÓN EN VIVO:

Reductions(sobj_filt)

[1] "pca"  "umap"

Layers(sobj_filt[["RNA"]])

[1] "counts"     "data"       "scale.data"

# Recompute and compare:
sobj_check <- NormalizeData(sobj_filt, verbose = FALSE)
sobj_check <- FindVariableFeatures(sobj_check, verbose = FALSE)
sobj_check <- ScaleData(sobj_check, verbose = FALSE)
sobj_check <- RunPCA(sobj_check, npcs = 30, verbose = FALSE)

GRÁFICO / SALIDA: UMAPs lado a lado coloreados por cluster si el PCA recalculado está disponible.

El codo es donde agregar más PCs explica poca varianza adicional.

Usa esto para establecer dims en FindNeighbors y RunUMAP.

ElbowPlot(sobj_filt, ndims = 30) +
  geom_vline(xintercept = 20, linetype = "dashed", color = "#991b1b") +
  annotate("text", x = 21, y = 2.5, label = "~20 PCs", hjust = 0,
           color = "#991b1b") +
  ggtitle("Elbow Plot: variance explained per PC")

TipPREGUNTA 4.4

¿Cómo elegirías el número de PCs en un conjunto de datos real donde el gráfico de codo no tiene una curva pronunciada?

TipRESPUESTA

Combina cuatro señales: (1) varianza acumulada explicada, apuntando a 70-90%; (2) prueba de permutación JackStraw, seleccionar PCs significativos al alfa elegido; (3) estabilidad del clustering entre dimensiones (reejecutar FindClusters en dims=10, 15, 20, 25 y calcular el ARI entre particiones); (4) coherencia biológica: ¿se separan todas las poblaciones esperadas? Reporta la elección y el análisis de sensibilidad, no un número mágico.

DEMOSTRACIÓN EN VIVO:

# JackStraw is slow on full data; usable on subsamples for the same conclusion
sobj_filt <- JackStraw(sobj_filt, num.replicate = 50, dims = 25, verbose = FALSE)
sobj_filt <- ScoreJackStraw(sobj_filt, dims = 1:25)
JackStrawPlot(sobj_filt, dims = 1:25)

Warning: Removed 10928 rows containing missing values or values outside the scale range
(`geom_point()`).

GRÁFICO / SALIDA: JackStrawPlot: los PCs por encima de la diagonal son significativos

4.7.5 Paso 4.4a - Verificación de sanidad: ¿hicieron HVG y PCA lo que pediste?

• Fallo silencioso A. FindVariableFeatures(nfeatures = 2000) en un panel de 500 genes devuelve 500 sin ninguna advertencia. El análisis de HVG fue una operación nula (no-op).

n_features  <- nrow(sobj_filt[["RNA"]])
n_hvg_asked <- 2000
n_hvg_got   <- length(VariableFeatures(sobj_filt))
cat("Features in assay  :", n_features, "\n")

Features in assay  : 500 

cat("HVG requested      :", n_hvg_asked, "\n")

HVG requested      : 2000 

cat("HVG returned       :", n_hvg_got, "\n")

HVG returned       : 500 

if (n_hvg_got < n_hvg_asked) {
  cat(">> nfeatures > total features. HVG is the full panel; the call",
      "selected nothing.\n")
}

>> nfeatures > total features. HVG is the full panel; the call selected nothing.

• Fallo silencioso B. RunPCA devuelve el número de PCs que solicitas, incluso cuando solo una fracción lleva señal. Inspecciona la varianza de la cola para detectar PCs muertos.

pca_sdev <- sobj_filt@reductions$pca@stdev
pca_var  <- pca_sdev ^ 2
cat("\nVariance explained by last 5 PCs (out of", length(pca_var), "):\n")

Variance explained by last 5 PCs (out of 30 ):

print(round(tail(pca_var / sum(pca_var) * 100, 5), 3))

[1] 1.804 1.795 1.789 1.785 1.769

cat("If the tail is below ~0.5 percent each, those PCs are mostly noise.\n")

If the tail is below ~0.5 percent each, those PCs are mostly noise.

TipPREGUNTA 4.4a

¿Cuántos PCs llevan más del 1 por ciento de la varianza total?

La línea de abajo lo calcula; usa el resultado como una cota inferior para dims= más adelante.

n_pcs_above_1pct <- sum((pca_var / sum(pca_var)) > 0.01)
cat("\nNumber of PCs with > 1% variance explained:", n_pcs_above_1pct, "\n")

Number of PCs with > 1% variance explained: 30 

cat("Recommended lower bound for dims=: 1:", n_pcs_above_1pct, sep = "", "\n")

Recommended lower bound for dims=: 1:30

cat("(The elbow plot above gave ~20; this floor confirms or contradicts it.)\n")

(The elbow plot above gave ~20; this floor confirms or contradicts it.)

TipRESPUESTA

Calculado en el script. Típicamente 10-15 PCs llevan >1% de varianza en este conjunto de datos; el codo da aproximadamente el mismo número. Si difieren por más de un factor de 2, o bien el codo se está leyendo mal o el conjunto de datos tiene una estructura de varianza inusual.

DEMOSTRACIÓN EN VIVO:

pca_var <- sobj_filt@reductions$pca@stdev^2
pct <- pca_var / sum(pca_var) * 100
sum(pct > 1)

[1] 30

plot(pct, type="b", xlab="PC", ylab="% variance")

GRÁFICO / SALIDA: Gráfico de dispersión de varianza por PC; el punto de aplanamiento es el piso para dims

4.7.6 Paso 4.4b - Decidir si integrar (Harmony)

Antes de agrupar (clustering), pregunta: ¿se separan los donantes en el PCA? Si es así, el UMAP estará dominado por el donante, no por la biología, y el clustering codificará el lote (batch).

p_pca_donor <- DimPlot(sobj_filt, reduction = "pca", group.by = "donor_id") +
  ggtitle("PCA by donor")
p_pca_cond  <- DimPlot(sobj_filt, reduction = "pca", group.by = "condition") +
  ggtitle("PCA by condition")
p_pca_donor | p_pca_cond

• Diagnóstico rápido de R^2: ¿cuánto de PC1 y PC2 se explica por el donante?

pc_coords <- Embeddings(sobj_filt, reduction = "pca")[, 1:2]
for (pc in 1:2) {
  fit <- summary(lm(pc_coords[, pc] ~ sobj_filt$donor_id))
  cat(sprintf("PC%d ~ donor_id R^2: %.3f\n", pc, fit$r.squared))
}

PC1 ~ donor_id R^2: 0.065
PC2 ~ donor_id R^2: 0.007

TipConsejo

Regla práctica de decisión: - R^2 < 0.10 : efecto de donante mínimo, no se necesita integración - R^2 0.10 a 0.30: límite, considerar integración si los tipos celulares están mezclados entre donantes - R^2 > 0.30 : el donante domina, se recomienda integración

En este conjunto de datos, el efecto de donante está entre bajo y límite, así que el clustering y la anotación a través del Bloque 5 (Sección 4.8) continúan sobre el PCA/UMAP simple sin Harmony.

El Paso 5.4e (Sección 4.8.9) ejecuta Harmony de todas formas. Una vez que los tipos celulares están anotados, para construir un panel de comparación lado a lado: ver que Harmony apenas cambia la disposición en datos en los que ya confías es lo que te da confianza al leer la misma comparación en un nuevo conjunto de datos donde aún no conoces la respuesta.

Sintaxis de referencia (paquete harmony actual, método de objeto Seurat):

library(harmony)
# NOT RUN
sobj_harm <- RunHarmony(sobj_filt, group.by.vars = "donor_id",
                          reduction      = "pca",
                           reduction.save = "harmony",
                           verbose        = FALSE)
   DimPlot(sobj_harm, reduction = "harmony", group.by = "donor_id") +
     ggtitle("After Harmony (donor)")

En los pasos posteriores FindNeighbors / RunUMAP, cambia reduction = "pca" a reduction = "harmony" (dims permanece igual; Harmony devuelve la misma dimensionalidad que la reducción de entrada).

TipPREGUNTA 4.4b

¿A qué R^2 de donante cambiarías a Harmony en tus propios datos? El bloque anterior ya imprimió el R^2 para PC1 y PC2; aplica esta regla directamente:

r2_pc1 <- summary(lm(pc_coords[, 1] ~ sobj_filt$donor_id))$r.squared
verdict <- if (r2_pc1 < 0.10) {
  "NO integration needed"
} else if (r2_pc1 < 0.30) {
  "BORDERLINE: integrate only if cell types are also separated by donor"
} else {
  "INTEGRATE (Harmony or similar)"
}
cat(sprintf("\nDecision for this dataset: R^2(PC1)= %.3f -> %s\n", r2_pc1, verdict))

Decision for this dataset: R^2(PC1)= 0.065 -> NO integration needed

cat("\nTrade-off of integration:\n")

Trade-off of integration:

cat("  Pro: removes technical donor variance; cell types pool across donors.\n")

  Pro: removes technical donor variance; cell types pool across donors.

cat("  Con: may remove TRUE biological inter-donor variance (e.g., one donor\n")

  Con: may remove TRUE biological inter-donor variance (e.g., one donor

cat("       genuinely lacks a cell type). Always cross-check post-integration\n")

       genuinely lacks a cell type). Always cross-check post-integration

cat("       that cell type proportions per donor remain plausible.\n")

       that cell type proportions per donor remain plausible.

TipRESPUESTA

Umbrales (en el script): R^2 < 0.10 = sin integración; 0.10-0.30 = límite, integrar solo si los tipos celulares también están separados por donante; >0.30 = integrar. Compensación: la integración elimina la varianza técnica de donante para que los tipos celulares se agrupen entre donantes. También elimina la varianza biológica VERDADERA entre donantes: si un donante genuinamente carece de un tipo celular, la integración puede fusionar sus células con células de otros donantes que SÍ lo tienen, ocultando la diferencia. Después de integrar, siempre verifica que las proporciones de tipo celular por donante sigan siendo plausibles. El R^2 de este conjunto de datos cae en el rango bajo a límite, así que el clustering y la anotación continúan sin Harmony a través del Bloque 5 (Sección 4.8). El Paso 5.4e (Sección 4.8.9) ejecuta Harmony de todas formas y el Paso 5.4f (Sección 4.8.10) construye una comparación de 4 paneles para que los estudiantes vean directamente si la integración habría cambiado algo, en lugar de tomar el umbral de R^2 por fe.

DEMOSTRACIÓN EN VIVO:

pc_coords <- Embeddings(sobj_filt, "pca")[, 1:5]
sapply(1:5, function(i)
  summary(lm(pc_coords[, i] ~ sobj_filt$donor_id))$r.squared)

[1] 0.064676642 0.007479984 0.042866818 0.793805724 0.035593229

GRÁFICO / SALIDA: DimPlot(sobj_filt, reduction='pca', group.by='donor_id') lado a lado con group.by='condition'; ver también la comparación de 4 paneles en el Paso 5.4f (Sección 4.8.10).

4.7.7 Paso 4.4c - Encuentra el error

Lee el código a continuación cuidadosamente. Hay dos cosas incorrectas. Encuéntralas ambas antes de ejecutar nada.

sobj_filt <- FindNeighbors(sobj_filt, dims = 0:20, verbose = FALSE)
sobj_filt <- RunUMAP(sobj_filt, dims = c(1, 2, 3, 5, 7, 11, 13))

Tu respuesta:

REVELAR:

1. dims = 0:20 incluye el PC 0, que no existe. La indexación de PCA comienza en 1. FindNeighbors descartará silenciosamente el 0 o dará error dependiendo de la versión de Seurat. Usa siempre 1:N.
2. Al pasar dims como un vector no contiguo, se omiten PCs intermedios. UMAP se ejecutará pero el resultado representa solo las 7 dimensiones seleccionadas, no la estructura capturada por el codo. Usa siempre un rango contiguo comenzando desde 1.

4.7.8 Paso 4.5 - FindNeighbors: una discrepancia común de dimensiones

RunPCA se llamó con npcs = 30. Pasar dims = 1:50 a FindNeighbors solicita 50 PCs que no existen.

# >>> DELIBERATE ERROR - read the message, then continue. <<<
try({
  # npcs=30 was used in RunPCA. What happens if we request 50 dims?
  sobj_filt <- FindNeighbors(sobj_filt, dims = 1:50, verbose = FALSE)
})

Error in FindNeighbors.Seurat(sobj_filt, dims = 1:50, verbose = FALSE) : 
  More dimensions specified in dims than have been computed

El error ocurre porque el objeto PCA contiene solo 30 componentes. - dims = 1:50 solicita componentes que no existen. - Causa común: se cambia npcs en RunPCA pero las llamadas posteriores no se actualizan.

sobj_filt <- FindNeighbors(sobj_filt, dims = 1:20, verbose = FALSE)

4.7.9 Paso 4.6 - Sensibilidad de resolución y clustering final

¿Cuánto cambia el clustering entre la resolución 0.2 y 1.2?

res_values <- c(0.2, 0.5, 1.2)
n_clusters <- sapply(res_values, function(r) {
  tmp <- FindClusters(sobj_filt, resolution = r, verbose = FALSE)
  length(unique(tmp$seurat_clusters))
})

cat("Resolution vs number of clusters:\n")

Resolution vs number of clusters:

for (i in seq_along(res_values)) {
  cat("  resolution =", res_values[i], "->", n_clusters[i], "clusters\n")
}

  resolution = 0.2 -> 12 clusters
  resolution = 0.5 -> 12 clusters
  resolution = 1.2 -> 12 clusters

cat("\nUsing resolution = 0.5 for the rest of the analysis.\n")

Using resolution = 0.5 for the rest of the analysis.

sobj_filt <- FindClusters(sobj_filt, resolution = 0.5, verbose = FALSE)
cat("Clusters found:", length(unique(sobj_filt$seurat_clusters)), "\n")

Clusters found: 12 

print(table(sobj_filt$seurat_clusters))

  0   1   2   3   4   5   6   7   8   9  10  11 
527 460 347 221 198 197 165 135 118 116 107  82 

TipPREGUNTA 4.6

Una resolución más baja da menos clusters, más grandes; una resolución más alta da más clusters, más pequeños. Ninguna es intrínsecamente correcta.

• ¿Qué evidencia usas para defender una resolución elegida?

TipRESPUESTA

Tres cosas: (1) coherencia biológica: cada cluster tiene un perfil de marcadores distinguible, defendible frente a la literatura; (2) estabilidad entre resoluciones cercanas: los clusters no deberían fragmentarse dramáticamente ante una perturbación pequeña; el ARI entre resoluciones 0.4 y 0.6 debería ser alto; (3) sanidad de los pasos posteriores: el número de clusters coincide con las expectativas para el tejido (PBMC con 3,000 células: 8-14 clusters es razonable; 25 es demasiado; 4 es muy pocos). Muestra una visualización estilo clustree a través de un barrido de resoluciones y reporta qué resolución y por qué.

DEMOSTRACIÓN EN VIVO:

for (r in c(0.3, 0.5, 0.8, 1.2)) {
  sobj_filt <- FindClusters(sobj_filt, resolution = r, verbose = FALSE)
  cat(sprintf("res=%.1f -> %d clusters\n", r, length(unique(Idents(sobj_filt)))))
}

res=0.3 -> 12 clusters
res=0.5 -> 12 clusters
res=0.8 -> 12 clusters
res=1.2 -> 12 clusters

GRÁFICO / SALIDA: Si clustree está instalado: clustree(sobj_filt, prefix=‘RNA_snn_res.’)

4.7.10 Paso 4.7 - Visualización de UMAP

Tres preguntas para responder inmediatamente después de generar el UMAP:

1. ¿Las células se agrupan por biología o por donante? (verificación de efecto de lote)
2. ¿Los genes marcadores canónicos se mapean a los clusters esperados?
3. ¿La condición (Sano vs COVID-19) muestra alguna estructura espacial? set.seed(42)

• umap.method='uwot' es el valor predeterminado actual pero Seurat imprime un aviso único cuando no se indica explícitamente. Establecerlo silencia el aviso.

sobj_filt <- RunUMAP(
  sobj_filt,
  dims        = 1:20,
  umap.method = "uwot",
  metric      = "cosine",
  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 session

p1 <- DimPlot(sobj_filt, group.by = "seurat_clusters",
              label = TRUE, label.size = 4) +
  NoLegend() + ggtitle("Clusters (resolution = 0.5)")

p2 <- DimPlot(sobj_filt, group.by = "orig.ident") +
  ggtitle("By donor: batch effect?") +
  theme(legend.text = element_text(size = 7))

p1 | p2

p_cond <- DimPlot(sobj_filt, group.by = "condition",
                  cols = c("Healthy" = "#0d7377", "COVID19" = "#991b1b"),
                  pt.size = 0.5) +
  ggtitle("Healthy vs COVID-19")

p_sev <- DimPlot(sobj_filt, group.by = "severity",
                 pt.size = 0.5) +
  ggtitle("COVID-19 severity") +
  theme(legend.text = element_text(size = 8))

p_cond | p_sev

TipPREGUNTA 4.7

Si los donantes se separan visiblemente en el UMAP, ¿es efecto de lote o biología? ¿Qué información adicional te ayudaría a decidir?

TipRESPUESTA

No se puede saber solo con el UMAP. Dos verificaciones: (1) ¿los donantes se agrupan dentro de regiones de tipo celular compartidas o forman sus propias regiones? Los mismos tipos celulares en la misma región del UMAP entre donantes = biología esperada; los mismos tipos celulares en regiones DIFERENTES del UMAP por donante = lote (batch). Colorea el UMAP por tipo celular y por donante en el mismo gráfico. (2) Calcula el R^2 de PC1 contra donor_id (ver 4.4b); si es alto, el lote domina el embedding.

DEMOSTRACIÓN EN VIVO:

DimPlot(sobj_filt, reduction="umap", group.by="donor_id") 

DimPlot(sobj_filt, reduction="umap", group.by="singler_label")

Warning: The following requested variables were not found: singler_label

Error in `[.data.frame`:
! undefined columns selected

GRÁFICO / SALIDA: Los UMAPs lado a lado hacen la decisión obvia en la mayoría de los casos.

4.7.11 Paso 4.8 - Marcadores canónicos de PBMC (RNA)

Coteja la expresión de marcadores con las ubicaciones de los clusters antes de la anotación formal

FeaturePlot(sobj_filt,
            features   = c("CD3E", "CD19", "CD14", "NCAM1", "IL7R", "FCGR3A"),
            ncol       = 3,
            min.cutoff = "q05")

DotPlot(sobj_filt,
        features = c("CD3E", "CD4", "CD8A", "CD19", "MS4A1",
                     "CD14", "FCGR3A", "NCAM1", "PPBP"),
        group.by = "seurat_clusters") +
  RotatedAxis() +
  ggtitle("Marker expression by cluster")

4.7.12 Paso 4.9 - Proporciones de clusters entre condiciones

Este es un resultado preliminar. Debe validarse después de la anotación.

props <- sobj_filt@meta.data %>%
  group_by(condition, seurat_clusters) %>%
  summarise(n = n(), .groups = "drop") %>%
  group_by(condition) %>%
  mutate(prop = n / sum(n))

ggplot(props, aes(x = condition, y = prop, fill = seurat_clusters)) +
  geom_bar(stat = "identity", width = 0.6) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(title    = "Cluster proportions: Healthy vs COVID-19",
       subtitle = "Check per-donor before interpreting any group difference",
       x = "Condition", y = "Proportion", fill = "Cluster") +
  theme_classic(base_size = 12)

4.8 Bloque 5 - Anotación de tipos celulares (basada en RNA)

Objetivo: asignar etiquetas de tipo celular a partir de una referencia curada (SingleR con el Human Primary Cell Atlas), y luego verificar el objeto etiquetado antes de que comience el flujo de trabajo de CITE-seq. La anotación basada únicamente en RNA tiene debilidades conocidas para los subtipos de células T; el Bloque 6 (Sección 5.1) las revisará con proteína.

4.8.1 Paso 5.1 - Ejecutar SingleR contra el Human Primary Cell Atlas

Opción A — Ruta de instalación estándar

👉 Este flujo de trabajo usa la referencia oficial del Human Primary Cell Atlas y es el enfoque confirmatorio recomendado.

cat("Loading reference (Human Primary Cell Atlas)...\n")
ref     <- celldex::HumanPrimaryCellAtlasData()
rna_mat <- LayerData(sobj_filt, assay = "RNA", layer = "data")

cat("Running SingleR...\n")
singler_res <- SingleR(
  test    = rna_mat,
  test    = rna_mat,
  ref     = ref,
  labels  = ref$label.main,
  BPPARAM = BiocParallel::SerialParam()
)

Opción B — Alternativa con archivo RDS pre‑guardado

Si la instalación de celldex u otras dependencias falla:

1. Descarga el archivo pre‑generado Obtén singleR.rds de los materiales del curso proporcionados.
2. Carga el archivo en R:

library(here) # install.packages("here")

here() starts at /private/var/folders/zz/cdcfnwts145c0xjqs_c0y6mm0000gn/T/RtmppqY2dt/file177bf1a6b81e1

singler_res <- readRDS(
  here("checkpoints", "step5.1", "singleR.rds")
)

3. Usa las etiquetas directamente Estas etiquetas pueden integrarse en tu pipeline de análisis como un checkpoint, permitiéndote continuar sin instalar celldex.

👉 Este flujo de trabajo es exploratorio y asegura la reproducibilidad incluso en entornos restringidos (por ejemplo, macOS ARM64 o acceso limitado a internet).

cat("\nCell type distribution:\n")

Cell type distribution:

print(sort(table(singler_res$labels), decreasing = TRUE))

             T_cells             Monocyte              NK_cell 
                 786                  447                  354 
              B_cell                   DC    Endothelial_cells 
                 307                  156                  140 
         Neutrophils           Macrophage     Pre-B_cell_CD34- 
                 138                  103                   71 
    Epithelial_cells            Platelets          Gametocytes 
                  68                   49                   14 
                 CMP            Myelocyte          Hepatocytes 
                  12                   10                    8 
                 MEP              Neurons          Fibroblasts 
                   6                    6                    5 
                  BM     Pro-B_cell_CD34+                  GMP 
                   3                    3                    1 
          HSC_-G-CSF            HSC_CD34+ Neuroepithelial_cell 
                   1                    1                    1 
       Pro-Myelocyte    Tissue_stem_cells 
                   1                    1 

cat("\nLow-confidence cells (NA in pruned labels):",
    sum(is.na(singler_res$pruned.labels)), "\n")

Low-confidence cells (NA in pruned labels): 17 

4.8.2 Paso 5.2 - Graficar una etiqueta antes de que esté en el objeto

Olvidar transferir las etiquetas de SingleR al objeto Seurat antes de llamar a DimPlot(group.by = "singler_label") produce un error opaco. La etiqueta solo existe en el objeto de resultado de SingleR hasta que la copias.

# >>> DELIBERATE ERROR - read the message, then continue. <<<
try({
  # singler_res is a separate object. Labels must be transferred to sobj_filt first.
  DimPlot(sobj_filt, group.by = "singler_label")
})

Warning: The following requested variables were not found: singler_label

Error in `[.data.frame`(data, , group) : undefined columns selected

singler_label todavía no existe en sobj_filt@meta.data. singler_res es un DataFrame (clase de Bioconductor). Las etiquetas deben agregarse explícitamente a los metadatos de Seurat.

sobj_filt$singler_label  <- singler_res$labels
sobj_filt$singler_pruned <- singler_res$pruned.labels

Verificar

cat("singler_label in metadata:", "singler_label" %in% colnames(sobj_filt@meta.data), "\n")

singler_label in metadata: TRUE 

cat("Distribution:\n")

Distribution:

print(sort(table(sobj_filt$singler_label), decreasing = TRUE))

             T_cells             Monocyte              NK_cell 
                 782                  445                  351 
              B_cell                   DC    Endothelial_cells 
                 305                  154                  138 
         Neutrophils           Macrophage     Pre-B_cell_CD34- 
                 136                  102                   71 
    Epithelial_cells            Platelets          Gametocytes 
                  68                   49                   14 
                 CMP            Myelocyte          Hepatocytes 
                  12                   10                    8 
                 MEP              Neurons          Fibroblasts 
                   6                    6                    4 
                  BM     Pro-B_cell_CD34+                  GMP 
                   3                    3                    1 
          HSC_-G-CSF            HSC_CD34+ Neuroepithelial_cell 
                   1                    1                    1 
       Pro-Myelocyte    Tissue_stem_cells 
                   1                    1 

Las células con puntuaciones similares entre múltiples tipos son ambiguas. Una dispersión amplia = alta confianza. Dispersión estrecha = baja confianza. plotScoreHeatmap() dibuja una fila por cada etiqueta de referencia posible: el catálogo COMPLETO de HumanPrimaryCellAtlasData, ~37 filas, sin importar cuántas células se asignaron realmente a cada una. Observa esto antes de que ocurra cualquier limpieza. La mayoría de esas 37 filas mostrarán casi ninguna señal para cualquier célula en este conjunto de datos; ese desorden visual es la evidencia que motiva la consolidación en el siguiente paso.

plotScoreHeatmap(singler_res,
                 main = "SingleR annotation scores: full reference (37 possible labels)")

TipPREGUNTA 5.2

¿Cuántas de las filas en este mapa de calor muestran alguna señal real (un grupo de células con una celda claramente brillante)? ¿Qué te dice el resto de las filas sobre cuántas de las etiquetas impresas arriba son probablemente ruido?

TipRESPUESTA

En este conjunto de datos, típicamente solo entre 5 y 8 de las aproximadamente 37 filas muestran un bloque brillante claro: una separación limpia donde un grupo-columna de células se ilumina intensamente para esa fila y permanece oscuro para el resto. Esas filas corresponden a las poblaciones celulares reales realmente presentes en PBMC: células T, células B, monocitos, células NK, y algunas otras. Las filas restantes son casi uniformemente tenues en todas las células. Una fila tenue no significa que SingleR cometiera un error; significa que ninguna célula de este conjunto de datos puntuó alto contra esa etiqueta de referencia, lo cual es esperado para etiquetas como Hepatocitos o Neuronas en una muestra de sangre. La lectura práctica: si la fila de una etiqueta en este mapa de calor nunca se ilumina intensamente en ningún lugar, cualquier célula asignada a esa etiqueta por el clasificador crudo es sospechosa, y esa es exactamente la razón por la que aparecieron tantos valores distintos en la ‘distribución de tipos celulares’ impresa en el Paso 5.2 (Sección 4.8.2). Las filas tenues son la evidencia visual que motiva colapsar las etiquetas de baja frecuencia juntas en el Paso 5.2b (Sección 4.8.3), en lugar de confiar en todas ellas como poblaciones igualmente reales.

DEMOSTRACIÓN EN VIVO:

Contar etiquetas distintas con al menos una célula asignada con confianza:

label_tab <- sort(table(sobj_filt$singler_label), decreasing = TRUE)
print(label_tab)

             T_cells             Monocyte              NK_cell 
                 782                  445                  351 
              B_cell                   DC    Endothelial_cells 
                 305                  154                  138 
         Neutrophils           Macrophage     Pre-B_cell_CD34- 
                 136                  102                   71 
    Epithelial_cells            Platelets          Gametocytes 
                  68                   49                   14 
                 CMP            Myelocyte          Hepatocytes 
                  12                   10                    8 
                 MEP              Neurons          Fibroblasts 
                   6                    6                    4 
                  BM     Pro-B_cell_CD34+                  GMP 
                   3                    3                    1 
          HSC_-G-CSF            HSC_CD34+ Neuroepithelial_cell 
                   1                    1                    1 
       Pro-Myelocyte    Tissue_stem_cells 
                   1                    1 

Compara con cuántas filas se “iluminan” visualmente en el mapa de calor anterior.

GRÁFICO / SALIDA: El plotScoreHeatmap completo del Paso 5.2 (Sección 4.8.2) (37 filas); cuenta cuántas muestran un bloque brillante frente a atenuación uniforme.

4.8.3 Paso 5.2b - Consolidar a las etiquetas top-N y verificar la plausibilidad de la distribución

HumanPrimaryCellAtlasData cubre docenas de tipos celulares a través de muchos tejidos. En un panel PBMC de 500 genes, típicamente devuelve más de 20 etiquetas distintas, la mayoría respaldadas por solo un puñado de células (ruido de baja confianza, no poblaciones reales). Graficar todas ellas produce una leyenda ilegible y frustra el propósito de cada comparación posterior.

singler_label (detalle completo) se mantiene intacto en los metadatos para los casos en que la etiqueta exacta importa. singler_label_top colapsa cada etiqueta fuera de las N más frecuentes en “Other (n small labels)” y es lo que usa por defecto cada gráfico de aquí en adelante. El mismo paso que construye la consolidación también verifica si la distribución subyacente parece una anotación PBMC real, ya que ambas preguntas provienen de la misma tabla.

Lo que decide este paso, y lo que no: esta es una decisión basada en conteos sobre qué CATEGORÍAS merecen su propio color en un gráfico. Nunca observa un solo gen o la expresión de una sola célula. Una célula etiquetada “T_cell” permanece “T_cell” aquí sin importar si realmente expresa algún marcador de célula T; lo único que cambia es si esa etiqueta obtiene su propio espacio en la leyenda o se pliega en “Other” porque muy pocas células la comparten.

El Paso 5.4a (Sección 4.8.6), más adelante, plantea una pregunta diferente sobre esta: dada una etiqueta que sobrevivió a este filtro, ¿la célula individual que la lleva realmente muestra la evidencia de expresión que esa etiqueta implica? Esa es una verificación por célula, basada en marcadores, no por etiqueta, basada en conteos.

El Paso 5.4c (Sección 4.8.7) luego construye singler_label_clean directamente a partir de singler_label_top, agregando solo las células que fallaron la verificación de marcadores del Paso 5.4a (Sección 4.8.6) como una nueva categoría “Ambiguous”. Los dos pasos no son redundantes: este decide qué mostrar, Paso 5.4a/5.4c (Sección 4.8.6 / Sección 4.8.7) decide en quién confiar dentro de lo que se muestra.

TOP_N_LABELS <- 8

label_tab <- sort(table(sobj_filt$singler_label), decreasing = TRUE)
top_labels <- names(label_tab)[seq_len(min(TOP_N_LABELS, length(label_tab)))]
n_collapsed <- length(label_tab) - length(top_labels)

sobj_filt$singler_label_top <- ifelse(
  sobj_filt$singler_label %in% top_labels,
  sobj_filt$singler_label,
  sprintf("Other (%d small labels)", n_collapsed)
)
sobj_filt$singler_label_top <- factor(
  sobj_filt$singler_label_top,
  levels = c(top_labels, sprintf("Other (%d small labels)", n_collapsed))
)

cat(sprintf("Kept top %d labels, collapsed %d smaller labels into 'Other':\n",
            length(top_labels), n_collapsed))

Kept top 8 labels, collapsed 18 smaller labels into 'Other':

print(table(sobj_filt$singler_label_top))

                T_cells                Monocyte                 NK_cell 
                    782                     445                     351 
                 B_cell                      DC       Endothelial_cells 
                    305                     154                     138 
            Neutrophils              Macrophage Other (18 small labels) 
                    136                     102                     260 

Reglas prácticas de sanidad para PBMC, verificadas contra la tabla completa (sin colapsar) de etiquetas:

• Se esperan entre 3 y 7 etiquetas dominantes (células T, células B, monocitos, células NK, DC)
• Una etiqueta > 80% de todas las células: sospechoso (colapso de anotación)
• Más de 15 etiquetas distintas en 3,000 células: referencia demasiado granular
• Cualquier etiqueta no inmune individual > 5% (por ejemplo, hepatocitos, fibroblastos): referencia incorrecta

n_labels <- length(label_tab)
top_prop <- max(label_tab) / sum(label_tab)
cat(sprintf("\nDistinct labels (full table): %d  |  Dominant label proportion: %.1f%%\n",
            n_labels, top_prop * 100))

Distinct labels (full table): 26  |  Dominant label proportion: 29.3%

if (n_labels > 15) cat(">> Many distinct labels. Reference may be too granular for PBMC.\n")

>> Many distinct labels. Reference may be too granular for PBMC.

if (top_prop > 0.80) cat(">> One label dominates. Check whether reference matches the tissue.\n")

TipPREGUNTA 5.2b

Observa las etiquetas colapsadas en “Other” (las filas de label_tab más allá del rango 8). ¿Alguna de ellas es biológicamente plausible para PBMC (por ejemplo, “Platelets”, “DC”) o todas son implausibles (por ejemplo, “Hepatocytes”, “Neurons”)? ¿Qué harías diferente si una etiqueta plausible fuera colapsada?

TipRESPUESTA

En este conjunto de datos, HumanPrimaryCellAtlasData típicamente devuelve más de 20 etiquetas para un panel PBMC de 500 genes. La mayor parte de lo que cae fuera del top 8 es biológicamente implausible para sangre (Hepatocitos, Fibroblastos, Células_endoteliales, Neuronas, Gametocitos, Astrocitos) y refleja ruido de baja confianza de una referencia que cubre muchos tejidos, no una señal específica de PBMC. Algunas etiquetas colapsadas SÍ pueden ser biológicamente plausibles pero raras en esta cohorte (Plaquetas, DC, subconjuntos de Pro-B_cell) - esas son poblaciones minoritarias reales, solo que pequeñas. La consolidación top-N es una ayuda de visualización, no una corrección a la anotación en sí: singler_label (sin colapsar) se preserva en los metadatos exactamente por esta razón. Si una población rara plausible importa para tu análisis (por ejemplo, si estás estudiando específicamente plaquetas o DCs), aumenta TOP_N_LABELS, o sigue usando singler_label directamente para ese análisis en particular en lugar de singler_label_top.

DEMOSTRACIÓN EN VIVO:

Inspecciona lo que cayó en “Other”:

label_tab <- sort(table(sobj_filt$singler_label), decreasing = TRUE)
print(label_tab[-(1:8)])

    Pre-B_cell_CD34-     Epithelial_cells            Platelets 
                  71                   68                   49 
         Gametocytes                  CMP            Myelocyte 
                  14                   12                   10 
         Hepatocytes                  MEP              Neurons 
                   8                    6                    6 
         Fibroblasts                   BM     Pro-B_cell_CD34+ 
                   4                    3                    3 
                 GMP           HSC_-G-CSF            HSC_CD34+ 
                   1                    1                    1 
Neuroepithelial_cell        Pro-Myelocyte    Tissue_stem_cells 
                   1                    1                    1 

Aumenta TOP_N_LABELS si una población rara plausible está siendo colapsada:

TOP_N_LABELS <- 12

GRÁFICO / SALIDA: Salida de consola (tabla de etiquetas colapsadas y sus conteos).

4.8.4 Paso 5.3 - UMAP anotado y el mapa de calor depurado

p_ann <- DimPlot(sobj_filt, group.by = "singler_label_top",
                 label = TRUE, label.size = 3, repel = TRUE) +
  ggtitle("SingleR annotation (top labels)") + NoLegend()

p_cond <- DimPlot(sobj_filt, group.by = "condition",
                  cols = c("Healthy" = "#0d7377", "COVID19" = "#991b1b"),
                  pt.size = 0.4) +
  ggtitle("Condition")

p_ann | p_cond

Mismo mapa de calor que el Paso 5.2 (Sección 4.8.2), restringido a las etiquetas que sobrevivieron a la consolidación top-N. Compáralo directamente con la versión completa de 37 filas mostrada anteriormente: esto es lo que se ve “legible” una vez que las etiquetas que no llevaban ninguna señal real fueron dejadas de lado.

cat("Labels passed to labels.use:\n")

Labels passed to labels.use:

print(top_labels)

[1] "T_cells"           "Monocyte"          "NK_cell"          
[4] "B_cell"            "DC"                "Endothelial_cells"
[7] "Neutrophils"       "Macrophage"       

heatmap_result <- tryCatch({
  plotScoreHeatmap(singler_res,
                   labels.use = top_labels,
                   main = "SingleR annotation scores (top labels only; compare to Step 5.2)")
  "filtered"
}, error = function(e) {
  message("labels.use raised an error in this SingleR version: ", conditionMessage(e))
  message("Falling back to the full heatmap.")
  plotScoreHeatmap(singler_res,
                   main = "SingleR annotation scores (wider spread = more confident)")
  "full (fallback)"
})

cat("Heatmap actually drawn:", heatmap_result, "\n")

Heatmap actually drawn: filtered 

cat("If this says 'full (fallback)', the plot above is identical to Step 5.2\n")

If this says 'full (fallback)', the plot above is identical to Step 5.2

cat("by design: that is the fallback path, not the filtered version.\n")

by design: that is the fallback path, not the filtered version.

4.8.5 Paso 5.4 - Verificar el objeto anotado

Compara el objeto ahora con su estado en Sección 4.4 (Bloque 1) antes de pasar al flujo de trabajo de CITE-seq.

cat("Object after RNA workflow + annotation:\n")

Object after RNA workflow + annotation:

print(sobj_filt)

An object of class Seurat 
524 features across 2673 samples within 2 assays 
Active assay: RNA (500 features, 500 variable features)
 3 layers present: counts, data, scale.data
 1 other assay present: ADT
 2 dimensional reductions calculated: pca, umap

cat("\nReductions :", paste(SafeReductions(sobj_filt), collapse = ", "), "\n")

Reductions : pca, umap 

cat("Assays     :", paste(SafeAssays(sobj_filt), collapse = ", "), "\n")

Assays     : RNA, ADT 

cat("New metadata columns vs pre-QC:\n")

New metadata columns vs pre-QC:

print(setdiff(colnames(sobj_filt@meta.data), colnames(sobj@meta.data)))

[1] "scDblFinder.class" "scDblFinder.score" "is_dbl"           
[4] "RNA_snn_res.0.3"   "RNA_snn_res.0.8"   "RNA_snn_res.1.2"  
[7] "singler_label"     "singler_pruned"    "singler_label_top"

TipPREGUNTA 5.4

• ¿Qué slots siguen vacíos / sin cambios en el ensayo ADT?
• ¿Qué te dice eso sobre lo que el Bloque 6 (Sección 5.1) necesita hacer primero?

TipRESPUESTA

ADT aún tiene solo la capa de conteos. Sin capa data (sin normalización), sin scale.data, sin var.features, sin reducciones que referencien a ADT. El Bloque 6 (Sección 5.1) debe construir el lado de ADT desde cero: normalizar (CLR margin=2), y luego o bien usarlo directamente para gráficos y gating (no se necesita scale.data) o ejecutar RunPCA sobre los conteos de ADT antes de WNN.

DEMOSTRACIÓN EN VIVO:

Layers(sobj_filt[["ADT"]])

[1] "counts" "data"  

Reductions(sobj_filt)

[1] "pca"  "umap"

GRÁFICO / SALIDA: Salida de consola

4.8.6 Paso 5.4a - Validación de marcadores canónicos por tipo celular

Antes de confiar en una anotación, verifica que el marcador de RNA canónico de cada linaje principal se detecte en una fracción razonable de las células asignadas a ese linaje. La verificación de sanidad de distribución de etiquetas ya se ejecutó en el Paso 5.2b (Sección 4.8.3), justo después de que las etiquetas fueron consolidadas; este paso es un tipo diferente de verificación, sobre la evidencia por marcador en lugar de los conteos de etiquetas.

El Paso 5.2b (Sección 4.8.3) preguntó “¿cuántas células comparten esta etiqueta, y ese conteo es plausible para PBMC?” Esa es una pregunta sobre etiquetas como categorías: nunca inspeccionó un solo gen. Este paso hace una pregunta diferente: “para una célula que lleva esta etiqueta, ¿su RNA realmente muestra el marcador que esa etiqueta implica?” Esa es una verificación por célula, basada en expresión, ejecutada independientemente de qué tan común o rara fuera la etiqueta. Una etiqueta puede pasar la verificación de frecuencia del Paso 5.2b (Sección 4.8.3) (lo suficientemente común como para mantener su propia categoría) y aun así fallar esta (la mayoría de las células que la llevan no expresan el marcador esperado), que es exactamente lo que la tabla de tasas de detección a continuación está construida para detectar. La salida de este paso (detection_rates) alimenta el Paso 5.4c (Sección 4.8.7), que construye singler_label_clean a partir de singler_label_top y agrega “Ambiguous” solo para las células que fallaron esta verificación de marcadores, sobre las categorías que el Paso 5.2b (Sección 4.8.3) ya decidió que valía la pena mantener.

Mapea los patrones de etiqueta de SingleR a marcadores canónicos

check_markers <- list(
  "T_cell|T cell" = c("CD3E", "CD3D"),
  "B_cell|B cell" = c("CD79A", "MS4A1"),
  "Monocyte"      = c("CD14", "LYZ"),
  "NK"            = c("NKG7", "GNLY"),
  "DC|Dendritic"  = c("FCER1A")
)

cat(sprintf("\n%-25s %-10s %-12s %s\n", "Label", "Marker", "Cells (n)", "Detection rate"))

Label                     Marker     Cells (n)    Detection rate

cat(sprintf("%-25s %-10s %-12s %s\n",   "-----", "------", "---------", "--------------"))

-----                     ------     ---------    --------------

for (lbl in unique(sobj_filt$singler_label)) {
  for (lbl_pat in names(check_markers)) {
    if (grepl(lbl_pat, lbl, ignore.case = TRUE)) {
      cell_mask <- sobj_filt$singler_label == lbl
      n_cells   <- sum(cell_mask)
      for (mk in check_markers[[lbl_pat]]) {
        if (mk %in% rownames(sobj_filt[["RNA"]])) {
          rna_vec <- FetchData(sobj_filt, vars = mk)[cell_mask, 1]
          rate    <- mean(rna_vec > 0) * 100
          flag    <- if (rate < 25) "LOW (dropout or wrong label)" else "ok"
          cat(sprintf("%-25s %-10s %-12d %.1f%% %s\n",
                      substr(lbl, 1, 25), mk, n_cells, rate, flag))
        }
      }
    }
  }
}

DC                        FCER1A     154          8.4% LOW (dropout or wrong label)
B_cell                    CD79A      305          15.7% LOW (dropout or wrong label)
B_cell                    MS4A1      305          17.4% LOW (dropout or wrong label)
Monocyte                  CD14       445          25.8% ok
Monocyte                  LYZ        445          24.7% LOW (dropout or wrong label)
T_cells                   CD3E       782          35.4% ok
T_cells                   CD3D       782          35.3% ok
Pre-B_cell_CD34-          CD79A      71           11.3% LOW (dropout or wrong label)
Pre-B_cell_CD34-          MS4A1      71           11.3% LOW (dropout or wrong label)
NK_cell                   NKG7       351          16.8% LOW (dropout or wrong label)
NK_cell                   GNLY       351          16.0% LOW (dropout or wrong label)
Pro-B_cell_CD34+          CD79A      3            0.0% LOW (dropout or wrong label)
Pro-B_cell_CD34+          MS4A1      3            0.0% LOW (dropout or wrong label)

TipPREGUNTA 5.4a

¿Qué tipo celular anotado tiene la tasa de detección más baja de su marcador canónico? El código de abajo lo encuentra explícitamente.

detection_rates <- data.frame(label = character(0), marker = character(0),
                              rate = numeric(0), stringsAsFactors = FALSE)
for (lbl in unique(sobj_filt$singler_label)) {
  for (lbl_pat in names(check_markers)) {
    if (grepl(lbl_pat, lbl, ignore.case = TRUE)) {
      cell_mask <- sobj_filt$singler_label == lbl
      for (mk in check_markers[[lbl_pat]]) {
        if (mk %in% rownames(sobj_filt[["RNA"]])) {
          rna_vec <- FetchData(sobj_filt, vars = mk)[cell_mask, 1]
          detection_rates <- rbind(detection_rates,
            data.frame(label = lbl, marker = mk,
                       rate = mean(rna_vec > 0) * 100,
                       stringsAsFactors = FALSE))
        }
      }
    }
  }
}
if (nrow(detection_rates) > 0) {
  worst <- detection_rates[which.min(detection_rates$rate), ]
  cat(sprintf("\nLowest detection: %s expressing %s at %.1f%%\n",
              worst$label, worst$marker, worst$rate))
  cat("\nDecision rule:\n")
  cat("  rate >= 50%  -> annotation consistent with RNA (no action)\n")
  cat("  rate 25-50%  -> dropout likely (ADT in Block 6 should rescue)\n")
  cat("  rate < 25%   -> suspect misannotation; cross-check with ADT now\n")
}

Lowest detection: Pro-B_cell_CD34+ expressing CD79A at 0.0%

Decision rule:
  rate >= 50%  -> annotation consistent with RNA (no action)
  rate 25-50%  -> dropout likely (ADT in Block 6 should rescue)
  rate < 25%   -> suspect misannotation; cross-check with ADT now

TipRESPUESTA

Calculado en el script. La tasa más baja normalmente aparece para una etiqueta de célula T expresando CD3E o CD3D; el dropout de RNA de CD3E en PBMC suele ser del 30-50%, así que tasas de detección entre 25-50% indican dropout (recuperable en el Bloque 6 (Sección 5.1) vía ADT CD3). Una tasa de detección por debajo del 25% para cualquier marcador es más probable que sea una mala anotación. Umbrales de decisión: >=50% ok; 25-50% explicado por dropout, verificar con ADT; <25% sospechoso, revisar la anotación.

DEMOSTRACIÓN EN VIVO:

lbl_mask <- sobj_filt$singler_label == "T_cells"
mean(FetchData(sobj_filt, "CD3E")[lbl_mask, 1] > 0) * 100

[1] 35.42199

GRÁFICO / SALIDA: Salida de consola

4.8.7 Paso 5.4c - Marcar etiquetas inconsistentes como Ambiguous

Las células cuya tasa de detección de marcador canónico cayó por debajo del 25% en el Paso 5.4a (Sección 4.8.6) se marcan, no se eliminan. singler_label_clean lleva los mismos valores que singler_label_top excepto esas células marcadas, que se convierten en “Ambiguous”. Cada columna previa (singler_label, singler_label_top) permanece intacta.

low_conf_labels <- character(0)
if (nrow(detection_rates) > 0) {
  per_label_rate  <- aggregate(rate ~ label, data = detection_rates, FUN = min)
  low_conf_labels <- per_label_rate$label[per_label_rate$rate < 25]
}

sobj_filt$singler_label_clean <- as.character(sobj_filt$singler_label_top)
ambiguous_mask <- sobj_filt$singler_label %in% low_conf_labels
sobj_filt$singler_label_clean[ambiguous_mask] <- "Ambiguous"
sobj_filt$singler_label_clean <- factor(sobj_filt$singler_label_clean)

cat(sprintf("\nLabels flagged Ambiguous (canonical marker detection < 25%%): %d\n",
            length(low_conf_labels)))

Labels flagged Ambiguous (canonical marker detection < 25%): 6

if (length(low_conf_labels) > 0) cat(" ", paste(low_conf_labels, collapse = ", "), "\n")

  B_cell, DC, Monocyte, NK_cell, Pre-B_cell_CD34-, Pro-B_cell_CD34+ 

cat(sprintf("Cells marked Ambiguous: %d (%.1f%% of total)\n",
            sum(ambiguous_mask), 100 * mean(ambiguous_mask)))

Cells marked Ambiguous: 1329 (49.7% of total)

cat("\nFinal label distribution (singler_label_clean):\n")

Final label distribution (singler_label_clean):

print(table(sobj_filt$singler_label_clean))

              Ambiguous       Endothelial_cells              Macrophage 
                   1329                     138                     102 
            Neutrophils Other (18 small labels)                 T_cells 
                    136                     186                     782 

TipPREGUNTA 5.4c

Las células marcadas como Ambiguous siguen en sobj_filt y siguen contando para ncol(sobj_filt). ¿Por qué mantenerlas en lugar de eliminarlas directamente? ¿Qué cambiaría silenciosamente en los conteos de células reportados en cada gráfico de aquí en adelante si se eliminaran?

TipRESPUESTA

Mantenerlas preserva un registro honesto de lo que el pipeline de anotación realmente produjo: una fracción real de células no pudo ser tipificada con confianza solo con RNA, y esa fracción es en sí misma informativa (a menudo se reduce una vez que se agrega ADT en el Bloque 6 (Sección 5.1), que es todo el propósito del flujo de trabajo de CITE-seq). Eliminarlas en esta etapa reduciría silenciosamente ncol(sobj_filt), lo cual cambia cada denominador posterior: resúmenes de control de calidad, proporciones por condición, conteos de células por donante, y cualquier cálculo de porcentaje se calcularían sobre una población más pequeña y no documentada. Un lector de un gráfico posterior no tendría forma de saber que se descartaron células aquí a menos que la eliminación se declarara explícitamente cada vez. Marcar y filtrar solo en el paso de graficado (Paso 5.4d, Sección 4.8.8) mantiene el conteo de células del objeto significativo durante el resto del script, y la etiqueta Ambiguous en sí se convierte en un resultado digno de reportar (por ejemplo, en el Bloque 6 (Sección 5.1) puedes verificar si ADT resuelve algunas de estas células).

DEMOSTRACIÓN EN VIVO:

table(sobj_filt$singler_label_clean == "Ambiguous")

FALSE  TRUE 
 1344  1329 

ncol(sobj_filt)  # unchanged regardless of how many cells are Ambiguous

[1] 2673

GRÁFICO / SALIDA: Salida de consola

4.8.8 Paso 5.4d - Volver a graficar el UMAP solo con etiquetas confiables

Mismo UMAP que el Paso 5.3 (Sección 4.8.4), restringido a las células que NO fueron marcadas como Ambiguous. Las células no se eliminan del objeto; simplemente se excluyen de este gráfico específico mediante un filtrado estilo cells.highlight sobre una copia de los metadatos usada para graficar.

confident_cells <- colnames(sobj_filt)[sobj_filt$singler_label_clean != "Ambiguous"]
cat(sprintf("Plotting %d / %d cells (%.1f%%) with confident annotation.\n",
            length(confident_cells), ncol(sobj_filt),
            100 * length(confident_cells) / ncol(sobj_filt)))

Plotting 1344 / 2673 cells (50.3%) with confident annotation.

DimPlot(sobj_filt, cells = confident_cells,
        group.by = "singler_label_clean",
        label = TRUE, label.size = 3, repel = TRUE) +
  ggtitle("Confident annotation only (Ambiguous cells excluded from view)") +
  NoLegend()

4.8.9 Paso 5.4e - Integrar con Harmony

El Paso 4.4b (Sección 4.7.6) encontró que el efecto de donante estaba entre bajo y límite en este conjunto de datos (R^2 en PC1/PC2), así que la integración no fue necesaria para continuar. La ejecutamos de todas formas aquí para construir el panel de comparación en el Paso 5.4f (Sección 4.8.10): con vs sin Harmony es un hábito que vale la pena ver en datos que ya entiendes, antes de depender de él en datos que no.

library(harmony)

Loading required package: Rcpp

• This is Harmony2 version 2.0.5
• Read the guide: run vignette("quickstart", package="harmony")
• Get help: Visit the website at <https://korsunskylab.github.io/harmony2/> and
report issues on <https://github.com/immunogenomics/harmony/issues>

sobj_filt <- RunHarmony(
  sobj_filt,
  group.by.vars  = "donor_id",
  reduction      = "pca",
  dims.use       = 1:20,
  reduction.save = "harmony",
  verbose        = FALSE
)

sobj_filt <- RunUMAP(
  sobj_filt,
  reduction      = "harmony",
  dims           = 1:20,
  reduction.name = "umap.harmony",
  umap.method    = "uwot",
  metric         = "cosine",
  verbose        = FALSE
)

cat("Reductions now available:", paste(SafeReductions(sobj_filt), collapse = ", "), "\n")

Reductions now available: pca, umap, harmony, umap.harmony 

El Paso 4.4b (Sección 4.7.6) calculó el R^2 de donante sobre el PCA simple antes de que se tomara ninguna decisión. Ahora que Harmony se ha ejecutado, calcula el mismo R^2 sobre el embedding armonizado y compara directamente: esta es la ganancia o pérdida real de integrar, en las mismas unidades usadas para tomar la decisión original, no solo una impresión visual de un UMAP.

pc_coords_pca     <- Embeddings(sobj_filt, reduction = "pca")[, 1:2]
pc_coords_harmony <- Embeddings(sobj_filt, reduction = "harmony")[, 1:2]

r2_comparison <- data.frame(
  dimension = c("PC1", "PC2"),
  r2_before_harmony = sapply(1:2, function(pc)
    summary(lm(pc_coords_pca[, pc] ~ sobj_filt$donor_id))$r.squared),
  r2_after_harmony = sapply(1:2, function(pc)
    summary(lm(pc_coords_harmony[, pc] ~ sobj_filt$donor_id))$r.squared)
)
r2_comparison$change <- r2_comparison$r2_after_harmony - r2_comparison$r2_before_harmony

cat("\nDonor R^2 before vs after Harmony:\n")

Donor R^2 before vs after Harmony:

print(r2_comparison, row.names = FALSE)

 dimension r2_before_harmony r2_after_harmony      change
       PC1       0.064676642      0.029133436 -0.03554321
       PC2       0.007479984      0.003635754 -0.00384423

Interpretación: un cambio negativo grande significa que Harmony eliminó exitosamente la varianza impulsada por donante de esa dimensión. Un cambio cercano a cero confirma la lectura del Paso 4.4b (Sección 4.7.6): había poco efecto de donante que eliminar en primer lugar, así que la integración costó varianza biológica real por esencialmente ninguna ganancia en corrección de lote.

TipPREGUNTA 5.4e

¿El cambio en el R^2 de donante es grande o cercano a cero en este conjunto de datos? ¿Coincide eso con lo que predijo la decisión del Paso 4.4b (Sección 4.7.6)? Si ejecutaras esto en un conjunto de datos con R^2 por encima de 0.30, ¿qué cambio esperarías ver en esta misma tabla?

TipRESPUESTA

En este conjunto de datos, el R^2 de donante ya era bajo a límite antes de Harmony (Paso 4.4b, Sección 4.7.6), así que el cambio después de Harmony también debería ser pequeño: a Harmony le queda poca varianza impulsada por donante que eliminar, que es exactamente lo que predijo la decisión del Paso 4.4b (Sección 4.7.6) de omitir la integración. Este es el caso confirmatorio: si el cambio hubiera resultado grande aquí, significaría que la verificación original de R^2 solo en PC1/PC2 pasó por alto un efecto de donante que vivía en PCs posteriores, y la decisión de no integrar habría sido incorrecta. En un conjunto de datos donde el R^2 previo a Harmony estaba por encima de 0.30 (efecto de donante dominante), el patrón esperado es una caída grande en el R^2 después de Harmony (a menudo cayendo hacia 0.05-0.15) en las dimensiones que se le indicó a Harmony corregir, junto con una disposición de tipo celular visiblemente diferente en el panel de comparación: esa combinación es lo que realmente se ve el éxito de la integración en números, no solo en un UMAP que se ve más bonito.

DEMOSTRACIÓN EN VIVO:

print(r2_comparison, row.names = FALSE)

 dimension r2_before_harmony r2_after_harmony      change
       PC1       0.064676642      0.029133436 -0.03554321
       PC2       0.007479984      0.003635754 -0.00384423

• Un valor de columna ‘change’ negativo grande significa que Harmony eliminó señal de donante.
• Un cambio cercano a cero significa que había poca señal de donante que eliminar.

GRÁFICO / SALIDA: La tabla r2_comparison; combinar con la comparación de UMAP de 4 paneles en el Paso 5.4f (Sección 4.8.10).

4.8.10 Paso 5.4f - Panel de comparación: UMAP antes/después de la anotación, con/sin Harmony

Cuatro paneles responden cuatro preguntas diferentes sobre el mismo objeto:

1. Pre-anotación: ¿se ven razonables los clusters no supervisados? 2. Post-anotación (sin Harmony): ¿tiene sentido la anotación en el embedding que realmente usaste a través de Sección 4.8 (Bloque 5, Sección 4.8)? 3. Post-anotación, Harmony: ¿cambió la integración qué células están cerca de cuáles otras? 4. Color de donante en el UMAP de Harmony: ¿logró Harmony realmente mezclar a los donantes, o no lo necesitaba (consistente con la decisión del Paso 4.4b (Sección 4.7.6))?

p1_preannot <- DimPlot(sobj_filt, reduction = "umap", group.by = "seurat_clusters",
                        label = TRUE, label.size = 3) +
  ggtitle("1. Pre-annotation (unsupervised clusters)") + NoLegend()

p2_postannot <- DimPlot(sobj_filt, reduction = "umap", group.by = "singler_label_clean",
                         label = TRUE, label.size = 2.5, repel = TRUE) +
  ggtitle("2. Post-annotation, no Harmony") + NoLegend()

p3_harmony_annot <- DimPlot(sobj_filt, reduction = "umap.harmony", group.by = "singler_label_clean",
                             label = TRUE, label.size = 2.5, repel = TRUE) +
  ggtitle("3. Post-annotation, with Harmony") + NoLegend()

p4_harmony_donor <- DimPlot(sobj_filt, reduction = "umap.harmony", group.by = "donor_id") +
  ggtitle("4. Harmony UMAP by donor")

(p1_preannot | p2_postannot) / (p3_harmony_annot | p4_harmony_donor)

TipPREGUNTA 5.4f

Compara los paneles 2 y 3. Si la disposición de tipos celulares se ve casi idéntica antes y después de Harmony, ¿qué confirma eso sobre la decisión del Paso 4.4b (Sección 4.7.6) de no integrar? Si se ve diferente, ¿en qué panel confiarías para el resto del análisis, y por qué?

TipRESPUESTA

Si los paneles 2 y 3 se ven casi idénticos, eso confirma la lectura del Paso 4.4b (Sección 4.7.6) sobre los datos: el efecto de donante en PC1/PC2 ya era bajo (rango R^2 < 0.10-0.30), así que Harmony tiene poco trabajo que hacer y converge a esencialmente la misma disposición de tipos celulares. Ese es el resultado esperado y tranquilizador aquí, y el panel 4 (UMAP de Harmony por donante) debería seguir mostrando a los donantes razonablemente mezclados dentro de cada región de tipo celular, igual que antes de la integración. Si los paneles 2 y 3 se vieran claramente diferentes en cambio, eso significaría que el R^2 de PC1/PC2 del Paso 4.4b (Sección 4.7.6) subestimó un efecto de donante que vivía en PCs posteriores, y la decisión correcta sería confiar en el panel integrado con Harmony (3) para el análisis posterior, ya que corrige activamente el lote en lugar de solo medirlo en dos dimensiones. De cualquier manera, el panel de comparación es la verificación que valida (o revierte) una decisión de integración tomada anteriormente a partir de un diagnóstico parcial.

DEMOSTRACIÓN EN VIVO:

Cuantifica la comparación visual en lugar de estimarla a simple vista:

table(sobj_filt$singler_label_clean, useNA = "ifany")

              Ambiguous       Endothelial_cells              Macrophage 
                   1329                     138                     102 
            Neutrophils Other (18 small labels)                 T_cells 
                    136                     186                     782 

Compara la composición de clusters antes/después de Harmony usando el índice de Rand ajustado si quieres un número en lugar de un gráfico:

mclust::adjustedRandIndex(cluster_labels_no_harmony, cluster_labels_harmony)

Error:
! object 'cluster_labels_no_harmony' not found

GRÁFICO / SALIDA: El panel de comparación 2x2 en sí (Paso 5.4f, Sección 4.8.10) es la respuesta; no se necesita un gráfico separado.

4.8.11 Paso 5.4g - Marcadores de expresión diferencial por tipo celular anotado

Explícito sobre lo que realmente usa esta prueba, ya que el objeto ahora lleva tanto un PCA/UMAP simple como un PCA/UMAP corregido con Harmony del Paso 5.4e (Sección 4.8.9):

TipExplicación

• Assay: RNA, capa “data” (conteos normalizados logarítmicamente del Bloque 4, Sección 4.7).
  • Harmony solo toca los embeddings de PCA/UMAP (usados para visualización y clustering); nunca modifica los valores de expresión de RNA que lee FindAllMarkers. Ejecutar esto en el objeto armonizado o en el objeto previo a Harmony da resultados idénticos, porque la prueba de DE aquí compara la expresión directamente entre grupos de células, no su posición en ningún embedding.
• Agrupación: singler_label_clean (Paso 5.4c, Sección 4.8.7), células Ambiguous excluidas. No son un grupo biológico coherente, así que probarlas como su propio “cluster” no significaría nada.
• Filtros: only.pos = TRUE (marcadores, no anti-marcadores), min.pct = 0.25 (expresado en al menos el 25% de un grupo), logfc.threshold = 0.5 (al menos un cambio de 1.4 veces).

DefaultAssay(sobj_filt) <- "RNA"
Idents(sobj_filt) <- "singler_label_clean"

de_markers <- FindAllMarkers(
  subset(sobj_filt, idents = "Ambiguous", invert = TRUE),
  only.pos        = TRUE,
  min.pct         = 0.25,
  logfc.threshold = 0.5,
  verbose         = FALSE
)

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 session

Warning in FindMarkers.default(object = data.use, cells.1 = cells.1, cells.2 =
cells.2, : No features pass logfc.threshold threshold; returning empty
data.frame

TOP_N_MARKERS <- 5
top_markers <- de_markers %>%
  group_by(cluster) %>%
  slice_max(order_by = avg_log2FC, n = TOP_N_MARKERS) %>%
  ungroup()

cat(sprintf("Top %d markers per confident cell type:\n", TOP_N_MARKERS))

Top 5 markers per confident cell type:

print(top_markers[, c("cluster", "gene", "avg_log2FC", "pct.1", "pct.2", "p_val_adj")],
      n = nrow(top_markers))

# A tibble: 20 × 6
   cluster           gene  avg_log2FC pct.1 pct.2 p_val_adj
   <fct>             <chr>      <dbl> <dbl> <dbl>     <dbl>
 1 Endothelial_cells IFIT2      0.706 0.971 0.668  6.44e- 9
 2 Endothelial_cells OAS1       0.656 0.957 0.673  1.50e- 7
 3 Endothelial_cells RSAD2      0.569 0.978 0.663  6.39e- 6
 4 Endothelial_cells MX1        0.525 0.935 0.672  3.75e- 5
 5 Endothelial_cells ISG15      0.507 0.949 0.676  9.17e- 7
 6 Macrophage        RSAD2      0.851 0.961 0.673  2.38e-11
 7 Macrophage        MX1        0.805 0.971 0.677  2.45e- 6
 8 Macrophage        ISG15      0.682 0.961 0.683  5.42e- 7
 9 Macrophage        IFIT3      0.670 0.971 0.688  5.23e- 7
10 Macrophage        IFIT1      0.587 0.961 0.687  1.89e- 4
11 Neutrophils       MX1        0.853 0.963 0.67   6.24e-14
12 Neutrophils       IFIT2      0.851 1     0.666  6.55e-14
13 Neutrophils       OAS1       0.818 0.978 0.671  2.48e-13
14 Neutrophils       IFIT1      0.802 0.978 0.677  4.19e-14
15 Neutrophils       RSAD2      0.741 0.963 0.665  1.19e- 9
16 T_cells           FOXP3      0.754 0.347 0.283  6.50e- 1
17 T_cells           CD3G       0.665 0.345 0.283  1   e+ 0
18 T_cells           CCR7       0.557 0.353 0.281  8.39e- 1
19 T_cells           CD3D       0.546 0.353 0.279  5.88e- 1
20 T_cells           CD4        0.515 0.353 0.29   1   e+ 0

DotPlot(sobj_filt,
        features = unique(top_markers$gene),
        idents   = setdiff(levels(Idents(sobj_filt)), "Ambiguous"),
        group.by = "singler_label_clean") +
  RotatedAxis() +
  ggtitle(sprintf("Top %d markers per cell type (confident labels only)", TOP_N_MARKERS))

TipPREGUNTA 5.4g

Elige un tipo celular del DotPlot. ¿Su marcador principal coincide con un marcador canónico que ya conoces para ese linaje (Paso 1.8, Sección 4.4.9)? Si no, ¿es eso una señal de alerta sobre la anotación, o un hallazgo fuera de la lista canónica que vale la pena investigar más a fondo?

TipRESPUESTA

Para la mayoría de los linajes principales en este conjunto de datos, el marcador de DE principal debería coincidir con un gen canónico de la lista del Paso 1.8 (Sección 4.4.9): CD3D/CD3E para células T, MS4A1/CD79A para células B, CD14/LYZ para monocitos, NKG7/GNLY para células NK. Una coincidencia es tranquilizadora: la evidencia estadística independiente (prueba de DE) concuerda con el conocimiento biológico previo (marcadores canónicos), que es la forma más fuerte de validación disponible sin un ensayo ortogonal. Una discrepancia es más interesante que alarmante por sí sola: primero descarta explicaciones técnicas (¿el gen canónico siquiera está en el panel de 500 genes? Verifica rownames(sobj_filt[['RNA']])), luego considera si el cluster es un subtipo conocido pero menos de libro de texto (por ejemplo, un subconjunto de Treg cuyo marcador principal es FOXP3 en lugar de genes CD3 genéricos) antes de tratarlo como una señal de alerta sobre la anotación en sí.

DEMOSTRACIÓN EN VIVO:

Verifica si un marcador canónico siquiera existe en este panel de 500 genes:

canonical_check <- c("CD3D","CD3E","MS4A1","CD79A","CD14","LYZ","NKG7","GNLY")
canonical_check[!canonical_check %in% rownames(sobj_filt[["RNA"]])]

character(0)

Luego compara contra el marcador principal del DotPlot para el tipo celular en cuestión

top_markers[top_markers$cluster == "T_cells", ]

# A tibble: 5 × 7
    p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene 
    <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
1 0.00130      0.754 0.347 0.283     0.650 T_cells FOXP3
2 0.00452      0.665 0.345 0.283     1     T_cells CD3G 
3 0.00168      0.557 0.353 0.281     0.839 T_cells CCR7 
4 0.00118      0.546 0.353 0.279     0.588 T_cells CD3D 
5 0.00542      0.515 0.353 0.29      1     T_cells CD4  

GRÁFICO / SALIDA: El DotPlot del Paso 5.4g (Sección 4.8.11); comparar visualmente contra marcadores canónicos.

4.8.12 Paso 5.5 - Guardar el checkpoint anotado

saveRDS(top_markers, "checkpoints/outputs/top_markers.rds")
saveRDS(sobj_filt, "checkpoints/outputs/sobj_preprocessed.rds")

cat("Checkpoint saved: outputs/sobj_preprocessed.rds\n")

Checkpoint saved: outputs/sobj_preprocessed.rds

cat("To reload: sobj_filt <- readRDS('outputs/sobj_preprocessed.rds')\n")

To reload: sobj_filt <- readRDS('outputs/sobj_preprocessed.rds')

Mira el UMAP coloreado por condición. ¿Algún cluster aparece exclusiva o predominantemente en donantes con COVID-19? Mira el DotPlot: ¿qué clusters podrían ser monocitos según la expresión de CD14 y FCGR3A? Después del descanso, agregamos la capa de proteína para probar estas interpretaciones.

Tip📌 Nota final para los usuarios

Al final de este paso, puedes descargar el objeto Seurat preprocesado para continuar con la siguiente etapa del análisis de CITE‑seq:

👉 Este archivo contiene el objeto Seurat preprocesado (sobj_preprocessed.rds) y sirve como punto de partida para el análisis multimodal posterior. Al guardar y compartir este checkpoint, todos los usuarios pueden continuar de manera consistente sin repetir los pasos anteriores.

           used  (Mb) gc trigger   (Mb) limit (Mb) max used   (Mb)
Ncells 12097809 646.1   20894572 1115.9         NA 20894572 1115.9
Vcells 25254397 192.7   55391439  422.7      24576 55391439  422.7