# SummarizedExperiment subsetting: keep assays, rowData, and colData aligned

When you filter samples or features in a Bioconductor `SummarizedExperiment`, subset the
container with `[i, j]` -- never the assay matrix directly. The `[` method subsets every assay,
`rowData`/`rowRanges`, and `colData` in one operation, so the pieces cannot drift out of sync.
The official vignette's motivation: "Improperly accounting for metadata and observational data
has resulted in a number of incorrect results and retractions."

## The one rule

Rows are features (genes, transcripts, exons). Columns are samples. Always do
`se_subset <- se[i, j]` and only then extract data with `assay(se_subset)`, `colData(se_subset)`,
`rowData(se_subset)`.

The anti-pattern that corrupts analyses silently:

```r
library(SummarizedExperiment)

# WRONG: subsetting the bare matrix breaks the link to colData.
# Downstream code that pairs assay columns with colData rows by POSITION
# now attributes each sample's data to the wrong sample.
counts <- assay(se)
counts <- counts[, se$dex == "trt"]
stopifnot(ncol(counts) == nrow(colData(se)))  # fails -- desynced

# RIGHT: subset the container; assays and colData move together.
se_trt <- se[, se$dex == "trt"]
stopifnot(ncol(assay(se_trt)) == nrow(colData(se_trt)))  # always true
```

Positional desync is the worst failure mode because nothing errors: a design matrix built from
the full `colData` applied to a filtered assay silently swaps sample labels.

## `drop` is ignored -- the container never collapses

Per the method documentation, the `drop` argument is "ignored by these methods."
`x[i, j, drop=TRUE]` always returns a `SummarizedExperiment`, even for a single row or column.
This differs from base R, where `m[, 1]` becomes a vector:

```r
sub <- se[, 1]          # still a SummarizedExperiment, dim: NFEATURES x 1
class(sub)              # "SummarizedExperiment"
dim(sub)                # e.g. 200 1

# The bare matrix follows base R rules and drops:
vec <- assay(se)[, 1]   # a vector, dimnames lost
```

If you need a single sample's values as a vector, be explicit: `assay(se[, 1])[, 1]`.
If you need a colData column as a vector, use the `[[` accessor (note: `se[[1]]` is the
FIRST colData COLUMN, not the first assay row): `se[[TREATMENT_COL]]` or `se$treatment`.

## `[, j]` keeps all features; `[i, ]` keeps all samples

`i` and `j` may each be numeric, logical, character, or missing. Missing means "keep all":

```r
se_treated   <- se[, se$dex == "trt"]   # sample filter; every feature kept
se_sig       <- se[rowSums(assay(se)) > 10, ]  # feature filter; every sample kept
se_both      <- se[rowSums(assay(se)) > 10, se$dex == "trt"]
se_by_name   <- se[c("ENSG00000000003", "ENSG00000000005"), ]  # character i
```

## Logical indexing: NA creates rows, it does not drop them

Base R subsetting semantics apply inside the container: a logical vector containing NA
produces NA placeholder rows in assays, rowData, and colData -- consistently, so alignment
holds, but the object now contains all-NA rows/columns. Clean your index first:

```r
keep <- se$dex == "trt"
keep[is.na(keep)] <- FALSE   # decide: NA samples are out, explicitly
se <- se[, keep]
```

`subset(x, subset, select)` evaluates `subset` in the context of `rowData(x)` (or
`rowRanges(x)` for a RangedSummarizedExperiment) and `select` in the context of
`colData(x)`; missing values in the logical result are taken as false there, which is safer:

```r
se_trt <- subset(se, select = dex == "trt")
```

## Character indexing and duplicated names

Character `i`/`j` values are matched against `rownames(x)`/`colnames(x)` with bounds
checking -- an unknown name errors. But if dimnames contain DUPLICATES, matching takes the
first occurrence silently. Guarantee unique names before character subsetting, or prefer
logical/numeric indices:

```r
stopifnot(!anyDuplicated(rownames(se)), !anyDuplicated(colnames(se)))
se_sel <- se[SAMPLE_IDS_OF_INTEREST, ]
```

## Assay dimnames are overlaid on extraction

`assay(x)` copies the top-level dimnames of the SE onto the returned matrix by default.
The stored matrix may have different (or NULL) dimnames underneath. Two consequences:

1. `assay(se)` and the matrix you supplied at construction can LOOK identical while the
   stored object differs -- use `assay(se, withDimnames = FALSE)` to see (or pass to
   performance-critical code) the matrix as stored.
2. At construction, assay dimnames must be NULL or identical to the object's dimnames, and
   the constructor validates counts: "nb of cols in 'assay' must equal nb of rows in
   'colData'". If your counts and metadata came from different files, reorder first --
   order matters, the constructor does not reorder for you:

```r
counts <- counts[, rownames(coldata)]      # align columns to colData order
rowranges <- rowranges[rownames(counts), ] # align rows to assay order
stopifnot(identical(rownames(coldata), colnames(counts)))
se <- SummarizedExperiment(assays = list(counts = counts),
                           rowRanges = rowranges, colData = coldata)
```

## metadata() does not subset

`metadata()` is a plain list describing the whole experiment (design formulas, processing
parameters). It is NOT subset-aware: any per-sample or per-feature vector you stash there
goes stale the moment you subset. Keep per-feature data in `rowData`, per-sample data in
`colData`, experiment-wide data in `metadata`.

## Replacement and combining

- `x[i, j] <- value` requires `value` to be a SummarizedExperiment with dimensions,
  dimnames, and assays consistent with the replaced region. You cannot inject a bare matrix.
- `cbind(...)` combines samples: duplicate `rowData` columns across objects must contain
  identical data; assays combine by name matching. `rbind(...)` combines features with the
  mirror constraint on `colData`. `metadata` from all objects is concatenated into a list
  with no name checking -- conflicting entries are silently stacked.

## Quick alignment audit

After any chain of manipulations, run this before modeling:

```r
audit_se <- function(x) {
  stopifnot(
    all(vapply(assays(x), function(a) identical(dim(a), dim(x)), logical(1))),
    nrow(colData(x)) == ncol(x),
    nrow(rowData(x)) == nrow(x),
    is.null(colnames(x)) || identical(rownames(colData(x)), colnames(x)) || is.null(rownames(colData(x)))
  )
  invisible(TRUE)
}
audit_se(se)
```