Setup: getting the data onto your machine

This notebook introduces the dataset the tutorial uses, shows how to get the raw files onto your own machine, and sets out what running everything costs in time, memory, and disk. It assumes the array background from Background.

Where to type each block

Every code block in this chapter carries a grey label telling you where it runs. There are only two places, and mixing them up is the most common way to get stuck:

RStudio Console — the pane with the > prompt, where R lives. Every block of R code goes here.
Terminal / PowerShell — a system shell, for the handful of non-R commands such as git clone. On macOS that is Terminal.app, on Windows it is PowerShell, on Linux it is whatever shell you use.

One trap worth naming: RStudio has a Terminal tab sitting right next to the Console, and it is a system shell, not R. R code pasted there will fail with a string of syntax errors — so when a block says RStudio Console, make sure you have clicked the Console tab. If you are not using RStudio, R at a terminal gives you the same > prompt and works fine.

One thing that runs before each chapter, and what it is not allowed to do

Every chapter begins with a hidden chunk that does source("_setup.R"). That file is deliberately limited to two housekeeping jobs: it puts the library holding the EPIC annotation packages on R’s search path, and it registers the figure font and ggplot theme so the plots share the site’s look.

What it does not do is supply anything the code you read depends on. Each chapter loads the packages it uses with its own visible library() call, writes its file paths out in full as data/<file>, and names colors by hex code. So if you copy a block out of a chapter into a fresh R session, it runs — you are never silently relying on something defined in a file you cannot see. If a chapter’s code ever fails for you with object not found or could not find function, that is a bug in the chapter, not something you are missing.

Building the software environment

Everything except the pipeline chapter needs only R. If you are new to computing, start there and ignore conda entirely — you will still be able to work through preprocessing, QC, cell composition, batch effects, the EWAS itself, and annotation.

You need R 4.3 or newer (download from CRAN) — install the current release — and, if you want, RStudio. Everything below is typed into the Console pane — the one with the > prompt.

One note on versions: the published numbers and the Zenodo checkpoints were produced on R 4.2.3 (Bioconductor 3.16), so a value may differ in its last digits on a newer stack. The code itself is unchanged. See SESSIONINFO.md.

You need four CRAN and eighteen Bioconductor packages. The whole installation is short enough to read, so option 1 is to select the block below, copy it, and paste it into the Console. Nothing to download first, nothing to configure, and you can see exactly what it does before you run it. Press Enter and wait.

RStudio Console

repos <- "https://cloud.r-project.org"
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager", repos = repos)

BiocManager::install(c(
  "data.table", "ggplot2", "knitr", "DT",
  "minfi", "illuminaio", "wateRmelon", "GEOquery", "Biobase",
  "IlluminaHumanMethylationEPICmanifest",
  "IlluminaHumanMethylationEPICanno.ilm10b4.hg19",
  "FlowSorted.Blood.EPIC", "genefilter",
  "sva", "limma", "bacon",
  "sesame", "sesameData", "GenomicRanges", "rtracklayer",
  "missMethyl", "methylGSA"
), ask = FALSE, update = FALSE, upgrade = "never")

Expect 20-40 minutes the first time. On Windows and macOS most of these arrive as pre-built binaries — there are simply a lot of them — while on Linux they build from source, which is slower. You will see a great deal of output scroll past; that is normal.

Other ways to run the same installation

The block above is all you need. These alternatives do exactly the same thing, and exist only if you prefer them.

Option 2 — run the shipped script without downloading it. The same commands live in install_packages.R, which also checks your R version first. One line in the Console fetches and runs it:

RStudio Console

source("https://raw.githubusercontent.com/krferrier/Methylation-EWAS-tutorial/main/tutorial/install_packages.R")

Option 3 — run your local copy of that script. If you cloned the repository (see Section 3), set the working directory to tutorial/ first — Session → Set Working Directory → Choose Directory…, or open tutorial/ as an RStudio Project so it is set for you every time — then:

RStudio Console

source("install_packages.R")

Option 4 — from a terminal. The same script runs outside RStudio, from tutorial/:

Terminal

Rscript install_packages.R

The non-R sections: chapter 07 runs the EWAS through a Snakemake pipeline, and that genuinely needs conda. Chapter 06 already runs the same association test directly in R, so nothing analytical is lost — chapter 07 shows how to scale it up. The comb-p region-level section of chapter 08 also needs conda; the rest of that chapter does not.

If you already use conda, or you want to run chapter 07, four environments are defined in the repository’s envs/ directory. Three of them exist because their toolchains cannot coexist: Snakemake needs Python 3, comb-p runs on Python 2.7, and methylGSA is only packaged for an older R than the rest of the stack.

Terminal

conda env create -f envs/methyl.yml       # chapters 01-06, 08   -- R 4.5
conda env create -f envs/smk.yml          # chapter 07           -- Snakemake
conda env create -f envs/combp.yml        # chapter 08 DMRs      -- comb-p
conda env create -f envs/methylgsa.yml    # chapter 08 methylGSA -- R 4.2

Environment	Chapters	Core tool
`ewas-methyl`	01-06, 08	R 4.5, Bioconductor 3.22
`ewas-smk`	07	Snakemake 9.26.1
`ewas-combp`	08 (DMR section)	comb-p on Python 2.7
`ewas-methylgsa`	08 (methylGSA section)	methylGSA on R 4.2

If you only want the R chapters, envs/methyl.yml is the only one you need. methylGSA is one of three probe-bias corrections chapter 08 compares, and the other two — gometh and the sesame-based enrichment — run in the main environment.

Each file lists only the packages the tutorial uses directly and lets conda resolve the rest. envs/smk.yml deliberately holds nothing but the workflow engine: the pipeline declares its own per-rule R environment and Snakemake builds it on first run, which is also where METAL is compiled from source. envs/methylgsa.yml is split out for a packaging reason rather than a scientific one: bioconda’s newest methylGSA build is for R 4.2, so keeping it in envs/methyl.yml would hold that whole environment a Bioconductor generation back. Readers on the R-only route never see this — Bioconductor 3.23 ships a current methylGSA, so BiocManager installs it alongside everything else.

If you do not have conda, Miniforge is a minimal installer that defaults to the conda-forge channel; its README covers macOS, Windows, and Linux.

The exact software state that produced every number in this tutorial — full sessionInfo() with all chapters’ packages loaded, plus troubleshooting for two version-compatibility problems you may hit — is recorded in SESSIONINFO.md.

Raw Methylation data files

Each sample on an Illumina array is stored as two IDAT files — one per color channel:

GSM3853168_200932680028_R01C01_Grn.idat.gz   (green channel)
GSM3853168_200932680028_R01C01_Red.idat.gz   (red channel)

The filename encodes the physical location of the sample: 200932680028 is the Sentrix barcode (the BeadChip / “slide”), and R01C01 is the array position (row 1, column 1) on that chip.

For this tutorial, we are going to use a subset of 96 samples from the Grady Trauma Project (GSE132203), which is available from the NCBI GEO database. GEO stores each sample’s IDAT pair in its own per-sample supplementary FTP directory, e.g.:

https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3853nnn/GSM3853168/suppl/

Note the GSM3853nnn part: GEO buckets samples into directories of a thousand by replacing the last three digits of the accession with nnn. GSM3853168 therefore lives under GSM3853nnn. The download scripts below construct this for you.

The 96-sample subset is 192 files, 1.42 GB — two IDATs per sample. Do not download these by hand from the GEO website; use one of the three options below.

First, get the sample sheet

All three download options below read data/sample_sheet.csv, so you need that file before you can fetch anything. It is a 96-row table, one row per array, that pairs each GEO accession with the phenotype data for that sample. The download scripts use two of its columns: sample_id (the GSMnnnnnnn accession) and Basename (the filename stem of the array, without the _Grn.idat.gz / _Red.idat.gz suffix).

The sample sheet ships with this tutorial’s repository, so there are two ways to get it.

This gets you the sample sheet, the notebooks, and the helper scripts in one step, which is what you want if you plan to work through the tutorial yourself. Run it in a terminal — not in the R Console:

Terminal

git clone https://github.com/krferrier/Methylation-EWAS-tutorial.git

That creates a folder called Methylation-EWAS-tutorial, and the sample sheet lands at Methylation-EWAS-tutorial/tutorial/data/sample_sheet.csv.

Now point R at the tutorial folder inside it. Every path in the chapters is written relative to that folder, so this one step decides whether the rest of the tutorial works. In RStudio: Session → Set Working Directory → Choose Directory…, then pick Methylation-EWAS-tutorial/tutorial. Better still, open that folder as an RStudio Project (File → Open Project…) and it is set for you every time you come back.

Check it before moving on:

RStudio Console

getwd()                               # should end in .../Methylation-EWAS-tutorial/tutorial
file.exists("data/sample_sheet.csv")  # must print TRUE

If you only want the sample sheet, download it straight from the repository: sample_sheet.csv (in your browser: right-click → Save As).

Taking this route means you have to build the folder structure yourself, because the chapters look for data/sample_sheet.csv — the file inside a folder called data, not loose in your Downloads folder. Make a folder for the tutorial wherever you keep your work, put a data folder inside it, and save the file in there. You can do that from the Console:

RStudio Console

dir.create("~/ewas-tutorial/data", recursive = TRUE)   # make the folders
setwd("~/ewas-tutorial")                               # work from the top one
download.file(
  "https://raw.githubusercontent.com/krferrier/Methylation-EWAS-tutorial/main/tutorial/data/sample_sheet.csv",
  "data/sample_sheet.csv")
file.exists("data/sample_sheet.csv")                   # must print TRUE

Note that this gets you the sample sheet only — not the notebooks or the helper scripts.

Either way, check that you have 96 samples plus a header row before continuing.

The layout every chapter assumes

No chapter uses a path like C:/Users/you/Downloads/.... They all use short relative paths — data/sample_sheet.csv, data/idats/ — which R resolves against your working directory. So the working directory has to be the folder that contains data:

<your working directory>          <- getwd() prints this; "tutorial/" if you cloned
├── data/
│   ├── sample_sheet.csv          <- you have this now
│   └── idats/                    <- created for you by the download step below
├── 00_setup.qmd                  <- these only exist if you cloned
├── 01_qc.qmd
└── ...

Two lines confirm you are in the right place, and they are worth running before every download step:

RStudio Console

getwd()                               # where R currently is
file.exists("data/sample_sheet.csv")  # TRUE means the paths below will resolve

If you see cannot open the connection

This is the most common error in the whole tutorial, and it is worth recognising on sight:

Error in file(file, "rt") : cannot open the connection
In addition: Warning message:
In file(file, "rt") :
  cannot open file 'data/sample_sheet.csv': No such file or directory

You may see that short path or a longer one naming a folder that does not exist — either way the fix is the same. Nothing is wrong with the file and nothing is wrong with the code: R is looking in the wrong folder, because your working directory is not the one that contains data.

Step 1 — find the folder the file is actually in. Search for sample_sheet.csv (or for the repository folder, Methylation-EWAS-tutorial) in Finder on macOS or File Explorer on Windows. What you need is the folder containing data, not the file itself — if you cloned the repository, that is the tutorial folder. Then copy that folder’s path:

macOS — click the folder, then press Cmd + Option + C. (Or right-click the folder, hold Option, and choose “Copy … as Pathname”.)
Windows — right-click the folder and choose “Copy as path”.

Step 2 — point R at it. Paste the path into setwd():

RStudio Console

setwd("<paste the folder path here>")

On Windows the path you copied will use backslashes, which R reads as escape characters. Swap them for forward slashes — C:/Users/you/Methylation-EWAS-tutorial/tutorial — or double them. Or skip the typing altogether: Session → Set Working Directory → Choose Directory… in RStudio does the same thing and cannot get the slashes wrong.

Step 3 — check that it worked.

RStudio Console

getwd()                               # confirm where R thinks it is
file.exists("data/sample_sheet.csv")  # should print TRUE now

Once that last line prints TRUE, re-run whatever failed and it will work. For a fuller walkthrough of working directories, see section 6 of Jacinta Kong’s Setup R page.

The Basename column is a relative path

Basename reads idats/GSM3853168_200932680028_R01C01 — a path relative to data/, not an absolute one. That is why the download options below all write into data/idats/, and why read.metharray() further down composes the path as file.path("data/idats", basename(ss$Basename)). If you put the IDATs somewhere else, adjust that path rather than editing the sample sheet.

Downloading the 96-sample subset

Which option should I pick?

Start with option 1: it is R code you paste into the Console, it needs nothing installed beyond the packages above, and it behaves the same on macOS, Windows, and Linux. Option 2 pulls one archive from Zenodo instead of 192 separate files, so it is faster — but it needs either a terminal or a manual download in your browser. Option 3 fetches the same GEO files as option 1, from a terminal instead of the Console.

Option 1 — from GEO, using R (paste this into the Console)

This reads the sample sheet, builds each GEO URL, and downloads all 192 files. It is safe to re-run: files already present are skipped, so an interrupted download resumes where it stopped.

RStudio Console

ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE)
dir.create("data/idats", recursive = TRUE, showWarnings = FALSE)

# R's default download timeout is 60 seconds, which is not enough for a 7 MB
# file on a slow connection. Raise it or downloads will fail part-way through.
options(timeout = max(600, getOption("timeout")))

for (k in seq_len(nrow(ss))) {
  gsm    <- ss$sample_id[k]                       # e.g. GSM3853168
  stub   <- basename(ss$Basename[k])              # GSM3853168_200932680028_R01C01
  bucket <- paste0(substr(gsm, 1, nchar(gsm) - 3), "nnn")   # GSM3853nnn

  for (channel in c("Grn", "Red")) {
    fname <- sprintf("%s_%s.idat.gz", stub, channel)
    dest  <- file.path("data/idats", fname)
    if (file.exists(dest)) next                   # already have it
    url <- sprintf("https://ftp.ncbi.nlm.nih.gov/geo/samples/%s/%s/suppl/%s",
                   bucket, gsm, fname)
    # mode = "wb" is required on Windows: without it R corrupts binary files.
    download.file(url, dest, mode = "wb", quiet = TRUE)
  }
  cat(sprintf("[%3d/%d] %s\n", k, nrow(ss), gsm))
}

cat("files downloaded:", length(list.files("data/idats")), "of 192\n")

Option 2 — one archive from Zenodo (one file instead of 192)

The subset is deposited as a single archive alongside the other tutorial data, so this is one download and one extraction, with no per-sample loop and no GEO path construction. You can either access and download the archive from Zenodo manually from the tutorial data record — that link always resolves to the newest version — or you can use the provided bash script get_data.sh as follows:

Terminal

./get_data.sh A_idats

If using the bash script, the 192 IDATs will be located in data/idats/ relative to wherever the script was run. So, make sure to run the script from the location you want the dataset to live. If you are downloading the archive manually, I recommend following the same filepath structure such that the samples are all located in data/idats/, which is where the sample sheet’s Basename column expects them.

Option 3 — from GEO, using the command line

curl is preinstalled on macOS. -C - resumes a partial file, so re-running is safe.

Terminal

mkdir -p data/idats

# Pull the accession (column 1) and Basename (column 11) out of the sample
# sheet, stripping the quotes the CSV uses, then fetch both channels for each.
awk -F, 'NR > 1 { gsub(/"/, "", $1); gsub(/"/, "", $11); print $1, $11 }' \
    data/sample_sheet.csv |
while read -r gsm basename_path; do
  stub="$(basename "$basename_path")"        # GSM3853168_200932680028_R01C01
  bucket="${gsm%???}nnn"                     # GSM3853nnn
  for channel in Grn Red; do
    f="${stub}_${channel}.idat.gz"
    [ -s "data/idats/$f" ] && continue       # already downloaded
    url="https://ftp.ncbi.nlm.nih.gov/geo/samples/${bucket}/${gsm}/suppl/${f}"
    if curl -C - -sfL -o "data/idats/$f" "$url"; then
      echo "ok    $f"
    else
      echo "FAIL  $f"
    fi
  done
done

ls data/idats/*.idat.gz | wc -l      # expect 192

Open PowerShell (not Command Prompt) and run:

PowerShell

New-Item -ItemType Directory -Force -Path data\idats | Out-Null
$ss = Import-Csv data\sample_sheet.csv

foreach ($row in $ss) {
    $gsm    = $row.sample_id
    $stub   = Split-Path $row.Basename -Leaf
    $bucket = $gsm.Substring(0, $gsm.Length - 3) + "nnn"

    foreach ($channel in @("Grn", "Red")) {
        $fname = "${stub}_${channel}.idat.gz"
        $dest  = "data\idats\$fname"
        if (Test-Path $dest) { continue }
        $url = "https://ftp.ncbi.nlm.nih.gov/geo/samples/$bucket/$gsm/suppl/$fname"
        Invoke-WebRequest -Uri $url -OutFile $dest
    }
    Write-Host "$($gsm) done"
}

(Get-ChildItem data\idats\*.idat.gz).Count    # expect 192

You do not need to unzip the IDATs

minfi reads gzipped IDATs (.idat.gz) directly. Leave them compressed.

Reading the IDATs into R

The sample sheet’s Basename column holds the path stem of each array — everything except the _Grn.idat.gz / _Red.idat.gz suffix. read.metharray() takes those stems and pairs the two channels for you:

Code

RStudio Console

library(minfi)

# The sample sheet does double duty: it lists the files to read, and it holds
# the phenotype row for each array that we attach to the object below.
ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE)

# Type the variables we intend to model, now, once. Level order matters: R uses
# the FIRST level as the reference category, so listing "Control" first is what
# makes chapter 06 report the effect of PTSD relative to controls. Left to
# itself, factor() sorts alphabetically, puts "Case" first, and silently
# reverses the sign of every effect estimate downstream.
ss$ptsd            <- factor(ss$ptsd,            levels = c("Control", "Case"))
ss$sex             <- factor(ss$sex,             levels = c("F", "M"))
ss$childhood_abuse <- factor(ss$childhood_abuse, levels = c("neg", "pos"))

# Read the 192 IDATs (two channels x 96 arrays) into one RGChannelSet.
RGset <- read.metharray(
  basenames = file.path("data/idats", basename(ss$Basename)),
  force     = TRUE          # see the callout below -- required for this dataset
)

# read.metharray() reads intensities and nothing else -- it does not carry the
# sample sheet along with it. Attach the phenotype table now so that every later
# chapter can recover it with pData(RGset) rather than re-reading the CSV and
# hoping the row order still matches. It does match here, because the files were
# read in sample-sheet order.
colnames(RGset) <- ss$sample_id                       # label arrays by GEO accession
pData(RGset)    <- DataFrame(ss, row.names = ss$sample_id)   # DataFrame() comes from
                                                             # S4Vectors, attached by minfi

RGset

class: RGChannelSet 
dim: 1051943 96 
metadata(0):
assays(2): Green Red
rownames(1051943): 1600101 1600111 ... 99810990 99810992
rowData names(0):
colnames(96): GSM3853168 GSM3853169 ... GSM3853880 GSM3853881
colData names(15): sample_id series ... childhood_abuse sex
Annotation
  array: IlluminaHumanMethylationEPIC
  annotation: ilm10b4.hg19

Code

dim(RGset)                  # 1051943 addresses x 96 arrays

[1] 1051943      96

Code

# Save it. Chapters 01 and 02 both start from this object, and re-reading 192
# IDATs costs about two and a half minutes every time. The file is roughly
# 460 MB, so keep it on a disk you have room on.
saveRDS(RGset, "data/01_RGset.rds")

The result is an RGChannelSet: raw red and green intensities per address, with no normalization and no probe filtering yet. Saved as data/01_RGset.rds, it is the input to Quality control and to normalization — both chapters open by reading that exact file. It is also what the B_qc tier gives you if you skip this step, so if reading the IDATs is impractical on your machine you can download the object instead and carry on from chapter 01 unchanged.

force = TRUE is required here, and it is not a formality

Without it this dataset fails with:

[read.metharray] Trying to parse IDAT files with different array size
but seemingly all of the same type.

Ten of the twelve slides in this series were scanned with 1,052,641 addresses; the remaining two (201220980045 and 201228780138, 16 samples) carry 1,051,943. read.metharray() refuses to parse the different array sizes unless you set force = TRUE, which restricts the object to the addresses common to every array. This is why dim(RGset) reports the smaller number.

This is a scanner-decoding difference, not a data-quality problem, and taking the intersection is the correct way to handle the situation. A probe missing from some arrays cannot be compared across all of them anyway. But do not automatically use force = TRUE on other data — if the array types genuinely differ (EPIC v1 mixed with v2, say), the same error is telling you something you should not override. Check what the sizes actually are before you force it:

Code

RStudio Console

library(illuminaio)
sizes <- sapply(file.path("data/idats", paste0(basename(ss$Basename), "_Grn.idat.gz")),
                function(f) nrow(readIDAT(f)$Quants))
table(sizes)

sizes
1051943 1052641 
     16      80

Computational cost: time, memory, and disk

Everything below was measured on the machine that produced this tutorial — 32 logical cores (24 physical), 126 GB RAM, local NVMe — running the final filtered 87-sample cohort and 756,251 CpGs. Your times will differ, but this should provide a general estimate to guide what sample size you can analyse with your computing resources. The expensive steps are cell-type deconvolution, normalization, and SVA; The EWAS itself actually uses less resources depending on how it’s run.

Benchmarks per chapter

Each computationally intensive step of the tutorial was benchmarked for runtime and max memory needed. The table below shows one row per step. Memory is the peak total across every R process, which is what has to fit in RAM.

Chapter	Step	Wall time	Peak memory
01 QC	Read 96 IDAT pairs (`read.metharray`)	2 min 21 s	5.8 GB
01 QC	Detection p-values (`detectionP`)	48 s	3.8 GB
02 Normalization	Functional normalization (`preprocessFunnorm`)	2 min 4 s	12.5 GB
03 Probe filtering	Mask + detection-p + sex chromosomes	26 s	5.6 GB
04 Cell composition	Deconvolution (`estimateCellCounts2`, IDOL)	1 min 54 s	15.0 GB
05 Batch effects	Stratified `ComBat` on slide	57 s	4.1 GB
05 Batch effects	`num.sv` + `sva` at k = 6	5 min 1 s	8.9 GB
06 EWAS	EWAS using `lmFit` + `eBayes` + BACON	18 s	4.1 GB
Chapter 07	Non-stratified EWAS on 87 samples, 8 workers	6 min 12 s	12.0 GB
Chapter 07	Stratified EWAS: 2 strata (45 and 42) run simultaneously, 4 workers each	11 min 37 s	6.3 GB
08 Annotation	GENCODE v41 join over 756,251 CpGs	11 s	1.3 GB

On average, the whole tutorial should take you about 30 minutes of compute. The potential limiting resource is memory. The normalization and cell-type deconvolution takes about 12.5GB and 15GB of memory to run, respectively. Those using a laptop with 8GB of RAM may not be able to perform these steps on their local computer. Saved checkpoints of these results are available in the Zenodo record so that you can still run the rest of the steps in these chapters.

For the Chapter 07 EWAS run using a snakemake pipeline, the first row is the whole cohort in one EWAS. The second splits it into two — here by sex, giving 45 and 42 samples — and runs those two in parallel. Note that the split run needs less memory, not more: peak memory is driven by samples × workers, and each stratum is running half the samples with half the workers. The amount of memory required for the EWAS using snakemake can be reduced by reducing the number of workers. The trade-off for reducing the number of workers is that the wall time will increase.

Time versus Memory tradeoffs

You may have noticed that the EWAS method shown in Chapter 6 runs significantly faster than the EWAS performed using snakemake in Chapter 07, and at this sample size it also uses less memory — 4.1 GB against 12.0 GB for the same 87 samples at 8 workers. Those two numbers are not fixed in the same way, though. The pipeline’s memory is a function of how many workers you give it, so it is a dial you control; the Chapter 06 approach holds the whole matrix in one process, so its memory is set by the data and grows with the cohort. 87 samples is quite small, and a full-scale EWAS will have several hundreds to thousands of samples. The faster method shown in Chapter 06 will very quickly require more memory resources than you have available as the number of samples you try to run at once increase. The Snakemake pipleine in Chapter 07 takes a bit longer, especially for smaller datasets, because part of what it’s doing under the hood is breaking the big dataset into smaller, more reasonably sized chunks. The Snakemake EWAS method was designed to minimize memory cost for full-scale analyses; the consequence of this is longer runtimes.

Disk space can also be a constraint. The input files require the largest amount of space, but some of the intermediate files produced can also be large. For the 87 sample subset, I would make sure you have at least 6 GB of free disk space.

Checking what resources your own machine has

Before you plan a run, find out what you are working with.

Cores

We can use R to determine the number of cores available on your machine. This uses base R and should work for Windows, Mac, or Linux OS.

Code

RStudio Console

# --- cores -----------------------------------------------------------------
parallel::detectCores()                      # logical cores, incl. hyperthreading

[1] 32

Code

parallel::detectCores(logical = FALSE)       # physical cores

[1] 32

detectCores() reports logical cores. With hyperthreading that is often twice the physical count, but do not count on the factor of two: Intel’s hybrid designs (12th generation and later) pair hyperthreaded performance cores with single-threaded efficiency cores, so the machine these benchmarks were run on — a Core i9-14900K — reports 32 logical against 24 physical. Set workers from the physical number.

If you would rather not open R, or want to confirm what R is reporting:

Terminal (macOS)

sysctl -n hw.logicalcpu      # logical cores, incl. hyperthreading
sysctl -n hw.physicalcpu     # physical cores -- use this for `workers`

Terminal (Linux)

nproc     # logical cores, incl. hyperthreading

# physical cores = "Core(s) per socket" x "Socket(s)"
lscpu | awk -F: '/Core\(s\) per socket/{c=$2} /Socket\(s\)/{s=$2} END{print c*s}'

Recent lscpu versions indent their output under group headings, so anchoring a pattern to the start of the line (grep '^Core(s) per socket') silently matches nothing. The awk above does not anchor, and works either way.

Press Ctrl + Shift + Esc to open Task Manager.
Click on the Performance tab (the speedometer icon on the left or top).
Select CPU in the left sidebar.
Read Cores and Logical processors in the lower-right block of stats. Set workers from Cores.

RAM

You can check how much RAM is available to use with Task Manager on Windows, Activity Monitor on Mac, and System Monitor on Linux.

Press Cmd + Space to open Spotlight, type Activity Monitor, and press Enter.
Click on the Memory tab near the top of the window.
Look at the bottom of the window at Memory Pressure and Physical Memory usage graphs to see how much RAM is free or used.

From a terminal:

Terminal (macOS)

sysctl -n hw.memsize | awk '{ printf "%.1f GB installed\n", $1/1024/1024/1024 }'
vm_stat                      # live breakdown of free / active / wired pages

Open your applications menu and search for System Monitor (or KSysGuard / Plasma System Monitor on KDE desktops).
Click on the Resources tab (or Performance tab depending on your distribution).
Look at the Memory section to see a live graph and numbers showing used and free RAM.

From a terminal:

Terminal (Linux)

free -h                      # 'total' is installed, 'available' is what you can use

Press Ctrl + Shift + Esc to open Task Manager.
Click on the Performance tab (the speedometer icon on the left or top).
Select Memory in the left sidebar.
View your In use, Available, and Free RAM amounts at the bottom of the screen.

From PowerShell:

PowerShell

Get-CimInstance Win32_ComputerSystem |
  Select-Object @{n='InstalledGB';e={[math]::Round($_.TotalPhysicalMemory/1GB,1)}}

“Available” is the number that matters

Installed RAM is the headline figure, but the operating system and everything else you have open are already using some of it. Compare the peak-memory column in the table above against available memory, not installed.

Disk Space

You can check free disk space using a graphical user interface (GUI) by opening your system’s built-in file manager or storage settings panel on Windows, Mac, and Linux.

Click the Apple logo in the top-left corner of your screen.
Select System Settings (or System Preferences on older versions).
Click General in the sidebar, then select Storage.
View a colorful breakdown of your used and available disk space.

From a terminal, checking the disk holding the current folder:

Terminal (macOS)

df -h .                      # look at the 'Avail' column

Open your default file manager (such as Nautilus in GNOME or Dolphin in KDE).
Click on Other Locations or look at the left sidebar under Devices or Hard Disk to see available storage.

From a terminal, checking the disk holding the current folder:

Terminal (Linux)

df -h .                      # look at the 'Avail' column

Open File Explorer (press Windows Key + E).
Click on This PC in the left sidebar.
Look at the Devices and drives section to see a visual bar and remaining free gigabytes for your drives (like the C: drive).
Alternatively, open the Settings app, go to System, and select Storage for a detailed breakdown.

From PowerShell:

PowerShell

Get-PSDrive -PSProvider FileSystem |
  Select-Object Name, @{n='FreeGB';e={[math]::Round($_.Free/1GB,1)}}

Optional: the full 795-sample dataset

If you want to try a full-scale model after finishing the tutorial, GEO also serves the entire series as one tar archive. It is 11.9 GB compressed and expands to 1,590 IDAT files, so budget disk space and time accordingly — and expect the EWAS itself to take longer than the 87-sample walkthrough.

Terminal

mkdir -p data/idats_full && cd data/idats_full
curl -C - -O "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE132nnn/GSE132203/suppl/GSE132203_RAW.tar"
tar -xvf GSE132203_RAW.tar
ls *.idat.gz | wc -l      # expect 1590

PowerShell

New-Item -ItemType Directory -Force -Path data\idats_full | Out-Null
Set-Location data\idats_full
$u = "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE132nnn/GSE132203/suppl/GSE132203_RAW.tar"
Invoke-WebRequest -Uri $u -OutFile GSE132203_RAW.tar
tar -xvf GSE132203_RAW.tar      # tar ships with Windows 10 and later
(Get-ChildItem *.idat.gz).Count

You can also manually download the full sample set by selecting the link for the GSE132203_RAW.tar from the Supplementary Files section at the bottom of the GEO Accession page for the Grady Trauma Project.

The archive also contains the platform manifest and a filelist.txt inventory.

Building a sample sheet for the full series

The 96-sample sheet was built for you, but for the full series you need to make your own. The phenotype data lives in the series matrix file, which GEOquery downloads and parses for you — you do not need to touch the tar archive to get it.

Two details make this less tidy than it sounds, and both are shown in the code below. First, GEO packs the phenotype into columns named characteristics_ch1, characteristics_ch1.1, and so on, each holding a "label: value" string rather than a typed value; you have to split them and recover the label. Second, the Basename is not a column at all — it has to be recovered from the supplementary-file URL.

Code

RStudio Console

library(GEOquery)

# The series matrix is ~30 MB; the default 60-second timeout is not enough.
options(timeout = max(1200, getOption("timeout")))

# getGPL = FALSE skips the platform manifest, which we do not need here
# and which is large.
gse <- getGEO("GSE132203", GSEMatrix = TRUE, getGPL = FALSE)[[1]]
pd  <- Biobase::pData(gse)
nrow(pd)                     # 795 samples

[1] 795

Code

# Every "characteristics_ch1*" column holds a "label: value" string. Split on the
# first colon, using the label GEO supplies as the column name.
ch <- grep("^characteristics_ch1", names(pd), value = TRUE)
pheno <- lapply(ch, function(k) {
  v     <- as.character(pd[[k]])
  label <- sub(":.*$", "", v[1])
  value <- trimws(sub("^[^:]*:", "", v))
  setNames(data.frame(value, stringsAsFactors = FALSE), make.names(label))
})
pheno <- do.call(cbind, pheno)

# The Grn supplementary URL ends in "<GSM>_<slide>_<pos>_Grn.idat.gz". Strip that
# suffix to recover the Basename stem minfi wants.
stem <- sub("_Grn\\.idat\\.gz$", "", basename(as.character(pd$supplementary_file)))

full <- data.frame(
  sample_id  = pd$geo_accession,
  series     = "GSE132203",
  sentrix_id = pd$title,                       # e.g. 200928190033_R01C01
  slide      = sub("_.*$", "", pd$title),
  array_pos  = sub("^.*_", "", pd$title),
  Basename   = file.path("idats_full", stem),  # matches the download path above
  stringsAsFactors = FALSE
)
full <- cbind(full, pheno)

# Everything out of the series matrix is text, including the numbers.
for (k in c("age", "mergedcapsandpsswinthin30days",
            "childabphyssexemot_ctq_01modandsev")) {
  full[[k]] <- suppressWarnings(as.numeric(full[[k]]))
}

# Derive the two variables the tutorial models. PTSD is the CAPS/PSS indicator;
# 50 samples have it missing and will drop out of any model that uses it.
full$ptsd <- factor(ifelse(full$mergedcapsandpsswinthin30days == 1, "Case", "Control"),
                    levels = c("Control", "Case"))
full$sex  <- ifelse(full$gender == "Female", "F", "M")

write.csv(full, "data/sample_sheet_full.csv", row.names = FALSE)
table(full$ptsd, full$sex, useNA = "ifany")

         
            F   M
  Control 375 182
  Case    160  28
  <NA>     36  14

Running that gives 795 rows and 24 columns, with this breakdown:

	Female	Male
Control	375	182
Case	160	28
missing	36	14

Three things to notice before you model any of it. PTSD status is more imbalanced than in the teaching subset — 188 cases against 557 controls, and among males only 28 cases — where the 96-sample subset was deliberately balanced on case/control and sex. And 50 samples have no PTSD value at all, so they will be excluded from any model using that variable, exactly as 9 of the 96 do in this tutorial.

Code

RStudio Console

ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE)
i  <- match(ss$sample_id, full$sample_id)
sum(!is.na(i))                                   # 96

[1] 96

Code

identical(basename(ss$Basename), basename(full$Basename[i]))   # TRUE

[1] TRUE

Note that Basename above points at idats_full/, matching the download commands in this section. If you put the full series somewhere else, change that one line rather than editing the sheet afterwards.

The sample sheet

minfi reads IDATs given a Basename (the path up to _Grn/_Red). We pair each Basename with the phenotype row for that GSM. Here is the first 8 rows of the sample sheet:

Code

# The same CSV you downloaded above -- this is the file the whole tutorial reads.
ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE)
ss <- ss[order(ss$slide, ss$array_pos), ]   # chip, then position within the chip
cat("Samples:", nrow(ss), " | Chips:", length(unique(ss$slide)), "\n")

Samples: 96  | Chips: 12

Code

knitr::kable(head(ss[, c("sample_id","slide","array_pos","gender","age",
                         "ptsd","race")], 8),
             row.names = FALSE,
             caption = "First 8 rows of the sample sheet")

First 8 rows of the sample sheet
sample_id	slide	array_pos	gender	age	ptsd	race
GSM3853168	200932680028	R01C01	Male	25	Control	AA
GSM3853169	200932680028	R02C01	Male	59	Control	AA
GSM3853170	200932680028	R03C01	Male	43	Case	AA
GSM3853171	200932680028	R04C01	Male	21	Control	AA
GSM3853172	200932680028	R05C01	Male	45	Control	AA
GSM3853173	200932680028	R06C01	Male	19	Case	AA
GSM3853174	200932680028	R07C01	Male	47	Case
GSM3853175	200932680028	R08C01	Male	46	Control	AA

What the next notebook does

With the sample sheet provided, the QC notebook reads all 96 IDAT pairs into a single RGChannelSet and runs the quality checks that decide which samples are trustworthy before we begin normalization.

--- title: "Setup: getting the data onto your machine" --- ```{r} #| label: setup #| include: false source("_setup.R") ``` This notebook introduces the dataset the tutorial uses, shows how to get the raw files onto your own machine, and sets out what running everything costs in time, memory, and disk. It assumes the array background from [Background](00a_background.qmd). ::: {.callout-important} ## Where to type each block Every code block in this chapter carries a grey label telling you where it runs. There are only two places, and mixing them up is the most common way to get stuck: - **RStudio Console** — the pane with the `>` prompt, where R lives. Every block of R code goes here. - **Terminal** / **PowerShell** — a system shell, for the handful of non-R commands such as `git clone`. On macOS that is Terminal.app, on Windows it is PowerShell, on Linux it is whatever shell you use. One trap worth naming: RStudio has a **Terminal** tab sitting right next to the Console, and it is a system shell, not R. R code pasted there will fail with a string of syntax errors — so when a block says *RStudio Console*, make sure you have clicked the Console tab. If you are not using RStudio, `R` at a terminal gives you the same `>` prompt and works fine. ::: ::: {.callout-note collapse="true"} ## One thing that runs before each chapter, and what it is not allowed to do Every chapter begins with a hidden chunk that does `source("_setup.R")`. That file is deliberately limited to two housekeeping jobs: it puts the library holding the EPIC annotation packages on R's search path, and it registers the figure font and ggplot theme so the plots share the site's look. What it does **not** do is supply anything the code you read depends on. Each chapter loads the packages it uses with its own visible `library()` call, writes its file paths out in full as `data/<file>`, and names colors by hex code. So if you copy a block out of a chapter into a fresh R session, it runs — you are never silently relying on something defined in a file you cannot see. If a chapter's code ever fails for you with `object not found` or `could not find function`, that is a bug in the chapter, not something you are missing. ::: ## Building the software environment {#sec-env} Everything except the pipeline chapter needs **only R**. If you are new to computing, start there and ignore conda entirely — you will still be able to work through preprocessing, QC, cell composition, batch effects, the EWAS itself, and annotation. ::: {.panel-tabset} ### R only (recommended if you are getting started) You need **R 4.3 or newer** ([download from CRAN](https://cran.r-project.org/)) — install the current release — and, if you want, [RStudio](https://posit.co/download/rstudio-desktop/). Everything below is typed into the **Console** pane — the one with the `>` prompt. One note on versions: the published numbers and the Zenodo checkpoints were produced on R 4.2.3 (Bioconductor 3.16), so a value may differ in its last digits on a newer stack. The code itself is unchanged. See [`SESSIONINFO.md`](https://github.com/krferrier/Methylation-EWAS-tutorial/blob/main/tutorial/SESSIONINFO.md). You need four CRAN and eighteen Bioconductor packages. The whole installation is short enough to read, so **option 1 is to select the block below, copy it, and paste it into the Console.** Nothing to download first, nothing to configure, and you can see exactly what it does before you run it. Press Enter and wait. ```{.r filename="RStudio Console"} repos <- "https://cloud.r-project.org" if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager", repos = repos) BiocManager::install(c( "data.table", "ggplot2", "knitr", "DT", "minfi", "illuminaio", "wateRmelon", "GEOquery", "Biobase", "IlluminaHumanMethylationEPICmanifest", "IlluminaHumanMethylationEPICanno.ilm10b4.hg19", "FlowSorted.Blood.EPIC", "genefilter", "sva", "limma", "bacon", "sesame", "sesameData", "GenomicRanges", "rtracklayer", "missMethyl", "methylGSA" ), ask = FALSE, update = FALSE, upgrade = "never") ``` Expect 20-40 minutes the first time. On Windows and macOS most of these arrive as pre-built binaries — there are simply a lot of them — while on Linux they build from source, which is slower. You will see a great deal of output scroll past; that is normal. ::: {.callout-note collapse="true"} ## Other ways to run the same installation The block above is all you need. These alternatives do exactly the same thing, and exist only if you prefer them. **Option 2 — run the shipped script without downloading it.** The same commands live in [`install_packages.R`](https://github.com/krferrier/Methylation-EWAS-tutorial/blob/main/tutorial/install_packages.R), which also checks your R version first. One line in the Console fetches and runs it: ```{.r filename="RStudio Console"} source("https://raw.githubusercontent.com/krferrier/Methylation-EWAS-tutorial/main/tutorial/install_packages.R") ``` **Option 3 — run your local copy of that script.** If you cloned the repository (see @sec-sample-sheet-get), set the working directory to `tutorial/` first — **Session → Set Working Directory → Choose Directory…**, or open `tutorial/` as an RStudio Project so it is set for you every time — then: ```{.r filename="RStudio Console"} source("install_packages.R") ``` **Option 4 — from a terminal.** The same script runs outside RStudio, from `tutorial/`: ```{.bash filename="Terminal"} Rscript install_packages.R ``` ::: The non-R sections: **chapter 07** runs the EWAS through a Snakemake pipeline, and that genuinely needs conda. Chapter 06 already runs the same association test directly in R, so nothing analytical is lost — chapter 07 shows how to scale it up. The `comb-p` region-level section of chapter 08 also needs conda; the rest of that chapter does not. ### conda If you already use conda, or you want to run chapter 07, four environments are defined in the repository's `envs/` directory. Three of them exist because their toolchains cannot coexist: Snakemake needs Python 3, `comb-p` runs on Python 2.7, and `methylGSA` is only packaged for an older R than the rest of the stack. ```{.bash filename="Terminal"} conda env create -f envs/methyl.yml # chapters 01-06, 08 -- R 4.5 conda env create -f envs/smk.yml # chapter 07 -- Snakemake conda env create -f envs/combp.yml # chapter 08 DMRs -- comb-p conda env create -f envs/methylgsa.yml # chapter 08 methylGSA -- R 4.2 ``` | Environment | Chapters | Core tool | |---|---|---| | `ewas-methyl` | 01-06, 08 | R 4.5, Bioconductor 3.22 | | `ewas-smk` | 07 | Snakemake 9.26.1 | | `ewas-combp` | 08 (DMR section) | comb-p on Python 2.7 | | `ewas-methylgsa` | 08 (methylGSA section) | methylGSA on R 4.2 | If you only want the R chapters, `envs/methyl.yml` is the only one you need. `methylGSA` is one of three probe-bias corrections chapter 08 compares, and the other two — `gometh` and the sesame-based enrichment — run in the main environment. Each file lists only the packages the tutorial uses directly and lets conda resolve the rest. `envs/smk.yml` deliberately holds nothing but the workflow engine: the pipeline declares its own per-rule R environment and Snakemake builds it on first run, which is also where METAL is compiled from source. `envs/methylgsa.yml` is split out for a packaging reason rather than a scientific one: bioconda's newest `methylGSA` build is for R 4.2, so keeping it in `envs/methyl.yml` would hold that whole environment a Bioconductor generation back. Readers on the R-only route never see this — Bioconductor 3.23 ships a current `methylGSA`, so `BiocManager` installs it alongside everything else. If you do not have conda, [Miniforge](https://github.com/conda-forge/miniforge) is a minimal installer that defaults to the conda-forge channel; its README covers macOS, Windows, and Linux. ::: The exact software state that produced every number in this tutorial — full `sessionInfo()` with all chapters' packages loaded, plus troubleshooting for two version-compatibility problems you may hit — is recorded in [`SESSIONINFO.md`](https://github.com/krferrier/Methylation-EWAS-tutorial/blob/main/tutorial/SESSIONINFO.md). ## Raw Methylation data files {#sec-idat-files} Each sample on an Illumina array is stored as **two IDAT files** — one per color channel: ``` GSM3853168_200932680028_R01C01_Grn.idat.gz (green channel) GSM3853168_200932680028_R01C01_Red.idat.gz (red channel) ``` The filename encodes the physical location of the sample: `200932680028` is the **Sentrix barcode** (the BeadChip / "slide"), and `R01C01` is the **array position** (row 1, column 1) on that chip. For this tutorial, we are going to use a subset of 96 samples from the Grady Trauma Project (GSE132203), which is available from the NCBI GEO database. GEO stores each sample's IDAT pair in its own per-sample supplementary FTP directory, e.g.: ``` https://ftp.ncbi.nlm.nih.gov/geo/samples/GSM3853nnn/GSM3853168/suppl/ ``` Note the `GSM3853nnn` part: GEO buckets samples into directories of a thousand by replacing the last three digits of the accession with `nnn`. `GSM3853168` therefore lives under `GSM3853nnn`. The download scripts below construct this for you. The 96-sample subset is **192 files, 1.42 GB** — two IDATs per sample. Do not download these by hand from the GEO website; use one of the three options below. ## First, get the sample sheet {#sec-sample-sheet-get} All three download options below read `data/sample_sheet.csv`, so you need that file before you can fetch anything. It is a 96-row table, one row per array, that pairs each GEO accession with the phenotype data for that sample. The download scripts use two of its columns: `sample_id` (the `GSMnnnnnnn` accession) and `Basename` (the filename stem of the array, without the `_Grn.idat.gz` / `_Red.idat.gz` suffix). The sample sheet ships with this tutorial's repository, so there are two ways to get it. ::: {.panel-tabset} #### Clone the whole repository This gets you the sample sheet, the notebooks, and the helper scripts in one step, which is what you want if you plan to work through the tutorial yourself. Run it in a **terminal** — not in the R Console: ```{.bash filename="Terminal"} git clone https://github.com/krferrier/Methylation-EWAS-tutorial.git ``` That creates a folder called `Methylation-EWAS-tutorial`, and the sample sheet lands at `Methylation-EWAS-tutorial/tutorial/data/sample_sheet.csv`. **Now point R at the `tutorial` folder inside it.** Every path in the chapters is written relative to that folder, so this one step decides whether the rest of the tutorial works. In RStudio: **Session → Set Working Directory → Choose Directory…**, then pick `Methylation-EWAS-tutorial/tutorial`. Better still, open that folder as an RStudio Project (**File → Open Project…**) and it is set for you every time you come back. Check it before moving on: ```{.r filename="RStudio Console"} getwd() # should end in .../Methylation-EWAS-tutorial/tutorial file.exists("data/sample_sheet.csv") # must print TRUE ``` #### Download just the one file If you only want the sample sheet, download it straight from the repository: [sample_sheet.csv](https://raw.githubusercontent.com/krferrier/Methylation-EWAS-tutorial/main/tutorial/data/sample_sheet.csv) (in your browser: right-click → Save As). Taking this route means **you have to build the folder structure yourself**, because the chapters look for `data/sample_sheet.csv` — the file inside a folder called `data`, not loose in your Downloads folder. Make a folder for the tutorial wherever you keep your work, put a `data` folder inside it, and save the file in there. You can do that from the Console: ```{.r filename="RStudio Console"} dir.create("~/ewas-tutorial/data", recursive = TRUE) # make the folders setwd("~/ewas-tutorial") # work from the top one download.file( "https://raw.githubusercontent.com/krferrier/Methylation-EWAS-tutorial/main/tutorial/data/sample_sheet.csv", "data/sample_sheet.csv") file.exists("data/sample_sheet.csv") # must print TRUE ``` Note that this gets you the sample sheet only — not the notebooks or the helper scripts. ::: Either way, check that you have 96 samples plus a header row before continuing. ### The layout every chapter assumes {#sec-layout} No chapter uses a path like `C:/Users/you/Downloads/...`. They all use short relative paths — `data/sample_sheet.csv`, `data/idats/` — which R resolves against your **working directory**. So the working directory has to be the folder that contains `data`: ``` <your working directory> <- getwd() prints this; "tutorial/" if you cloned ├── data/ │ ├── sample_sheet.csv <- you have this now │ └── idats/ <- created for you by the download step below ├── 00_setup.qmd <- these only exist if you cloned ├── 01_qc.qmd └── ... ``` Two lines confirm you are in the right place, and they are worth running before every download step: ```{.r filename="RStudio Console"} getwd() # where R currently is file.exists("data/sample_sheet.csv") # TRUE means the paths below will resolve ``` ::: {.callout-warning} ## If you see `cannot open the connection` This is the most common error in the whole tutorial, and it is worth recognising on sight: ``` Error in file(file, "rt") : cannot open the connection In addition: Warning message: In file(file, "rt") : cannot open file 'data/sample_sheet.csv': No such file or directory ``` You may see that short path or a longer one naming a folder that does not exist — either way the fix is the same. Nothing is wrong with the file and nothing is wrong with the code: R is looking in the wrong folder, because your working directory is not the one that contains `data`. **Step 1 — find the folder the file is actually in.** Search for `sample_sheet.csv` (or for the repository folder, `Methylation-EWAS-tutorial`) in Finder on macOS or File Explorer on Windows. What you need is the folder *containing* `data`, not the file itself — if you cloned the repository, that is the `tutorial` folder. Then copy that folder's path: - **macOS** — click the folder, then press Cmd + Option + C. (Or right-click the folder, hold Option, and choose "Copy … as Pathname".) - **Windows** — right-click the folder and choose "Copy as path". **Step 2 — point R at it.** Paste the path into `setwd()`: ```{.r filename="RStudio Console"} setwd("<paste the folder path here>") ``` On Windows the path you copied will use backslashes, which R reads as escape characters. Swap them for forward slashes — `C:/Users/you/Methylation-EWAS-tutorial/tutorial` — or double them. Or skip the typing altogether: **Session → Set Working Directory → Choose Directory…** in RStudio does the same thing and cannot get the slashes wrong. **Step 3 — check that it worked.** ```{.r filename="RStudio Console"} getwd() # confirm where R thinks it is file.exists("data/sample_sheet.csv") # should print TRUE now ``` Once that last line prints `TRUE`, re-run whatever failed and it will work. For a fuller walkthrough of working directories, see section 6 of [Jacinta Kong's Setup R page](https://jacintak.github.io/teaching/introR/#working-directories-in-r). ::: ::: {.callout-note} ## The `Basename` column is a relative path `Basename` reads `idats/GSM3853168_200932680028_R01C01` — a path relative to `data/`, not an absolute one. That is why the download options below all write into `data/idats/`, and why `read.metharray()` further down composes the path as `file.path("data/idats", basename(ss$Basename))`. If you put the IDATs somewhere else, adjust that path rather than editing the sample sheet. ::: ## Downloading the 96-sample subset {#sec-download} ::: {.callout-tip} ## Which option should I pick? Start with **option 1**: it is R code you paste into the Console, it needs nothing installed beyond the packages above, and it behaves the same on macOS, Windows, and Linux. Option 2 pulls one archive from Zenodo instead of 192 separate files, so it is faster — but it needs either a terminal or a manual download in your browser. Option 3 fetches the same GEO files as option 1, from a terminal instead of the Console. ::: ### Option 1 — from GEO, using R (paste this into the Console) This reads the sample sheet, builds each GEO URL, and downloads all 192 files. It is safe to re-run: files already present are skipped, so an interrupted download resumes where it stopped. ```{.r filename="RStudio Console"} ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE) dir.create("data/idats", recursive = TRUE, showWarnings = FALSE) # R's default download timeout is 60 seconds, which is not enough for a 7 MB # file on a slow connection. Raise it or downloads will fail part-way through. options(timeout = max(600, getOption("timeout"))) for (k in seq_len(nrow(ss))) { gsm <- ss$sample_id[k] # e.g. GSM3853168 stub <- basename(ss$Basename[k]) # GSM3853168_200932680028_R01C01 bucket <- paste0(substr(gsm, 1, nchar(gsm) - 3), "nnn") # GSM3853nnn for (channel in c("Grn", "Red")) { fname <- sprintf("%s_%s.idat.gz", stub, channel) dest <- file.path("data/idats", fname) if (file.exists(dest)) next # already have it url <- sprintf("https://ftp.ncbi.nlm.nih.gov/geo/samples/%s/%s/suppl/%s", bucket, gsm, fname) # mode = "wb" is required on Windows: without it R corrupts binary files. download.file(url, dest, mode = "wb", quiet = TRUE) } cat(sprintf("[%3d/%d] %s\n", k, nrow(ss), gsm)) } cat("files downloaded:", length(list.files("data/idats")), "of 192\n") ``` ### Option 2 — one archive from Zenodo (one file instead of 192) The subset is deposited as a single archive alongside the other tutorial data, so this is one download and one extraction, with no per-sample loop and no GEO path construction. You can either access and download the archive from Zenodo manually from the [tutorial data record](https://doi.org/10.5281/zenodo.22135215) — that link always resolves to the newest version — or you can use the provided bash script `get_data.sh` as follows: ```{.bash filename="Terminal"} ./get_data.sh A_idats ``` If using the bash script, the 192 IDATs will be located in `data/idats/` relative to wherever the script was run. So, make sure to run the script from the location you want the dataset to live. If you are downloading the archive manually, I recommend following the same filepath structure such that the samples are all located in `data/idats/`, which is where the sample sheet's `Basename` column expects them. ### Option 3 — from GEO, using the command line ::: {.panel-tabset} #### macOS / Linux `curl` is preinstalled on macOS. `-C -` resumes a partial file, so re-running is safe. ```{.bash filename="Terminal"} mkdir -p data/idats # Pull the accession (column 1) and Basename (column 11) out of the sample # sheet, stripping the quotes the CSV uses, then fetch both channels for each. awk -F, 'NR > 1 { gsub(/"/, "", $1); gsub(/"/, "", $11); print $1, $11 }' \ data/sample_sheet.csv | while read -r gsm basename_path; do stub="$(basename "$basename_path")" # GSM3853168_200932680028_R01C01 bucket="${gsm%???}nnn" # GSM3853nnn for channel in Grn Red; do f="${stub}_${channel}.idat.gz" [ -s "data/idats/$f" ] && continue # already downloaded url="https://ftp.ncbi.nlm.nih.gov/geo/samples/${bucket}/${gsm}/suppl/${f}" if curl -C - -sfL -o "data/idats/$f" "$url"; then echo "ok $f" else echo "FAIL $f" fi done done ls data/idats/*.idat.gz | wc -l # expect 192 ``` #### Windows (PowerShell) Open **PowerShell** (not Command Prompt) and run: ```{.powershell filename="PowerShell"} New-Item -ItemType Directory -Force -Path data\idats | Out-Null $ss = Import-Csv data\sample_sheet.csv foreach ($row in $ss) { $gsm = $row.sample_id $stub = Split-Path $row.Basename -Leaf $bucket = $gsm.Substring(0, $gsm.Length - 3) + "nnn" foreach ($channel in @("Grn", "Red")) { $fname = "${stub}_${channel}.idat.gz" $dest = "data\idats\$fname" if (Test-Path $dest) { continue } $url = "https://ftp.ncbi.nlm.nih.gov/geo/samples/$bucket/$gsm/suppl/$fname" Invoke-WebRequest -Uri $url -OutFile $dest } Write-Host "$($gsm) done" } (Get-ChildItem data\idats\*.idat.gz).Count # expect 192 ``` ::: ::: {.callout-note} ## You do not need to unzip the IDATs `minfi` reads gzipped IDATs (`.idat.gz`) directly. Leave them compressed. ::: ## Reading the IDATs into R {#sec-read-idats} The sample sheet's `Basename` column holds the path stem of each array — everything except the `_Grn.idat.gz` / `_Red.idat.gz` suffix. `read.metharray()` takes those stems and pairs the two channels for you: ```{r} #| label: read-idats #| filename: "RStudio Console" library(minfi) # The sample sheet does double duty: it lists the files to read, and it holds # the phenotype row for each array that we attach to the object below. ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE) # Type the variables we intend to model, now, once. Level order matters: R uses # the FIRST level as the reference category, so listing "Control" first is what # makes chapter 06 report the effect of PTSD relative to controls. Left to # itself, factor() sorts alphabetically, puts "Case" first, and silently # reverses the sign of every effect estimate downstream. ss$ptsd <- factor(ss$ptsd, levels = c("Control", "Case")) ss$sex <- factor(ss$sex, levels = c("F", "M")) ss$childhood_abuse <- factor(ss$childhood_abuse, levels = c("neg", "pos")) # Read the 192 IDATs (two channels x 96 arrays) into one RGChannelSet. RGset <- read.metharray( basenames = file.path("data/idats", basename(ss$Basename)), force = TRUE # see the callout below -- required for this dataset ) # read.metharray() reads intensities and nothing else -- it does not carry the # sample sheet along with it. Attach the phenotype table now so that every later # chapter can recover it with pData(RGset) rather than re-reading the CSV and # hoping the row order still matches. It does match here, because the files were # read in sample-sheet order. colnames(RGset) <- ss$sample_id # label arrays by GEO accession pData(RGset) <- DataFrame(ss, row.names = ss$sample_id) # DataFrame() comes from # S4Vectors, attached by minfi RGset dim(RGset) # 1051943 addresses x 96 arrays # Save it. Chapters 01 and 02 both start from this object, and re-reading 192 # IDATs costs about two and a half minutes every time. The file is roughly # 460 MB, so keep it on a disk you have room on. saveRDS(RGset, "data/01_RGset.rds") ``` The result is an `RGChannelSet`: raw red and green intensities per address, with no normalization and no probe filtering yet. Saved as `data/01_RGset.rds`, it is the input to [Quality control](01_qc.qmd) and to [normalization](02_normalization.qmd) — both chapters open by reading that exact file. It is also what the `B_qc` tier gives you if you skip this step, so if reading the IDATs is impractical on your machine you can download the object instead and carry on from chapter 01 unchanged. ::: {.callout-warning} ## `force = TRUE` is required here, and it is not a formality Without it this dataset fails with: ``` [read.metharray] Trying to parse IDAT files with different array size but seemingly all of the same type. ``` Ten of the twelve slides in this series were scanned with 1,052,641 addresses; the remaining two (`201220980045` and `201228780138`, 16 samples) carry 1,051,943. `read.metharray()` refuses to parse the different array sizes unless you set `force = TRUE`, which restricts the object to the addresses common to every array. This is why `dim(RGset)` reports the smaller number. This is a scanner-decoding difference, not a data-quality problem, and taking the intersection is the correct way to handle the situation. A probe missing from some arrays cannot be compared across all of them anyway. But do not automatically use `force = TRUE` on other data — if the array *types* genuinely differ (EPIC v1 mixed with v2, say), the same error is telling you something you should not override. Check what the sizes actually are before you force it: ```{r} #| label: idat-sizes #| filename: "RStudio Console" library(illuminaio) sizes <- sapply(file.path("data/idats", paste0(basename(ss$Basename), "_Grn.idat.gz")), function(f) nrow(readIDAT(f)$Quants)) table(sizes) ``` ::: ## Computational cost: time, memory, and disk {#sec-cost} Everything below was measured on the machine that produced this tutorial — 32 logical cores (24 physical), 126 GB RAM, local NVMe — running the final filtered 87-sample cohort and 756,251 CpGs. Your times will differ, but this should provide a general estimate to guide what sample size you can analyse with your computing resources. The expensive steps are cell-type deconvolution, normalization, and SVA; The EWAS itself actually uses less resources depending on how it's run. ### Benchmarks per chapter Each computationally intensive step of the tutorial was benchmarked for runtime and max memory needed. The table below shows one row per step. Memory is the peak **total across every R process**, which is what has to fit in RAM. | Chapter | Step | Wall time | Peak memory | |---|----------------------------------------------|---:|---:| | [01 QC](01_qc.qmd) | Read 96 IDAT pairs (`read.metharray`) | 2 min 21 s | 5.8 GB | | [01 QC](01_qc.qmd) | Detection *p*-values (`detectionP`) | 48 s | 3.8 GB | | [02 Normalization](02_normalization.qmd) | Functional normalization (`preprocessFunnorm`) | 2 min 4 s | 12.5 GB | | [03 Probe filtering](03_probe_filtering.qmd) | Mask + detection-*p* + sex chromosomes | 26 s | 5.6 GB | | [04 Cell composition](04_cell_composition.qmd) | Deconvolution (`estimateCellCounts2`, IDOL) | 1 min 54 s | **15.0 GB** | | [05 Batch effects](05_batch_effects.qmd) | Stratified `ComBat` on slide | 57 s | 4.1 GB | | [05 Batch effects](05_batch_effects.qmd) | `num.sv` + `sva` at *k* = 6 | 5 min 1 s | 8.9 GB | | [06 EWAS](06_ewas.qmd) | EWAS using `lmFit` + `eBayes` + BACON | 18 s | 4.1 GB | | [Chapter 07](07_pipeline.qmd) | Non-stratified EWAS on 87 samples, 8 workers | 6 min 12 s | **12.0 GB** | | [Chapter 07](07_pipeline.qmd) | Stratified EWAS: 2 strata (45 and 42) run simultaneously, 4 workers each | 11 min 37 s | 6.3 GB | | [08 Annotation](08_annotation.qmd) | GENCODE v41 join over 756,251 CpGs | 11 s | 1.3 GB | On average, the whole tutorial should take you about 30 minutes of compute. The potential limiting resource is memory. The normalization and cell-type deconvolution takes about 12.5GB and 15GB of memory to run, respectively. Those using a laptop with 8GB of RAM may not be able to perform these steps on their local computer. Saved checkpoints of these results are available in the Zenodo record so that you can still run the rest of the steps in these chapters. For the Chapter 07 EWAS run using a snakemake pipeline, the first row is the whole cohort in one EWAS. The second splits it into two — here by sex, giving 45 and 42 samples — and runs those two in parallel. Note that the split run needs *less* memory, not more: peak memory is driven by **samples × workers**, and each stratum is running half the samples with half the workers. The amount of memory required for the EWAS using snakemake can be reduced by reducing the number of workers. The trade-off for reducing the number of workers is that the wall time will increase. ::: {.callout-warning} ## Time versus Memory tradeoffs You may have noticed that the EWAS method shown in Chapter 6 runs significantly faster than the EWAS performed using snakemake in Chapter 07, and at this sample size it also uses less memory — 4.1 GB against 12.0 GB for the same 87 samples at 8 workers. Those two numbers are not fixed in the same way, though. The pipeline's memory is a function of how many workers you give it, so it is a dial you control; the Chapter 06 approach holds the whole matrix in one process, so its memory is set by the data and grows with the cohort. 87 samples is quite small, and a full-scale EWAS will have several hundreds to thousands of samples. The faster method shown in Chapter 06 will very quickly require more memory resources than you have available as the number of samples you try to run at once increase. The Snakemake pipleine in Chapter 07 takes a bit longer, especially for smaller datasets, because part of what it's doing under the hood is breaking the big dataset into smaller, more reasonably sized chunks. The Snakemake EWAS method was designed to minimize memory cost for full-scale analyses; the consequence of this is longer runtimes. ::: Disk space can also be a constraint. The input files require the largest amount of space, but some of the intermediate files produced can also be large. For the 87 sample subset, I would make sure you have at least 6 GB of free disk space. ### Checking what resources your own machine has Before you plan a run, find out what you are working with. #### Cores We can use R to determine the number of cores available on your machine. This uses base R and should work for Windows, Mac, or Linux OS. ```{r} #| label: check-resources #| filename: "RStudio Console" # --- cores ----------------------------------------------------------------- parallel::detectCores() # logical cores, incl. hyperthreading parallel::detectCores(logical = FALSE) # physical cores ``` `detectCores()` reports *logical* cores. With hyperthreading that is often twice the physical count, but do not count on the factor of two: Intel's hybrid designs (12th generation and later) pair hyperthreaded performance cores with single-threaded efficiency cores, so the machine these benchmarks were run on — a Core i9-14900K — reports **32 logical against 24 physical**. Set `workers` from the physical number. If you would rather not open R, or want to confirm what R is reporting: ::: {.panel-tabset} ##### macOS ```{.bash filename="Terminal (macOS)"} sysctl -n hw.logicalcpu # logical cores, incl. hyperthreading sysctl -n hw.physicalcpu # physical cores -- use this for `workers` ``` ##### Linux ```{.bash filename="Terminal (Linux)"} nproc # logical cores, incl. hyperthreading # physical cores = "Core(s) per socket" x "Socket(s)" lscpu | awk -F: '/Core$s$ per socket/{c=$2} /Socket$s$/{s=$2} END{print c*s}' ``` Recent `lscpu` versions indent their output under group headings, so anchoring a pattern to the start of the line (`grep '^Core(s) per socket'`) silently matches nothing. The `awk` above does not anchor, and works either way. ##### Windows - Press Ctrl + Shift + Esc to open Task Manager. - Click on the Performance tab (the speedometer icon on the left or top). - Select CPU in the left sidebar. - Read **Cores** and **Logical processors** in the lower-right block of stats. Set `workers` from **Cores**. ::: #### RAM You can check how much RAM is available to use with Task Manager on Windows, Activity Monitor on Mac, and System Monitor on Linux. ::: {.panel-tabset} ##### macOS - Press Cmd + Space to open Spotlight, type Activity Monitor, and press Enter. - Click on the Memory tab near the top of the window. - Look at the bottom of the window at Memory Pressure and Physical Memory usage graphs to see how much RAM is free or used. From a terminal: ```{.bash filename="Terminal (macOS)"} sysctl -n hw.memsize | awk '{ printf "%.1f GB installed\n", $1/1024/1024/1024 }' vm_stat # live breakdown of free / active / wired pages ``` ##### Linux - Open your applications menu and search for System Monitor (or KSysGuard / Plasma System Monitor on KDE desktops). - Click on the Resources tab (or Performance tab depending on your distribution). - Look at the Memory section to see a live graph and numbers showing used and free RAM. From a terminal: ```{.bash filename="Terminal (Linux)"} free -h # 'total' is installed, 'available' is what you can use ``` ##### Windows - Press Ctrl + Shift + Esc to open Task Manager. - Click on the Performance tab (the speedometer icon on the left or top). - Select Memory in the left sidebar. - View your In use, Available, and Free RAM amounts at the bottom of the screen. From PowerShell: ```{.powershell filename="PowerShell"} Get-CimInstance Win32_ComputerSystem | Select-Object @{n='InstalledGB';e={[math]::Round($_.TotalPhysicalMemory/1GB,1)}} ``` ::: ::: {.callout-note} ## "Available" is the number that matters Installed RAM is the headline figure, but the operating system and everything else you have open are already using some of it. Compare the peak-memory column in the table above against **available** memory, not installed. ::: #### Disk Space You can check free disk space using a graphical user interface (GUI) by opening your system's built-in file manager or storage settings panel on Windows, Mac, and Linux. ::: {.panel-tabset} ##### macOS - Click the Apple logo in the top-left corner of your screen. - Select System Settings (or System Preferences on older versions). - Click General in the sidebar, then select Storage. - View a colorful breakdown of your used and available disk space. From a terminal, checking the disk holding the current folder: ```{.bash filename="Terminal (macOS)"} df -h . # look at the 'Avail' column ``` ##### Linux - Open your default file manager (such as Nautilus in GNOME or Dolphin in KDE). - Click on Other Locations or look at the left sidebar under Devices or Hard Disk to see available storage. From a terminal, checking the disk holding the current folder: ```{.bash filename="Terminal (Linux)"} df -h . # look at the 'Avail' column ``` ##### Windows - Open File Explorer (press Windows Key + E). - Click on This PC in the left sidebar. - Look at the Devices and drives section to see a visual bar and remaining free gigabytes for your drives (like the C: drive). - Alternatively, open the Settings app, go to System, and select Storage for a detailed breakdown. From PowerShell: ```{.powershell filename="PowerShell"} Get-PSDrive -PSProvider FileSystem | Select-Object Name, @{n='FreeGB';e={[math]::Round($_.Free/1GB,1)}} ``` ::: ## Optional: the full 795-sample dataset {#sec-full-series} If you want to try a full-scale model after finishing the tutorial, GEO also serves the entire series as one tar archive. It is **11.9 GB** compressed and expands to 1,590 IDAT files, so budget disk space and time accordingly — and expect the EWAS itself to take longer than the 87-sample walkthrough. ::: {.panel-tabset} #### macOS / Linux command line ```{.bash filename="Terminal"} mkdir -p data/idats_full && cd data/idats_full curl -C - -O "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE132nnn/GSE132203/suppl/GSE132203_RAW.tar" tar -xvf GSE132203_RAW.tar ls *.idat.gz | wc -l # expect 1590 ``` #### Windows (PowerShell) ```{.powershell filename="PowerShell"} New-Item -ItemType Directory -Force -Path data\idats_full | Out-Null Set-Location data\idats_full $u = "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE132nnn/GSE132203/suppl/GSE132203_RAW.tar" Invoke-WebRequest -Uri $u -OutFile GSE132203_RAW.tar tar -xvf GSE132203_RAW.tar # tar ships with Windows 10 and later (Get-ChildItem *.idat.gz).Count ``` #### Web browser You can also manually download the full sample set by selecting the link for the `GSE132203_RAW.tar` from the Supplementary Files section at the bottom of the [GEO Accession page for the Grady Trauma Project](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE132203). ::: The archive also contains the platform manifest and a `filelist.txt` inventory. ### Building a sample sheet for the full series The 96-sample sheet was built for you, but for the full series you need to make your own. The phenotype data lives in the **series matrix file**, which `GEOquery` downloads and parses for you — you do not need to touch the tar archive to get it. Two details make this less tidy than it sounds, and both are shown in the code below. First, GEO packs the phenotype into columns named `characteristics_ch1`, `characteristics_ch1.1`, and so on, each holding a `"label: value"` *string* rather than a typed value; you have to split them and recover the label. Second, the Basename is not a column at all — it has to be recovered from the supplementary-file URL. ```{r} #| label: full-sample-sheet #| filename: "RStudio Console" library(GEOquery) # The series matrix is ~30 MB; the default 60-second timeout is not enough. options(timeout = max(1200, getOption("timeout"))) # getGPL = FALSE skips the platform manifest, which we do not need here # and which is large. gse <- getGEO("GSE132203", GSEMatrix = TRUE, getGPL = FALSE)[[1]] pd <- Biobase::pData(gse) nrow(pd) # 795 samples # Every "characteristics_ch1*" column holds a "label: value" string. Split on the # first colon, using the label GEO supplies as the column name. ch <- grep("^characteristics_ch1", names(pd), value = TRUE) pheno <- lapply(ch, function(k) { v <- as.character(pd[[k]]) label <- sub(":.*$", "", v[1]) value <- trimws(sub("^[^:]*:", "", v)) setNames(data.frame(value, stringsAsFactors = FALSE), make.names(label)) }) pheno <- do.call(cbind, pheno) # The Grn supplementary URL ends in "<GSM>_<slide>_<pos>_Grn.idat.gz". Strip that # suffix to recover the Basename stem minfi wants. stem <- sub("_Grn\\.idat\\.gz$", "", basename(as.character(pd$supplementary_file))) full <- data.frame( sample_id = pd$geo_accession, series = "GSE132203", sentrix_id = pd$title, # e.g. 200928190033_R01C01 slide = sub("_.*$", "", pd$title), array_pos = sub("^.*_", "", pd$title), Basename = file.path("idats_full", stem), # matches the download path above stringsAsFactors = FALSE ) full <- cbind(full, pheno) # Everything out of the series matrix is text, including the numbers. for (k in c("age", "mergedcapsandpsswinthin30days", "childabphyssexemot_ctq_01modandsev")) { full[[k]] <- suppressWarnings(as.numeric(full[[k]])) } # Derive the two variables the tutorial models. PTSD is the CAPS/PSS indicator; # 50 samples have it missing and will drop out of any model that uses it. full$ptsd <- factor(ifelse(full$mergedcapsandpsswinthin30days == 1, "Case", "Control"), levels = c("Control", "Case")) full$sex <- ifelse(full$gender == "Female", "F", "M") write.csv(full, "data/sample_sheet_full.csv", row.names = FALSE) table(full$ptsd, full$sex, useNA = "ifany") ``` Running that gives 795 rows and 24 columns, with this breakdown: | | Female | Male | |---------|---:|---:| | Control | 375 | 182 | | Case | 160 | 28 | | missing | 36 | 14 | Three things to notice before you model any of it. **PTSD status is more imbalanced than in the teaching subset** — 188 cases against 557 controls, and among males only 28 cases — where the 96-sample subset was deliberately balanced on case/control and sex. And **50 samples have no PTSD value at all**, so they will be excluded from any model using that variable, exactly as 9 of the 96 do in this tutorial. ```{r} #| label: verify-full-sheet #| filename: "RStudio Console" ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE) i <- match(ss$sample_id, full$sample_id) sum(!is.na(i)) # 96 identical(basename(ss$Basename), basename(full$Basename[i])) # TRUE ``` Note that `Basename` above points at `idats_full/`, matching the download commands in this section. If you put the full series somewhere else, change that one line rather than editing the sheet afterwards. ## The sample sheet {#sec-sample-sheet} `minfi` reads IDATs given a **Basename** (the path up to `_Grn/_Red`). We pair each Basename with the phenotype row for that GSM. Here is the first 8 rows of the sample sheet: ```{r} #| label: sample-sheet # The same CSV you downloaded above -- this is the file the whole tutorial reads. ss <- read.csv("data/sample_sheet.csv", stringsAsFactors = FALSE) ss <- ss[order(ss$slide, ss$array_pos), ] # chip, then position within the chip cat("Samples:", nrow(ss), " | Chips:", length(unique(ss$slide)), "\n") knitr::kable(head(ss[, c("sample_id","slide","array_pos","gender","age", "ptsd","race")], 8), row.names = FALSE, caption = "First 8 rows of the sample sheet") ``` ## What the next notebook does {#sec-next} With the sample sheet provided, the [QC notebook](01_qc.qmd) reads all 96 IDAT pairs into a single `RGChannelSet` and runs the quality checks that decide which samples are trustworthy before we begin normalization.