library(knitr)
library(reg2gene)
library(InteractionSet)
library(GenomicRanges)
library(rmarkdown)

opts_chunk$set(warning = FALSE,
               message= FALSE,
               fig.align='center',
               fig.path='Figures', 
               dev='png',
               fig.show='hold', 
               cache=FALSE)

Introduction

reg2gene R package was build to perform two main categories of tasks:

1. to associate target genes to regulatory elements genome-wide based on gene expression ~ enhancer activity models

2. to annotate user provided genomic regions to genes -

TASK 1: Building models & predicting target genes set of reg2gene functions allows users to perform:

Schematic representation of quantification and data integration:

#knitr::include_graphics("https://github.com/BIMSBbioinfo/reg2gene/blob/master/vignettes/Figures/QuantificationDataIntegrationSimplified.png")

knitr::include_graphics("/data/akalin/Projects/AAkalin_reg2gene/reg2gene/vignettes/Figures/QuantificationDataIntegrationSimplified.png")

TASK 2: Annotate user provided genomic regions to the genes

One can annotate regions of interest (ChIP-Seq peaks or similarily defined genomic regions) to the genes they are associated with, based on the result of modelling enhancer-gene associations.

Optionally, one can as well associate genes with diseases reported in the disease-gene databases.


This vignette describes necessary input data for reg2gene package and demonstrates following functionalities :


Building models with reg2gene

Quantifying gene expression using bwToGeneExp()

Before modelling, one needs to either have info about prequantified gene expression results (eg RPKMs), or he/she can quantify (and normalize) gene expression for genes of interest using reg2gene function bwToGeneExp() as follows:

# toy bigwig files
test.bw <- system.file("extdata", "test.bw",package = "reg2gene")
test2.bw <- system.file("extdata", "test2.bw",package = "reg2gene")

# toy GRanges object
regTSS_toy <- GRanges(c(rep("chr1",2),"chr2",rep("chr1",3)),
                      IRanges(c(1,7,9,15,1,15),c(4,8,14,20,4,20)),
                                            c(rep("+",3),rep("-",3)))
regTSS_toy$reg <-  regTSS_toy[c(1,1,3,5,5,5)]
regTSS_toy$name2 <- regTSS_toy$name <- paste0("TEST_Reg",
                                        c(1,1,3,5,5,5))

# run quntification f()                                        
 bwToGeneExp(exons = regTSS_toy,target = c(test.bw,test2.bw))

In the regTSS_toy toy example there are 3 genes: TEST_Reg1, TEST_Reg3 and TEST_Reg5 with two, one, and five exons. This function firstly quantifies exon expressions over pre-defined exon regions (regTSS_toy) using signal from RNA-Seq tracks (bigwig files), then it sums over all exons of a gene of interest to obtain levels of gene expression. With this function, gene expression levels can be quantified over a set of samples, cell types or conditions (list of .bigWig files).

Quantifying enhancer activity using regActivity()

Additionaly to the gene expression, one needs prequantified enhancer activity for running gene expression models. To achieve that one can use regActivity() from this package to quantify signal from ChIP-Seq, DNase-seq, WGBS coverage tracks as follows:

# toy bigwig files
test.bw <- system.file("extdata", "test.bw",package = "reg2gene")
test2.bw <- system.file("extdata", "test2.bw",package = "reg2gene")

# toy GRanges object
regTSS_toy <- GRanges(c(rep("chr1",4),rep("chr2",2)),
                      IRanges(c(1,7,9,15,1,15),c(4,8,14,20,4,20)),
                                            c(rep("+",3),rep("-",3)))
regTSS_toy$reg <-  regTSS_toy[c(1,1,3:6)]
regTSS_toy$name2 <- regTSS_toy$name <- paste0("TEST_Reg",
                                        c(1,1,3:length(regTSS_toy)))
# run quantification of enhance activity
regActivity(windows = regTSS_toy,target = c(test.bw,test2.bw))   

Preparing data prior modelling procedure using regActivityAroundTSS()

After quantification has been done, quantified enhancer activities and gene expressions need to be stored in one GRangesList object with per-gene GRanges which store info about gene expression and enhancer activities of
all enhancers within a certain range from the TSS (transcriptional start site) of a gene (by default it is +/-1Mb). This is done using regActivityAroundTSS() as follows:

# create toy GRanges objects for gene expression (regTSS_toy) and enhancer 
# activity (regReg_toy) with scores stored in bw columns

regTSS_toy <- GRReg1_toy
  regTSS_toy$bw1 <- rep(1,length(GRReg1_toy))
  regTSS_toy$bw2 <- rep(2,length(GRReg1_toy))
  regTSS_toy$bw3 <- rep(3,length(GRReg1_toy))
regReg_toy <- GRReg2_toy
   regReg_toy$bw1 <- rep(3,length(regReg_toy))
   regReg_toy$bw2 <- rep(4,length(regReg_toy))

# combine these object into per gene GRangesList

regActivityAroundTSS(regActivity = regReg_toy,geneExpression = regTSS_toy, 
                     upstream=1,downstream=1)

Individual per gene GRanges objects have the following metadata columns:

 1. featureType: either "gene" or "regulatory" 
 2. name: name/id for gene and enhancers. Gene name could 
 be id from a database enhancer name should be in the format as follows 
 "chr:start-end" 
 3. name2: a secondary name for the feature, such as gene symbol "PAX6" etc.
 not necessary for enhancers could be NA 
 4. other columns: numeric values for gene expression or regulatory actvity.
 Column names represent sample names/ids

Importantly, only enhancers located within predefined (+/-) upstream/downstream regions of TSS are identified, extracted and reported in output (together with info about gene expression).

!!! Sample id's (corresponding to the cell types or conditions) are included in output object only if both, 1) gene expression values and 2) quantified regulatory activity are available in TSS and regActivity objects. Non-overlapping cell types are excluded.

Modelling gene expression ~ enhancer activity using associateReg2Gene()

To link enhancer activity and gene expression reg2gene function utilizes correlation methods (Spearman and Pearson corr coeficients) and/or use other statistical approaches, elastic net and random forests, to model gene expression ~ enhancer activity across different cell types (This approach is based on the observation that enhancers show very high tissue specificity and the level of their activity correlates with gene expression (Visel et al. 2009., Ernst et al. 2011).

To run models one uses regActivityAroundTSS() object

###############################
#STEP 1.  Getting random and predefined .8 correlation

 require(GenomicRanges)
 require(doMC)
 require(glmnet)
 require(foreach)
 require(stringr)
 require(qvalue)

 ####################################
 # create example

 x <- c(2.000346,2.166255,0.7372374,0.9380581,2.423209, 
      2.599857,4.216959,2.589133,1.848172,3.039659)
 y <- c(2.866875,2.817145,2.1434456,2.9039771,3.819091,5.009990,
      5.048476,2.884551,2.780067,4.053136)
 corrM <- rbind(x,y)

 # define Granges object
  gr0 <- GRanges(seqnames=rep("chr1",2),IRanges(1:2,3:4))

    GeneInfo <- as.data.frame(matrix(rep(c("gene","regulatory"),each=3),
                ncol = 3,byrow = TRUE),stringsAsFactors=FALSE)

        colnames(GeneInfo) <- c("featureType","name","name2")

       mcols(gr0) <- DataFrame(cbind(GeneInfo,corrM))


       gr0

And run associateReg2Gene() as follows ( gr0):

###############################
#STEP 1.  Getting random and predefined .8 correlation

 require(GenomicRanges)
 require(doMC)
 require(glmnet)
 require(foreach)
 require(stringr)
 require(qvalue)

 ####################################
 # create example

 x <- c(2.000346,2.166255,0.7372374,0.9380581,2.423209, 
      2.599857,4.216959,2.589133,1.848172,3.039659)
 y <- c(2.866875,2.817145,2.1434456,2.9039771,3.819091,5.009990,
      5.048476,2.884551,2.780067,4.053136)
 corrM <- rbind(x,y)

 # define Granges object
  gr0 <- GRanges(seqnames=rep("chr1",2),IRanges(1:2,3:4))

    GeneInfo <- as.data.frame(matrix(rep(c("gene","regulatory"),each=3),
                ncol = 3,byrow = TRUE),stringsAsFactors=FALSE)

        colnames(GeneInfo) <- c("featureType","name","name2")

       mcols(gr0) <- DataFrame(cbind(GeneInfo,corrM))

    print("associateReg2Gene(input=gr0,cores = 1,B=100,method=\"spearman\")")
    associateReg2Gene(input = gr0,cores = 1,B=100)     

This example was created to have correlation between gene expression and enhancer activity equal to 0.8. This estimated association statistic is reported as coefs and is recalculated as 0.8. Associated P-values are calculated in and reported as pval. Corresponding qval is calculated. However, q-values are set to be NA since, this toy example contains only one gene~enhancer pair, thus one p-value is calculated but q-values cannot be calculated based on only one p-value.

Importantly, this function does not report statistically significant gene-enhancer associations, whereas it reports all input gene-enhancer associations and corresponding test statistics. It is up to researcher to decide which gene-enhancer associations they consider to be statistically significant pairs.

MODELLING ADDITIONS: meta-analysis & voting

In cases when more that one modelling procedure has been performed, one can try to combine these results. This package can perform meta-analysis and voting procedure.

Model voting by voteInteractions()

You can perform voting analysis to combine models that differ in algorithm, or method, or cohort (just as meta-analysis) used when modelling by running voteInteractions() as follows:

voteInteractions( interactions=interactions, cutoff.stat="pval", cutoff.val=0.05, vote.threshold=0.5)

require(GenomicRanges)
require(InteractionSet)

gr2 <- gr <- GRanges(seqnames=rep("chr1",3),IRanges(1:3,3:5))
   x <- 1:5
   y <- 2:6
   z <- 10:14
   a <- rep(0,length(x))


   GeneInfo <- as.data.frame(matrix(c(rep("gene",3),rep("regulatory",6)),
               ncol = 3,byrow = TRUE),stringsAsFactors=FALSE)
               colnames(GeneInfo) <- c("featureType","name","name2")

 mcols(gr) <- DataFrame(cbind(GeneInfo,rbind(x,y,z)))
 mcols(gr2) <- DataFrame(cbind(GeneInfo,rbind(x,y,a)))

 # create associateReg2Gene output objects, GInteractions will all 
 # output results

 AssocObject <- reg2gene::associateReg2Gene(gr)
 AssocObject2 <- reg2gene::associateReg2Gene(gr2)

 # input for meta-analysis is list of such objects

 interactions <- list(AssocObject,AssocObject2)
 names(interactions) <- c("H3K4me1","H327ac")

 # Run voteInteractions


 voteInteractions( interactions, 
                  cutoff.stat="pval",
                  cutoff.val=0.05,
                  vote.threshold=0.5)

This f() selects POSITIVES (statistically associated gene~enhancer pairs) for each result of associateReg2Gene() analysis that one wants to combine by majority voting, thus the result of the voting procedure is a list of gene~enhancer associations which have been confirmed by at least two or more models.

Schematic representation of majority voting:

#knitr::include_graphics("https://github.com/BIMSBbioinfo/reg2gene/master/vignettes/Figures/VOTING.png")

knitr::include_graphics("/data/akalin/Projects/AAkalin_reg2gene/reg2gene/vignettes/Figures/VOTING.png")

Additionally, one can perform meta-analysis:

Meta-analysis using metaInteractions()

When one have information coming from different cohorts, or produced by different organizations/data centers,and thus there is no overlap between samples, one can run meta-analysis to combine p-values from individual models which should improve reprodicibility of modelling results. But before running meta-analysis, please read more when and where to use it. One runs meta-analysis using metaInteractions() to combine association P-values and coefficients from different data sets using Fisher's method and weigthed averaging respectively as follows:

metaInteractions(interactions = interactions)

# creating datasets

require(GenomicRanges)
require(InteractionSet)

gr2 <- gr <- GRanges(seqnames=rep("chr1",3),IRanges(1:3,3:5))
   x <- 1:5
   y <- 2:6
   z <- 10:14
   a <- rep(0,length(x))


   GeneInfo <- as.data.frame(matrix(c(rep("gene",3),rep("regulatory",6)),
               ncol = 3,byrow = TRUE),stringsAsFactors=FALSE)
               colnames(GeneInfo) <- c("featureType","name","name2")

 mcols(gr) <- DataFrame(cbind(GeneInfo,rbind(x,y,z)))
 mcols(gr2) <- DataFrame(cbind(GeneInfo,rbind(x,y,a)))

 # create associateReg2Gene output objects, GInteractions will all 
 # output results

 AssocObject <- reg2gene::associateReg2Gene(gr)
 AssocObject2 <- reg2gene::associateReg2Gene(gr2)

 # input for meta-analysis is list of such objects

 interactions <- list(AssocObject,AssocObject2)
 names(interactions) <- c("Roadmap","Blueprint")



 # Run metaA
 metaInteractions( interactions=interactions)

Meta-analysis - simplified pictorial representation:

#knitr::include_graphics("https://github.com/BIMSBbioinfo/reg2gene/blob/master/vignettes/Figures/Meta-Analysis_Simplified.png")


knitr::include_graphics("/data/akalin/Projects/AAkalin_reg2gene/reg2gene/vignettes/Figures/Meta-Analysis_Simplified.png")

Benchmarking enhancer-gene associations using benchmarkInteractions()

BACKGROUND

Associated gene-enhacer pairs (an output of associateReg2Gene()) can be easily benchmarked using the second set of linked genes and enhancers (for example, information about the 4D genome architecture coordinates obtained using chromatin conformation capture and related methods (CHiA-PET, 4C, 5C, HiC, PC-HiC).

Schematic representation of possible benchmarking procedure

#knitr::include_graphics("https://github.com/BIMSBbioinfo/reg2gene/blob/master/vignettes/Figures/BenchSimpleE.png")

knitr::include_graphics("/data/akalin/Projects/AAkalin_reg2gene/reg2gene/vignettes/Figures/BenchSimpleE.png")

benchmarkInteractions() takes two GInteractions objects as an input, and the first object (reg2Gene) is benchmarked with respect to the second one (benchData) as follows:

benchmarkInteractions(interactions, benchInteractions, preFilter = T, binary=TRUE)

   interactions <- GInteractions(GRReg1_toy,GRReg1_toy$reg)[2]
   benchInteractions <- GInteractions(GRReg2_toy,GRReg2_toy$reg)

   # removing confusing meta-data
   mcols(interactions) <- NULL

benchmarkInteractions(interactions = interactions,
              benchInteractions = benchInteractions, 
              preFilter = T,
              binary=TRUE) 

Criss-cross overlap of interacting regions is performed: if anchor1 from the benchmark dataset is overlapping anchor2 from the tested set, than anchor2 from the benchmark dataset needs to overlap anchor1 from the tested set, or vice-versa as shown:

#knitr::include_graphics("https://github.com/BIMSBbioinfo/reg2gene/master/vignettes/Figures/benchmarkF.png")

knitr::include_graphics("/data/akalin/Projects/AAkalin_reg2gene/reg2gene/vignettes/Figures/benchmarkF.png")

Assessing modelling performance using confusionMatrix()

confusion matrix can be calculated to assess performance of the modelling procedure based on the external benchmark dataset and the resuls of benchmarkInteractions()

interactionsBench <- GInteractions(GRReg1_toy,GRReg1_toy$reg)

Bench <- interactionsBench$anchor1.Bench1Exp
Filter <- interactionsBench$anchor1.Filter1Exp
      mcols(interactionsBench) <- NULL

   interactionsBench$Pval <- seq(0, 1, length.out = length(GRReg1_toy))
   interactionsBench$Bench <- Bench
   interactionsBench$Filter <- Filter

  interactionsBench

using confusionMatrix() as follows:

confusionMatrix(interactionsBench, thresholdID = "Pval", thresholdValue = 0.05, benchCol = "Bench", prefilterCol = "Filter", statistics = "ConfusionMatrix")

interactionsBench <- GInteractions(GRReg1_toy,GRReg1_toy$reg)

Bench <- interactionsBench$anchor1.Bench1Exp
Filter <- interactionsBench$anchor1.Filter1Exp
      mcols(interactionsBench) <- NULL

   interactionsBench$Pval <- seq(0, 1, length.out = length(GRReg1_toy))
   interactionsBench$Bench <- Bench
   interactionsBench$Filter <- Filter

confusionMatrix(interactionsBench,
                thresholdID = "Pval",
                thresholdValue = 0.05,
                benchCol = "Bench",
                prefilterCol = "Filter",
                statistics = "ConfusionMatrix")

Visualizing enhancer-gene associations using plotInteractions()

For user provided enhancer~gene associations (imported as GInteractions) object plot associations as loops in the context of the genome using plotInteractions() as follows (FTO/IRX3/IRX5 region example)

plotInteractions(interactions = GenomeInteractions, statistics ="pval", coloring = "color")

enhancers <- GRanges(rep("chr16",6),
                      IRanges(c(53112601,55531601,53777201,
                                53778801,54084001,53946467),
                              c(53114200,55533250, 53778800, 
                                53780400, 54084400 ,53947933)))

genes <- GRanges(rep("chr16",6),
                     IRanges(c(53737874, 54964773, 54320676,
                               53737874, 54964773, 54320676),
                             c(53737874, 54964773, 54320676,
                               53737874, 54964773, 54320676)))

GenomeInteractions <- GInteractions(enhancers,genes)

GenomeInteractions$name2 <- c("FTO","IRX5","IRX3")

GenomeInteractions$pval <- c(0.20857403, 0.72856090, 0.03586015,
                             0.32663439, 0.32534945, 0.03994488)

GenomeInteractions$color <- c("red","blue","grey")


 plotInteractions(interactions = GenomeInteractions,
                  statistics ="pval",
                  coloring = "color")

You can choose specific gene to plot: FTO

 plotInteractions(interactions = GenomeInteractions,
                        selectGene="FTO")

Or choose specific region to plot: chr16:53112601-53114200

NOTE! This region does not necessarilly need to be equal to the regulatory regions reported in the interactions input objects, whereas it only needs to overlap some regulatory regions.

  plotInteractions(interactions = GenomeInteractions,
               selectRegulatoryRegion = "chr16:53112601-53114200")

Additionaly, benchmark datasets can be plotted as well

 benchInteractions = list(GenomeInteractions[1:3])

  plotInteractions(interactions = GenomeInteractions,
                         coloring = "color",
                         statistics = "pval",
                         benchInteractions = benchInteractions)

Annotating genomic regions to genes

Annotating regulatory regions to genes using reg2gene()

User can annotate regions of interest (ChIP-Seq peaks or any other genomic regions according to the provided geneAnnotations GRanges object and/or enhancer~gene associations object (e.g. interactions argument) using reg2gene() as follows:

reg2gene(windows, geneAnnotations=geneAnnotations , interactions =interactions, identified=T)

Creating toy example:

# creating toy example
# 1. windows of interest
   windows <- GRanges(c("chr1:1-2", # 1. overlap prom
                                 "chr2:1-2",  # 2. overlap enh
                                 "chr3:1-2", # 3. overlap tss +/- 1,000,000
                                 "chr4:1-2")) # 4. do not overlap tss +/- 1Mb

     annotationsEnh <- GRanges(c("chr1:1-2",
                                 "chr2:1-2",
                                 "chr3:100000-100002",
                                 "chr4:10000001-10000002"))

     annotationsGenes <- GRanges(c("chr1:1-2",
                                   "chr2:100000-100002",
                                   "chr3:99999-100002",
                                   "chr4:10000001-10000002"))

     seqlengths(annotationsEnh) <- seqlengths(annotationsGenes) <- rep(10000002,
                                                                       4)
  # example of interactions
         interactions = GInteractions(annotationsEnh,annotationsGenes,
                                 name=c("gen1","gen2","gen3","gen4"))

  # example of geneAnnotations   
     geneAnnotations=second(interactions)
     mcols(geneAnnotations) <- mcols(interactions) 

     print("windows of interest")
     windows
      print("toy example of interactions object")
     interactions
      print("toy example of geneAnnotations object")
     geneAnnotations

Running reg2gene annotation function with created toy example:

  # run annotation function: 
 reg2gene(windows=windows,
          interactions=interactions, 
          geneAnnotations = geneAnnotations)

 # which regions are not identified

 reg2gene(windows=windows,
          interactions=interactions, 
          geneAnnotations = geneAnnotations,
          identified=FALSE)

 # if interactions are not available, assign interactions based solely on the 
 # proximity to promoters
 reg2gene(windows = windows,
                 interactions=NULL,
                 geneAnnotations = geneAnnotations)

If interactions object is missing, function annotates the input regions to their nearby genes.

When GIneraction object - interactions - is used as an input, then hierarchical association procedure is runned as follows: promoters,enhancers, nearby genes, eg: 1) genomic regions of interest are first considered to be promoters and associated with nearby genes if they are located within a certain distance from TSS of nerbay gene (default +/-1000bp); otherwise 2) remaning genomic regions are overlapped with enhancer regions, and genes associated to that enhancer regions are reported, 3) if no overlap with either promoters nor enhancers is identified, then closest gene is reported if it is located within 1Mb 4) if no gene located within 1Mb has been identified then, this region is filtered out.

IMPORTANT! interactions can store info about enhancer~gene interactions, and in that case Anchor1 in GInteractions object needs to be regulatory region, whereas anchor2 is the location of gene/TSS.

However, it can store info about interactions regardless whether one anchor of these interactions correspond to promoters, for example HiC interactions. In that case set annotateInteractions=FALSE and interactions are firstly annotated to the genes (unannotated interactions are removed from the dataset), and windows are subsequently annotated to genes using annotated interactions object.

Schematic representation of reg2gene:

#knitr::include_graphics("https://github.com/BIMSBbioinfo/reg2gene/master/vignettes/Figures/annotatExample.png")

knitr::include_graphics("/data/akalin/Projects/AAkalin_reg2gene/reg2gene/vignettes/Figures/annotatExample.png")

In addition, if argument identified is set to be FALSE reg2gene() will report genomic regions for which corresponding genes were not identified.

 reg2gene(windows,
          geneAnnotations=geneAnnotations,
          interactions,
          identified=F)

sessionInfo

sessionInfo()

Short background of the problem

Enhancers

Enhancers represent distal regulatory regions in the genome that can be located up to 1Mb from the transcription start sites of genes, increase gene expression regardless of their position, orientation and distance to the promoter.

Enhancer-like chromatin marks

Different enhancer-like chromatin marks have been previously used to map these regulatory regions and assess their activity: level of histone modifications (especially H3K4me1 and H3K27ac), nucleosome depletion, open chromatin accessibility, DNA methylation, nucleotide conservation, etc. (Mathelier et al., 2015).
To improve reproducibility of mapping efforts, different chromatin marks were integrated in chromatin features (multiple modifications and more complex elements linked together, Stricker et al. 2017), combined with an unsupervised machine-learning approaches and higher number of cell types included in the mapping analysis (for example Ernst and Kellis 2012 used ChromHMM Core 15-state model and reported more than 900,000 potential enhancer regions across 127 epigenomes from the Roadmap Epigenomics Project (Kundaje et al. 2015).

Mapping potential gene targets of an enhancer regions

1) Mapping of regulatory regions and their targeted genes can be done computationally by correlating gene expression with levels of enhancer- associated chromatin features. Ernst et al. 2011 correlated gene expression with different histone modification marks, including enhancer associated marks H3K27ac and H3K4me1 (ENCODE Project Consortium 2012), Sheffield et al. 2013 correlated DNase I Hypersensitivity and gene expression, whereas Varley et al. 2013 was focused on DNA methylation. This package enables easy integration of all these datasets.

1) A long-range interactions can be mediated by chromatin looping - a mechanism by which enhancers and promoters are brought together in the 3D space of a nucleus, which, at least partially, enables precise regulation of gene expression (Gorkin et al. 2014, Rennie et al. 2017). Thus information about the 3D genome architecture, generated by chromosome conformation capture and related techniques, is can be used as a proxy of enhancer-mediated regulation of gene expression or in these case to benchmark results of reg2gene modelling approach (Dekker et al. 2002, Dosie et al. 2006, Simonis et al. 2006, Fullwood et al. 2009, Lieberman-Aiden et al. 2009).


Literature

Akalin, Altuna, et al. "Genomation: a toolkit to summarize, annotate and visualize genomic intervals." Bioinformatics 31.7 (2014): 1127-1129.

Bolstad, Benjamin Milo. "preprocessCore: A collection of pre-processing functions." R package version 1.0 (2013).

Dabney, Alan, John D. Storey, and G. R. Warnes. "qvalue: Q-value estimation for false discovery rate control." R package version 1.0 (2010).

Dekker, Job, et al. "Capturing chromosome conformation." science 295.5558 (2002): 1306-1311.

Dostie, Josée, et al. "Chromosome Conformation Capture Carbon Copy (5C): a massively parallel solution for mapping interactions between genomic elements. " Genome research 16.10 (2006): 1299-1309.

ENCODE Project Consortium. "An integrated encyclopedia of DNA elements in the human genome." Nature 489.7414 (2012): 57-74.

Ernst, Jason, et al. "Mapping and analysis of chromatin state dynamics in nine human cell types." Nature 473.7345 (2011): 43-49.

Ernst, Jason, and Manolis Kellis. "ChromHMM: automating chromatin-state discovery and characterization." Nature methods 9.3 (2012): 215-216.

Fullwood, Melissa J., et al. "An oestrogen-receptor-α-bound human chromatin interactome." Nature 462.7269 (2009): 58-64.

Gorkin, David U., Danny Leung, and Bing Ren. "The 3D genome in transcriptional regulation and pluripotency." Cell stem cell 14.6 (2014): 762-775.

Friedman, Jerome, Trevor Hastie, and Rob Tibshirani. "glmnet: Lasso and elastic-net regularized generalized linear models." R package version 1.4 (2009).

Kundaje, Anshul, et al. "Integrative analysis of 111 reference human epigenomes. " Nature 518.7539 (2015): 317-330.

Lieberman-Aiden, Erez, et al. "Comprehensive mapping of long-range interactions reveals folding principles of the human genome." science 326.5950 (2009): 289-293.

Lonsdale, John, et al. "The genotype-tissue expression (GTEx) project." Nature genetics 45.6 (2013): 580-585.

Love, Michael I., Wolfgang Huber, and Simon Anders. "Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2." Genome biology 15.12 (2014): 550.

Mathelier, Anthony, Wenqiang Shi, and Wyeth W. Wasserman. "Identification of altered cis-regulatory elements in human disease." Trends in Genetics 31.2 (2015): 67-76.

Sheffield, Nathan C., et al. "Patterns of regulatory activity across diverse human cell types predict tissue identity, transcription factor binding, and long-range interactions." Genome research 23.5 (2013): 777-788.

Simonis, Marieke, et al. "Nuclear organization of active and inactive chromatin domains uncovered by chromosome conformation capture–on-chip (4C)." Nature genetics 38.11 (2006): 1348-1354.

Stricker, Stefan H., Anna Köferle, and Stephan Beck. "From profiles to function in epigenomics." Nature Reviews Genetics (2016).

Varley, Katherine E., et al. "Dynamic DNA methylation across diverse human cell lines and tissues." Genome research 23.3 (2013): 555-567.

Visel, Axel, et al. "VISTA Enhancer Browser—a database of tissue-specific human enhancers." Nucleic acids research 35.suppl 1 (2007): D88-D92.

Wright, Marvin N., and Andreas Ziegler. "ranger: A fast implementation of random forests for high dimensional data in C++ and R." arXiv preprint arXiv:1508.04409 (2015).



BIMSBbioinfo/reg2gene documentation built on May 3, 2019, 6:42 p.m.