Reference Data Standardization - FunGen-xQTL Consortium

Downloads and standardises the genome, annotation and variant reference files that the rest of the protocol reads.

Overview¶

Every study that uses this protocol should start from the same reference files: the same genome build, the same gene annotation, the same variant lists. When those differ between steps or between cohorts, the failures are quiet ones - coordinates that shift by a base, genes that exist in one annotation and not the other, meta-analyses that silently drop variants.

This notebook downloads and standardises that reference set once, so the rest of the protocol reads from a consistent source.

When to run it. Run this module once, before downstream workflows that require a reference genome, gene annotation, transcript annotation, or aligner index.

This module follows the TOPMed workflow from Broad for the aquisition and preparation of reference files for use throughout the pipeline. The processed reference data can be found in this folder on Synapse. The PDF document compiled by Data Standardization Working Group is on Synapse as well as on ADSP Dashboard. It contains the reference data to use for the project.

Notes on the reference files¶

Why the GTF is modified. Three reasons. RSEM requires the GTF to use the same chromosome naming as the FASTA (a chr prefix); for STAR this could be handled with --sjdbGTFchrPrefix "chr", but RSEM needs it in the file. The GTEx gene-model collapsing script collapse_annotation.py expects transcript_type rather than transcript_biotype, so that attribute is renamed here -- it could instead be patched in collapse_annotation.py when building the container, but since the GTF is already being customised for the chr prefix, this is folded in. ERCC records are appended at the same time. The first two could be dropped if the RSEM requirement goes away.

The collapsed gene model. Gene-level expression and eQTLs in GTEx are computed on a collapsed gene model -- all isoforms of a gene combined into one transcript -- under three rules: transcripts annotated retained_intron or read_through are excluded, and transcripts overlapping annotated read-through transcripts may be blacklisted (blacklists exist for GENCODE v19, v24 and v25; none for v26); the union of each gene’s exon intervals is taken; and intervals overlapping between genes are excluded from all of them. That last rule mainly removes genes annotated on both strands, which cannot be quantified unambiguously from unstranded RNA-seq -- for a stranded protocol it can be skipped with --collapse_only. Rules quoted from the GTEx gene_model documentation; see also the GTEx Portal.

STAR and RSEM indices. Both take the reference fasta and a gtf that must be unzipped and must be the pre-collapsing GTF. sjdbOverhang sets how much genomic sequence around each annotated junction goes into the splice-junction database; it should ideally be the read length minus one. The value 100 used here follows the TOPMed pipeline (discussion). STAR needs at least 40 GB of memory.

When to run it. First, before any data processing, and again only when you change genome build or annotation release.

Input¶

The four download_* workflows take no options -- they fetch the reference files from their public sources. Everything else consumes files those downloads produce, or that this module derives:

--hg-reference (required by hg_reference_1, hg_gtf_1, gene_annotation): the reference genome FASTA, e.g. input/reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy.fasta.
--hg-gtf (required by hg_gtf_1, gene_annotation): the gene model GTF.
--gff3-file (required by gff3_to_gtf): a GFF3 to convert to GTF.
--ercc-reference (required by hg_reference_2): the ERCC spike-in FASTA, merged into the genome.
--ercc-gtf (required by ercc_gtf): the ERCC annotation to patch and append.
--stranded / --no-stranded (required by hg_gtf_2): whether the library is stranded.
--TAD-list, --region-list, --TAD-genotype (all required by Tad_annotateion): the TAD definitions, the regions to annotate and the genotypes to annotate them against.
--genoFile (used by download_dbsnp): the dbSNP VCF to index, e.g. output/reference_data/00-All.vcf.gz.
--cwd (default output) and --modular-script-dir (default code/script).

Cluster resources: --numThreads (default 8), --job-size (default 1), --walltime (default 5h), --mem (default 16G), --container.

Reference files shipped with the protocol, under input/reference_data/:

GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy.fasta -- the genome before ERCC merge.
GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta (+ .fai) -- after merge.
Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf -- the collapsed gene model.
00-All.add_chr.variants.gz (+ .tbi) -- chr-renamed dbSNP variants.
look_up_gene_id.tsv, STAR_Index/, TAD/ -- gene-ID map, STAR index and TAD definitions.

Output¶

<genome>.fasta / .fa (+ .fai) -- from hg_reference_1..3: the reference genome after HLA/ALT and decoy removal, merged with the ERCC spike-in reference and indexed.
<genome>.dict -- the sequence dictionary Picard tools require alongside the FASTA.
<gtf>.collapse_only.gene.gtf -- from gene_annotation: the collapsed gene model, one record per gene, used wherever a gene-level annotation is needed.
ERCC92.gtf and ERCC92.genes.patched.gtf -- from ercc_gtf: the ERCC annotation, and the patched form appended to the gene model.
<name>.ref.flat -- from RefFlat_generation: the RefFlat annotation Picard uses for RNA metrics.
00-All.vcf.gz (+ .tbi) -- from download_dbsnp: the dbSNP VCF, bgzipped and tabix-indexed.
RSEM_Index/rsem_reference.grp and rsem_reference.n2g.idx.fa -- from RSEM_index.
look_up_gene_id.txt -- the gene-ID lookup table written alongside the gene annotation.
STAR_Index/ and the RSEM index -- from STAR_index and RSEM_index.
the TAD annotations -- from Tad_annotateion.

All paths are relative to --cwd.

tests/fixtures/reference_data_preparation/ERCC92.gtf:

ERCC-00002	ERCC	exon	1	1061	0.000000	+	.	gene_id "ERCC-00002"; transcript_id "DQ459430";
ERCC-00003	ERCC	exon	1	1023	0.000000	+	.	gene_id "ERCC-00003"; transcript_id "DQ516784";

output/reference_data/ERCC92.ref.flat:

	DQ459430	ERCC-00002	+	0	1061	1061	1061	1	0,	1061,
	DQ516784	ERCC-00003	+	0	1023	1023	1023	1	0,	1023,

Only the ERCC products have committed example output. output/reference_data/ holds ERCC92.gtf, ERCC92.genes.patched.gtf and ERCC92.ref.flat; the genome merge, collapsed gene model, STAR and RSEM indices and the TAD annotations were not run here, so no dump is shown for them.

Minimal Working Example¶

Preparing the reference data¶

These are not one sequential pipeline, so the blocks below are entry points rather than numbered steps. The four downloads are independent of each other and of everything else. Three do have to run in order -- format the reference genome, then process the transcript and gene model, then assemble the gene annotation -- because each reads what the one before it wrote. The remaining generators (STAR, RSEM, RefFlat, rsID and TAD) each read the formatted references and are independent of one another, so run only the ones the analysis needs.