library(diffloop)
library(diffloopdata)
library(ggplot2)
library(GenomicRanges)
library(ggrepel)
library(DESeq2)
library(digest)
library(edgeR)
library(limma)#Read in loop file ##If not a mango file ## Control_HDAC_SMC3 <- loopsMake("~/Desktop/HDAC_HiChIP/")
Control_HDAC_SMC3 <- loopsMake.mango("~/Desktop/HDAC_HiChIP")
celltypes <- c("Control_SMC3", "HDAC_SMC3")
Control_HDAC_SMC3 <- updateLDGroups(Control_HDAC_SMC3, celltypes)
dim(Control_HDAC_SMC3)## anchors interactions samples colData rowData
## 1 40390 134649 2 2 1#Subset loops and perform mango correction
Control_HDAC_SMC3 <- subsetLoops(Control_HDAC_SMC3, Control_HDAC_SMC3@rowData$loopWidth >= 5000)
Control_HDAC_SMC3_Mango <- mangoCorrection(Control_HDAC_SMC3, FDR = 0.01)
dim(Control_HDAC_SMC3_Mango)## anchors interactions samples colData rowData
## 1 24242 35855 2 2 3#Filter loops based on width and minumun number of loops per sample
Control_HDAC_SMC3_Filtered <- filterLoops(Control_HDAC_SMC3_Mango, width = 5000, nreplicates = 2, nsamples = 1)
dim(Control_HDAC_SMC3_Filtered)## anchors interactions samples colData rowData
## 1 14663 17810 2 2 3loopMetrics(Control_HDAC_SMC3_Filtered)## Control_SMC3 HDAC_SMC3
## unique 18947 78190#Call differential loops for single replicate
Control_HDAC_SMC3_Filtered@rowData$diffCount <- Control_HDAC_SMC3_Filtered@counts[,2] - Control_HDAC_SMC3_Filtered@counts[,1]
Control_HDAC_SMC3_Filtered@rowData$diffCount <- Control_HDAC_SMC3_Filtered@counts[,1] - Control_HDAC_SMC3_Filtered@counts[,2]
#Annotating Differential Loops to Peaks - Here using SMC3 peaks
Control_SMC3_pad <- rmchr(padGRanges(bedToGRanges("~/Desktop/SKCM-WNEG_M_vs_SKCM-WNEG_G_macs1_peaks.bed"), pad = 5000))
HDAC_SMC3_pad <- rmchr(padGRanges(bedToGRanges("~/Desktop/SKCM-WH81_M_vs_SKCM-WH81_G_macs1_peaks.bed"), pad = 5000))
#Pad enhancers using SMC3 peaks and promoters using human TSS
enhancer <- union(Control_SMC3_pad, HDAC_SMC3_pad)
promoter <- padGRanges(getHumanTSS(), pad = 5000)
#Annotate Differential Loops
Control_HDAC_SMC3_annotated_loops <- annotateLoops(Control_HDAC_SMC3_Filtered, enhancer = enhancer, promoter = promoter)
head(Control_HDAC_SMC3_annotated_loops)## An object of class "loops"
## Slot "anchors":
## GRanges object with 6 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] 1 6843325-6847408 *
## [2] 1 7265853-7269456 *
## [3] 1 7372237-7374913 *
## [4] 1 7722547-7731877 *
## [5] 1 7809412-7814422 *
## [6] 1 8019895-8023080 *
## -------
## seqinfo: 23 sequences from an unspecified genome; no seqlengths
##
## Slot "interactions":
## left right
## [1,] 1 2
## [2,] 1 3
## [3,] 4 5
## [4,] 6 20
## [5,] 7 24
## [6,] 8 16
##
## Slot "counts":
## Control_SMC3 HDAC_SMC3
## [1,] 0 2
## [2,] 0 2
## [3,] 4 3
## [4,] 0 2
## [5,] 0 2
## [6,] 0 2
##
## Slot "colData":
## sizeFactor groups
## Control_SMC3 1 Control_SMC3
## HDAC_SMC3 1 HDAC_SMC3
##
## Slot "rowData":
## loopWidth mango.FDR mango.P diffCount loop.type
## 1 422288 0.0001708118 1.510107e-05 -2 e-p
## 2 528208 0.0001056135 8.316585e-06 -2 e-p
## 3 84705 0.0031685219 4.464202e-04 1 e-e
## 4 340927 0.0002607523 2.548722e-05 -2 e-e
## 5 369299 0.0002064391 1.909710e-05 -2 e-e
## 6 183103 0.0015446761 2.030407e-04 -2 e-e#Summary of Differential Enhancer-Promoter and Enhancer-Enhancer SMC3-mediated loops
Control_HDAC_SMC3_annotated_loops_summary <- summary(Control_HDAC_SMC3_annotated_loops)
head(Control_HDAC_SMC3_annotated_loops_summary)## chr_1 start_1 end_1 chr_2 start_2 end_2 Control_SMC3 HDAC_SMC3 loopWidth
## 1 1 6843325 6847408 1 7265853 7269456 0 2 422288
## 2 1 6843325 6847408 1 7372237 7374913 0 2 528208
## 3 1 7722547 7731877 1 7809412 7814422 4 3 84705
## 4 1 8019895 8023080 1 8361130 8363700 0 2 340927
## 5 1 8063135 8066846 1 8432315 8436265 0 2 369299
## 6 1 8072825 8076859 1 8256794 8259096 0 2 183103
## mango.FDR mango.P diffCount loop.type region
## 1 0.0001708118 1.510107e-05 -2 e-p chr1:6818325-7294456
## 2 0.0001056135 8.316585e-06 -2 e-p chr1:6818325-7399913
## 3 0.0031685219 4.464202e-04 1 e-e chr1:7697547-7839422
## 4 0.0002607523 2.548722e-05 -2 e-e chr1:7994895-8388700
## 5 0.0002064391 1.909710e-05 -2 e-e chr1:8038135-8461265
## 6 0.0015446761 2.030407e-04 -2 e-e chr1:8047825-8284096Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.