In jamesdalg/GS3DKViz: Extra functions and visualization for the GenoScan3DKnock package

---

title: "GS3DKViz_Arc_Plots" author: "James Dalgleish" date: "7/9/2019" output: rmarkdown::html_vignette vignette: > %\VignetteEngine{knitr::rmarkdown} %\VignetteIndexEntry{GS3DKViz_Arc_Plots} %\VignetteEncoding{UTF-8}

---

---

title: "GS3DKViz_Arc_Plots" author: "James Dalgleish" date: "12/2/2020" output: html_document

---

knitr::opts_chunk$set(echo = TRUE,message=F,error=F,warning=F)
knitr::opts_chunk$set(root.dir="./inst/extdata/")
library(Gviz)
library(GenomicInteractions)
library(magrittr)
library(GS3DKViz)

We'll start by loading a few libraries that will enable us to generate plots and loading the data.

For simple plots, we've created a convenient function for plotting.

 gs3dk_gint=readGint(system.file(
   "extdata","signficant_genes_genoscan_3D_knock.csv",
   package = "GS3DKViz"))
   filtered_gs3dk_gint=gs3dk_gint[(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::pull(trait))=="FEV1"]
   plotInteraction(gint = filtered_gs3dk_gint,
  chr=unique(janitor::clean_names(tibble::as_tibble(mcols(filtered_gs3dk_gint)))$chr),
  bounds_offset=1.5e4,
  main="Promoter–enhancer interactions, FEV1"
  )

A details plot is below, with unconnected promoters in orange. The code behind this plot will be explained towards the end of the tutorial.

extra_prom_gr_tm=data.table::fread(system.file(
"extdata","more_promoters_topmed.csv",
  package = "GS3DKViz")
) %>% janitor::clean_names( ) %>%  GRanges()
GS3DKViz::plotInteractionDetail(gint = gs3dk_gint,
  chr="chr19", #unique(janitor::clean_names(tibble::as_tibble(mcols(gs3dk_gint)))$chr)
  bounds_offset=0,#1.5e4
  main="Promoter–enhancer interactions, AD",
  extra_prom_gr = extra_prom_gr_tm )

I have also included a custom function to convert data frames to GenomicInteractions objects, but in any way that you create such an object is fine. The counts column will contain the gene-enhancer score (from either ABC or GeneCards). I implement a scaling algorithm here so that the larger scaled GeneHancer does not have vastly larger arc height. GenomicInteractions Objects just require 6 columns in the following order: chr1-chromosome of the first range start1-start positions of the first range end1-end positions of the first range chr2-Chromosomes for the second set of ranges start1-start position of the second range start2-start position of the second range Thus, we'll need to split our data into start and end positions for both promoters (range 1, which can be accessed by the standard anchorOne function), and enhancers (range 2, which can be accessed with the anchorTwo function). I will provide the detailed code here on how our results were generated, but if desired, one can skip to the next section and just pull in our sample data. If executing this code, be sure to install data.table and tidyverse.

Essentially, GenomicInteractions object will work to reproduce the plots, so long as the counts have a meaningful measure of interaction score and all the interactions are on the same chromosome.

#read in data
gs3dk_gint=data.table::fread(system.file(
"extdata","signficant_genes_genoscan_3D_knock.csv",
  package = "GS3DKViz")) %>% 
  #automatically clean column names, converting to standard format.
janitor::clean_names() %>%
  #split the gene position into start and end
  tidyr::separate(gene_position,c("genestart","geneend"),  "-") %>% 
    #split the promoter position into start and end
  tidyr::separate(promoter_position,c("promoterstart","promoterend"),  "-") %>% 
      #split the enhancer position into start and end
  tidyr::separate(best_enhancer_position,c("enhstart","enhend"),  "-") %>%
  #convert the promoter starts to numeric, assign replacements for missing data.
  dplyr::mutate(promoterstart=as.numeric(promoterstart),chr1=chr,chr2=chr,
                #replace missing promoter starts with TSS
  promoterstart=ifelse(is.na(promoterstart),genestart,promoterstart),
  #replace missing promoter ends with TSS
  promoterend=ifelse(is.na(promoterend),genestart,promoterend)) %>%  
  #create metadata column that determines whether an interaction is from
  #the ABC model or from GeneHancer.
  dplyr::mutate(enhancer_type=ifelse(grepl(pattern="GH",x=enhancer_id),"GH","ABC")) %>% 
  #rearrange data columns for GenomicInteractions Objects
  dplyr::select(chr1,promoterstart,promoterend,chr2,enhstart,enhend,dplyr::everything()) %>% 
  #rename columns to standard GenomicInteractions format.
  dplyr::rename(start1=promoterstart,end1=promoterend,start2=enhstart,end2=enhend) %>%
  #remove chromosomes from second range (ranges should only be positions).
  dplyr::mutate(start2=gsub("chr19: ","",start2)) %>% 
  #create the GI.
  makeGenomicInteractionsFromDataFrame() %>%
  #remove arcs where the promoter and enhancers overlap.
  .[which(calculateDistances(.)!=0),]
#with the GenomicInteractions object created, we'll normalize the interaction score
#and store it in the counts column for plotting.
mcols(gs3dk_gint)=gs3dk_gint %>% tibble::as_tibble() %>% janitor::clean_names() %>% 
  dplyr::group_by(enhancer_type) %>% dplyr::mutate(interaction_score=interaction_score/mean(interaction_score,na.rm=T))
mcols(gs3dk_gint)$counts<-ifelse(!is.na(mcols(gs3dk_gint)$interaction_score),
                                 mcols(gs3dk_gint)$interaction_score,1) 
#save the object

save(list = c("gs3dk_gint"),file = "gs3dk_sample_data.rda",compress = "xz",
     compression_level = 9)

Now with the data loaded, we can create tracks with Gviz and plot them. First we'll load the data, extract promoter and enhancer granges objects We'll use PSMA4 to get the promoter, but let the bounds be defined by the min and max ranges in the available data.

load(system.file(
"extdata","gs3dk_sample_data.rda",
  package = "GS3DKViz"))
#extract promoter and enhancer data from GenomicInteractions object.
promoter_gr=anchorOneWithMetadata(gs3dk_gint)
enhancer_gr=anchorTwoWithMetadata(gs3dk_gint)
geneSymbol="PSMA4"
#set the chromosome containing the PSMA4 gene to be the plotted chromosome.
selected_chromosome=as.character(gs3dk_gint %>% as.data.frame() %>%
  janitor::clean_names() %>% dplyr::filter(gene_id==geneSymbol))

We will now construct the interaction, promoter, and enhancer tracks.

interaction_track <- InteractionTrack(gs3dk_gint, name = "Interaction", chromosome = "chr15")

bounds=c(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(start1) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(end1) %>% max(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(start2) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(end2) %>% max())
promoterTrack <- AnnotationTrack(promoter_gr, genome="hg38", name="Promoters",
                             id=mcols(promoter_gr)$gene_id,  featureAnnotation="id")
enhancerTrack <- AnnotationTrack(enhancer_gr, genome="hg38", name="Enhancers",
                             id=mcols(enhancer_gr)$gene_id,  featureAnnotation="id")
feature(enhancerTrack)<-mcols(gs3dk_gint)$enhancer_type

displayPars(promoterTrack) <- list(fill = "olivedrab1", col = NA, 
                                   fontcolor.feature = "black", fontsize=8,
                                   just.group="below",rotation=90,rotation.group=90,rotation.item=90)
displayPars(enhancerTrack) <- list(fill = "mediumpurple1", col = NA, 
                                   fontcolor.feature = "black", fontsize=10,
                                   just.group="below",rotation.item=90,
                                   collapse=T,mergeGroups=T,showOverplotting=T,groupAnnotation="group",group=mcols(gs3dk_gint)$enhancer_type)
displayPars(interaction_track) <- list(fill = "deepskyblue", col = NA, 
                                   fontcolor.feature = "black", fontsize=8,
                                   just.group="below",plot.anchors=T,plot.outside=T,col.outside="lightblue",                                   interaction.measure="counts",
                                   interaction.dimension="height",
                                   col.interactions="black",
                                   plot.trans=T,
                                              fontsize.legend=200
                                   )
itrack <- IdeogramTrack(genome = "hg38", chromosome =as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(gene_id==geneSymbol) %>% dplyr::pull(seqnames1)) )
gtrack <- GenomeAxisTrack()
displayPars(enhancerTrack)=list(rotation.item=0)
plotTracks(list(itrack,gtrack,interaction_track,promoterTrack,enhancerTrack),
           chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(gene_id==geneSymbol) %>% dplyr::pull(seqnames1)),
           from = (min(bounds))-1.5e4,
           to = (max(bounds))+1.5e4,
           type = c("b"),
           showSampleNames = TRUE, 
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=12,collapse=T,min.width=5,mergeGroups=T
           ,stacking="squish" 
           ,main="Promoter–enhancer interactions, FEV1/FVC",
           background.title = "black"
)

Now, we'll add more promoters from a separate csv file that do not have enhancer links in a separate track (which will be vertically combined later). Track properties are copied from the previous promoter track and color is changed to dark orange. We previously did the color change in Illustrator and combined the two tracks vertically using a separate plot that used stacking=dense. We exported to postscript from this point and embedded the fonts

extra_prom_gr_tm=data.table::fread(system.file(
"extdata","more_promoters_topmed.csv",
  package = "GS3DKViz")
) %>% janitor::clean_names( ) %>%  GRanges()
extraPromoterTrackTm <- AnnotationTrack(extra_prom_gr_tm, genome="hg38", name=" Promoters",
                             id=mcols(extra_prom_gr_tm)$gene_id,  featureAnnotation="id")
displayPars(extraPromoterTrackTm)=displayPars(promoterTrack)
displayPars(extraPromoterTrackTm) <- list(fill = "darkorange")
#postscript("fev1-fvc-psma4-dense-extra-prom.eps",width=27,height=9,family="sans")
plotTracks(list(itrack,gtrack,interaction_track,promoterTrack,extraPromoterTrackTm,enhancerTrack),
           chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(gene_id==geneSymbol) %>% dplyr::pull(seqnames1)),
           from = (min(bounds))-1.5e4,
           to = (max(bounds))+1.5e4,
           type = c("b"),
           showSampleNames = TRUE, 
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=12,collapse=T,min.width=5,mergeGroups=T
           ,stacking="squish" #comment out to restore stacking
           ,main="Promoter–enhancer interactions, FEV1/FVC",
           background.title = "black"
)
#dev.off()

Now we'll do FEV1 and FVC by changing the bounds and the chromosome to those matching detected FEV1 and FVC promoters.

interaction_track <- InteractionTrack(gs3dk_gint, name = "Interaction", chromosome = as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1") %>% dplyr::pull(seqnames1)))
displayPars(interaction_track)=list(col.interactions="black")
bounds=c(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1") %>% dplyr::pull(start1) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1") %>% dplyr::pull(end1) %>% max(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1") %>% dplyr::pull(start2) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1") %>% dplyr::pull(end2) %>% max())

plotTracks(list(itrack,gtrack,interaction_track,promoterTrack,enhancerTrack),
           chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1") %>% dplyr::pull(seqnames1)),
           from = (min(bounds))-1.5e4,
           to = (max(bounds))+1.5e4,
           type = c("b"),
           showSampleNames = TRUE, 
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=12,collapse=T,min.width=5,mergeGroups=T
           ,stacking="dense" #comment out to restore stacking
           ,main="Promoter–enhancer interactions, FEV1",
           background.title = "black"
)

interaction_track <- InteractionTrack(gs3dk_gint, name = "Interaction", chromosome = as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FVC") %>% dplyr::pull(seqnames1)))
displayPars(interaction_track)=list(col.interactions="black")
bounds=c(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FVC") %>% dplyr::pull(start1) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FVC") %>% dplyr::pull(end1) %>% max(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FVC") %>% dplyr::pull(start2) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FVC") %>% dplyr::pull(end2) %>% max())
plotTracks(list(itrack,gtrack,interaction_track,promoterTrack,enhancerTrack),
           chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FVC") %>% dplyr::pull(seqnames1)),
           from = (min(bounds))-1.5e4,
           to = (max(bounds))+1.5e4,
           type = c("b"),
           showSampleNames = TRUE, 
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=12,collapse=T,min.width=5,mergeGroups=T
           ,stacking="dense" #comment out to restore stacking
           ,main="Promoter–enhancer interactions, FVC",
           background.title = "black"
)

Another example using alternative colors for the interactions:

annotateInteractions(GIObject = gs3dk_gint,annotations =  list(promoter=promoter_gr, enhancer=enhancer_gr),id.col="gene_id")

knockoff_genes=mcols(gs3dk_gint) %>% tibble::as_tibble() %>% dplyr::filter(identified_by_knockoff_or_not=="Yes") %>% dplyr::pull(gene_id)

gs3dk_gint@regions@elementMetadata$node.class=(mcols(regions(gs3dk_gint)) %>% tibble::as_tibble() %>% dplyr::mutate(knockoff=ifelse((promoter.id %in% knockoff_genes | enhancer.id %in% knockoff_genes),"knockoff_detected","knockoff_removed"),node.class=knockoff) %>% dplyr::pull(knockoff))

interaction_track <- InteractionTrack(gs3dk_gint, name = "Interaction", chromosome = as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(seqnames1)))

displayPars(interaction_track)=list(col.interactions="black",col.interaction.types=c('knockoff_detected-knockoff_detected'='blue', 'knockoff_detected-knockoff_removed'='blue','knockoff_removed-knockoff_removed'='black'))
bounds=c(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(start1) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(end1) %>% max(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(start2) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(end2) %>% max())
plotTracks(list(itrack,gtrack,interaction_track,promoterTrack,enhancerTrack),
           chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="FEV1/FVC") %>% dplyr::pull(seqnames1)),
           from = (min(bounds))-1.5e4,
           to = (max(bounds))+1.5e4,
           type = c("b"),
           showSampleNames = TRUE, 
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=12,collapse=T,min.width=5,mergeGroups=T
           ,stacking="squish" #comment out to restore stacking
           ,main="Promoter–enhancer interactions, FEV1/FVC",
           background.title = "black"
)

Lastly, we'll do a details plot, which requires a details function and selection function. It's a bit complex, but the code is below.

annotateInteractions(GIObject = gs3dk_gint,annotations =  list(promoter=promoter_gr, enhancer=enhancer_gr),id.col="gene_id")
knockoff_genes=mcols(gs3dk_gint) %>% tibble::as_tibble() %>% dplyr::filter(identified_by_knockoff_or_not=="Yes") %>% dplyr::pull(gene_id)
gs3dk_gint@regions@elementMetadata$node.class=(mcols(regions(gs3dk_gint)) %>% tibble::as_tibble() %>% dplyr::mutate(knockoff=ifelse((promoter.id %in% knockoff_genes | enhancer.id %in% knockoff_genes),"knockoff_detected","knockoff_removed"),node.class=knockoff) %>% dplyr::pull(knockoff))
interaction_track <- InteractionTrack(gs3dk_gint, name = "Interaction", chromosome = as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(seqnames1)))
displayPars(interaction_track)=list(col.interactions="black",col.interaction.types=c('knockoff_detected-knockoff_detected'='blue', 'knockoff_detected-knockoff_removed'='blue','knockoff_removed-knockoff_removed'='black'))
bounds=c(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(start1) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(end1) %>% max(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(start2) %>% min(),
gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(end2) %>% max())
selFun <- function(identifier, start, end, track, GdObject, ...){
    gcount <- table(group(GdObject))
    pxRange <- Gviz:::.pxResolution(min.width = 50, coord = "x")
    return((end - start) < pxRange && gcount[identifier] == 1 && (distanceToNearest(GdObject[group(GdObject) == identifier]@range,GdObject[group(GdObject) != identifier]@range)@elementMetadata$distance)<50
           )
}
detFun <- function(identifier, GdObject.original, ...){
  plotTracks(list(GenomeAxisTrack(scale = 0.3, size = 0.2, cex = 0.7), 
        GdObject.original[group(GdObject.original) == identifier]),         add = TRUE, showTitle = FALSE)
}
deTrackEnh2 <- AnnotationTrack(name = "Enhancers",enhancer_gr, fun = detFun, 
                            selectFun = selFun,
                            groupDetails = TRUE, details.size = 0.5, 
                            detailsConnector.cex = 0.5, 
                            detailsConnector.lty = "dotted",
                            shape = c("smallArrow", "arrow"), 
                            groupAnnotation = "group",
                            id=mcols(enhancer_gr)$gene_id, 
                            group=subjectHits(findOverlaps(enhancer_gr, reduce(enhancer_gr))),featureAnnotation="id",stacking="hide",
                      )
displayPars(deTrackEnh2) <- list(fill = "mediumpurple1", col = NA, 
                                   fontcolor.feature = "black", fontsize=10,
                                   just.group="below",rotation.item=90,
                                   collapse=T,mergeGroups=T,showOverplotting=T,groupAnnotation="group",    group=mcols(gs3dk_gint)$enhancer_type,feature=mcols(gs3dk_gint)$enhancer_type,
                              groupDetails=T,detailsConnector.lty="solid",
                              detailsConnector.col="black",
                              detailsBorder.col="black",
                              detailsBorder.lty="solid",showId=F
                             ,min.distance=5,min.width=5,title="Enhancers") 
extra_prom_gr=data.table::fread(system.file("extdata","more_promoters.csv",
  package = "GS3DKViz")) %>% janitor::clean_names( ) %>%  dplyr::filter(!(gene_id %in%  janitor::clean_names(tibble::as_tibble(mcols(gs3dk_gint)))$gene_id)) %>% GRanges()
itrack <- IdeogramTrack(genome = "hg38", chromosome =unique(as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(seqnames1)) ))
combined_prom_gr=c(promoter_gr,extra_prom_gr)
CombinedPromoterTrack <- AnnotationTrack(combined_prom_gr, genome="hg38", name="Promoters",
                             id=mcols(combined_prom_gr)$gene_id,  featureAnnotation="id",groupAnnotation="feature")
displayPars(CombinedPromoterTrack)=list(fill = "olivedrab1", col = NA, 
                                   fontcolor.feature = "black", fontsize=8,                 just.group="below",rotation=90,rotation.group=90,rotation.item=90,
                                   min.width=1,min.distance=5)
feature(CombinedPromoterTrack)=c(rep("connected",length(promoter_gr)),rep("unconnected",length(extra_prom_gr)))
plotTracks(list(itrack,gtrack,interaction_track,CombinedPromoterTrack,deTrackEnh2),chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(seqnames1)),
           from = (min(bounds)-1.5e4),
           to = (max(bounds)+1.5e4),
           type = c("b"),
           showSampleNames = TRUE,
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=7,collapse=T,
           mergeGroups=T
           ,stacking="dense"
           ,main="Promoter–enhancer interactions, AD",
           background.title = "black",
           unconnected="orange",connected="green"
)
extra_prom_gr=data.table::fread(system.file("extdata","more_promoters.csv",
  package = "GS3DKViz")) %>% janitor::clean_names( ) %>%  dplyr::filter(!(gene_id %in%  janitor::clean_names(tibble::as_tibble(mcols(gs3dk_gint)))$gene_id)) %>% GRanges()

To export it to EPS, for further edits, do the following:

filename=gsub(":","_",gsub(" ","_",paste(date(),"AD-dense-details-extra_color_arcs.eps")))
postscript(filename,width=50,height=500,family="sans")
plotTracks(list(itrack,gtrack,interaction_track,promoterTrack,extraPromoterTrack,deTrackEnh2), 
           chromosome=as.character(gs3dk_gint %>% as.data.frame() %>% janitor::clean_names() %>% dplyr::filter(trait=="AD") %>% dplyr::pull(seqnames1)),
           from = (min(bounds)-1.5e4),
           to = (max(bounds)+1.5e4),
           type = c("b"),
           showSampleNames = TRUE,
           cex.sampleNames = 0.6,
           cex.main=3,
           fontsize=18,fontsize.item=7,collapse=T,
           #min.width=1,min.distance=5,
           mergeGroups=T
           ,stacking="dense" #comment out to restore stacking
           ,main="Promoter–enhancer interactions, AD",
           background.title = "black"
)
dev.off()
extrafont::embed_fonts(filename,outfile=paste0(filename,"-embed.eps"))

#Appendix:
#cross checks

highlighted_genes=c("QPCTL", "OPA3", "ZNF45", "ZNF155", "ZNF230","ZNF223","ZNF284","ZNF224","ZNF225","ZNF234","ZNF180","IGSF23","BCAM","NECTIN2","RELB")
setdiff(mcols(gs3dk_gint[is.ko.detected(gs3dk_gint) & gs3dk_gint@elementMetadata$seqnames1=="chr19"])$gene_id,highlighted_genes)
intersect(mcols(gs3dk_gint[is.ko.detected(gs3dk_gint) & gs3dk_gint@elementMetadata$seqnames1=="chr19"])$gene_id,highlighted_genes)
setdiff(gs3dk_gint[gs3dk_gint@elementMetadata$identified_by_knockoff_or_not=="Yes" & gs3dk_gint@elementMetadata$seqnames1=="chr19"]$gene_id,highlighted_genes)
intersect(gs3dk_gint[gs3dk_gint@elementMetadata$identified_by_knockoff_or_not=="Yes" & gs3dk_gint@elementMetadata$seqnames1=="chr19"]$gene_id,highlighted_genes)
length(intersect(gs3dk_gint[gs3dk_gint@elementMetadata$identified_by_knockoff_or_not=="Yes" & gs3dk_gint@elementMetadata$seqnames1=="chr19"]$gene_id,highlighted_genes))
length(highlighted_genes)
old=GS3DKViz::readGint("../inst/extdata/signficant_genes_genoscan_3D_knock.csv")
new=GS3DKViz::readGint("../inst/extdata/signficant_genes_genoscan_3D_knock_12-20.csv")
table(old==new)

.collapseAnnotation=function (grange, minXDist, elements, GdObject, offset = 0) 
{

    needsRestacking <- TRUE
    annoSplit <- merged <- NULL
    anno <- as.data.frame(grange)
    for (i in colnames(anno)) if (is.factor(anno[, i])) 
        anno[, i] <- as.character(anno[, i])
    cols <- c("strand", "density", "gdensity", 
        "feature", "id", "start", "end", 
        if (is(GdObject, "GeneRegionTrack")) c("gene", 
            "exon", "transcript", "symbol", 
            "rank") else "group")
    missing <- which(!cols %in% colnames(anno))
    for (i in missing) anno[, cols[missing]] <- if (cols[i] == 
        "density") 
        1
    else NA
    rRed <- if (length(grange) > 1) 
        reduce(grange, min.gapwidth = minXDist, with.revmap = TRUE)
    else grange
    if (length(rRed) < length(grange)) {
        needsRestacking <- TRUE
        mapping <- rep(seq_along(rRed$revmap), elementNROWS(rRed$revmap))
        identical <- mapping %in% which(table(mapping) == 1)
        newVals <- anno[identical, cols]
        if (nrow(newVals)) {
            newVals$seqnames <- elements[as.character(anno[identical, 
                "seqnames"])] == 1
            newVals$gdensity <- ifelse(elements[as.character(anno[identical, 
                "seqnames"])] == 1, 1, NA)
        }
        grange <- grange[!identical]
        rRed <- rRed[-(mapping[identical])]
        index <- mapping[!identical]
        annoSplit <- split(anno[!identical, ], index)
        cid <- function(j) sprintf("[Cluster_%i]  ", j + 
            offset)
        newVals <- rbind(newVals, as.data.frame(t(sapply(seq_along(annoSplit), 
            function(i) {
                x <- annoSplit[[i]]
                if (is(GdObject, "GeneRegionTrack")) {
                  c(strand = ifelse(length(unique(x[, "strand"])) == 
                    1, as.character(x[1, "strand"]), "*"), 
                    density = sum(as.integer(x[, "density"])), 
                    gdensity = ifelse(is.na(head(x[, "gdensity"], 
                      1)), 1, sum(as.integer(x[, "gdensity"]))), 
                    feature = ifelse(length(unique(x[, "feature"])) == 
                      1, as.character(x[1, "feature"]), 
                      "composite"), id = ifelse(length(unique(x[, 
                      "id"])) == 1, as.character(x[1, "id"]), 
                      cid(i)), start = min(x[, "start"]), 
                    end = max(x[, "end"]), gene = ifelse(length(unique(x[, 
                      "gene"])) == 1, as.character(x[1, 
                      "gene"]), cid(i)), exon = ifelse(length(unique(x[, 
                      "exon"])) == 1, as.character(x[1, 
                      "exon"]), cid(i)), transcript = ifelse(length(unique(x[, 
                      "transcript"])) == 1, as.character(x[1, 
                      "transcript"]), cid(i)), symbol = ifelse(length(unique(x[, 
                      "symbol"])) == 1, as.character(x[1, 
                      "symbol"]), cid(i)), rank = min(as.integer(x[, 
                      "rank"])), seqnames = as.vector(nrow(x) == 
                      elements[x[1, "seqnames"]]))
                }
                else {
                  c(strand = ifelse(length(unique(x[, "strand"])) == 
                    1, as.character(x[1, "strand"]), "*"), 
                    density = sum(as.integer(x[, "density"])), 
                    gdensity = ifelse(is.na(head(x[, "gdensity"], 
                      1)), 1, sum(as.integer(x[, "gdensity"]))), 
                    feature = ifelse(length(unique(x[, "feature"])) == 
                      1, as.character(x[1, "feature"]), 
                      "composite"), id = ifelse(length(unique(x[, 
                      "id"])) == 1, as.character(x[1, "id"]), 
                      cid(i)), start = min(x[, "start"]), 
                    end = max(x[, "end"]), group = ifelse(length(unique(x[, 
                      "group"])) == 1, as.character(x[1, 
                      "group"]), cid(i)), seqnames = as.vector(nrow(x) == 
                      elements[x[1, "seqnames"]]))
                }
            })), stringsAsFactors = FALSE))
        merged <- as.logical(newVals$seqnames)
        grange <- GRanges(seqnames = chromosome(GdObject), strand = newVals[, 
            "strand"], ranges = IRanges(start = as.integer(newVals[, 
            "start"]), end = as.integer(newVals[, "end"])))
        cnMatch <- match(c(colnames(values(GdObject)), "gdensity"), 
            colnames(newVals))
        mcols(grange) <- if (any(is.na(cnMatch))) 
            newVals[, setdiff(colnames(newVals), c("strand", 
                "start", "end", "seqnames"))]
        else newVals[, cnMatch]
    }
    else {
        grange2 <- GRanges(seqnames = chromosome(GdObject), strand = strand(grange), 
            ranges = ranges(grange))
        mcols(grange2) <- mcols(grange)
        grange <- grange2
    }
    browser()
    return(list(range = grange, needsRestacking = needsRestacking, 
        split = annoSplit, merged = merged, offset = length(annoSplit)))
}

jamesdalg/GS3DKViz documentation built on Jan. 1, 2021, 4:27 a.m.