--- a
+++ b/ATAC/Functions/user.CoveragePlot.R
@@ -0,0 +1,56 @@
+
+user.CoveragePlot <- function(scATAC.object, interest.genes, interested.conns, ccans, gene_model, idents = NULL, extend.upstream = 1000, extend.downstream = 1000, coaccess_cutoff = 0.2, signac = T, cicero = T, heights = NULL){
+    require(GenomicRanges)
+    res <- sapply(interest.genes, function(x){
+        cat("Plot coverage: ", x, "\n")
+        idx1 <- which(interested.conns$peak1_nearestGene == x)
+        idx2 <- which(interested.conns$peak2_nearestGene == x)
+        idx <- unique(idx1, idx2)
+        conns <- data.frame(Peak1 = gsub("_", "-", interested.conns$Peak1[idx]),
+                            Peak2 = gsub("_", "-", interested.conns$Peak2[idx]),
+                            coaccess = as.numeric(interested.conns$coaccess[idx]))
+        region.highlights <- c(conns$Peak1, conns$Peak2)
+        region.highlights <- StringToGRanges(region.highlights)
+        start.loc <- min(start(region.highlights))
+        end.loc <- max(end(region.highlights))
+        region.highlights$color <- rep(c("#ffa500", "#c71585"), each = nrow(conns)) #第一个颜色标记promoter
+
+        links <- ConnectionsToLinks(conns = conns, ccans = ccans)
+        Links(scATAC.object) <- links
+
+        plot.regions <- paste0(as.character(seqnames(region.highlights))[1], "-", start.loc, "-", end.loc)
+        if(signac){
+            if(is.null(idents)){
+                p <- CoveragePlot(scATAC.object, 
+                    region = plot.regions,
+                    region.highlight = region.highlights, 
+                    extend.upstream = extend.upstream, extend.downstream = extend.downstream,
+                    heights = heights)
+            }else{
+                p <- CoveragePlot(scATAC.object, 
+                    region = plot.regions,
+                    idents = idents,
+                    region.highlight = region.highlights, 
+                    extend.upstream = extend.upstream, extend.downstream = extend.downstream,
+                    heights = heights)                
+            }
+            print(p)
+        }
+
+        #cicero plot
+        if(cicero){
+            require(cicero)
+            require(monocle3)
+            p <- plot_connections(conns,
+                as.character(seqnames(region.highlights))[1], start.loc, end.loc,
+                gene_model = gene_model,
+                viewpoint = interested.conns$Peak1[idx[1]],
+                coaccess_cutoff = coaccess_cutoff, 
+                connection_width = .5, include_axis_track = F,
+                collapseTranscripts = "longest")
+            print(p)
+        }
+        return(conns)
+    })
+    return(res)
+}
\ No newline at end of file