lenhancer: LASSO based Enhancer identification

Introduction

Combining a novel strategy to identify communities of related control elements with a penalized regression approach, we developed individual gene-by-gene models to identify the potential control elements predictive of the expression of the gene.

The package requires the following data:

• expression (E): gene expression data in say n cell types
• regions (R): A set of putative enhancer CREs with a given distance of the genes
• regulation signals (S): regulation signals such as active chromatin marks, DNaseI-seq or ATAC-seq in the n cell types. You can normalise the way you prefer for library sizes, and region sizes.
• tfbs (T): Transcription factor binding signals in these regions. You can use as many TF binding data as you want, but TF ChIP-seq in the n cell types is highly recommended.

For each gene, the method consists of the following steps:

• m = consider peaks that are with 100kB (can be changed)
• Using correlation of T and S between these m regions, a community clustering algorithm is used to group the CREs into clusters called coCREs.
• singleton CREs within 10Kb and coCREs are then considered as predictors, and the gene's expression across cell types as response variable.
• Penalised regression (L1 regularisation, LASSO) along with cross-validation is used to identify enhancers for the gene.
• covariance testing is used to compute the p-value.

Installation

install.packages("devtools")
devtools::install_github("hadley/devtools")
devtools::install_github(c("vjbaskar/covTest", "vjbaskar/lenhancer"))

Example run

The below example run is for gene "Runx1". Please not the format of the four input data below and keep them the same for your input as well.

#### Input variables
library(lenhancer)
gene="Runx1"coCRE_cutoff=100000# 100Kb upper limit for considering regions as being mapped to the genesingleton_cutoff=20000# 20Kb upper limit for considering regions that are not coCREs in the modelcoCRE_corr_cutoff=0.5# The higher the value the tighter the community of CREsmin_tfb_events=2# lower limit for considering a region for walktrap community clusteringalphaVal=1# 1 is the lasso penalty, and 0 the ridge penalty.scale.predictors=TRUE# scale predictors True or Falsefamily="gaussian"# "gaussian","binomial","poisson","multinomial","cox","mgaussian"nfold=10# nfold cross-validation: >= 3 and can be as large as the sample (# of cell types) size (leave-one-out CV) 
data("expression")
data("regulation_signal")
data("regulation_tfbs")
data("region_gene_mapping")
gene_preds= findEnhancer(gene, expression, regulation_signal, regulation_tfbs, region_gene_mapping, min_tfb_events=2, coCRE_corr_cutoff=0.5, singleton_cutoff=20000, alphaVal=1, scale.predictors=TRUE, family="gaussian", nfoldxval=10)

Input data types

expression

expression data consists of gene expression values (such as log₂(FPKM)). gene as rows, cell types as columns.

> data(expression)
> head(expression)
ESCCESCMESCMESCMCPHBHEHPSergef1.1379911.0441181.8957152.2176831.6662600.6467382.0080901.2447081.802633Bcl7a2.2066942.5507762.8443593.4960803.1162373.4789032.1299271.4035242.419382Lnx22.0538562.5002701.8631021.8029511.9306531.3328762.2895952.3771882.195265Ppia10.07471410.3420729.82683710.1604309.5995579.5099541.9163611.5339362.086379Gkn1-9.965784-9.965784-9.965784-9.965784-9.965784-9.965784-9.965784-9.965784-9.965784Lrat-6.122783-6.090128-6.229214-7.160255-6.348240-7.412639-8.432611-5.641031-6.950677MACSergef2.6563687Bcl7a0.2551161Lnx22.8721921Ppia0.8198570Gkn1-9.9657843Lrat-5.6209541

regions

region_gene_mapping data frame consists of all the regulatory regions that are mapped to the genes in expression. There are many software such as ChIPseeker, ChIPpeakAnno and HOMER that can map regions to genes. If you map peaks to genes, and the peak lies within the gene you can mark their distance as any value less than singleton_cutoff (default, 20Kb)

> data("region_gene_mapping")
> head(region_gene_mapping)
chrstartendpeakgenedistance7786chr169251427692515568merge.105629Runx1858987792chr169251590792516390merge.105630Runx1850767797chr169252286392525592merge.105632Runx1758747815chr169257252892574113merge.105643Runx1273537820chr169257926592579735merge.105644Runx1217317824chr169260546492606424merge.105650Runx13000

NOTE: Please keep the colnames the same. findEnhancer function looks for c("peak", "gene", "distance") columns specifically in its codes and is case-senstive. NOTE that the region_gene_mapping$gene and row.names(expression) must overlap.

regulation signal

regulation_signal data consists of regions as rows and chromatin activity signal (such as log₂(FPKM)) for each cell type as columns.

NOTE that the colnames for expression and regulation_signal are the same. NOTE that the rownames for regulation_signal and region_gene_mapping$peak should be same.

> data("regulation_signal")
> head(regulation_signal)
ESCCESCMESCMESCMCPHBHEmerge.1056292.568202712.58062790.91683541.52567410.41299320.53309591.6784691-0.1319497merge.1056302.440974911.64298560.57976112.16399651.55150231.74881393.10671320.7631582merge.1056321.308105960.73718040.92147901.53224011.15761571.10915391.20703160.7889895merge.1056430.02400262-0.30863070.67541630.29509021.38505480.78988720.70383053.2788514merge.1056441.242603381.27217411.91497122.24363352.25431882.69060413.72038292.8612361merge.1056502.548227473.03953003.88577613.09257262.13488281.51246024.09953852.2658194HPMACmerge.1056290.99045660.8144609merge.1056301.64980182.3616716merge.1056321.20599284.0246107merge.1056431.51716225.4947875merge.1056443.39010370.6582151merge.1056501.3444846-3.3652031

TF binding data

regulation_tfbs data consists of TF binding data as 0 = non-binding and 1 = binding. The rownames are regions and colnames are TF ChIP-seq data. In theory, you can use any ChIP-seq data for computing overlap, but the most meaningful is to use the TFBS data for the cell types you are studying.

NOTE that rownames(regulation_tfbs) is same as rownames(regulation_signal)

> data("regulation_tfbs")
> head(regulation_tfbs)
ESC_ESRRBESC_NANOGESC_OCT4ESC_SOX2MES_CEBPBMES_ELK4MES_OCT4HB_CEBPBHB_GATA2merge.105677000000000merge.105727000000000merge.105746000000000merge.105753000000000merge.105769000000000merge.105629000000000HB_LMO2HB_TAL1HB_TEAD4HE_CEBPBHE_FLI1HE_LMO2HE_MEIS1HE_TAL1HP_CEBPBHP_FLI1merge.1056771100000000merge.1057270000000000merge.1057461110000000merge.1057530000000000merge.1057690000000000merge.1056290000000000HP_GATA1HP_GATA2HP_GFI1BHP_GFI1HP_LMO2HP_PU1HP_RUNX1HP_TAL1MAC_CEBPBMAC_FLI1merge.1056770000000000merge.1057270000000010merge.1057460000000000merge.1057530000000010merge.1057690000000010merge.1056290000000000MAC_LMO2MAC_PU1MAC_RUNX1MAC_TAL1merge.1056770100merge.1057270100merge.1057460000merge.1057530100merge.1057690100merge.1056290000

Output

The output gene_pred is a list and the main output are lambda_min = the predictors (regions) obtained when thresholding for lambda at min(MSE), lambda_1se = same as lambda_min but with lambda at min(MSE) + 1stdev and more conservative and p = p value

> summary(gene_preds$Runx1)
Length Class Mode expression 10 -none- numeric
meta 6 data.frame list reg_signal 10 data.frame list reg_tfbs 33 data.frame list lambda_1se 2 data.frame list lambda_min 6 data.frame list p 1 -none- numeric