Apéndice A — Extra

A.1 Bloque 8 - Proyecto final (opcional, después de clase): un objeto desconocido y averiado

Objetivo: Aplicar los hábitos de diagnóstico de los bloques anteriores a un único objeto que contiene cuatro problemas independientes. Para cada escenario, ejecuta el código de diagnóstico, escribe tu hipótesis como un comentario en el bloque de abajo, y luego ejecuta la corrección para confirmarla.

Este es el proyecto final: el resto del curso te enseñó a detectar fallos específicos de uno en uno. Aquí todos llegan en el mismo objeto, sin un orden particular, y sin etiquetas.

Este bloque va después del Cierre a propósito. La sesión en vivo termina en el Bloque 7; esto es trabajo para hacer en casa para quien quiera más práctica antes de la siguiente sesión, no algo para apresurarse en el salón. Ejecútalo en tu propio tiempo, idealmente uno o dos días después del curso, una vez que los hábitos de diagnóstico de los Bloques 1-6 hayan tenido oportunidad de asentarse.

El instructor publicará sobj_broken.rds por separado. Si el archivo no está presente, este bloque se omite (have_broken = FALSE) y el resto del script sigue ejecutándose.

Intenta primero con el objeto inyectado (necesario para el Escenario 4), y luego con el base. Si ninguno existe, el Bloque 8 se omite de forma controlada (have_broken = FALSE) para que el resto del informe siga renderizándose.

Código
broken_candidates <- c("checkpoints/sobj_broken_errors.rds", "checkpoints/sobj_broken.rds")
broken_path <- broken_candidates[file.exists(broken_candidates)][1]

have_broken <- !is.na(broken_path)

if (have_broken) {
  sobj_broken <- readRDS(broken_path)
  cat("Loaded broken object from:", broken_path, "\n")
  print(sobj_broken)
  head(sobj_broken@meta.data, 3)
} else {
  cat("Broken object not found in data/. Block 8 will be skipped.\n")
  cat("To enable Block 8, place one of these files (paths relative to this .R script):\n")
  cat("  ", paste(broken_candidates, collapse = "  or  "), "\n")
  cat("Scenario 4 (ADT batch / isotype) requires the *_errors.rds from\n")
  cat("inject_course_errors.R.\n")
}
Loaded broken object from: checkpoints/sobj_broken_errors.rds 
An object of class Seurat 
526 features across 3000 samples within 2 assays 
Active assay: ADT (26 features, 0 variable features)
 1 layer present: counts
 1 other assay present: RNA
 2 dimensional reductions calculated: pca, umap
           orig.ident nCount_RNA nFeature_RNA percent.mt donor_id condition
CELL000001    Donor01        336           36   18.15476  Donor01   Healthy
CELL000011    Donor01        344           40   15.11628  DONOR01   Healthy
CELL000021    Donor01        426           39    0.00000  DONOR01   Healthy
           severity age cell_type nCount_ADT nFeature_ADT percent.ribo
CELL000001           35    NKcell        267           25     46.13095
CELL000011           35        DC        300           25     64.53488
CELL000021           35     Bcell        472           26     66.66667
           RNA_snn_res.0.5 seurat_clusters adt_batch donor_num
CELL000001               8               8    batchA         1
CELL000011               9               9    batchB         1
CELL000021              10              10    batchB         1

A.1.1 Escenario 1 - Assay por defecto incorrecto

Contexto: Recibiste este objeto de un colaborador. Ejecutas el flujo de trabajo estándar de preprocesamiento de RNA. Todo se ejecuta sin errores, pero los resultados son completamente incorrectos: FindVariableFeatures devuelve menos características de las esperadas, el PCA explica casi toda la varianza en PC1, y el UMAP es una única mancha (blob).

Código
if (have_broken) {
  cat("Default assay        :", DefaultAssay(sobj_broken), "\n")
  cat("Features active assay:", nrow(sobj_broken), "\n")
  cat("All assays           :", paste(SafeAssays(sobj_broken), collapse = ", "), "\n")
}
Default assay        : ADT 
Features active assay: 26 
All assays           : RNA, ADT
Código
if (have_broken) {
  # What does FindVariableFeatures return on this object?
  test_hvg <- FindVariableFeatures(sobj_broken, nfeatures = 2000, verbose = FALSE)
  cat("Variable features found:", length(VariableFeatures(test_hvg)), "\n")
  cat("Feature names:\n")
  print(head(VariableFeatures(test_hvg), 15))
}
Variable features found: 26 
Feature names:
 [1] "CD27"          "CD4"           "CD86"          "CD8a"         
 [5] "CD56"          "CD38"          "CD20"          "CD25"         
 [9] "PD1"           "CD62L"         "LAG3"          "IgG1-isotype" 
[13] "IgG2a-isotype" "TIGIT"         "CD69"
Código
if (have_broken) {
  # Compare RNA and ADT feature counts
  cat("Features in RNA assay:", nrow(sobj_broken[["RNA"]]), "\n")
  cat("Features in ADT assay:", nrow(sobj_broken[["ADT"]]), "\n")
}
Features in RNA assay: 500 
Features in ADT assay: 26

📝 Escribe tu hipótesis como un comentario abajo, luego ejecuta la corrección para confirmarla. ¿Cuál es el DefaultAssay? ¿Por qué FindVariableFeatures devuelve tan pocas características? ¿Qué pasaría si ejecutaras RunPCA y RunUMAP en este objeto sin corregirlo? ¿Cómo detectarías este problema en un objeto Seurat publicado que descargaste?

Código
if (have_broken) {
  cat("Before fix:", DefaultAssay(sobj_broken), "\n")

  # The DefaultAssay is set to "ADT" (24 proteins).
  # All RNA functions are running silently on 24 protein features instead
  # of 500 RNA genes. FindVariableFeatures returns 24 features because that
  # is all that exists in the active assay.

  DefaultAssay(sobj_broken) <- "RNA"
  cat("After fix :", DefaultAssay(sobj_broken), "\n")
  cat("Features  :", nrow(sobj_broken), "\n")

  test_fixed <- FindVariableFeatures(sobj_broken, nfeatures = 2000, verbose = FALSE)
  cat("Variable features now:", length(VariableFeatures(test_fixed)), "\n")
}
Before fix: ADT 
After fix : RNA 
Features  : 500 
Variable features now: 500

Lecciones aprendidas:

  1. DefaultAssay() es la primera línea que se debe ejecutar en cualquier objeto recibido.
  2. Este error es completamente silencioso. Sin advertencia. Sin error. Solo resultados incorrectos.
  3. FindVariableFeatures, ScaleData, RunPCA, FindMarkers operan todos sobre el assay activo. Un PCA ejecutado sobre 24 características de ADT no es un PCA transcriptómico.
  4. Siempre verifica DefaultAssay() después de cualquier cambio de assay para confirmar que se restableció.

A.1.2 Escenario 2 - Inconsistencias en los metadatos

Contexto: Este objeto se ensambló a partir de muestras procesadas en múltiples sitios. Agrupar por condition y donor_id produce resultados inesperados. Faltan algunas muestras en los gráficos y ciertos donantes aparecen duplicados.

Código
if (have_broken) {
  cat("Metadata columns:\n")
  print(colnames(sobj_broken@meta.data))

  cat("\nNAs per column:\n")
  print(colSums(is.na(sobj_broken@meta.data)))
}
Metadata columns:
 [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" "adt_batch"       "donor_num"      

NAs per column:
     orig.ident      nCount_RNA    nFeature_RNA      percent.mt        donor_id 
              0               0               0               0               0 
      condition        severity             age       cell_type      nCount_ADT 
             30               0               0             200               0 
   nFeature_ADT    percent.ribo RNA_snn_res.0.5 seurat_clusters       adt_batch 
              0               0               0               0               0 
      donor_num 
              0
Código
if (have_broken) {
  # Inspect the problematic columns
  cat("Unique donor_id values:\n")
  print(sort(unique(sobj_broken$donor_id)))

  cat("\nCondition table (including NAs):\n")
  print(table(sobj_broken$condition, useNA = "always"))

  cat("\nCell type NAs:", sum(is.na(sobj_broken$cell_type)), "\n")
}
Unique donor_id values:
 [1] "Donor01" "DONOR01" "Donor02" "DONOR02" "Donor03" "DONOR03" "Donor04"
 [8] "DONOR04" "Donor05" "DONOR05" "Donor06" "DONOR06" "Donor07" "DONOR07"
[15] "donor1"  "donor2"  "donor3"  "donor4"  "donor5"  "donor6"  "donor7" 

Condition table (including NAs):

Covid-19  covid19  COVID19       HC  Healthy     <NA> 
      60       79     1840       50      941       30 

Cell type NAs: 200
Código
if (have_broken) {
  # What does a DimPlot by condition look like?
  DimPlot(sobj_broken, group.by = "condition") +
    ggtitle("Condition plot: how many groups appear?")
}

📝 Escribe tu hipótesis como un comentario abajo, luego ejecuta la corrección para confirmarla. ¿Cuántos problemas distintos encontraste? ¿Cuál sería la consecuencia biológica si cada uno quedara sin corregir? Por ejemplo: si las etiquetas de condición son inconsistentes, ¿qué le sucede a una comparación entre Sano y COVID-19?

Código
if (have_broken) {
  sobj_meta <- sobj_broken

  # PROBLEM 1: donor_id has three different formats for the same donors
  # e.g. "Donor01", "donor1", "DONOR01" all mean the same donor
  cat("=== FIX 1: donor_id ===\n")
  sobj_meta$donor_id <- toupper(gsub("[^A-Z0-9a-z]", "", sobj_meta$donor_id))
  print(table(sobj_meta$donor_id))

  # PROBLEM 2: condition has 4 variants of 2 values
  # "COVID19", "covid19", "Covid-19", "HC", NA
  cat("\n=== FIX 2: condition ===\n")
  sobj_meta$condition_clean <- dplyr::case_when(
    grepl("covid", sobj_meta$condition, ignore.case = TRUE) ~ "COVID-19",
    grepl("healthy|HC|control", sobj_meta$condition,
          ignore.case = TRUE)                               ~ "Healthy",
    TRUE                                                    ~ NA_character_
  )
  print(table(sobj_meta$condition_clean, useNA = "always"))

  # PROBLEM 3: 200 cells have NA in cell_type
  cat("\n=== FIX 3: cell_type NAs ===\n")
  cat("NAs before:", sum(is.na(sobj_meta$cell_type)), "\n")
  sobj_meta$cell_type[is.na(sobj_meta$cell_type)] <- "Unknown"
  cat("NAs after :", sum(is.na(sobj_meta$cell_type)), "\n")
}
=== FIX 1: donor_id ===

DONOR01 DONOR02 DONOR03 DONOR04 DONOR05 DONOR06 DONOR07  DONOR1  DONOR2  DONOR3 
    247     218     244     306     304     392     389     103      82     106 
 DONOR4  DONOR5  DONOR6  DONOR7 
    144     146     158     161 

=== FIX 2: condition ===

COVID-19  Healthy     <NA> 
    1979      991       30 

=== FIX 3: cell_type NAs ===
NAs before: 200 
NAs after : 0

Lecciones aprendidas:

  1. table(col, useNA = "always") y str(meta.data) antes de cada análisis.
  2. Los datos de múltiples sitios casi siempre tienen inconsistencias de formato.
  3. NA en condición: esa célula queda excluida de cualquier comparación entre Sano y COVID-19. Con 30 NAs, estás descartando silenciosamente el 1% de las células de todas las comparaciones de grupo.
  4. Las mayúsculas y minúsculas importan: "COVID19" == "covid19" es FALSE en R.

A.1.3 Escenario 3 - Fallo de normalización y efecto de lote

Contexto: El UMAP muestra una separación clara por donante en lugar de por tipo celular. Los datos de ADT también parecen distorsionados. Identifica ambos problemas y corrígelos.

Código
if (have_broken) {
  DimPlot(sobj_broken, group.by = "orig.ident") +
    ggtitle("Donors separate in UMAP: batch effect or biology?")
}

Código
if (have_broken) {
  DefaultAssay(sobj_broken) <- "ADT"
  adt_data <- LayerData(sobj_broken, layer = "data")

  cat("Row means (per protein); near 0 means margin=1 was used (wrong):\n")
  print(round(rowMeans(adt_data), 4))

  cat("\nColumn means (per cell, first 10); near 0 means margin=2 (correct):\n")
  print(round(colMeans(adt_data)[1:10], 4))

  DefaultAssay(sobj_broken) <- "RNA"
}
Warning: Layer 'data' is empty
Row means (per protein); near 0 means margin=1 was used (wrong):
numeric(0)

Column means (per cell, first 10); near 0 means margin=2 (correct):
 [1] NA NA NA NA NA NA NA NA NA NA
Código
if (have_broken) {
  rna_data <- LayerData(sobj_broken, assay = "RNA", layer = "data")
  cat("Any negative values in RNA data layer?",
      any(rna_data < 0), "\n")
  cat("(TRUE would indicate SCTransform residuals stored as 'data')\n")
  cat("\nRange of non-zero values in RNA data:\n")
  print(summary(rna_data@x))
}
Any negative values in RNA data layer? FALSE 
(TRUE would indicate SCTransform residuals stored as 'data')

Range of non-zero values in RNA data:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  2.553   3.392   4.530   4.620   5.715   8.034

📝 Escribe tu hipótesis como un comentario abajo, luego ejecuta la corrección para confirmarla. ¿Cómo determinaste qué margin se usó? ¿Cuál es la consecuencia del efecto de lote para cualquier comparación a nivel de condición? Si publicaras el UMAP original como una figura, ¿qué afirmación científica quedaría invalidada?

Código
if (have_broken) {
  sobj_renorm <- sobj_broken

  # STEP 1: Re-normalize RNA from counts (counts layer is untouched)
  DefaultAssay(sobj_renorm) <- "RNA"
  sobj_renorm <- NormalizeData(sobj_renorm,
    normalization.method = "LogNormalize",
    scale.factor         = 10000,
    verbose              = FALSE)
  cat("RNA re-normalized from counts.\n")

  # STEP 2: Re-normalize ADT with correct margin
  DefaultAssay(sobj_renorm) <- "ADT"
  sobj_renorm <- NormalizeData(sobj_renorm,
    normalization.method = "CLR",
    margin               = 2,
    verbose              = FALSE)
  cat("ADT re-normalized (CLR, margin=2).\n")

  # STEP 3: Rerun preprocessing
  DefaultAssay(sobj_renorm) <- "RNA"
  sobj_renorm <- FindVariableFeatures(sobj_renorm, verbose = FALSE)
  sobj_renorm <- ScaleData(sobj_renorm, verbose = FALSE)
  sobj_renorm <- RunPCA(sobj_renorm, npcs = 30, verbose = FALSE)

  # STEP 4: Harmony batch correction
  sobj_renorm <- RunHarmony(
    sobj_renorm,
    group.by.vars  = "orig.ident",
    reduction      = "pca",
    reduction.save = "harmony",
    verbose        = FALSE
  )

  sobj_renorm <- RunUMAP(sobj_renorm,
    reduction      = "harmony",
    dims           = 1:20,
    reduction.name = "umap.harmony",
    verbose        = FALSE)

  cat("Harmony applied.\n")
}
RNA re-normalized from counts.
ADT re-normalized (CLR, margin=2).
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
Harmony applied.
Código
if (have_broken) {
  p_before <- DimPlot(sobj_renorm, reduction = "umap",
                      group.by = "orig.ident") +
    ggtitle("Before Harmony") + theme(legend.position = "none")

  p_after  <- DimPlot(sobj_renorm, reduction = "umap.harmony",
                      group.by = "orig.ident") +
    ggtitle("After Harmony") + theme(legend.position = "none")

  p_before | p_after
}

Lecciones aprendidas:

  1. Medias de fila ~0 en los datos de ADT es la firma diagnóstica de margin=1. Incorrecto.
  2. Medias de columna ~0 en los datos de ADT es la firma de margin=2. Correcto.
  3. La capa de conteos (counts) es la verdad fundamental (ground truth). Vuelve a normalizar a partir de ella; nunca la sobrescribas.
  4. La separación por donante en el UMAP = efecto de lote. Ejecuta DimPlot(group.by="orig.ident") antes de interpretar cualquier UMAP biológicamente.
  5. Harmony corrige los efectos de lote a nivel de donante mientras preserva la variación biológica. Requiere al menos 3 donantes por grupo para ser confiable.

A.1.4 Escenario 4 | Lote de tinción de ADT y fondo de isotipo (CITE-seq, 25 min)

Contexto: Los datos de proteína parecen técnicamente correctos: se normalizan, se grafican, el gating funciona. Pero el PCA de ADT separa las células por lote de procesamiento, no por tipo celular, y un donante muestra una señal inflada en todas las proteínas. Hay dos problemas distintos de CITE-seq presentes: un efecto de lote de tinción y un alto fondo inespecífico. Ninguno es un problema de RNA y ninguno produce un error.

Código
if (have_broken) {
  if (!"adt_batch" %in% colnames(sobj_broken@meta.data)) {
    cat("This broken object has no 'adt_batch' column.\n")
    cat("Scenario 4 needs the injected object (data/sobj_broken_errors.rds\n")
    cat("from inject_course_errors.R). Skipping Scenario 4 diagnostics.\n")
  } else {
    DefaultAssay(sobj_broken) <- "ADT"
    adt_feats <- rownames(sobj_broken[["ADT"]])
    sobj_broken <- ScaleData(sobj_broken, features = adt_feats, verbose = FALSE)
    sobj_broken <- RunPCA(sobj_broken, features = adt_feats,
                          npcs = 15, reduction.name = "pca.adt.broken",
                          reduction.key = "pcaADTb_", verbose = FALSE)
    print(
      DimPlot(sobj_broken, reduction = "pca.adt.broken", group.by = "adt_batch") +
        ggtitle("ADT PCA colored by staining batch: should NOT separate")
    )
    DefaultAssay(sobj_broken) <- "RNA"
  }
}
Warning: No layers found matching search pattern provided
Error in `ScaleData()`:
! No layer matching pattern 'data' found. Please run NormalizeData and retry
Código
if (have_broken) {
  # Isotype antibodies bind nothing specific. High isotype = high background
  # (sticky/dying cells, over-staining). Per-cell total ADT correlated with
  # isotype signal is the fingerprint.
  DefaultAssay(sobj_broken) <- "ADT"
  adt_counts <- LayerData(sobj_broken, layer = "counts")

  isotypes <- grep("[Ii]sotype|IgG", rownames(adt_counts), value = TRUE)
  cat("Isotype control channels found:", paste(isotypes, collapse = ", "), "\n")

  if (length(isotypes) > 0) {
    iso_total  <- Matrix::colSums(adt_counts[isotypes, , drop = FALSE])
    adt_total  <- Matrix::colSums(adt_counts)
    cat("Spearman(total ADT, isotype signal):",
        round(cor(adt_total, iso_total, method = "spearman"), 3), "\n")
    # donor_id is inconsistent in this object; use a clean donor number
    dn <- sobj_broken$donor_num
    if (is.null(dn))
      dn <- suppressWarnings(as.integer(gsub("\\D", "",
                             as.character(sobj_broken$donor_id))))
    cat("Median isotype signal by donor number:\n")
    print(round(tapply(iso_total, dn, median), 2))
  }
  DefaultAssay(sobj_broken) <- "RNA"
}
Isotype control channels found: IgG1-isotype, IgG2a-isotype 
Spearman(total ADT, isotype signal): 0.215 
Median isotype signal by donor number:
   1    2    3    4    5    6    7 
16.5  0.0  0.0  0.0  0.0  0.0  0.0

📝 Escribe tu hipótesis como un comentario abajo, luego ejecuta la corrección para confirmarla. ¿Qué donante tiene el fondo más alto? ¿Eliminaría CLR por sí solo un efecto de lote de tinción? Si aplicaras el gating de células CD4+ con un único umbral fijo a través de ambos lotes, ¿qué le sucedería a los conteos por lote?

Código
if (have_broken) {
  if (!"adt_batch" %in% colnames(sobj_broken@meta.data)) {
    cat("No 'adt_batch' column; Scenario 4 solution requires the injected object.\n")
  } else {
  # STEP 1: Flag and optionally remove high-background cells using isotypes.
  DefaultAssay(sobj_broken) <- "ADT"
  adt_counts <- LayerData(sobj_broken, layer = "counts")
  isotypes   <- grep("[Ii]sotype|IgG", rownames(adt_counts), value = TRUE)

  sobj_clean <- sobj_broken
  if (length(isotypes) > 0) {
    iso_total <- Matrix::colSums(adt_counts[isotypes, , drop = FALSE])
    cut_hi    <- quantile(iso_total, 0.95)
    keep      <- iso_total <= cut_hi
    cat("High-background cells removed (top 5% isotype):", sum(!keep), "\n")
    sobj_clean <- subset(sobj_clean, cells = colnames(sobj_clean)[keep])
  }

  # STEP 2: Re-normalize ADT (CLR margin=2) from counts.
  sobj_clean <- NormalizeData(sobj_clean, normalization.method = "CLR",
                              margin = 2, verbose = FALSE)

  # STEP 3: Correct the ADT staining batch with Harmony in protein space.
  adt_feats  <- rownames(sobj_clean[["ADT"]])
  sobj_clean <- ScaleData(sobj_clean, features = adt_feats, verbose = FALSE)
  sobj_clean <- RunPCA(sobj_clean, features = adt_feats, npcs = 15,
                       reduction.name = "pca.adt", reduction.key = "pcaADT_",
                       verbose = FALSE)
  sobj_clean <- RunHarmony(sobj_clean, group.by.vars = "adt_batch",
                           reduction = "pca.adt",
                           reduction.save = "harmony.adt", verbose = FALSE)
  cat("ADT batch corrected in protein space (harmony.adt).\n")
  DefaultAssay(sobj_clean) <- "RNA"
  }
}
High-background cells removed (top 5% isotype): 138 
Warning in svd.function(A = t(x = object), nv = npcs, ...): You're computing
too large a percentage of total singular values, use a standard svd instead.
ADT batch corrected in protein space (harmony.adt).

Lecciones aprendidas:

  1. Una capa de ADT que se ve limpia aún puede estar dominada por estructura técnica. Siempre ejecuta un PCA de ADT y coloréalo por lote antes de confiar en el espacio de proteínas.
  2. Los controles de isotipo son la verdad fundamental (ground truth) para el fondo. Si el ADT total sigue la señal de isotipo, las células con valores altos son fondo, no biología.
  3. CLR por célula no elimina un efecto de tinción entre lotes. Se necesita una corrección de lote (Harmony sobre el PCA de ADT), separada del lote de RNA.
  4. Un único umbral de gating a través de los lotes clasifica erróneamente las células. Aplica el gating por lote o corrige el lote primero.