Home

Preliminary Steps

To build CRHs, you need to use the Activity-by-contact Score (Fulco et al. 2019). Please refer to: https://github.com/broadinstitute/ABC-Enhancer-Gene-Prediction to have deeper details on how to use the ABC-Score.

Preliminaries

Before introducing the code used to create CRHs

Package dependencies

Before executing any code, please be sure that all necessary packages are installed either through CRAN repository or Bioconductor

library(igraph)
library(GenomicRanges)
library(rtracklayer)
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(ensembldb)
library(EnsDb.Hsapiens.v86)

CRH construction

Since you obtained the output file from the ABC-Score, you should use process.ABC from code_CRHs script. Basically, the function adds new columns (typeOf, startProm, endProm) and removes X chromosome contacts, since we focus on autosomal chromosomes in our analysis. The data EnhancerPredictions_*.txt from the repo can be directly used here.

output_from_ABC = process.ABC("EnhancerPredictions_NEU_iPSC.txt")

Now we have file in the right format, we are able to create CRHs.

#If you need the pairs of promoters and distal elements in Pairs format 
#Useful when you want to compare CRHs directly with ABC-Score predictions or between different tissues
pairs = Pairs.from.ABC(output_from_ABC)

#CRHs can be built from the following code line
graph_from_ABC = create.CRHs.ABC(output_from_ABC)
#the return object is a igraph object with bipartite networks
#You can use all the functions provided by igraph to analyze graph_from_ABC

You need to add membership (the CRH), where the promoter or distal element belongs to analyze CRHs in the right way.

#Here we directly modify the data.frame, the column membership will be added
output_from_ABC = add.membership(graph_from_ABC, output_from_ABC)

CRH analysis

Components can be generated with igraph components function.

components = components(graph_from_ABC)
#returns:
#no : the number of clusters
#csize: the size of clusters
#membership: the cluster to which each element belongs

If you want to plot one specific CRH

decompose.ABC = decompose(graph_from_ABC)
plot(decompose.ABC[[1]], vertex.label=ifelse(names(V(decompose.ABC[[1]]))%in%unique(output_from_ABC$TargetGene), names(V(decompose.ABC[[1]])), NA))

(https://github.com/lmangnier/Hi-C_analysis/main/Rplot01.png)

Now if you need the number of unique promoters or distal elements per CRHs.

#Number of unique promoters
N_prom_per_CRH = aggregate(TagetGene~membership, output_from_ABC , function(x) length(unique(x)))
#Number of unique distal elements
N_distal_per_CRH = aggregate(name~membership, output_from_ABC , function(x) length(unique(x)))

The mean number of contacts per promoter or distal element

#Mean number of contacts by promoter at CRH-level
N_connect_prom_per_CRH = aggregate(TargetGene~membership,output_from_ABC,function(x) mean(table(x)))
#Mean number of contacts by distal element at CRH-level
N_connect_dist_per_CRH = aggregate(name~membership,output_from_ABC, function(x) mean(table(x)))