3  Data, references and public archives

An analysis needs three things: your reads, a reference genome, and a description of each sample. This chapter explains the file formats you will meet, how to choose a reference genome, and how to use public data safely.

3.1 The question

Do I have the right files, the right reference, and a correct description of each sample? Small mistakes here spread through the whole analysis. A sample labelled with the wrong condition, or reads aligned to the wrong genome version, give results that look normal but are wrong.

3.2 File formats

Table 3.1: Common file formats.

Format	What it holds	Made by
FASTQ (.fastq.gz)	Raw reads: the sequence of each read and a quality score for each base	The sequencer
BAM (.bam)	Aligned reads: where each read sits in the genome, and how sure the aligner is	The aligner (Li et al. 2009)
BED, narrowPeak, broadPeak	Regions, such as peaks: chromosome, start, end, and scores	The peak caller
bigWig (.bw)	Signal along the genome, for viewing in a browser	Coverage tools (Kent et al. 2010)
Bismark coverage (.cov)	Methylation per CpG: position, percent methylated, and counts of methylated and unmethylated reads	Bismark

FASTQ files are where the analysis starts. For paired-end sequencing, there are two files per sample: one for each end of the fragments (often called R1 and R2).

A BAM file stores each aligned read with its position and a mapping quality. The mapping quality says how sure the aligner is that the read belongs at that position (Chapter 5).

A peak file lists regions. The narrowPeak format adds a signal value, a p-value, a q-value and the position of the summit to each region.

A bigWig file stores the signal along the genome in a compact form. Genome browsers can read just the part they need, so even large files load quickly (Kent et al. 2010).

3.3 Reference genomes

Assemblies

A reference genome is the DNA sequence of a species, assembled into chromosomes. Each species has several versions, called assemblies. For example, human has hg19 (GRCh37) and hg38 (GRCh38), and mouse has mm10 (GRCm38) and mm39 (GRCm39). The same gene has different coordinates in different assemblies.

Use one assembly for the whole analysis. If you compare your results with published ones, use the same assembly as the publication, or convert the coordinates with a tool made for this.

Analysis sets

Some assemblies contain extra sequences: alternative versions of some regions, and unplaced pieces. If a region appears twice in the reference, a read from that region aligns equally well to both copies. The aligner then gives the read a low mapping quality, and the read is removed later (Chapter 5). To avoid this, analysis pipelines use an analysis set: a version of the assembly without the alternative copies.

Annotation

The annotation lists genes, transcripts and their start sites. For human and mouse, GENCODE is the most used annotation (Frankish et al. 2023). Annotations are updated often. Use the same version throughout one analysis, and report it.

Blacklists

Some regions of the genome give a strong signal in almost every experiment, whatever the target. Amemiya and colleagues describe them as regions “that have anomalous, unstructured, or high signal in next-generation sequencing experiments independent of cell line or experiment” (Amemiya et al. 2019). They are often repeats or regions near centromeres. Reads in these regions “should be identified and filtered before application of any threshold, normalization, or peak calling as they can dramatically bias the results” (Amemiya et al. 2019). A list of such regions is called a blacklist (or exclusion list).

NoteIn EpiNexus 1.2

EpiNexus supports 17 genomes. Only installed genomes are shown when you create a project.

• Human (hg38): the UCSC analysis set, with GENCODE v46 genes and the ENCODE unified exclusion list.
• Human (hg19): the UCSC analysis set, with GENCODE v46lift37 genes and the ENCODE Blacklist v2.
• Mouse (mm10, mm39): GENCODE genes (vM25, vM35), with the ENCODE Blacklist v2 (mm10) or the excluderanges list (mm39).
• Other animals: UCSC RefSeq genes. Blacklists exist for ce11 and dm6.
• Arabidopsis (tair10): Ensembl Plants, with the excluderanges list.

For genomes without a blacklist, the run tells you that blacklisted regions were not removed. Each run records the exact reference it used: the sequence file, the annotation and the blacklist, with their checksums.

3.4 Public data

Public archives hold data from thousands of published studies. You can use them to compare with your own data, to test a method, or to reanalyse a study.

• GEO (Gene Expression Omnibus) holds the description of each study and sample, and often the processed results (Barrett et al. 2013). GEO numbers start with GSE (a study, or “series”) or GSM (a sample).
• SRA (Sequence Read Archive) and ENA (European Nucleotide Archive) hold the raw reads (Katz et al. 2022; Burgin et al. 2023). They share their data with each other. Their numbers start with SRR or ERR (a sequencing run), SRX or ERX (an experiment) and PRJNA or PRJEB (a project).

One GEO sample (GSM) usually links to one experiment (SRX), which may have several runs (SRR).

Check what the data really are

Archive descriptions are written by people, and they are not always clear. Before you analyse public data, check each sample:

1. The assay. GEO and SRA list a library strategy for each sample, such as ChIP-Seq, ATAC-seq or DNase-Hypersensitivity. CUT&RUN and CUT&Tag are often listed as “ChIP-Seq” or “OTHER”.
2. The target. Find the antibody in the sample description. Inputs and IgG controls are often named only in the sample title.
3. The condition and the replicate. Read the sample characteristics, not only the title.
4. The genome. Check the species. Check if the authors’ processed files use an older assembly.

WarningPitfall: the wrong part of a study

Large studies are often split into several GEO series, one per assay, linked by a “SuperSeries”. While preparing this book, we imported GSE132292 expecting CUT&RUN data. It was the study’s DNase-seq series. The CUT&RUN data are in GSE132290. The analysis ran without error, and the results looked plausible. Only the library strategy in the archive showed the mistake. Always check the assay of each sample before you start.

NoteIn EpiNexus 1.2

On the Upload Data page, open the tab From a public repository (GEO / SRA / ENA).

• Enter up to 20 accessions: GSE, GSM, SRP, SRX, SRR, PRJNA or SAMN numbers (and their ENA and DDBJ forms). Studies with controlled access (dbGaP, EGA) cannot be imported.
• Click Look up. EpiNexus shows the samples and the files it would import. It fills in the target, the condition and the replicate from the archive text. Check every field. These are guesses from free text. Mark inputs as “input” and IgG controls as “IgG”.
• Click Import N sample(s). The files go straight from ENA to your project. Each file is checked against ENA’s MD5 checksum, so a damaged file is refused. If ENA has no FASTQ for a sample, EpiNexus downloads it from NCBI SRA and converts it.
• One import can hold up to 200 samples. Imported files count towards your storage quota of 100 GB.

When a sample’s library strategy in the archive does not match the assay of your run, the Review step of the run wizard shows a note.

Your own files

NoteIn EpiNexus 1.2

On the Upload Data page, the tab From my computer accepts the files listed in Table C.1: FASTQ files, BAM files, peak files, Bismark coverage files (.cov) and a few others. Files are sent in 50 MB pieces. If the connection breaks, the upload continues where it stopped. After the upload, EpiNexus checks that the transfer was complete and that the file looks like the format its name says. A file that passes shows Transfer verified · format check passed.

3.5 Pitfalls

WarningPitfall: mixing assemblies

Peaks called on hg19 and genes from hg38 do not match. The error is silent: genes are simply assigned to the wrong regions. Use one assembly for everything.

WarningPitfall: trusting sample titles

A title like “H3K27ac_Control_rep1” can mean an H3K27ac sample from the control condition, not a control sample. Read the sample characteristics, and set the target and condition yourself.

3.6 Further reading

• Amemiya et al. (2019): why blacklists are needed and how they were made.
• Frankish et al. (2023): GENCODE, the main gene annotation for human and mouse.
• Barrett et al. (2013), Katz et al. (2022) and Burgin et al. (2023): GEO, SRA and ENA.

Amemiya, Haley M., Anshul Kundaje, and Alan P. Boyle. 2019. ‘The ENCODE Blacklist: Identification of Problematic Regions of the Genome’. Scientific Reports 9: 9354. https://doi.org/10.1038/s41598-019-45839-z.

Barrett, Tanya, Stephen E. Wilhite, Pierre Ledoux, Carlos Evangelista, Irene F. Kim, et al. 2013. ‘NCBI GEO: Archive for Functional Genomics Data Sets—Update’. Nucleic Acids Research 41 (D1): D991–95. https://doi.org/10.1093/nar/gks1193.

Burgin, Josephine, Alisha Ahamed, Carla Cummins, Rajkumar Devraj, Khadim Gueye, et al. 2023. ‘The European Nucleotide Archive in 2022’. Nucleic Acids Research 51 (D1): D121–25. https://doi.org/10.1093/nar/gkac1051.

Frankish, Adam, Sílvia Carbonell-Sala, Mark Diekhans, Irwin Jungreis, Jane E. Loveland, et al. 2023. ‘GENCODE: Reference Annotation for the Human and Mouse Genomes in 2023’. Nucleic Acids Research 51 (D1): D942–49. https://doi.org/10.1093/nar/gkac1071.

Katz, Kenneth, Oleg Shutov, Richard Lapoint, Michael Kimelman, J. Rodney Brister, and Christopher O’Sullivan. 2022. ‘The Sequence Read Archive: A Decade More of Explosive Growth’. Nucleic Acids Research 50 (D1): D387–90. https://doi.org/10.1093/nar/gkab1053.

Kent, W. James, Ann S. Zweig, Galt Barber, Angie S. Hinrichs, and Donna Karolchik. 2010. ‘BigWig and BigBed: Enabling Browsing of Large Distributed Datasets’. Bioinformatics 26 (17): 2204–7. https://doi.org/10.1093/bioinformatics/btq351.

Li, Heng, Bob Handsaker, Alec Wysoker, Tim Fennell, Jue Ruan, et al. 2009. ‘The Sequence Alignment/Map Format and SAMtools’. Bioinformatics 25 (16): 2078–79. https://doi.org/10.1093/bioinformatics/btp352.