title: "BUSCO plots"
author: "zagor"
date: " 09 10 2019"
output: html_document
editor_options:
chunk_output_type: console
knitr::opts_chunk$set(echo = TRUE, fig.width=12, fig.height=9)
if (!require("data.table")) install.packages("data.table")
library(data.table)
if (!require("RColorBrewer")) install.packages("RColorBrewer")
library(RColorBrewer)
if (!require("grr")) install.packages("grr")
library(grr)
if (!require("plyr")) install.packages("plyr")
library(plyr)
if (!require("ggplot2")) install.packages("ggplot2")
library(ggplot2)
if (!require("stringr")) install.packages("stringr")
library(stringr)
ca <- brewer.pal(4, "BrBG")
cb <- brewer.pal(4, "Dark2")
par(mar=c(6, 4, 4, 8) + 1)
cnt = 1
dir.create("../output/BUSCOimg")
myTable1 = read.table(file = "../output/BUSCO_stPanTr_embryophyta_odb9.tsv",
header = FALSE,
sep = "\t",
quote = "",
stringsAsFactors = FALSE,
fill=TRUE)
myTable2 = read.table(file = "../output/BUSCO_stPanTr_solanaceae_odb10.tsv",
header = FALSE,
sep = "\t",
quote = "",
stringsAsFactors = FALSE,
fill=TRUE)
myTable1 = data.table(myTable1)
myTable2 = data.table(myTable2)
myTable1$V1[myTable1$V1 == ""] = NA
myTable2$V1[myTable2$V1 == ""] = NA
goodIdx <- !is.na(myTable1$V1)
goodVals <- c(NA, myTable1$V1[goodIdx])
fillIdx <- cumsum(goodIdx)+1
myTable1$V1 = goodVals[fillIdx]
myTable1 = myTable1[-which(is.na(myTable1$V2)),]
goodIdx <- !is.na(myTable2$V1)
goodVals <- c(NA, myTable2$V1[goodIdx])
fillIdx <- cumsum(goodIdx)+1
myTable2$V1 = goodVals[fillIdx]
myTable2 = myTable2[-which(is.na(myTable2$V2)),]
data = myTable1
colnames(data) <-c("DataSet", "BUSCOs", "BUSCOclass")
data$BUSCOclass = ordered(data$BUSCOclass, levels = unique(data$BUSCOclass))
data$DataSet = ordered(data$DataSet, levels = c("Solanum_tuberosum-ITAG_DM_v1_pep",
"Solanum_tuberosum-PGSC_DM_v34_pep",
"Solanum_tuberosum_PGSC_DM_v4.04_sequences_chrUn_pep",
"Solanum_tuberosum_PGSC_DM_v4.04_sequences_merged",
"StPGSC4.04n_seq_4_PGSC-p-rep_ITAG-pep_gff-2019-04-23",
"stPanTr.aa",
"stPanTr_rep.aa",
"stPanTr_alt.aa"))
data$BUSCOs = as.numeric(data$BUSCOs)
m = sum(data$BUSCOs[1:4])
data$perc = format(as.numeric(data$BUSCOs)/m*100, digits = 0, nsmall = 1)
n = nrow(data)
yy = NULL
for (i in seq(1,n,4)) {
yy = c(yy,(rev(c(i:(i+3)))))
}
data$BUSCOs_y = data$BUSCOs[yy]
data$BUSCOclass_y = data$BUSCOclass[yy]
data <- ddply(data, "DataSet",
transform, label_ypos=cumsum(BUSCOs_y))
data$perc.rev = data$perc[yy]
data$BUSCOclass
data$DataSet
mycaption = paste0("DB: embryophyta_odb9, size: ", m)
ggplot(data = data, aes(x = DataSet, y = BUSCOs, fill = BUSCOclass), font.axis=2) +
theme_bw() +
scale_fill_manual(values=rev(ca)) +
xlab("") +
ylab("BUSCO count") +
scale_x_discrete(labels= c("ITAG_DM_v1","PGSC_DM_v34","PGSC_v4.04_chrUn","DM_v34+v4.04_chrUn",
"ITAG/PGSC merged", "stPanTr all", "stPanTr rep", "stPanTr alt")) +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank()) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=-0.3, size=10)) +
theme(legend.position = "right") +
guides(fill=guide_legend(title="BUSCO classes")) +
labs(caption = mycaption) # geom_text(aes(y=label_ypos, label=perc), vjust=1.0, color="white", size=3.5)
myplot = ggplot(data = data, aes(x = DataSet, y = as.numeric(perc), fill = BUSCOclass), font.axis=2) +
theme_bw() +
scale_fill_manual(values=rev(ca)) +
xlab("") +
ylab("BUSCO %") +
scale_x_discrete(labels= c("ITAG_DM_v1","PGSC_DM_v34","PGSC_v4.04_chrUn","DM_v34+v4.04_chrUn",
"ITAG/PGSC merged", "stPanTr all", "stPanTr rep", "stPanTr alt")) +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank()) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=-0.3, size=10)) +
theme(legend.position = "right") +
guides(fill=guide_legend(title="BUSCOclass")) +
labs(caption = mycaption)
pdf(paste0("../output/BUSCOimg/myPlot", str_pad(cnt, 3, pad = "0"), ".pdf"), width=9, height=8)
print(myplot)
dev.off()
cnt = cnt + 1
data = myTable2
colnames(data) <-c("DataSet", "BUSCOs", "BUSCOclass")
data$BUSCOclass = ordered(data$BUSCOclass, levels = unique(data$BUSCOclass))
data$DataSet = ordered(data$DataSet, levels = c("Solanum_tuberosum-ITAG_DM_v1_pep",
"Solanum_tuberosum-PGSC_DM_v34_pep",
"Solanum_tuberosum_PGSC_DM_v4.04_sequences_chrUn_pep",
"Solanum_tuberosum_PGSC_DM_v4.04_sequences_merged",
"StPGSC4.04n_seq_4_PGSC-p-rep_ITAG-pep_gff-2019-04-23",
"stPanTr.aa",
"stPanTr_rep.aa",
"stPanTr_alt.aa"))
data$BUSCOs = as.numeric(data$BUSCOs)
m = sum(data$BUSCOs[1:4])
data$perc = format(as.numeric(data$BUSCOs)/m*100, digits = 0, nsmall = 1)
n = nrow(data)
yy = NULL
for (i in seq(1,n,4)) {
yy = c(yy,(rev(c(i:(i+3)))))
}
data$BUSCOs_y = data$BUSCOs[yy]
data$BUSCOclass_y = data$BUSCOclass[yy]
data <- ddply(data, "DataSet",
transform, label_ypos=cumsum(BUSCOs_y))
data$perc.rev = data$perc[yy]
data$BUSCOclass
data$DataSet
mycaption = paste0("DB: solanaceae_odb10, size: ", m)
ggplot(data = data, aes(x = DataSet, y = BUSCOs, fill = BUSCOclass), font.axis=2) +
theme_bw() +
scale_fill_manual(values=rev(ca)) +
xlab("") +
ylab("BUSCO count") +
scale_x_discrete(labels= c("ITAG_DM_v1","PGSC_DM_v34","PGSC_v4.04_chrUn","DM_v34+v4.04_chrUn",
"ITAG/PGSC merged", "stPanTr all", "stPanTr rep", "stPanTr alt")) +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank()) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=-0.3, size=10)) +
theme(legend.position = "right") +
guides(fill=guide_legend(title="BUSCO classes")) +
labs(caption = mycaption) # geom_text(aes(y=label_ypos, label=perc), vjust=1.0, color="white", size=3.5)
myplot = ggplot(data = data, aes(x = DataSet, y = as.numeric(perc), fill = BUSCOclass), font.axis=2) +
theme_bw() +
scale_fill_manual(values=rev(ca)) +
xlab("") +
ylab("BUSCO %") +
scale_x_discrete(labels= c("ITAG_DM_v1","PGSC_DM_v34","PGSC_v4.04_chrUn","DM_v34+v4.04_chrUn",
"ITAG/PGSC merged", "stPanTr all", "stPanTr rep", "stPanTr alt")) +
theme(axis.line = element_line(colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.border = element_blank(),
panel.background = element_blank()) +
geom_bar(stat = "identity") +
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust=-0.3, size=10)) +
theme(legend.position = "right") +
guides(fill=guide_legend(title="BUSCOclass")) +
labs(caption = mycaption)
pdf(paste0("../output/BUSCOimg/myPlot", str_pad(cnt, 3, pad = "0"), ".pdf"), width=9, height=8)
print(myplot)
dev.off()
cnt = cnt + 1