Gene expression data, for people who have never seen any

What is in the matrix, where it comes from, how to get it, and what to check before you compute anything.

Before reading: Basic R; no molecular-biology background is assumed.

1 The measurement question

This is the gene-expression branch of the Omics data and measurement overview. The transferable task is to identify the biological unit, the feature, the scale, and the measurement process. RNA counts and fluorescence intensities are not interchangeable; neither directly reports an absolute molecule count per cell without further calibration.

For a later feature comparison, write the target on the scale actually analyzed, for example E[Y_g\mid A=1]-E[Y_g\mid A=0] for log-expression Y_g in two specified populations. For a network edge, specify which population and adjustment define the association. Different filtering, normalization, or cell-composition choices can change that interpretation. Several specimens from one person remain clustered measurements even if they occupy separate columns.

The other assay primers develop parallel measurement questions. The biological and lupus examples here illustrate expression data; they are not prerequisites for the statistical guides. This page builds toward two conclusions: an expression value describes a particular assay representation, and a bulk association can reflect measurement noise or cell mixtures as well as within-cell regulation.

2 RNA and the central dogma

Every analysis of gene expression starts from a table of numbers, one per gene per sample. This section assumes no biology and explains where those numbers come from.

Think of a cell as a workshop with one very long instruction manual. The workshop never works from the manual itself. It photocopies the pages it needs right now, and the machinery that builds things reads the photocopies. A gene expression experiment measures the pile of photocopies.

DNA is the manual. DNA is a very long string over a four-letter alphabet (A, C, G, and T). A typical nucleated human cell carries two sets of 23 chromosomes. Different cell types share broadly the same genome but use different parts of it. A gene is a stretch that can be transcribed into RNA. About 20,000 human genes code for proteins; many other genes produce RNAs that are not translated.

RNA is the photocopy. RNA is a chain of letters like DNA, with U in place of T, and is usually single stranded. A cell can make many RNA copies from a gene. Their abundance at a moment depends on both production and breakdown; RNA lifetimes vary widely. Slower degradation can increase RNA abundance even if the production rate stays fixed. An abundance measurement therefore does not by itself measure a transcription rate.

The central dogma says how information moves: from DNA to RNA to protein (Figure 1).

Figure 1: The familiar DNA–RNA–protein pathway: transcription produces messenger RNA, and translation reads it three letters at a time to build a protein. The toy coding sequence ATG GCA TTC becomes AUG GCA UUC; these triplets specify methionine, alanine, and phenylalanine. The dashed RNA-to-DNA arrow is reverse transcription, used in RNA-seq and also occurring naturally, for example in retroviruses. The diagram illustrates a protein-coding transcript, not every RNA’s fate.

1. Transcription copies a gene from DNA into RNA. An enzyme, RNA polymerase, reads the gene and builds the RNA copy. In human cells the copy is then edited: splicing cuts out stretches called introns and joins the remaining pieces (exons), and most finished copies receive a tail of repeated A’s, the poly-A tail.
2. Translation reads the finished RNA three letters at a time. A ribosome links together the amino acids those triplets specify, and the chain folds into a protein. Proteins do most of the cell’s work: they are the enzymes, the structural parts, and the signaling machinery.

The arrow can run backward in the lab (the dashed arrow in Figure 1). Retroviruses carry an enzyme, reverse transcriptase, that copies RNA into DNA, and RNA sequencing borrows it.

Not all RNA is the same. The RNA that carries a protein’s instructions is messenger RNA (mRNA). Most of the RNA in a cell is something else, mainly ribosomal RNA (rRNA), a structural part of the ribosome, along with transfer RNAs, microRNAs, and long non-coding RNAs. Sequencing protocols either capture the mRNA by its poly-A tail or remove the rRNA, because otherwise most of the reads would be spent on rRNA.

Why measure RNA? An expression assay provides a broad view of transcript abundance across many genes in one experiment. It complements measurements of proteins and other molecular layers. RNA abundance is an imperfect proxy for protein abundance because translation, protein turnover, and other regulation intervene. A transcript contrast and a protein contrast are related scientific questions, not the same measurement.

Same DNA, different pages. A neutrophil and a T cell express different sets of genes. Transcription factors and signaling pathways can coordinate many genes. Interferon signaling, often prominent in lupus, can raise a set of interferon-stimulated transcripts together. Such shared activity is one source of co-expression, statistical co-variation across the sampled people. Cell composition and technical effects can produce co-variation too. A co-expression edge therefore does not establish that one gene directly regulates another.

2.1 How RNA-seq turns RNA into numbers

In a common short-read RNA-seq workflow, RNA is converted to a DNA library, fragments are sequenced, and their origins are inferred from reference sequences. The resulting count summarizes assigned reads or fragments. It is an indirect measurement shaped by capture, library preparation, and mapping. Figure 2 illustrates this workflow; direct-RNA and other protocols differ.

Figure 2: How RNA-seq turns a blood sample into one column of the count matrix. Colors mark which gene a molecule came from (genes A, B, and C). The reads in step 4 are gray because the sequencer does not know where they came from; they take on a gene’s color only when they are aligned in step 5. The counts are toy numbers.

An aligner such as STAR or HISAT2 can map reads, and a tool such as featureCounts can summarize assignments. Transcript quantifiers such as Salmon and kallisto estimate abundances using their own mapping and statistical procedures; estimated counts can be fractional. Paired-end reads may be counted as one fragment, and UMI-based protocols can count deduplicated captured molecules. Record what the provided entries count, rather than inferring it from the file’s name.

• A count is a sample, not a census. Conditional on a fixed read budget and assignment probabilities, a gene count has binomial sampling variation; a rare-feature Poisson approximation is useful in some settings. Variation across people and preparation can add dispersion. The negative-binomial model below summarizes extra variation rather than counting all cellular molecules.
• Read depth is not total cellular RNA. The library total (12 and 8 in the toy table) reflects sequencing allocation, library yield, and processing. A larger total alone does not establish more RNA per cell. Between-specimen comparisons need exposure and normalization assumptions.
• Counts share a finite budget. If one feature takes a larger share, other shares can fall at unchanged absolute amounts. TMM and related methods estimate useful reference-based corrections under stability assumptions; they do not identify every global absolute change.

Every step in the picture is also a place where technical differences creep in: RNA quality, the day the library was made, the sequencing lane. Those differences are the batch effects described below.

3 The matrix

A gene-expression matrix usually puts features in rows and specimens in columns. A feature can be a gene, transcript, or array probe, and entry (g,i) records its assay value in specimen i. Annotation and filtering determine the number of rows; study design determines the number of columns. The operational quantity is a count, intensity, or derived representation, with an additional measurement argument needed for absolute RNA abundance.

Two more pieces of context will matter throughout: what a sample is, and the disease used as the running example.

• A bulk specimen, such as blood or tissue, contains many cells. Its RNA combines their contributions. The weights for a relative RNA measurement depend on RNA amount and recovery, not only cell counts: equal numbers of two cell types need not contribute equal RNA. A change in mixture can change many measured genes without changing expression within either cell type.
• Lupus (SLE) is the running example for several expression datasets. Blood studies often report elevated interferon-stimulated transcripts such as IFI44L, RSAD2, and IFI27, with substantial patient heterogeneity. The pattern depends on tissue, disease state, treatment, and assay; it is useful biological context rather than a compulsory pattern for every patient or dataset.

For an ideal relative bulk measurement, write q_g^{\mathrm{bulk}}=\sum_c w_cq_{gc}, where q_{gc} is gene g’s RNA fraction within cell type c and w_c is that cell type’s share of the assayed RNA. If two cell types contribute equal cell numbers but one contributes four times as much RNA per cell, their RNA weights are 1/5 and 4/5 under equal recovery. This is why a cell-count proportion is not automatically the right RNA mixture weight. Deconvolution needs assumptions connecting the reference and bulk measurement.

4 Two technologies, two kinds of numbers

Microarrays. Probes bind labeled target material, and a scanner records fluorescence. Intensities are continuous; preprocessing such as RMA on suitable Affymetrix platforms commonly produces log2 summaries. The protocol and processing record determine the scale, not an assumed numerical range. Several probes may map to one gene, and a probe may have ambiguous or outdated mapping. Check the platform-specific annotation and the summary supplied.

RNA sequencing. Direct read/fragment counts are nonnegative integers; abundance estimators can produce noninteger estimated counts. Low-abundance features often have zeros. A count model can use the original counts while accounting for effective library exposure through an offset; it need not first log-transform them. Transformed matrices serve different tasks, such as exploratory plots or explicitly defined correlations.

A common working model has conditional mean \mu_{gi} and dispersion \alpha_g:

\operatorname{Var}(C_{gi})=\mu_{gi}+\alpha_g\mu_{gi}^2.

The first term resembles Poisson sampling variance; the second allows extra variation. Dividing by \mu_{gi}^2 gives squared relative SD 1/\mu_{gi}+\alpha_g. More counts reduce the first component, but a dispersion floor remains. An appropriate design still needs independent biological replication to estimate population differences.

# What RNA-seq counts look like: a negative binomial with mean-dependent overdispersion.
set.seed(1)
mu <- 2^seq(0, 14, length.out = 2000)                 # gene means, from 1 to ~16000
counts <- rnbinom(2000, mu = mu, size = 1 / 0.1)      # dispersion 0.1 (moderate; human cohorts are often higher)
plot(log2(mu), log2(counts + 1), pch = ".", xlab = "log2 true mean", ylab = "log2(count + 1)")
abline(0, 1, col = "grey")

Figure 3: One simulated negative-binomial count at each of 2,000 specified means, with dispersion 0.1. The reference line is y = x on the displayed log2 scales; it is not a fitted regression. Low means produce zeros and discrete counts, while relative variation remains at high means.

Each dot is one negative-binomial draw at its specified mean, not a repeated-measures estimate of a gene’s variance. At low means, integer sampling and zeros strongly affect the log plot. At high means the relative dispersion remains; the points do not collapse onto the reference line. The plot illustrates the measurement model, not an empirical claim about every cohort.

A dataset’s technology and supplied processing stage determine the appropriate next steps. Public archives contain both array and sequencing cohorts; publication date or the word “expression” does not establish the data scale.

5 From counts to something you can correlate

For an RNA-seq correlation analysis on a defined transformed scale, work through four decisions. A count-model differential-expression analysis instead uses its count likelihood, design, and offsets; the transformed values below are not replacements for its count input.

Filter with the design in mind. Very sparse features can provide little information for the intended contrast and unstable transformed correlations. A feature expressed mainly in one group is a different case from one barely measured anywhere. edgeR’s filterByExpr() uses library sizes and the supplied group/design to determine adequate expression across relevant samples, with additional total-count and large-group rules. Its thresholds are parameters, not a universal number of retained genes. Record the rule and resulting family. The edgeR reference manual gives its exact definition.

Define the exposure and reference. Plain counts per million use \mathrm{CPM}_{gi}=10^6C_{gi}/L_i, where L_i is the column’s library total. TMM estimates a factor f_i and uses effective library size L_i^{\mathrm{eff}}=L_if_i. It can enter a count model as an offset or a log-CPM calculation without overwriting the raw matrix. The factor is estimated from relative feature behavior under reference/stability assumptions. DESeq2’s median-of-ratios factors also use a reference-based construction; neither alone identifies a shared absolute increase in all RNA.

TPM and FPKM/RPKM additionally account for feature length and use different denominators. TPM is a length-adjusted composition whose within-sample sum is one million. Between-sample comparisons of such values concern that relative representation and its assumptions. They are not interchangeable with integer counts in a count-model likelihood, and length correction alone does not provide molecules per cell. Importing transcript-quantifier output needs the relevant method’s handling of estimated counts and effective lengths.

Specify a transformation. \log_2(\mathrm{CPM}+c) gives an observed log-scale variable when the offset c>0 is fixed. It keeps zeros finite and changes the low end most strongly; it does not turn a zero into evidence of biological absence. The next chunk adds c in CPM units. edgeR’s prior.count uses a library-scaled count addition with a corresponding library-size adjustment, so it is not the same operation as adding the same integer to every cell. DESeq2’s vst and rlog provide other transformations aimed at reducing mean-dependent variance under their models.

A transformation need not make all genes identically noisy or all observations normally distributed. For correlation, the target is the correlation of the specified representation. Recovering latent biological correlation from noisy measurements needs additional assumptions. The simulation below isolates one way noise affects that comparison.

Check the mean-variance relationship after transforming. Plot per-gene standard deviation against per-gene mean on the transformed scale. A strong trend means the noise level depends on how highly a gene is expressed, and your correlations will inherit that. The next chunk runs steps 2 to 4 by hand in base R on a simulated count matrix, so you can see the plot before you have real data. The edgeR version for real data follows it.

# Steps 2 to 4 by hand: 20,000 genes x 20 samples, with sample depths that differ several-fold.
set.seed(2)
n <- 20
mu <- 2^runif(20000, -3, 13)                                  # gene means: many low, a few high
depth <- exp(rnorm(n, 0, 0.4))                                # per-sample depth factors
counts <- matrix(rnbinom(20000 * n, mu = outer(mu, depth), size = 1 / 0.16), nrow = 20000)
lib_size <- colSums(counts)                                   # library sizes, in reads
cpm <- t(t(counts) / lib_size * 1e6)                          # normalize: counts per million
offset <- 1                                                   # try 0.1 and 5
log_cpm <- log2(cpm + offset)                                 # transform
gene_mean <- rowMeans(log_cpm)
gene_sd <- apply(log_cpm, 1, sd)
plot(gene_mean, gene_sd, pch = ".",
     xlab = "mean of log2(CPM + offset)", ylab = "SD of log2(CPM + offset)")
lines(lowess(gene_mean, gene_sd), col = "red")

Figure 4: Across-specimen SD versus mean of log2(CPM + 1) for 20,000 simulated genes in 20 specimens, with a red smooth trend. At the lowest abundances, the offset compresses both the mean and SD toward zero. This is a simulation of the chosen representation, not a universal mean–variance curve.

The plot summarizes simulated across-specimen SD for each gene after normalization and transformation. With offset 1 CPM, nearly unobserved genes are often mapped close to \log_2(1)=0, so their SD can be small because the offset dominates. A small SD here does not imply precise measurement of a biological abundance. Changing the offset can alter both this pattern and correlations at the low end. Compare offsets in optional practice rather than treating the displayed curve as a universal property.

Now consider correlations. The next chunk gives two latent log-expression variables correlation 0.6, exponentiates them to obtain count means, and generates conditionally independent negative-binomial noise. It repeats this experiment at five mean-count levels. The reported correlation is between \log_2(\mathrm{count}+1) values, so its target differs from correlation of the latent log variables.

# Observed correlation between two genes with the same underlying correlation, by mean count.
ObsCor <- function(m, n = 200, rho = 0.6, reps = 200) {
  mean(replicate(reps, {
    z1 <- rnorm(n)
    z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n)               # correlated latent expression
    x <- rnbinom(n, mu = m * exp(0.4 * z1 - 0.08), size = 1 / 0.16)
    y <- rnbinom(n, mu = m * exp(0.4 * z2 - 0.08), size = 1 / 0.16)
    cor(log2(x + 1), log2(y + 1))
  }))
}
set.seed(3)
means <- c(1, 5, 20, 100, 1000)
round(setNames(sapply(means, ObsCor), means), 2)

   1    5   20  100 1000 
0.07 0.16 0.24 0.27 0.29 

The latent log correlation stays fixed while the average observed correlation rises from roughly 0.07 at mean count 1 to roughly 0.29 at 1,000. In this experiment, weaker observed edges at low counts arise without a change in the latent dependence. At high means, negative-binomial dispersion still adds noise, so more counts do not recover the latent log correlation exactly. Exponentiating also changes Pearson correlation; the latent abundances themselves need not have correlation 0.6.

A simpler additive-error model explains the attenuation idea. If X=U+e and Y=V+d, with centered errors independent of the latent variables and each other, then \operatorname{Cov}(X,Y)=\operatorname{Cov}(U,V) while both variances increase. Consequently

\operatorname{Corr}(X,Y)=\frac{\operatorname{Corr}(U,V)} {\sqrt{(1+\sigma_e^2/\sigma_U^2)(1+\sigma_d^2/\sigma_V^2)}}.

Here \sigma_U^2,\sigma_V^2>0 are latent variances and \sigma_e^2,\sigma_d^2 are error variances. This exact identity uses the additive independent-error assumptions; it is an analogy, not the formula for the transformed negative-binomial experiment. Correlated technical errors or normalization can create associations as well as weaken them.

On real data, the standard path uses edgeR, which normalizes with TMM instead of plain library-size scaling. Here count_matrix and sample_group stand for your data.

# The standard path with edgeR (install from Bioconductor: BiocManager::install("edgeR"))
library(edgeR)
y <- DGEList(counts = count_matrix, group = sample_group)
keep <- filterByExpr(y, group = sample_group)
y <- y[keep, , keep.lib.sizes = FALSE]
y <- calcNormFactors(y, method = "TMM")
logcpm <- cpm(y, log = TRUE, prior.count = 2)     # genes x samples, log2 scale

For arrays, verify scale and preprocessing from the platform and processing record before transforming again. Numerical ranges and boxplots are diagnostic clues, not proof: matching distributions can be imposed by quantile normalization, and differing distributions can have biological causes. A gene-level summary needs a prespecified probe-mapping and aggregation rule, including treatment of ambiguous or discordant probes. Keeping probe-level features is another defined analysis. Selecting whichever probe gives the largest group difference is not an innocuous annotation step.

6 Things that will bite you

• Batch effects. Inspect processing dates, sites, chips, and runs, and cross-tabulate them with group. PCA colored by batch can reveal structure, but association with a principal component does not identify its cause. Batch covariates, ComBat, and removeBatchEffect() serve different modeling/preprocessing roles; check the intended design and preserve the contrasts of interest. Perfect group–batch confounding prevents their direct separation without additional information. A correction can also change covariance, so a network analysis needs to assess the resulting representation.
• Cell composition. Blood cell mixtures vary with disease and treatment. In an ideal mixture, two genes enriched in the same cell type can co-vary as that type’s contribution changes, even with fixed within-type expression. Reference-based deconvolution estimates composition under measurement and reference assumptions. Measured blood counts provide other information but need specimen matching and an argument connecting cell counts to RNA contributions. Estimated composition is not an exact confounder or a cell-type-specific outcome.
• Sex and age. Check their distributions by group and their role in the intended contrast. XIST and Y-chromosome transcript patterns can help investigate sample identity, but expression depends on tissue, assay, and biology; they are not infallible classifiers. Investigate discrepancies using the original records instead of automatically overwriting metadata.
• Outliers and failed samples. PCA or clustering can reveal unusual specimens. Determine whether they reflect assay failure, a metadata error, or genuine heterogeneity. Prespecify defensible exclusion rules where possible and report their consequences for the analyzed population. Maximizing group separation is not a QC criterion.
• Scale. Always ask: counts, CPM, log-CPM, TPM, or array intensities, and on which base of logarithm? GEO deposits are inconsistent, and “normalized” in the submitter’s description can mean any of several things.
• Repeated samples. Longitudinal cohorts have several samples per subject. Choose the independent unit using Units of analysis, replication, and study design; retain all visits from one person together when people are the independent draws. Check the metadata for subject identifiers; they are sometimes only in the sample title string. GSE65391 has a subject field; in GSE49454 the patient appears only in the title (C2PL28V2 is patient C2PL28 at visit 2).

7 Where the data come from

Public repositories supply real gene-expression datasets for research and practice. The examples below explain how to obtain and inspect them. Elsewhere, controlled methods examples use synthetic data so that population truth is known; keep that distinction when interpreting a result.

GEO (Gene Expression Omnibus, NCBI). Three identifiers matter: GSE identifies a series, GSM a deposited sample record, and GPL a platform. A GSM can represent a specimen, technical measurement, or multiplexed pool; it is not automatically one independent person. A series can provide processed matrices, supplementary files, and links to raw data. Read the file descriptions and processing fields before importing. RNA-seq supplementary files can hold counts, RPKM, or another representation, so downloading the right study is not enough to establish the scale.

# Pull a GEO series in R. install.packages("BiocManager"); BiocManager::install("GEOquery")
library(GEOquery)
dir.create("data/raw", recursive = TRUE, showWarnings = FALSE)
gse <- getGEO("GSE65391", GSEMatrix = TRUE, destdir = "data/raw")   # a list; inspect its platforms before selecting
eset <- gse[[1]]
expr <- exprs(eset)                                  # genes/probes x samples
meta <- pData(eset)                                  # one row per sample: title, characteristics, etc.
anno <- fData(eset)                                  # probe annotation, including gene symbols
dim(expr); head(meta[, head(grep("characteristics", names(meta)), 5)])   # the first five metadata columns
stopifnot(identical(colnames(expr), rownames(meta)))   # verify specimen alignment
summary(as.vector(expr[, 1]))                        # inspect range; establish scale from processing metadata

Metadata in GEO are free text in “characteristics” columns and need cleaning every time. Each cell holds text of the form key: value, so strip the key first (for example with sub("^[^:]*: ", "", x)). Write the cleaning as a script, never by hand, so that anyone can rerun it.

For this RNA-seq example, the deposit supplies a specifically named RPKM file. Select it explicitly rather than reading whichever supplementary file happens to appear first; another file could be an archive or another assay. This example inspects a processed representation, not input for the count-model path above.

supp <- getGEOSuppFiles("GSE72509", baseDir = "data/raw")   # returns the paths of the files it saved
paths <- rownames(supp)
chosen <- paths[basename(paths) == "GSE72509_SLE_RPKMs.txt.gz"]
stopifnot(length(chosen) == 1L)
supp_table <- read.delim(gzfile(chosen), check.names = FALSE)

SRA (Sequence Read Archive). Raw sequencing resources use identifiers such as SRR (run), SRX (experiment), and SRP (study). Reprocessing can be useful when available summaries do not match the needed feature definition or pipeline. Several runs can come from one specimen; preserve that mapping instead of treating runs as biological replicates.

recount3 and ARCHS4 provide large reprocessed expression resources. A common processing pipeline can reduce one source of heterogeneity, but studies can still differ in collection, assay design, population, and metadata. Inspect each resource’s actual units and construction before substituting it for a deposited matrix.

ArrayExpress / BioStudies (EMBL-EBI). The European counterpart to GEO; many GEO array studies are mirrored there (E-GEOD-xxxxx), and some European studies are only there.

ImmPort. NIH’s immunology data repository: clinical trials and cohort studies in immunology with expression, cytometry, and clinical data, including lupus studies. Registration required; worth it for immunology and autoimmune disease.

GTEx is a multi-tissue donor resource. Tissue and donor selection matter when using it as a reference; it is not a universally exchangeable healthy-control group for another study. Several tissues from the same donor also retain a shared donor identity.

TCGA, through the GDC portal, provides cancer molecular and clinical data. Tumor, adjacent tissue, and external control samples define different comparisons. Check specimen types, person IDs, and processing compatibility before combining them.

Single-cell resources, including CELLxGENE and the Human Cell Atlas, provide cell-level data and annotations. Millions of cells can still come from a modest number of independently sampled donors. For a donor-population contrast, preserve donor and specimen IDs; a cell-type label alone is not enough. The Perez et al. lupus study below illustrates that distinction.

Access conditions. Some molecular and clinical records require controlled access through resources such as dbGaP or the European Genome-phenome Archive. A public study record does not mean every raw file or covariate is openly downloadable. Establish which files and identifiers are available before choosing a research question around them.

8 Lupus cohorts you will meet

These examples illustrate different observation units and deposited scales. Counts below describe the cited study/deposit, not a promise about the analysis-ready subset. Reconstruct people, visits, and technical repeats from the actual imported metadata before fitting a model.

• GSE65391. Banchereau et al. (2016, Cell), longitudinal pediatric whole-blood arrays. GEO lists 996 arrays: 924 SLE measurements from 158 patients and 72 control measurements, including 24 technical control replicates. The deposited total is therefore neither the number of people nor the number of distinct biological specimens. The linked clinical metadata make repeated visits, disease activity, and blood composition useful audit questions.
• GSE49454. Chiche et al. (2014, Arthritis and Rheumatology), adult whole-blood arrays. The study reports 157 SLE samples from 62 patients, alongside 20 controls. Compare the publication’s patient count with the deposit’s sample count; visit labels need explicit parsing and validation before person-level inference.
• GSE88884. Baseline arrays from the ILLUMINATE-1 and ILLUMINATE-2 tabalumab trials, with 1,760 SLE patients and 60 controls. The parent study includes later visits, but this accession’s stated comparison is baseline. Two trial sources plus separately obtained controls require careful technical and population comparisons; a trial label is not automatically a pure batch variable.
• GSE72509. Hung et al. (2015, Science), whole-blood RNA-seq with 117 sample records. The named supplementary matrix contains RPKM, not raw counts. Verify group labels, available covariates, and person identifiers from the deposit; a count-based likelihood cannot be applied merely because the original assay was sequencing.
• GSE174188. Perez et al. (2022, Science), more than 1.2 million PBMCs from 162 cases and 99 controls. GEO lists 88 multiplexed pool records and directs users to controlled-access data. Neither 88 pools nor 1.2 million cells equals the 261 donor units. Identify an accessible data product with the necessary donor, specimen, and cell-type mappings before defining a donor-level analysis.

Record whatever you find on first contact in the project README for the next person.

9 A first-contact checklist

When you open a new dataset, before anything else:

1. Technology and scale. Array or RNA-seq; counts, CPM, log, or intensities; which log base.
2. Dimensions. Genes by samples; how many genes survive filtering.
3. Metadata. Group labels, subject identifiers, batch, sex, age, time point. How many samples per group, per subject, per batch. Some deposits, such as GSE72509, have almost none of these.
4. Sample distributions and PCA. Compare technical and biological labels, remembering that a pattern is a diagnostic clue, not proof of its cause.
5. Identity checks. Compare available metadata, replicate relationships, and informative expression patterns; investigate discrepancies using the source records.
6. Biological context. Inspect relevant gene programs, including interferon-related transcripts in suitable lupus blood cohorts. Unexpected patterns prompt investigation of tissue, treatment, assay, labels, and heterogeneity; they do not by themselves prove the data are wrong.
7. Write everything in steps 1 to 6 into a script that produces a short report, and save the cleaned matrix and metadata together with the script that made them.

10 Exercises

Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.

Checkpoint · OM-01 — Identify measurements and biological replication

Task type: Interpretation.

A raw RNA-seq matrix has 15,000 genes and 30 specimen columns. Ten people each contributed three visits; five people are controls and five are patients. A metadata table gives person ID, visit, group, sequencing run, and library size.

Identify the feature, specimen, and independent biological unit under an independent-person sampling model. Explain why 30 columns are not 30 independent people, why counts require a scale/preprocessing decision before comparison, and what identifiers must be checked before joining the matrix to metadata.

Ready to continue when: you can distinguish measurement units from replication and describe the information needed to interpret the matrix.

Transfer: What changes if the three columns are technical repeats of the same specimen rather than three visits?

Hint

One matrix column records one specimen; independence is a statement about the sampling design.

Solution outline — open after trying

The features are measured genes, the columns are visit specimens, and the independent biological units in this design are ten people. The three visits within each person can be dependent. Library sizes and other measurement differences affect raw counts, so the target scale and normalization assumptions must be specified. Match each specimen column to one metadata record using specimen IDs; retain person IDs for clustering and check visit/group consistency, duplicated IDs, missing joins, and run–group confounding.

Additional practice (optional)

These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.

1. Pull GSE72509 with GEOquery and inspect which components contain measurements. Find the named supplementary file, establish its scale, identify controls in the metadata, and list the checklist items the available record does not support. Do not generalize one deposit’s file layout to all RNA-seq series.
2. Using the simulated count matrix above or any RNA-seq count matrix, compute log-CPM with offsets of 0.1, 1, and 5 CPM and plot the per-gene SD against the mean for each. Describe what the offset does to low-count genes and to correlations among them.
3. Simulate a two-cell-type mixture: 200 genes specific to type A, 200 to type B, mixing proportions varying across 100 samples. Compute the gene-gene correlation matrix and describe its block structure. Nothing in this simulation involves gene regulation; where did the network come from?
4. In GSE65391, identify the subject identifier and tabulate visits per subject. How unequal is it? This is the design the permutation guide is about. Repeat for GSE49454, where you must build the subject identifier from the sample title.

11 Reading

• Conesa, A., et al. (2016). A survey of best practices for RNA-seq data analysis. Genome Biology, 17, 13.
• Law, C. W., et al. (2016). RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR. F1000Research, 5, 1408. A worked pipeline you can run.
• Robinson, M. D. and Oshlack, A. (2010). A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology, 11, R25. The TMM paper.
• Leek, J. T., et al. (2010). Tackling the widespread and critical impact of batch effects in high-throughput data. Nature Reviews Genetics, 11(10), 733 to 739.
• Banchereau, R., et al. (2016). Personalized immunomonitoring uncovers molecular networks that stratify lupus patients. Cell, 165(3), 551 to 565. Read the first three figures to see what a lupus cohort looks like.
• Davis, S. and Meltzer, P. S. (2007). GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinformatics, 23(14), 1846 to 1847.