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])
