Apêndice A — Extra

A.1 Bloco 8 - Projeto final (opcional, pós-aula): um objeto desconhecido e quebrado

Objetivo: Aplicar os hábitos de diagnóstico dos blocos anteriores a um único objeto que contém quatro problemas independentes. Para cada cenário, execute o código de diagnóstico, escreva sua hipótese como um comentário no bloco abaixo, e então execute a correção para confirmá-la.

Este é o projeto final: o restante do curso ensinou você a identificar falhas específicas, uma de cada vez. Aqui todas chegam no mesmo objeto, sem uma ordem particular, e sem rótulos.

Este bloco vem depois do Encerramento de propósito. A sessão ao vivo termina no Bloco 7; isto é trabalho para casa para quem quiser mais prática antes da próxima sessão, não algo para ser feito às pressas na sala. Execute-o no seu próprio tempo, idealmente um ou dois dias depois do curso, uma vez que os hábitos de diagnóstico dos Blocos 1-6 tenham tido a chance de se consolidar.

O instrutor disponibilizará o sobj_broken.rds separadamente. Se o arquivo não estiver presente, este bloco é pulado (have_broken = FALSE) e o resto do script continua rodando normalmente.

Tenta primeiro o objeto injetado (necessário para o Cenário 4), depois o base. Se nenhum dos dois existir, o Bloco 8 é pulado de forma controlada (have_broken = FALSE) para que o restante do relatório ainda seja renderizado.

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 Cenário 1 - Assay padrão incorreto

Contexto: Você recebeu este objeto de um colaborador. Você executa o fluxo de trabalho padrão de pré-processamento de RNA. Tudo é executado sem erros, mas os resultados estão completamente errados: FindVariableFeatures retorna menos features do que o esperado, o PCA explica quase toda a variância em PC1, e o UMAP é uma ú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

📝 Escreva sua hipótese como um comentário abaixo, depois execute a correção para confirmar. Qual é o DefaultAssay? Por que o FindVariableFeatures retorna tão poucas features? O que aconteceria se você executasse RunPCA e RunUMAP neste objeto sem corrigi-lo? Como você detectaria esse problema em um objeto Seurat publicado que você baixou?

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

Lições aprendidas:

  1. DefaultAssay() é a primeira linha a ser executada em qualquer objeto recebido.
  2. Este erro é completamente silencioso. Sem aviso. Sem erro. Apenas resultados errados.
  3. FindVariableFeatures, ScaleData, RunPCA, FindMarkers operam todos sobre o assay ativo. Um PCA executado em 24 features de ADT não é um PCA transcriptômico.
  4. Sempre verifique DefaultAssay() após qualquer troca de assay para confirmar que ele foi redefinido.

A.1.2 Cenário 2 - Inconsistências nos metadados

Contexto: Este objeto foi montado a partir de amostras processadas em múltiplos locais. Agrupar por condition e donor_id produz resultados inesperados. Algumas amostras estão faltando nos gráficos e certos doadores aparecem 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?")
}

📝 Escreva sua hipótese como um comentário abaixo, depois execute a correção para confirmar. Quantos problemas distintos você encontrou? Qual seria a consequência biológica se cada um deles não fosse corrigido? Por exemplo: se os rótulos de condição são inconsistentes, o que acontece com uma comparação entre Saudável e 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

Lições aprendidas:

  1. table(col, useNA = "always") e str(meta.data) antes de cada análise.
  2. Dados de múltiplos locais quase sempre têm inconsistências de formatação.
  3. NA em condition: essa célula é excluída de qualquer comparação entre Saudável e COVID-19. Com 30 NAs, você está descartando silenciosamente 1% das células de todas as comparações de grupo.
  4. Diferenciação entre maiúsculas e minúsculas importa: "COVID19" == "covid19" é FALSE em R.

A.1.3 Cenário 3 - Falha de normalização e efeito de lote

Contexto: O UMAP mostra uma separação clara por doador em vez de por tipo celular. Os dados de ADT também parecem distorcidos. Identifique ambos os problemas e corrija-os.

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

📝 Escreva sua hipótese como um comentário abaixo, depois execute a correção para confirmar. Como você determinou qual margin foi usado? Qual é a consequência do efeito de lote para qualquer comparação em nível de condição? Se você publicasse o UMAP original como uma figura, qual afirmação científica seria 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
}

Lições aprendidas:

  1. Médias de linha ~0 nos dados de ADT é a assinatura diagnóstica do margin=1. Errado.
  2. Médias de coluna ~0 nos dados de ADT é a assinatura do margin=2. Correto.
  3. A camada de contagens (counts) é a verdade fundamental (ground truth). Renormalize a partir dela; nunca a sobrescreva.
  4. Separação por doador no UMAP = efeito de lote. Execute DimPlot(group.by="orig.ident") antes de interpretar qualquer UMAP biologicamente.
  5. O Harmony corrige efeitos de lote em nível de doador enquanto preserva a variação biológica. Requer pelo menos 3 doadores por grupo para ser confiável.

A.1.4 Cenário 4 | Lote de coloração de ADT e fundo de isotipo (CITE-seq, 25 min)

Contexto: Os dados de proteína parecem tecnicamente corretos: normalizam, plotam, o gating funciona. Mas o PCA de ADT separa as células por lote de processamento, não por tipo celular, e um doador mostra sinal inflado em todas as proteínas. Dois problemas distintos de CITE-seq estão presentes: um efeito de lote de coloração e um alto fundo inespecífico. Nenhum dos dois é um problema de RNA e nenhum gera um erro.

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

📝 Escreva sua hipótese como um comentário abaixo, depois execute a correção para confirmar. Qual doador tem o maior fundo? O CLR sozinho removeria um efeito de lote de coloração? Se você fizesse o gating de células CD4+ com um único limiar fixo através de ambos os lotes, o que aconteceria com as contagens 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).

Lições aprendidas:

  1. Uma camada de ADT que parece limpa ainda pode estar dominada por estrutura técnica. Sempre execute um PCA de ADT e colora-o por lote antes de confiar no espaço de proteínas.
  2. Os controles de isotipo são a verdade fundamental (ground truth) para o fundo. Se o ADT total acompanha o sinal de isotipo, as células com valores altos são fundo, não biologia.
  3. O CLR por célula não remove um efeito de coloração entre lotes. É necessária uma correção de lote (Harmony no PCA de ADT), separada do lote de RNA.
  4. Um único limiar de gating através dos lotes classifica erroneamente as células. Faça o gating por lote ou corrija o lote primeiro.