Choosing a dataset

In this tutorial we use one of the datasets available from the NCBI GEO database. You are encouraged to follow along with the tutorial using the Grady Trauma Project, and once you’ve gotten some experience following along, then practice with another dataset. The NCBI GEO database can be a bit tricky to navigate, so for ease, I’ve curated a catalog of 70 Illumina MethylationEPIC GEO Series (including the dataset used in the tutorial) that carry enough phenotype and covariate metadata to support a real association analysis. Browse it to understand the landscape of public EPIC data and to pick your own dataset for a project.

The full table is downloadable as a CSV at the bottom of the page so you can filter and load it in R or Python yourself.

The catalog

Each row is one study. A few were deposited under several GEO accessions — a SuperSeries alongside its SubSeries — and those appear as a single row with every accession listed, so the 70 Series below collapse to 60 studies. The table is searchable and sortable — click a column header to sort, or type in the box to filter (e.g. type Strong to see only the top-tier studies, or a gene/phenotype term to find studies on a topic). The Grady Trauma Project (GSE132203) the dataset used throughout this tutorial, is highlighted.

What the columns mean

Column	Meaning
`GEO accession`	GEO Series accession — the identifier you give to `getGEO()` or use in the FTP path to download IDATs. Where a study was deposited as a SuperSeries plus SubSeries, every accession is listed here and the row’s other columns describe the first one.
`Title`	The study’s deposited title.
`Suitability tier` / `Suitability score`	Our teaching-suitability grade (Strong / Moderate / Limited) and the underlying score, based on metadata completeness (sex, age, phenotype, raw data). Higher is better for a first EWAS.
`Methylation array`	Array generation(s) present. Most are EPIC v1 (850K); a few mix in 450K or EPIC v2, which matters because probe IDs differ across generations.
`N samples`	Total samples across all matrix files for the first accession listed (super-series are summed).
`Has sex`, `Has age`	Whether sex and age are in the sample characteristics — both are strong methylation covariates you will almost always model.
`Has ancestry`	Whether genetic ancestry / race is reported (relevant for ancestry-aware probe masking, as in Probe filtering).
`Has cell composition`	Whether cell-type proportions are deposited — not necessary, but saves a step in data pre-processing if it is available.
`Phenotypes`, `N phenotypes`	The variables you could actually test — every phenotype and exposure field the submitter deposited, named but without their values. A cell reading none deposited means exactly that: the study’s characteristics carry nothing beyond sex, age, tissue and identifiers, so there is no exposure or outcome here to model. Sample and array identifiers, estimated cell proportions, principal components, derived age-acceleration metrics and technical flags are excluded, so this is a shortlist of candidate EWAS outcomes rather than a dump of the characteristics. Names are the submitter’s own except where this tutorial defines a plainer one: Grady’s `mergedcapsandpsswinthin30days` and `childabphyssexemot_ctq_01modandsev` appear as PTSD and Childhood abuse, matching the derived columns in Setup.
`PubMed`	Linked publication, where one exists.

What makes a good first EWAS dataset

Sort by Suitability_tier, then look for: raw IDATs deposited (so you can start from raw data), both age and sex present, at least one genuine phenotype you find interesting, and ideally whole blood or PBMC so you get to practice cell-type deconvolution. Deposited cell proportions values are a bonus — they give you an answer key to check your own pipeline against, or you can skip this step and just use the deposited cell proportions.

How the catalog was assembled

The 70 Series are the result of a deliberate screen:

Exported all GEO Series for the EPIC platform (GPL21145) from the GEO browser.
Restricted to Homo sapiens.
Kept studies with ≥ 150 samples (enough for a stable EWAS).
Manually annotated the phenotype of interest from each title.
Manually removed obvious non-blood tissue and cell-differentiation studies.

Metadata was then read directly from the Series Matrix file header (GSE*_series_matrix.txt.gz) — source name, characteristics tag–value pairs, platform, and sample count. This is a quick and easy way to profile GEO studies without committing to downloading all the data.

Figure 1: The 70-Series catalog at a glance. (a) Teaching-suitability tiers. (b) How many studies carry each metadata field. (c) Study sizes span 152 to ~3,100 samples. (d) Almost all are EPIC v1 (850K).

The tutorial dataset

From the shortlist we chose GSE132203, the Grady Trauma Project, because it is unusually complete for teaching: raw IDATs are deposited, it has explicit age and sex plus several trauma/PTSD exposures, and it provides cell-type proportions, giving us answer keys to validate the deconvolution step against. There are also published manuscripts that we can compare our results with (Smith et al. 2011) (Katrinli et al. 2020).

Not every study gives you IDATs

Some GEO methylation studies deposit only a processed beta-value matrix, not raw IDATs. You can still use these for your rotation — the Normalization notebook includes notes on how to enter the workflow from a processed matrix and which steps you can and cannot still do.

Download the catalog

The catalog is available as a CSV. It holds exactly the 60 studies in the table above, in the same order, with 30 columns rather than the 13 shown:

methylation_geo_catalog.csv — ready to load with read_csv() (R) or pd.read_csv() (Python).

GSE and N_samples are comma-separated and aligned with each other, so the Congo study reads GSE224365, GSE224363, GSE224364 against 712, 356, 356. Every other column describes the first accession listed. The file is written by the same chunk that draws the table, so the download cannot drift from what you see above.

You can also grab the current view directly from the table using the CSV and Copy buttons above it.

Next: with a dataset in hand, head to Setup to meet the platform and see how the raw IDATs reach R — or jump straight to Quality control if you already know the array.

References

Katrinli, Seyma, Adriana Lori, Varun Kilaru, et al. 2020. “Association of HLA Locus Alleles with Posttraumatic Stress Disorder.” Brain, Behavior, and Immunity 87: 37–45. https://doi.org/10.1016/j.bbi.2020.04.038.

Smith, Alicia K, Karen N Conneely, Varun Kilaru, et al. 2011. “Differential Immune System DNA Methylation and Cytokine Regulation in Post-Traumatic Stress Disorder.” American Journal of Medical Genetics Part B 156 (6): 700–708. https://doi.org/10.1002/ajmg.b.31212.

--- title: "Choosing a dataset" resources: - methylation_geo_catalog.csv --- ```{r} #| label: setup #| include: false source("_setup.R") library(DT) ``` In this tutorial we use one of the datasets available from the NCBI GEO database. You are encouraged to follow along with the tutorial using the Grady Trauma Project, and once you've gotten some experience following along, then practice with another dataset. The NCBI GEO database can be a bit tricky to navigate, so for ease, I've curated a catalog of 70 Illumina MethylationEPIC GEO Series (including the dataset used in the tutorial) that carry enough phenotype and covariate metadata to support a real association analysis. Browse it to understand the landscape of public EPIC data and to pick your own dataset for a project. The full table is downloadable as a CSV at the bottom of the page so you can filter and load it in R or Python yourself. ## The catalog Each row is one study. A few were deposited under several GEO accessions — a SuperSeries alongside its SubSeries — and those appear as a single row with every accession listed, so the 70 Series below collapse to 60 studies. The table is searchable and sortable — click a column header to sort, or type in the box to filter (e.g. type `Strong` to see only the top-tier studies, or a gene/phenotype term to find studies on a topic). The Grady Trauma Project (GSE132203) the dataset used throughout this tutorial, is highlighted. ```{r} #| label: catalog-table #| echo: false # 00b_geo_catalog_source.csv is the raw scrape: one row per GEO Series, values # and all. The single published CSV (methylation_geo_catalog.csv, linked at the # foot of this page) is written from it below, so the download can never drift # from the table above it. cat70 <- read.csv("00b_geo_catalog_source.csv", stringsAsFactors = FALSE, check.names = FALSE) # --- Phenotypes: the variables you could actually test ----------------------- # Source is Phenotype_fields, which lists EVERY deposited phenotype field. # (Phenotype_examples caps at four fields and also carries their values, so # building this column from it under-reported 15 studies.) Field names are kept # as deposited; anything that is not a candidate EWAS outcome or exposure is # dropped: sample/array identifiers and plate positions, pedigree keys, batch and # provenance fields, estimated cell proportions, principal components (including # ancestry PCs), derived clock metrics, assay/library QC read-outs, inclusion # flags and dates. drop_pat <- paste0("^(", paste(c( "sample", "sample[ _-]?(id|name|title|number|group|plate)", "methylationid", "dyad", "title", "id", "geo.*", "barcode", "sentrix.*", "position", "well", "plate", "meth-plate", "slide", "chip", "array.*", "serial number", "batch", "enrollment batch", "scan.*date", "description", "supplementary.*", "center", "recruite?_?site", "(subject|donor|participant|patient|individual|pair|session|generation)([ _-]?id)?", ".*family", "replicate.*", "is_technical_replicate", "cd4t?", "cd8t?", "nk", "bcell", "mono", "neu", "gran", "eos", "baso", "plasmablast", "cd4\\+.*", "cd8\\+.*", "cell proportion.*", ".*cell proportion", "pc[0-9]+", "anc\\.?pca?[0-9]+", ".*[._]pca?[0-9]+", "age.?acceleration", "predicted.*age", "dnam.*age", ".*methylation age", "bioanalyzer.*", ".*concentration.*", ".*_pcr", "nm_pcr", "included in.*", "passed.*quality.*", "qc.*", "\$untagged\$" ), collapse = "|"), ")$") drop_pat2 <- "\\bid\\b|identifier|_id$|^pmid|untagged|(^|[^a-z])date($|[^a-z])|date$" # Two fields this tutorial gives a plain-language name to: 00_setup.qmd derives # `ptsd` from the first and `childhood_abuse` from the second. rename_map <- c( "mergedcapsandpsswinthin30days" = "PTSD", "childabphyssexemot_ctq_01modandsev" = "Childhood abuse" ) clean_pheno <- function(x) { if (is.na(x) || !nzchar(trimws(x))) return("") f <- trimws(strsplit(x, "|", fixed = TRUE)[[1]]) # submitters sometimes encode the coding in the name itself, e.g. # "art treatment status (0 = spontaneous; 1 = assisted reproduction)"; # strip parentheticals containing "=", keep units such as (kg/m^2) f <- trimws(gsub("\$[^)]*=[^)]*\$", "", f)) f <- f[nzchar(f)] f <- f[!grepl(drop_pat, tolower(f))] f <- f[!grepl(drop_pat2, tolower(f))] hit <- f %in% names(rename_map) f[hit] <- rename_map[f[hit]] paste(unique(f), collapse = ", ") } cat70$Phenotypes <- vapply(cat70$Phenotype_fields, clean_pheno, character(1)) cat70$N_phenotypes <- vapply(strsplit(cat70$Phenotypes, ", "), function(v) sum(nzchar(v)), integer(1)) # --- Collapse GEO Series that belong to one publication ---------------------- # A SuperSeries and its SubSeries are separate accessions for the same study. # Group on PubMed id, keep the highest-scoring Series as the displayed row # (ties broken by sample count, which selects the SuperSeries), and list the # sibling accessions alongside it. grp <- ifelse(is.na(cat70$PubMed), paste0("row:", seq_len(nrow(cat70))), paste0("pmid:", cat70$PubMed)) # compute the ordering ONCE and apply it to the frame and the key together -- # deriving it twice reorders the key by an index taken from the already-sorted # frame, which silently mixes unrelated studies into the same group ord <- order(-cat70$Suitability_score, -cat70$N_samples) cat70 <- cat70[ord, ] grp <- grp[ord] acc <- tapply(cat70$GSE, grp, paste, collapse = ", ") nsamp <- tapply(cat70$N_samples, grp, paste, collapse = ", ") nser <- table(grp) cat70 <- cat70[!duplicated(grp), ] cat70$Accessions <- acc[grp[!duplicated(grp)]] cat70$N_series <- as.integer(nser[grp[!duplicated(grp)]]) # per-accession sizes, aligned with Accessions; the table keeps the numeric # N_samples of the first accession so the column stays sortable cat70$N_samples_all <- nsamp[grp[!duplicated(grp)]] # Title first, then the phenotypes you could test, then everything else show_cols <- c("Accessions", "Title", "Phenotypes", "Suitability_tier", "Suitability_score", "Meth_array", "N_samples", "Has_sex", "Has_age", "Has_ancestry", "Has_cell_composition", "N_phenotypes", "PubMed") show_cols <- show_cols[show_cols %in% colnames(cat70)] # order the frame once so the table and the exported CSV share a row order tier_rank <- c(Strong = 1, Moderate = 2, Limited = 3) cat70 <- cat70[order(tier_rank[cat70$Suitability_tier], -cat70$Suitability_score), ] tab <- cat70[, show_cols] # --- the single published CSV ----------------------------------------------- # Same 60 rows in the same order as the table, with more columns. GSE and # N_samples are comma-separated and aligned with each other; every other column # describes the first accession listed. out <- cat70 out$GSE <- out$Accessions out$N_samples <- out$N_samples_all csv_cols <- c("GSE", "N_series", "Title", "Phenotypes", "N_phenotypes", "Suitability_tier", "Suitability_score", "Meth_array", "N_samples", "Has_sex", "Has_age", "Has_ancestry", "Has_cell_composition", "PubMed", "Suitability_notes", "Other_platforms", "Has_tissue", "N_phenotype_fields", "Phenotype_fields", "N_char_fields", "All_char_tags", "Sex_field", "Age_field", "Tissue_field", "Tissue_examples", "First_sample_source_name", "Submission_date", "Last_update", "Platform_GPLs", "N_matrix_files") csv_cols <- csv_cols[csv_cols %in% colnames(out)] write.csv(out[, csv_cols], "methylation_geo_catalog.csv", row.names = FALSE, na = "") pretty <- c(Accessions = "GEO accession", Title = "Title", Phenotypes = "Phenotypes", Suitability_tier = "Suitability tier", Suitability_score = "Suitability score", Meth_array = "Methylation array", N_samples = "N samples", Has_sex = "Has sex", Has_age = "Has age", Has_ancestry = "Has ancestry", Has_cell_composition = "Has cell composition", N_phenotypes = "N phenotypes", PubMed = "PubMed") idx <- function(nm) which(colnames(tab) %in% nm) - 1L # DT is 0-indexed truncate_js <- function(n) DT::JS( "function(data, type, row, meta){", sprintf("return type === 'display' && data != null && data.length > %d ?", n), sprintf("'' + data.substr(0, %d) + '…' : data;", n), "}") htmltools::tagList( htmltools::tags$style(htmltools::HTML(" /* wrap the header labels instead of forcing every column as wide as its title, which is what left the short numeric columns full of empty space */ #catalog-tbl thead th { white-space: normal; vertical-align: bottom; text-align: center; line-height: 1.2; } #catalog-tbl thead th:nth-child(2) { text-align: left; } #catalog-tbl td, #catalog-tbl th { padding-left: 6px; padding-right: 6px; } ")), datatable( tab, elementId = "catalog-tbl", rownames = FALSE, colnames = unname(pretty[colnames(tab)]), filter = "top", extensions = "Buttons", options = list( pageLength = 5, lengthMenu = c(5, 10, 25, nrow(tab)), scrollX = TRUE, autoWidth = TRUE, dom = "Blfrtip", buttons = c("copy", "csv"), columnDefs = list( # Title is the one column that earns width; everything else is centred list(targets = idx("Title"), className = "dt-left", width = "300px", render = truncate_js(140)), list(targets = idx("Phenotypes"), className = "dt-left", width = "200px", # an empty cell here is a fact about the deposit, not a gap in the # table, so say so rather than leaving it blank render = DT::JS( "function(data, type, row, meta){", "if (type !== 'display') return data;", "if (data == null || data === '')", "return 'none deposited';", "return data.length > 100 ?", "'' + data.substr(0, 100) + '…' : data;", "}")), list(targets = idx("Accessions"), className = "dt-center", width = "110px"), list(targets = idx(c("Suitability_tier", "Suitability_score", "N_samples", "N_phenotypes", "Has_sex", "Has_age", "Has_ancestry", "Has_cell_composition", "PubMed")), className = "dt-center", width = "70px"), list(targets = idx("Meth_array"), className = "dt-center", width = "100px") ) ), caption = htmltools::tags$caption( style = "caption-side: bottom; text-align: left; font-size: 90%;", sprintf(paste("%d studies (%d GEO Series) on the EPIC array, ranked by teaching", "suitability."), nrow(tab), sum(cat70$N_series)), "Series sharing a publication — a SuperSeries and its SubSeries — are shown", "as one row, with every accession listed and the numbers describing the", "first one. Hover a truncated Title or Phenotypes cell for its full text." ) ) |> formatStyle( "Accessions", target = "row", backgroundColor = styleEqual("GSE132203", "#FFF3CD") ) ) ``` ## What the columns mean | Column | Meaning | |---|---| | `GEO accession` | GEO Series accession — the identifier you give to `getGEO()` or use in the FTP path to download IDATs. Where a study was deposited as a SuperSeries plus SubSeries, every accession is listed here and the row's other columns describe the first one. | | `Title` | The study's deposited title. | | `Suitability tier` / `Suitability score` | Our teaching-suitability grade (**Strong / Moderate / Limited**) and the underlying score, based on metadata completeness (sex, age, phenotype, raw data). Higher is better for a first EWAS. | | `Methylation array` | Array generation(s) present. Most are EPIC v1 (850K); a few mix in 450K or EPIC v2, which matters because probe IDs differ across generations. | | `N samples` | Total samples across all matrix files for the first accession listed (super-series are summed). | | `Has sex`, `Has age` | Whether sex and age are in the sample characteristics — both are strong methylation covariates you will almost always model. | | `Has ancestry` | Whether genetic ancestry / race is reported (relevant for ancestry-aware probe masking, as in [Probe filtering](03_probe_filtering.qmd)). | | `Has cell composition` | Whether cell-type proportions are deposited — not necessary, but saves a step in data pre-processing if it is available. | | `Phenotypes`, `N phenotypes` | The variables you could actually test — every phenotype and exposure field the submitter deposited, named but without their values. A cell reading *none deposited* means exactly that: the study's characteristics carry nothing beyond sex, age, tissue and identifiers, so there is no exposure or outcome here to model. Sample and array identifiers, estimated cell proportions, principal components, derived age-acceleration metrics and technical flags are excluded, so this is a shortlist of candidate EWAS outcomes rather than a dump of the characteristics. Names are the submitter's own except where this tutorial defines a plainer one: Grady's `mergedcapsandpsswinthin30days` and `childabphyssexemot_ctq_01modandsev` appear as **PTSD** and **Childhood abuse**, matching the derived columns in [Setup](00_setup.qmd). | | `PubMed` | Linked publication, where one exists. | ::: {.callout-tip} ## What makes a good *first* EWAS dataset Sort by `Suitability_tier`, then look for: **raw IDATs deposited** (so you can start from raw data), **both age and sex present**, at least one **genuine phenotype** you find interesting, and ideally **whole blood or PBMC** so you get to practice cell-type deconvolution. Deposited cell proportions values are a bonus — they give you an answer key to check your own pipeline against, or you can skip this step and just use the deposited cell proportions. ::: ## How the catalog was assembled The 70 Series are the result of a deliberate screen: 1. Exported **all GEO Series for the EPIC platform** (GPL21145) from the GEO browser. 2. Restricted to *Homo sapiens*. 3. Kept studies with **≥ 150 samples** (enough for a stable EWAS). 4. Manually annotated the phenotype of interest from each title. 5. Manually removed obvious non-blood tissue and cell-differentiation studies. Metadata was then read directly from the **Series Matrix file header** (`GSE*_series_matrix.txt.gz`) — source name, `characteristics` tag–value pairs, platform, and sample count. This is a quick and easy way to profile GEO studies without committing to downloading all the data. ```{r} #| label: fig-catalog-overview #| echo: false #| fig-cap: "The 70-Series catalog at a glance. (a) Teaching-suitability tiers. (b) How many studies carry each metadata field. (c) Study sizes span 152 to ~3,100 samples. (d) Almost all are EPIC v1 (850K)." knitr::include_graphics("data/00b_catalog_overview.png") ``` ## The tutorial dataset From the shortlist we chose **GSE132203**, the Grady Trauma Project, because it is unusually complete for teaching: raw IDATs are deposited, it has explicit age and sex plus several trauma/PTSD exposures, and it provides cell-type proportions, giving us answer keys to validate the deconvolution step against. There are also published manuscripts that we can compare our results with [@smith2011grady] [@katrinli2020ptsd]. ::: {.callout-note} ## Not every study gives you IDATs Some GEO methylation studies deposit only a processed beta-value matrix, not raw IDATs. You can still use these for your rotation — the [Normalization](02_normalization.qmd) notebook includes notes on how to enter the workflow from a processed matrix and which steps you can and cannot still do. ::: ## Download the catalog The catalog is available as a CSV. It holds exactly the 60 studies in the table above, in the same order, with 30 columns rather than the 13 shown: - [**methylation_geo_catalog.csv**](methylation_geo_catalog.csv) — ready to load with `read_csv()` (R) or `pd.read_csv()` (Python). `GSE` and `N_samples` are comma-separated and aligned with each other, so the Congo study reads `GSE224365, GSE224363, GSE224364` against `712, 356, 356`. Every other column describes the **first** accession listed. The file is written by the same chunk that draws the table, so the download cannot drift from what you see above. You can also grab the current view directly from the table using the **CSV** and **Copy** buttons above it. --- **Next:** with a dataset in hand, head to [Setup](00_setup.qmd) to meet the platform and see how the raw IDATs reach R — or jump straight to [Quality control](01_qc.qmd) if you already know the array.