Cis regulatory Elements and gene pair

[maxglenwall]

In an article by Trevino et al., 2020 they extracted data related to links between gene-distal CRE accessibility to gene expression. The identified CRE-gene pairs that represent potential enhancer-gene interactions (Table S2

Peak chromosome	Peak start	Peak end	Gene ID	Gene symbol	Gene TSS position	Gene strand	Distance peak to TSS	correlation	p-value	FDR adjusted p-value
chr1	191217	191718	ENSG00000243485	MIR1302-2HG	29553	+	-161913	0,588541157	1,90E-202	6,66E-200
chr1	981547	982048	ENSG00000188976	NOC2L	959308	-	22488	-0,415042811	7,02E-163	4,46E-161
chr1	981547	982048	ENSG00000272512	AL645608.8	998050	-	-16252	0,400419562	2,85E-159	1,59E-157
chr1	981547	982048	ENSG00000188290	HES4	1000171	-	-18373	0,505822625	4,59E-184	6,82E-182

I extracted this information Peak2GeneLinks function (p2g[[1]]) but it is only outputting information on chromosome start and end location, length of loop or length of the gene (i guess), correlation, and FDR.

How can we extract this information shown in Table2 using ArchR?

Thanks in advance.

[rcorces]

I'm not sure I understand your question but perhaps this excerpt from an upcoming version of the documentation will help:

---

To identify peak-to-gene links in ArchR, we use the addPeak2GeneLinks() function.

projHeme5 <- addPeak2GeneLinks(
	ArchRProj = projHeme5,
	reducedDims = "IterativeLSI"
)

We can then retrieve these peak-to-gene links in a similar fashion to how we retrieved co-accessibility interactions by using the getPeak2GeneLinks() function. As we saw previously, this function allows for a user-specified cutoff for correlation and resolution for linkages.

p2g <- getPeak2GeneLinks(
	ArchRProj = projHeme5,
	corCutOff = 0.45,
	resolution = 1,
	returnLoops = FALSE
)

When returnLoops is set to false, this function returns a DataFrame object anaolgous to the DataFrame object returned by getCoAccessibility(). The primary difference is that the indexes for the scATAC-seq peaks are stored in a column called idxATAC and the indexes for the scRNA-seq gene are stored in a column called idxRNA.

p2g

Here, idxATAC and idxRNA refer to the row indices of the peak or gene in the corresponding geneSet or peakSet which can be accessed via the metadata component of the p2g object.

metadata(p2g)

It is important not to confuse idxRNA with the value in the idx column of the geneSet. Within the geneSet, idx corresponds to the chronological position of the given gene across each chromosome. As such, idx is only unique across an individual chromosome and is not relevant for mapping idxRNA to a gene name.

So, if we wanted to add the gene name and peak coordinates to our p2g DataFrame object, we would do the following:

p2g$geneName <- mcols(metadata(p2g)$geneSet)$name[p2g$idxRNA]
p2g$peakName <- (metadata(p2g)$peakSet %>% {paste0(seqnames(.), "_", start(.), "_", end(.))})[p2g$idxATAC]
p2g

You may also notice that there is other information stored within the metadata() of p2g including pointers to files called seATAC and seRNA which represent SummarizedExperiment objects for the ATAC-seq and RNA-seq data used to identify peak-to-gene linkages and these each include both the raw and normalized data matrices used in the analysis.

metadata(p2g)$seATAC
metadata(p2g)$seRNA

If we set returnLoops = TRUE, then getPeak2GeneLinks() will return a loop track GRanges object that connects the peak and gene. As for co-accessibility, the start and end of the IRanges object represent the position of the peak and gene being linked. When resolution = 1, this links the center of the peak to the single-base TSS of the gene.

p2g <- getPeak2GeneLinks(
	ArchRProj = projHeme5,
	corCutOff = 0.45,
	resolution = 1,
	returnLoops = TRUE
)

p2g[[1]]

[ankushs0128]

Thanks, Ryan for the detailed explanation. I was wondering how can we merge these two data frames together as an only column that we can utilize is a correlation, which is not very optimal to merge to dataframes.

[rcorces]

I'm not sure what two data frames you are referring to. What are you trying to merge?

[ankushs0128]

p2g DataFrame with gene name and then getPeak2GeneLinks() that will return a loop track GRanges object that connects the peak and gene

[rcorces]

No thats not relevant here. The "loops" do not have a one-to-one mapping with entries in the DataFrame. Depending on the resolution used, multiple p2g links might get merged into a single loop. You cant annotate the "loops" in the same way you annotate the data frame. You would be better off creating your own loops object from the dataframe if thats what you want to do.

[maxglenwall]

Thanks, Ryan and Ankush,
I have both the data frames containing this information and I am now wondering How to create the required loop objects from the data frame? I can get data frame 500 basepairs of the peak the e.g chr3_110886497_110886997 and in another data frame, Start of the peak linked with the TSS of the genes in one data frame with a width of 18000

seqnames	start	end	width
chr3	110886747	111070071	183325

In this data frame, the start of the peak region is 250 base pairs plus/minus the peak start base pairing of the peak region.
In other loops , where the TSS (chr14 | 100725705 | 10084605) are lets say 120000 bases apart from the detected peak (chr14_100845801_100846301).

In such cases, How can we create our own loop objects from the dataframe?

Apologies

Thanks in advance

Max

[aa20g217]

Thanks, Ryan, for the explanation. In the p2g data frame (returnLoops=FALSE), we have information about the peaks to which particular genes are linked. We are interested in finding the length of each loop (start and end position of a loop not the peak position) for a particular gene. This info is present in the loop track GRanges object, but there is no way to extract it based on gene names.
Is there any straight forward way to do this?

[rcorces]

In my previous post, I showed how to use the p2g dataframe and extract informtation on the corresponding peaks and genes given the idxRNA and idxATAC columns from the p2g dataframe. If you know the peak and you know the gene then you know the corresponding positions and you know what the loop is.

> p2g
DataFrame with 38758 rows and 6 columns
        idxATAC    idxRNA Correlation         FDR  VarQATAC   VarQRNA
      <integer> <integer>   <numeric>   <numeric> <numeric> <numeric>
1            52         2    0.513157 2.72469e-33  0.904505  0.389495
2            67         4    0.561850 4.60984e-41  0.573333  0.274663
3            34         5    0.505633 3.36115e-32  0.755000  0.790495
4             7         7    0.657990 6.27030e-61  0.613340  0.452126
5            12         7    0.521996 1.31288e-34  0.927350  0.452126
...         ...       ...         ...         ...       ...       ...
38754    143363     18589    0.467486 4.49743e-27  0.702483  0.652331
38755    143352     18590    0.523039 9.12450e-35  0.663340  0.972743
38756    143356     18590    0.506386 2.62159e-32  0.532524  0.972743
38757    143368     18590    0.579947 2.65930e-44  0.741444  0.972743
38758    143393     18593    0.535505 1.06460e-36  0.815397  0.334176

[ankushs0128]

Hopefully its a correct way to get loop object from peak2genelinks

p2g <- getPeak2GeneLinks(
    ArchRProj = ArchRproject,
    corCutOff = 0.45,
    resolution = 1,
    returnLoops = FALSE
)
p2g
metadata(p2g)
metadata(p2g)$seATAC




p2g$geneName <- mcols(metadata(p2g)$geneSet)$name[p2g$idxRNA]
tssgenes<- as.data.frame(metadata(p2g)$geneSet)
colnames(tssgenes)[2] <- "StartTSS"
colnames(tssgenes)[3] <- "EndTSS"
p2g$peakName <- (metadata(p2g)$peakSet %>% {paste0(seqnames(.), "_", start(.), "_", end(.))})[p2g$idxATAC]
temp1=as.data.frame(p2g)
temp  <- strsplit(p2g$peakName, "_")
temp <- matrix(unlist(temp), ncol=3, byrow=TRUE)
temp <- as.data.frame(temp)
temp1$chr <- temp$V1
temp1$start <- temp$V2
temp1$end <- temp$V3
mergedp2g <- merge(temp1,tssgenes ,by.x="geneName", by.y= "name")
mergedp2g$peakCentre=as.numeric(mergedp2g$start)+250
mergedp2g$looplength=as.numeric(mergedp2g$StartTSS)-as.numeric(mergedp2g$peakCentre)

write.table(p2g,"__peak2genelinks_corr_0.45_res-1_loopobject.txt")