aryeelab/diffloop-paper

Diffloop Visualization of interactions

Opened this issue · 5 comments

koushik20 commented

Hi,

First, thank you for this wonderful package and documentation.

I am working on finding a differential loop calling for my 4 samples (8 reps) data using diffloop and below is my loopobject
`> dim(km_anno)

anchors interactions samples colData rowData
1 49597 96655 8 2 7`

but when I do visualization of the interaction by using plotTopLoops() or loopPlot() as I understand giving me an aggregate of interactions. Attached is an example pic for reference. Is there any other way to view each and every interaction rather than aggregate interactions that make this plot look almost the same across the cell lines?
Thank you!
Screen Shot 2022-04-19 at 3 37 24 PM

caleblareau commented
koushik20 commented

I didn't find any difference while zooming in on the plots. But thanks for the claratification on how the loops are plotted.

koushik20 commented

This is the code I am using to generate the loops. Am I doing something wrong?


## loading packages
library(diffloop)
library(diffloopdata)
library(ggplot2)
library(GenomicRanges)
library(ggrepel)
library(DESeq2)
library(dplyr)


## diffloop downstream analysis with bedpe files
bed_dir <- file.path("/Users/ayalurik/Documents/HiChIP_diffloop/diffloop")
bed_dir

samples <- c("kura_1","kura_2","uwb_1","uwb_2",
             "ft246_1","ft246_2","ft33_1","ft33_2")
full <- loopsMake(bed_dir, samples)
celltypes <- c("kura", "kura", "uwb", "uwb", "ft246", "ft246", "ft33", "ft33")
full <- updateLDGroups(full, celltypes)
head(full, 4)
dim(full)

## Quality control
# remove and loops that merged together from import
full1 <- subsetLoops(full, full@rowData$loopWidth >= 5000)
dim(full1)

# Filtering loops that are FDR > 0.01
noCNV_corrected <- mangoCorrection(full1, FDR = 0.01)
dim(noCNV_corrected)

# filtering loops that are present strongly (>= 5 PETs) in one replicate
# but absent (== 0 PETs) in the other replicate
cm <- noCNV_corrected@counts
k_dis <- ((cm[,1]>=5 &cm[,2]==0)|(cm[,2]>=5&cm[,1]==0))
m_dis <- ((cm[,3]>=5 &cm[,4]==0)|(cm[,4]>=5&cm[,3]==0))
qc_filt <- subsetLoops(noCNV_corrected, !(k_dis | m_dis))
dim(qc_filt)

p1 <- loopDistancePlot(qc_filt)
p1

## counts the number of loops for each sample and returns whether they are single,
## self, unique or none
loopMetrics(qc_filt)

# PC plot
pcp1dat <- qc_filt
pcp1dat@colData$sizeFactor <- 1
samples <- c("kura_1", "kura_2", "uwb_1", "uwb_2",
             "ft246_1", "ft246_2", "ft33_1", "ft33_2")
pcp1 <- pcaPlot(pcp1dat) + geom_text_repel(aes(label=samples)) + 
  ggtitle("PC Plot with no
  Size Factor Correction") + theme(legend.position="none")
pcp1

samples <- c("kura_1", "kura_2", "uwb_1", "uwb_2",
             "ft246_1", "ft246_2", "ft33_1", "ft33_2")
pcp2 <- pcaPlot(qc_filt) + geom_text_repel(aes(label=samples)) + 
  ggtitle("PC Plot with Size Factor Correction") +
  theme(legend.position = "none")
pcp2


## Differential Loop Calling
km_filt <- qc_filt
dim(km_filt)

# First model using edgeR over-dispersed Poisson regression
km_res <- loopAssoc(km_filt, method = "edgeR", coef = 2)
head(km_res@rowData)

## Epigenetic Annotation
h3dir <- file.path("/Users/ayalurik/Documents/HiChIP_diffloop")
kura <- paste0(h3dir, "/", "kura_H3k27ac_hg19_narrowpeak.bed")
uwb <- paste0(h3dir, "/", "UWB_H3k27ac_hg19_narrowpeak.bed")
ft246 <- paste0(h3dir, "/", "FT33_H3k27ac_hg19_narrowpeak.bed")
ft33 <- paste0(h3dir, "/", "FT246_H3k27ac_hg19_narrowpeak.bed")

h3k27ac.k <- rmchr(padGRanges(bedToGRanges(kura), pad = 1000))
h3k27ac.u <- rmchr(padGRanges(bedToGRanges(uwb), pad = 1000))
h3k27ac.f2 <- rmchr(padGRanges(bedToGRanges(ft246), pad = 1000))
h3k27ac.f3 <- rmchr(padGRanges(bedToGRanges(ft33), pad = 1000))
ku <- GenomicRanges::union(h3k27ac.k, h3k27ac.u)
ft <- GenomicRanges::union(h3k27ac.f2, h3k27ac.f3)
enhancer <- GenomicRanges::union(ku,ft)
promoter <- padGRanges(getHumanTSS(), pad = 1000)
km_anno <- annotateLoops(km_res, enhancer = enhancer, promoter = promoter)

# No.of differential loops when FDR <= 0.05
diflop <- subsetLoops(km_res, km_res@rowData$FDR <= 0.05)

# exploring out results
summary <- as.data.frame(km_anno@anchors)
stats <- as.data.frame(km_anno@rowData)
interactions <- as.data.frame(km_anno@interactions)

## saving a rda file
# save(km_anno, file = "diffloop_Kura_UWB_vs_FT246_FT33.rda")

## Visulaization
chr20reg <- GRanges(seqnames = c("19"), 
                    ranges = IRanges(start = c(42347500),end = c(49542500)))
p3 <- loopPlot(km_anno, chr20reg)

plotTopLoops(km_anno, organism = "h", FDR = 0.01, n = 10)

manyLoopPlots(km_anno, chr20reg)

Thank you!

caleblareau commented
koushik20 commented

Thank you for your quick response. Sure will run and check.