## ---- echo=FALSE, results="hide", message=FALSE-------------------------- require(knitr) opts_chunk$set(error=FALSE) ## ----style, echo=FALSE, results='asis'----------------------------------- BiocStyle::markdown() ## ----setup, echo=FALSE, message=FALSE------------------------------------ require(InteractionSet) ## ------------------------------------------------------------------------ all.regions <- GRanges("chrA", IRanges(0:9*10+1, 1:10*10)) ## ------------------------------------------------------------------------ index.1 <- c(1,5,10) index.2 <- c(3,2,6) region.1 <- all.regions[index.1] region.2 <- all.regions[index.2] ## ------------------------------------------------------------------------ gi <- GInteractions(region.1, region.2) ## ------------------------------------------------------------------------ gi <- GInteractions(index.1, index.2, all.regions) gi ## ------------------------------------------------------------------------ anchors(gi) ## ------------------------------------------------------------------------ anchors(gi, id=TRUE) ## ------------------------------------------------------------------------ regions(gi) ## ------------------------------------------------------------------------ temp.gi <- gi anchorIds(temp.gi) <- list(1:3, 5:7) ## ------------------------------------------------------------------------ temp.gi <- gi annotation <- rep(c("E", "P", "N"), length.out=length(all.regions)) regions(temp.gi)$anno <- annotation ## ------------------------------------------------------------------------ anchors(temp.gi) ## ------------------------------------------------------------------------ temp.gi <- gi super.regions <- GRanges("chrA", IRanges(0:19*10+1, 1:20*10)) replaceRegions(temp.gi) <- super.regions ## ------------------------------------------------------------------------ temp.gi <- gi extra.regions <- GRanges("chrA", IRanges(10:19*10+1, 11:20*10)) appendRegions(temp.gi) <- extra.regions ## ------------------------------------------------------------------------ metadata(gi)$description <- "I am a GInteractions object" metadata(gi) ## ------------------------------------------------------------------------ set.seed(1000) norm.freq <- rnorm(length(gi)) # obviously, these are not real frequencies. gi$norm.freq <- norm.freq mcols(gi) ## ------------------------------------------------------------------------ sub.gi <- gi[1:2] sub.gi anchors(sub.gi, id=TRUE) ## ------------------------------------------------------------------------ identical(regions(gi), regions(sub.gi)) ## ------------------------------------------------------------------------ c(gi, sub.gi) new.gi <- gi regions(new.gi) <- resize(regions(new.gi), width=20) c(gi, new.gi) ## ------------------------------------------------------------------------ new.gi <- gi anchorIds(new.gi) <- list(1:3, 5:7) combined <- c(gi, new.gi) ## ------------------------------------------------------------------------ combined <- swapAnchors(combined) ## ------------------------------------------------------------------------ order(combined) sorted <- sort(combined) anchors(sorted, id=TRUE) ## ------------------------------------------------------------------------ anchors(sorted, type="first") ## ------------------------------------------------------------------------ doubled <- c(combined, combined) duplicated(doubled) ## ------------------------------------------------------------------------ anchorIds(new.gi) <- list(4:6, 1:3) new.gi <- swapAnchors(new.gi) swap.gi <- swapAnchors(gi) match(new.gi, swap.gi) ## ------------------------------------------------------------------------ new.gi==swap.gi ## ------------------------------------------------------------------------ all.regions <- GRanges(rep(c("chrA", "chrB"), c(10, 5)), IRanges(c(0:9*10+1, 0:4*5+1), c(1:10*10, 1:5*5))) index.1 <- c(5, 15, 3, 12, 9, 10) index.2 <- c(1, 5, 11, 13, 7, 4) gi <- GInteractions(index.1, index.2, all.regions) ## ------------------------------------------------------------------------ pairdist(gi) ## ------------------------------------------------------------------------ intrachr(gi) ## ------------------------------------------------------------------------ of.interest <- GRanges("chrA", IRanges(30, 60)) olap <- findOverlaps(gi, of.interest) olap ## ------------------------------------------------------------------------ anchors(gi[queryHits(olap)]) ## ------------------------------------------------------------------------ paired.interest <- GInteractions(of.interest, GRanges("chrB", IRanges(10, 40))) olap <- findOverlaps(gi, paired.interest) olap ## ------------------------------------------------------------------------ all.genes <- GRanges("chrA", IRanges(0:9*10, 1:10*10)) all.enhancers <- GRanges("chrB", IRanges(0:9*10, 1:10*10)) out <- linkOverlaps(gi, all.genes, all.enhancers) head(out) ## ------------------------------------------------------------------------ out <- linkOverlaps(gi, all.genes) head(out) ## ------------------------------------------------------------------------ all.chrs <- as.character(seqnames(regions(gi))) group.by.chr <- paste0(all.chrs[anchors(gi, type="first", id=TRUE)], "+", all.chrs[anchors(gi, type="second", id=TRUE)]) ## ------------------------------------------------------------------------ swap.gi <- swapAnchors(gi) boundingBox(swap.gi, group.by.chr) ## ------------------------------------------------------------------------ sgi <- GInteractions(index.1, index.2, all.regions, mode="strict") class(sgi) ## ------------------------------------------------------------------------ sgi <- as(gi, "StrictGInteractions") ## ------------------------------------------------------------------------ anchorIds(sgi) <- list(7:12, 1:6) anchors(sgi, id=TRUE) ## ------------------------------------------------------------------------ set.seed(1000) Nlibs <- 4 counts <- matrix(rpois(Nlibs*length(gi), lambda=10), ncol=Nlibs) lib.data <- DataFrame(lib.size=seq_len(Nlibs)*1e6) iset <- InteractionSet(counts, gi, colData=lib.data) iset ## ------------------------------------------------------------------------ norm.freq <- matrix(rnorm(Nlibs*length(gi)), ncol=Nlibs) iset2 <- InteractionSet(list(counts=counts, norm.freq=norm.freq), gi, colData=lib.data) iset2 ## ------------------------------------------------------------------------ assay(iset) assay(iset2, 2) ## ------------------------------------------------------------------------ interactions(iset) ## ------------------------------------------------------------------------ anchors(iset, id=TRUE) # easier than anchors(interactions(iset)), id=TRUE) regions(iset) ## ------------------------------------------------------------------------ lib.size <- seq_len(Nlibs) * 1e6 colData(iset)$lib.size <- lib.size iset$lib.size <- lib.size # same result. ## ------------------------------------------------------------------------ new.gi <- interactions(iset) new.iset <- iset interactions(new.iset) <- new.gi ## ------------------------------------------------------------------------ regions(interactions(new.iset))$score <- 1 regions(new.iset)$score <- 1 # same as above. ## ------------------------------------------------------------------------ iset[1:3,] iset[,1:2] ## ------------------------------------------------------------------------ rbind(iset, iset[1,]) ## ------------------------------------------------------------------------ cbind(iset, iset[,3]) ## ------------------------------------------------------------------------ row.indices <- 1:10 col.indices <- 9:15 row.regions <- all.regions[row.indices] col.regions <- all.regions[col.indices] Nr <- length(row.indices) Nc <- length(col.indices) counts <- matrix(rpois(Nr*Nc, lambda=10), Nr, Nc) cm <- ContactMatrix(counts, row.regions, col.regions) ## ------------------------------------------------------------------------ cm <- ContactMatrix(counts, row.indices, col.indices, all.regions) cm ## ------------------------------------------------------------------------ as.matrix(cm) ## ------------------------------------------------------------------------ anchors(cm, type="row") ## ------------------------------------------------------------------------ regions(cm) ## ------------------------------------------------------------------------ temp.cm <- cm as.matrix(temp.cm[1,1]) <- 0 ## ------------------------------------------------------------------------ anchorIds(temp.cm) <- list(1:10, 1:7) ## ------------------------------------------------------------------------ regions(temp.cm)$GC <- runif(length(regions(temp.cm))) # not real values, obviously. regions(temp.cm) ## ------------------------------------------------------------------------ metadata(temp.cm)$description <- "I am a ContactMatrix" metadata(temp.cm) ## ------------------------------------------------------------------------ cm[1:5,] cm[,1:5] ## ------------------------------------------------------------------------ anchors(cm[,1:5], type="column", id=TRUE) ## ------------------------------------------------------------------------ t(cm) ## ------------------------------------------------------------------------ rbind(cm, cm[1:5,]) # extra rows cbind(cm, cm[,1:5]) # extra columns ## ------------------------------------------------------------------------ temp.cm <- cm anchorIds(temp.cm) <- list(10:1, 15:9) anchors(sort(temp.cm), id=TRUE) ## ------------------------------------------------------------------------ order(temp.cm) ## ------------------------------------------------------------------------ temp.cm <- rbind(cm, cm) duplicated(temp.cm) unique(temp.cm) ## ------------------------------------------------------------------------ pairdist(cm) ## ------------------------------------------------------------------------ intrachr(cm) ## ------------------------------------------------------------------------ of.interest <- GRanges("chrA", IRanges(50, 100)) olap <- overlapsAny(cm, of.interest) olap ## ------------------------------------------------------------------------ and.mask <- outer(olap$row, olap$column, "&") or.mask <- outer(olap$row, olap$column, "|") ## ------------------------------------------------------------------------ my.matrix <- as.matrix(cm) my.matrix[!and.mask] <- NA my.matrix ## ------------------------------------------------------------------------ olap <- overlapsAny(cm, GInteractions(of.interest, GRanges("chrB", IRanges(1, 20)))) olap ## ------------------------------------------------------------------------ counts <- rpois(length(gi), lambda=10) desired.rows <- 2:10 desired.cols <- 1:5 new.cm <- inflate(gi, desired.rows, desired.cols, fill=counts) new.cm anchors(new.cm, id=TRUE) as.matrix(new.cm) ## ------------------------------------------------------------------------ inflate(gi, GRanges("chrA", IRanges(50, 100)), GRanges("chrA", IRanges(10, 50)), fill=counts) ## ------------------------------------------------------------------------ inflate(gi, "chrA", "chrB", fill=counts) ## ------------------------------------------------------------------------ new.cm <- inflate(iset, desired.rows, desired.cols, sample=3) as.matrix(new.cm) ## ------------------------------------------------------------------------ new.iset <- deflate(cm) new.iset ## ------------------------------------------------------------------------ deflate(cm, collapse=FALSE) ## ------------------------------------------------------------------------ x <- GRanges("chrA", IRanges(42, 48)) rse <- linearize(iset, x) rse ## ------------------------------------------------------------------------ rowRanges(rse) ## ------------------------------------------------------------------------ sessionInfo()