Generation of Topologically Associated Domains and their Boundaries

Generation of Topologically Associated Domains and their Boundaries#

Merges tissue-specific TAD calls into one generalized set and extends it into per-gene fine-mapping regions.

Overview#

Analyses that work locus by locus need a definition of “locus”. Using a fixed window around each gene is simple but arbitrary: regulatory contacts do not respect a fixed distance, and two genes in the same regulatory neighbourhood get analysed as if independent. Topologically associating domain boundaries give a definition grounded in chromatin architecture instead, and this step generates the boundary files that downstream region-based analyses consume.

When to run it. Run this module before region-based association or fine-mapping when TAD-defined regions are preferred to fixed-width cis windows.

Managing TAD redundancy#

# merge the TADs to one file
cat ../../reference_data/TAD/Hippocampus_Schmitt2016-raw_TADs.txt ../../reference_data/TAD/Cortex_DLPFC_Schmitt2016-raw_TADs.txt > ../../reference_data/TAD/brain_TADs.txt

The merged TADs are processed in several passes. First, neighboring TABs are merged, including a step that manually corrects the start and end position of each chromosome. Preliminary merges then remove exact duplicates (same chromosome, start, and end). Next, the overlap between TADs is determined — the overlap with the preceding TAD, the overlap with the following TAD, and TADs entirely contained within another TAD — with two overlap values tracked per TAD (preceding and following). Finally, TADs are merged recursively: each TAD is extended to the right boundary of the last TAD and to the left boundary of the next, repeating until no further merging is possible, which yields the set of generalized TADs.

How to find boundaries and build TADB#

After the generalized TAD is built, the intervals between TADs are not simply the boundaries. Most generalized TADs have a structure like the one below.

The black line is the genome, and the red and blue shaded triangles are TAD results from hippocampus and cortex, respectively; in the code they are denoted as specific TADs because they are specific to one tissue. The two specific TAD lists (blue and red) do not have overlaps internally, but when pooled together red and blue can overlap. If their overlap exceeds 80%, they are merged into a larger TAD called a generalized TAD (green triangle), giving the generalized TAD list.

It is not necessary that specific TADs line up this way. For example, if two specific TADs do not overlap enough, or if a specific TAD does not overlap any other TAD, each forms a generalized TAD solely, as shown below.

At the boundary of a generalized TAD there are regions not covered by either of the two specific TADs, meaning not all information belongs to a TAD; these can belong to a TAB (boundary region). A generalized TAB is defined to include not only the intervals but also the ambiguous regions, shown as the yellow triangles below. To form TADBs, the generalized TAD is combined with its adjacent TABs — in the diagram, for generalized TAD2 the TADB is (generalized) TAB1 + TAD2 + TAB2.

Generalized TADs formed from a single specific TAD have no ambiguous regions, so their boundaries are also the boundaries of the generalized TAB; under this definition, adjacent generalized TADs such as TAD2 and TAD3 can share the same TADB.

The strategy is to derive a list of left and right boundaries. For a generalized TAD formed by two or more specific TADs, the left boundary is the start of the second specific TAD; for one formed from a single TAD, the left boundary is just that TAD’s start. Right boundaries are defined symmetrically using the second-to-last specific TAD’s end. Left boundaries are shown as red vertical lines and right boundaries as blue vertical lines, and the definition also covers more complicated cases.

Given these boundaries, a TADB is built by extending the left side of a generalized TAD to the last right boundary (including the start of the chromosome, 0) and the right side to the next left boundary (including the end of each chromosome).

Build TADBs and remove redundancy#

The procedure produces 1401 generalized TADs in total. Gene start and end coordinates are obtained from the gene-model GTF. Most generalized TADs are formed from 2 specific TADs, though some are formed from as many as 4; the same boundary strategy applies in every case, and generalized TADs formed from one versus more specific TADs have different boundary definitions.

Left boundaries should include the end of the chromosome, because the left boundary of a TAD is the right boundary of the former TADB; the last TADB takes the end of the chromosome as its right boundary. Each TAD is then mapped to the closest left/right boundary. After extending to TADBs the number of regions does not change, as expected, but after merging some TADBs become fully covered by other TADBs — for example, the TADB formed by generalized TAD1 can be entirely contained within the extension of generalized TAD2 (blue: right boundary, red: left boundary).

Two steps are therefore used to exclude TADBs that are fully covered by another single TADB (including the case where two TADBs are identical). After removing subsets and duplicates, 1381 TADBs remain.

Length distribution#

The average length of the generalized TADs is a little shorter than 2Mb. The generalized TADBs have an average length a little longer than 2Mb, and the number of extreme cases increases; as expected, a TADB is longer than its TAD because it also includes the boundaries.

Compare TADB and cis windows#

Each TADB covers many genes. For fine-mapping, each gene is assigned a cis window (the gene body ± 1Mb, i.e. gene start − 1Mb to gene end + 1Mb), so the next step checks whether the TADBs fully cover the cis windows of the genes inside them. Verifying that every gene is covered by a TADB confirms the TADBs are correct, since they should span the whole genome. The only genes in the full gene list not covered by a TADB are on chromosome Y, MT, or are ERCC spike-ins; these are expected to be absent because the TAD data does not include them. A value of 0 marks a cis window not fully covered by the TADB its gene lies in: even after constructing the TADBs, many genes still have a cis window (gene body ± 1Mb) that extends beyond its TADB region.

Construct TADB-enhanced cis windows and extended TADBs#

For each gene, its cis window (the gene body ± 1Mb, i.e. gene start − 1Mb to gene end + 1Mb) and the TADB it falls in are combined, extending the region outward to the smaller start and the larger end of the two. This yields one extended cis window per gene, most with a length close to or larger than the original cis window, so that every gene maps to an extended cis window and none is left uncovered. Compared with the TADB-enhanced cis window, the extended TADB allows a single gene to belong to multiple extended regions.

Why the regions are extended#

Fine mapping conventionally uses a cis window, the gene body extended 1 Mb upstream and downstream, as the region to search. That choice has no biological justification, so this module uses topologically associated domains instead, built by merging TAD calls from hippocampus and cortex brain tissue (McArthur et al. 2021, Schmitt et al. 2017).

Variants sitting within boundaries can affect gene expression too, which is why the regions are extended to the adjacent boundaries rather than stopping at the TAD edge, so that no possible causal variant is left out. The extension is deliberately conservative: the information a 1 Mb cis window provides should not be lost either, so both the TADB and the TADB-enhanced cis window are extended, each serving a different purpose.

An extended TADB starts from the generalized TADB. Every gene inside it is found, each gene’s start and end are extended by a 1 Mb cis window, and the outermost boundary of all those windows and of the TADB becomes the boundary of the extended TADB. On the reference data this yields 1,381 TADBs, which can then be used as functional units for epigenetic analysis.

A TADB-enhanced cis window starts instead from the cis window of each gene. The outermost boundary of the generalized TADB the gene sits in and of the gene’s own 1 Mb window becomes the boundary of the enhanced window; where a gene falls in two generalized TADBs, the outermost of both the TADB and the cis window is taken. The result is one cis window per gene.

When to run it. Once per genome build, as reference data preparation. The output is an input to the analyses that split work by region, rather than something rerun per study.

Input#

--tad-input: tab-separated TAD coordinates with columns chr, start, end (no header). Example: input/tadb/protocol_example.brain_TADs.txt.
--gene-coords: gene coordinate table with columns index, #chr, start, end, gene_id, gene_name (derived from the gene-model GTF). Example: input/tadb/protocol_example.gene_start_end.tsv.
--cwd: output directory (default output/tadb).
--overlap-cutoff: minimum percent overlap to merge two specific TADs into one generalized TAD (default 80).

For a genome-wide build, the raw inputs are Cortex_DLPFC_Schmitt2016-raw_TADs.txt and Hippocampus_Schmitt2016-raw_TADs.txt (hg38 cortex/hippocampus TAD regions, cf. McArthur et al 2021 and Schmitt et al 2017; available at 3dgenome.fsm.northwestern.edu) together with the gene model Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf.

Output#

{cwd}/generalized_TAD.tsv - the redundancy-managed merged TADs, chr, start, end, no header. 13 rows on the example:
```
chr22  10700000  12700000
chr22  13000000  15000000
chr22  15200000  17200000
```
{cwd}/generalized_TADB.tsv - the generalized TAD boundaries with a TADB index. Note the first boundary starts at 0, since a TADB spans the ambiguous region either side of a TAD:
```
#chr   start     end       index
chr22  0         13000000  TADB1
chr22  12650000  15300000  TADB2
```

{cwd}/TADB_enhanced_cis.bed - one TADB-enhanced cis window per gene, 1,387 rows on the example:

#chr   start  end       gene_id
chr22  0      13000000  ENSG00000226444
chr22  0      13000000  ENSG00000276871

{cwd}/extended_TADB.bed - the extended TAD boundaries, same count as the generalized TADBs but reaching further to take in the cis windows of the genes they contain:
```
#chr   start     end       index
chr22  0         13626642  TADB_1
chr22  12650000  16304556  TADB_2
```

Minimal Working Example#

One workflow, default, which builds the generalized TADs and then the boundary and cis-window products in a single pass.

Timing: TBD (on toy dataset)

sos run pipeline/generalized_TADB.ipynb default \
    --tad-input input/tadb/protocol_example.brain_TADs.txt \
    --gene-coords input/tadb/protocol_example.gene_start_end.tsv \
    --cwd output/tadb

Command Interface#

sos run pipeline/generalized_TADB.ipynb -h

usage: sos run pipeline/generalized_TADB.ipynb
               [workflow_name | -t targets] [options] [workflow_options]
  workflow_name:        Single or combined workflows defined in this script
  targets:              One or more targets to generate
  options:              Single-hyphen sos parameters (see "sos run -h" for details)
  workflow_options:     Double-hyphen workflow-specific parameters

Workflows:
  default

Global Workflow Options:
  --tad-input VAL (as path, required)
                        Brain/tissue TAD coordinates: a tab-separated file with
                        columns chr, start, end (no header)
  --gene-coords VAL (as path, required)
                        Gene coordinates table derived from the gene-model GTF
                        (columns: index, #chr, start, end, gene_id, gene_name)
  --cwd output/tadb (as path)
                        Output directory
  --modular-script-dir code/script (as path)
  --overlap-cutoff 80 (as int)
                        Minimum percent overlap for two specific TADs to be
                        merged into one generalized TAD
  --container ''
                        Container image (kept empty; we no longer use
                        containers)
  --job-size 1 (as int)
                        Resource options (kept for parity with the protocol)
  --walltime 20h
  --mem 40G
  --numThreads 1 (as int)

Sections
  default:

Workflow implementation#

[global]
# Brain/tissue TAD coordinates: a tab-separated file with columns chr, start, end (no header)
parameter: tad_input = path
# Gene coordinates table derived from the gene-model GTF (columns: index, #chr, start, end, gene_id, gene_name)
parameter: gene_coords = path
# Output directory
parameter: cwd = path("output/tadb")
parameter: modular_script_dir = path('code/script')  # override with --modular-script-dir
# Minimum percent overlap for two specific TADs to be merged into one generalized TAD
parameter: overlap_cutoff = 80
# Container image (kept empty; we no longer use containers)
parameter: container = ''
# Resource options (kept for parity with the protocol)
parameter: job_size = 1
parameter: walltime = "20h"
parameter: mem = "40G"
parameter: numThreads = 1

[default]
input: tad = tad_input, gene = gene_coords
output: gentad = f'{cwd}/generalized_TAD.tsv',
        gentadb = f'{cwd}/generalized_TADB.tsv',
        cis = f'{cwd}/TADB_enhanced_cis.bed',
        ext = f'{cwd}/extended_TADB.bed'
task: trunk_workers = 1, trunk_size = job_size, walltime = walltime, mem = mem, cores = numThreads, tags = f'{step_name}_{_output["cis"]:bn}'
bash: container=container, expand = "${ }", stderr = f'{_output["cis"]:n}.stderr', stdout = f'{_output["cis"]:n}.stdout'
    Rscript ${modular_script_dir}/reference_data/generalized_TADB.R \
        --tad-input ${_input["tad"]} \
        --gene-coords ${_input["gene"]} \
        --overlap-cutoff ${overlap_cutoff} \
        --output-gentad ${_output["gentad"]} \
        --output-gentadb ${_output["gentadb"]} \
        --output-cis ${_output["cis"]} \
        --output-ext ${_output["ext"]}