Plot custom gene density with R

Jit — Thu, 08 Mar 2018 16:30:34 -0600
library(karyoploteR)

pp <- getDefaultPlotParams(plot.type=2)
pp$data1outmargin <- 100
pp$data2outmargin <- 100
pp$topmargin <- 450

gff.file <- "http://plasmodb.org/common/downloads/Current_Release/PvivaxP01/gff/data/PlasmoDB-35_PvivaxP01.gff"
data.points.colors <- c("#FFBD07", "#00A6ED",  "#FF3366", "#8EA604", "#C200FB")

header.lines <- readLines(gff.file, n = 30)
## Error in file(con, "r"): cannot open the connection to 'http://plasmodb.org/common/downloads/Current_Release/PvivaxP01/gff/data/PlasmoDB-33_PvivaxP01.gff'
#The lines with the standard chromosomes start with "##sequence-region PvP01".
#Select them.
ll <- header.lines[grepl(header.lines, pattern = "##sequence-region PvP01")]
## Error in eval(expr, envir, enclos): object 'header.lines' not found
#split them by space, and create a data.frame
gg <- data.frame(do.call(rbind, strsplit(ll, split = " ")))
## Error in strsplit(ll, split = " "): object 'll' not found
gg[,3] <- as.numeric(as.character(gg[,3]))
## Error in eval(expr, envir, enclos): object 'gg' not found
gg[,4] <- as.numeric(as.character(gg[,4]))
## Error in eval(expr, envir, enclos): object 'gg' not found
#and create a GRanges with the information
PvP01.genome <- toGRanges(gg[,c(2,3,4)])
## Error in is(A, "GRanges"): object 'gg' not found
PvP01.genome
## Error in eval(expr, envir, enclos): object 'PvP01.genome' not found

#kp <- plotKaryotype(genome=PvP01.genome)
#kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL, plot.type=2, plot.params = pp)
kp <- plotKaryotype(genome=PvP01.genome, plot.type=2, plot.params = pp)

#kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL)

#kpAddCytobandsAsLine(kp)

features <- import(gff.file)

table(features$type)

genes <- features[features$type=="gene"]

#kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL)

#kpAddCytobandsAsLine(kp)

#kpPlotRegions(kp, data=genes)

#kp <- plotKaryotype(genome=PvP01.genome, ideogram.plotter = NULL, plot.type=2)
kpAddMainTitle(kp, "Plasmodium Vivax - PvP01 with genes", cex=2)


kpPlotRegions(kp, data=genes[strand(genes)=="+"], avoid.overlapping = FALSE, col="deepskyblue")

kpPlotRegions(kp, data=genes[strand(genes)=="-"], avoid.overlapping = FALSE, col="gold", data.panel=2)

kpAddLabels(kp, "strand +", cex=0.8, col="#888888")

kpAddLabels(kp, "strand -", data.panel=2, cex=0.8, col="#888888")



#plot all
pp <- getDefaultPlotParams(plot.type = 4)
pp$data1inmargin <- 0
pp$bottommargin <- 20

kp <- plotKaryotype(genome=PvP01.genome, plot.type=4, ideogram.plotter = NULL,
                    labels.plotter = NULL, plot.params = pp,
                    main="Gene Density")
kpAddCytobandsAsLine(kp)
kpAddChromosomeNames(kp, srt=45)
kpPlotDensity(kp, genes, window.size = 10e2, col="#ddaacc")



#plot all
pp <- getDefaultPlotParams(plot.type = 2)
pp$data1inmargin <- 0
pp$bottommargin <- 20

kp <- plotKaryotype(genome=PvP01.genome, plot.params = pp)
#kp <- kpPlotDensity(kp, genes)

kpAddChromosomeNames(kp, srt=45)
kpPlotDensity(kp, genes, col="#ddaacc", data.panel = 1)

kpAbline(kp, h=0.4, data.panel = 2, r0=0.2, r1=0, col=data.points.colors[3])
BOL: Plot custom gene density with R

Plot custom gene density with R
