R/noncoding_regions.R
define_regions_noncoding.RdBuilds the regions of one gene-centric non-coding mask on one chromosome:
one region per gene, whose candidate variants are the records that the
mask's annotation links to that gene, wherever they lie on the chromosome.
This is how STAARpipeline's gene-centric non-coding analysis forms its
sets, computed here once per chromosome and mask instead of once per gene.
The regions table carries the candidates as a list column
variant_rows that extract_variant_set takes in place
of its interval lookup; everything downstream (the filter's QC, type, MAF,
MAC and minimum-size tests, LD pruning, collapsing, the tests) is unchanged.
define_regions_noncoding(
gds,
chr,
category,
filter_spec,
variant_index = NULL,
context = NULL,
rule = c("glow", "staar"),
gene_list = NULL,
Annotation_dir = "annotation/info/FunctionalAnnotation",
Annotation_name_catalog = NULL,
verbose = 0
)An open SeqVarGDSClass of ONE chromosome, or a path to
such a file. A file holding several chromosomes stops with a message.
The chromosome, as the GDS names it (for example 22).
One of the eight mask names.
The run's variant_filter: its QC label,
pass value and variant type are applied here. Its MAF, MAC and
min_variants are NOT applied here; they belong to the per-set
extraction, as for every other region type. Its annotation clauses are
ignored.
A glow_variant_index of the GDS from
build_variant_index, or NULL (built here).
A context from an earlier call on the same file (the
attribute context of a returned table), or NULL (built here). A
scan builds one per chromosome and passes it to every mask, so the
chromosome-wide vectors are read once.
"glow" (default) or "staar".
NULL (the bundled table of the mask's family, on
chr) or a character vector of gene names to restrict the links to.
Character path prefix for annotations in the GDS.
Data frame mapping annotation names to GDS sub-paths. NULL uses the bundled catalog.
Integer. 1 prints the selection counts.
A data.frame with one row per gene with at least one linked record,
in gene-table order: region_id and label (the gene),
chr, start and end (the first and last position of
the linked records, a descriptive span that is wide for an enhancer far
from its gene), n_records (the linked records), category
(the mask), gene_rank (the gene's row in its table) and the list
column variant_rows (the linked records' row numbers in the
chromosome's variant index, unique and increasing). Attributes:
context (reusable for the next mask), n_selected_records
(the records the mask selected before linkage) and
n_linked_records (the (gene, record) links kept).
Every mask starts from the records that pass the filter's QC label and
variant type. The masks read five annotation fields of the aGDS, named by
their catalog names: GENCODE.Category, GENCODE.Info,
GeneHancer, CAGE and DHS; annotate_favor
writes all five by default.
| Mask | Gene table | Selection | Linkage field | TSS window |
upstream | protein-coding | category has upstream | GENCODE.Info | no |
downstream | protein-coding | category has downstream | GENCODE.Info | no |
UTR | protein-coding | category has UTR3 or UTR5 | GENCODE.Info | no |
promoter_CAGE | protein-coding | CAGE non-empty | GENCODE.Info | yes |
promoter_DHS | protein-coding | DHS non-empty | GENCODE.Info | yes |
enhancer_CAGE | protein-coding | CAGE and GeneHancer non-empty | GeneHancer | no |
enhancer_DHS | protein-coding | DHS and GeneHancer non-empty | GeneHancer | no |
ncRNA | ncRNA genes | category has ncRNA_exonic or ncRNA_splicing | GENCODE.Info | no |
"Has" means that the ";"-separated token list of the category contains the
value. The protein-coding table is the bundled genes_info, the ncRNA
table the bundled ncRNA_gene (21,104 names). The two promoter masks
keep, among the selected records, those inside the TSS \(\pm\) 3 kb
windows of GenomicFeatures::promoters(genes(TxDb.Hsapiens.UCSC.hg38.knownGene),
3000, 3000), genome-wide, as STAARpipeline does: the window is not tied to
the linked gene, so a variant can lie in gene A's window and be linked to
gene B by its annotation. The packages TxDb.Hsapiens.UCSC.hg38.knownGene,
GenomicFeatures, GenomicRanges and IRanges are needed
for those two masks only; the aGDS must be on hg38, which the code does not
verify.
Under the default rule = "glow" a record is linked by the meaning of
its annotation field, so that FAVOR v1 strings and the FAVOR 2.0 layout
written by annotate_favor give the same links wherever the
information is present in both. Every balanced parenthesized detail is
removed first; ";" separates the gene groups of a multi-category record
(one group per category token, in order) and "," the genes of a group;
NONE is a sentinel; hyphenated names stay whole.
upstream, downstream: the genes of the group at the
position of the mask's category token, so an upstream;downstream
record joins both masks with its own group in each.
UTR: all genes of the UTR groups (a record annotated to two
genes' UTRs joins both).
promoter_CAGE, promoter_DHS: the first gene of the
first group; a NONE first token links to no gene.
enhancer_CAGE, enhancer_DHS: the first listed
connected gene of the first GeneHancer element; a header without
targets links to no gene. Later targets are not used.
ncRNA: all genes of all groups, restricted to the ncRNA table.
In the FAVOR 2.0 layout the group boundaries of a two-category record are
not serialized, and such a record is then linked to every gene it names
under each of its masks. The link is then restricted to the mask's gene
table, without fallback to another gene. rule = "staar" reproduces
STAARpipeline 0.9.8's literal string operations mask by mask (the whole
category string must equal one of the values the function tests, so an
upstream;downstream record joins neither of those two masks, while
the UTR function accepts UTR5;UTR3 and the ncRNA function
ncRNA_exonic;splicing; the UTR text before the first "("; the
promoter first token after a split on "(",
",", ";" and "-"; the ncRNA first ";" group with one greedy removal of its
parenthesized text and its first three "," tokens) and exists for the
validation that enumerates the differences.
if (FALSE) { # \dontrun{
idx <- build_variant_index("chr22.gds")
spec <- variant_filter(rare_maf_cutoff = 0.01)
up <- define_regions_noncoding("chr22.gds", 22, "upstream", spec, idx)
ctx <- attr(up, "context") # reuse for the next mask
enh <- define_regions_noncoding("chr22.gds", 22, "enhancer_DHS", spec, idx,
context = ctx)
vset <- extract_variant_set("chr22.gds", up[1, ], spec, variant_index = idx)
} # }