#########################################
###### Author:Katharine Lee #############
############## 2017.07 ##################
#########################################
reads_compare = function(
input,
output,
outputFormat = "png",
group = NULL,
# themes = "pale","grey70",
# names = NULL,
thredhold = NULL,
# bar_colours = NULL,
width = 16,
height = 8
){
# library(ggplot2)
# library(scales)
# library(gtable)
# library(dplyr)
#
# input = system.file("example", "total_otu_stat.xls", package = "metaplot")
# input = "~/mix004.xls"
t = read.table(input,header = T,sep = "\t",check.names = F)[,-1]
if(nrow(t)==2){
df = data.frame(Samples = colnames(t),otu_num = t(t[1,]),count_num = t(t[2,]))[-1:-5,]
colnames(df) = c("Samples","OTUs","Tags")
df$OTUs = as.numeric(df$OTUs)
df$Tags = as.numeric(df$Tags)
df$name = as.character(1)
if(!is.null(group)) {
group <- unlist(strsplit(group,","))
row_names = dplyr::arrange(data.frame(sample = df$Samples, group = group),group)$sample
df$Samples = factor(df$Samples,levels = row_names)
}
AXIS1_MIN = min(df$Tags)
AXIS1_MAX = max(df$Tags)
AXIS2_MIN = min(df$OTUs)
AXIS2_MAX = max(df$OTUs)
scale_to_value1 <- function(values) rescale(values, to = c(AXIS1_MIN, AXIS1_MAX))
scale_to_value2 <- function(values) rescale(values, to = c(AXIS2_MIN, AXIS2_MAX))
p = ggplot() +
geom_bar(aes(x = Samples,y = Tags,fill =name) ,position = "dodge",data = df ,stat="identity",size=1) +
geom_line(aes(x = c(1:nrow(df)),y = scale_to_value1(OTUs),colour = name),data = df,size=1) +
geom_point(aes(x = c(1:nrow(df)),y = scale_to_value1(OTUs)),data = df,size=3,colour = 'olivedrab3')+
scale_y_continuous(
breaks=pretty_breaks(n=7),
# limits=c(AXIS1_MIN, AXIS1_MAX),
# expand = c(0,0),
sec.axis = sec_axis( ~ scale_to_value2(.),
breaks=pretty_breaks(n=7),
name = "Number of OTUs"))+
scale_colour_manual(name = NULL,
label = "OTUs",
values = "olivedrab3")+
scale_fill_manual(name = NULL,
label = "Tags",
values = "grey85") +
xlab('Samples')+ylab('Number of Tags')+
# theme(plot.title = element_text(lineheight=3, face="bold", color="black",size=24)) +
theme(legend.text=element_text(size=14),
legend.title=element_text(size=14))+
theme_bw()+
# theme(axis.text.x=element_text(angle=90))+
theme(legend.position = "bottom")+
theme(axis.text.x=element_text(hjust=.5, vjust=.5,angle=90,colour="grey20",size=12),
axis.text.y=element_text(hjust=.5, vjust=.5,colour="grey20",size=12),
axis.title.x=element_text(colour="grey20",size=12),
axis.title.y=element_text(colour="grey20",size=12),
legend.text = element_text(size = 12))
if(!is.null(thredhold)) p = p + geom_hline(yintercept = thredhold,colour="tomato", linetype="dashed")
# if(!is.null(colour)) {col <- unlist(strsplit(colour,","))} else col = c("grey50","grey70")
}
if(nrow(t)==3){
reads = data.frame(Samples = colnames(t),count_num = t(t[2,]))[-1:-5,]
colnames(reads) = c("Samples","num")
reads$Samples = unlist(strsplit(as.character(reads$Samples),":"))
reads$type = "clean"
raw = data.frame(Samples = colnames(t),raw = t(t[3,]))[-1:-5,]
colnames(raw) = c("Samples","num")
raw = raw[which(raw$Samples %in% reads$Samples),]
raw$type = "raw"
if(!is.null(group)) {
group <- unlist(strsplit(group,","))
row_names = dplyr::arrange(data.frame(sample = raw$Samples, group = group),group)$sample
raw$Samples = factor(raw$Samples,levels = row_names)
reads$Samples = factor(reads$Samples,levels = row_names)
}
df = data.frame(rbind(raw,reads))
df$type = factor(df$type,levels = c("raw","clean"))
otus = data.frame(Samples = colnames(t),otu_num = t(t[1,]))[-1:-5,]
colnames(otus) = c("Samples","otu_num")
otus$Samples = unlist(strsplit(as.character(otus$Samples),":"))
otus$otu_num = as.numeric(otus$otu_num)
rownames(otus) = otus$Samples
otus = otus[which(otus$Samples == unique(df$Samples)),]
otus$name = as.character(1)
otus = otus[levels(df$Samples),]
AXIS1_MIN = min(df$num)
AXIS1_MAX = max(df$num)
AXIS2_MIN = min(otus$otu_num)
AXIS2_MAX = max(otus$otu_num)
scale_to_value1 <- function(values) rescale(values, to = c(AXIS1_MIN, AXIS1_MAX))
scale_to_value2 <- function(values) rescale(values, to = c(AXIS2_MIN, AXIS2_MAX))
p = ggplot() +
geom_bar(aes(x = Samples,y = num,fill = type) ,position = "dodge",data = df ,stat="identity",size=1) +
geom_line(aes(x = c(1:nrow(otus)),y = scale_to_value1(otu_num),colour = name),data = otus,size=1) +
geom_point(aes(x = c(1:nrow(otus)),y = scale_to_value1(otu_num)),data = otus,size=3,colour = 'olivedrab3')+
scale_y_continuous(
breaks=pretty_breaks(n=7),
# limits=c(AXIS1_MIN, AXIS1_MAX),
# expand = c(0,0),
sec.axis = sec_axis( ~ scale_to_value2(.),
breaks=pretty_breaks(n=7),
name = "Number of OTUs"))+
scale_colour_manual(name = NULL,
label = "OTUs",
values = "olivedrab3")+
scale_fill_manual(name = NULL,
label = c("Valid Tags","Non-chimera Tags"),
values = c("grey65","grey85")) +
# theme(axis.text.x=element_text(angle=45,colour="grey20",face="bold",size=12),
# axis.text.y=element_text(colour="grey20",face="bold",hjust=1,vjust=0.8,size=15),
# axis.title.x=element_text(colour="grey20",face="bold",size=16),
# axis.title.y=element_text(colour="grey20",face="bold",size=16)) +
xlab('Samples')+ylab('Number of Tags')+
# theme(plot.title = element_text(lineheight=3, face="bold", color="black",size=24)) +
theme(legend.text=element_text(size=14),
legend.title=element_text(size=14))+
theme_bw()+
# theme(axis.text.x=element_text(angle=90))+
theme(legend.position = "bottom")+
theme(axis.text.x=element_text(hjust=.5, vjust=.5,angle=90,colour="grey20",size=12),
axis.text.y=element_text(colour="grey20",size=12),
axis.title.x=element_text(colour="grey20",size=12),
axis.title.y=element_text(colour="grey20",size=12),
legend.text = element_text(size = 12))
if(!is.null(thredhold)) p = p + geom_hline(yintercept = thredhold,colour="tomato", linetype="dashed")
# if(!is.null(colour)) {col <- unlist(strsplit(colour,","))} else col = c("grey50","grey70")
}
Format = strsplit(outputFormat,",")[[1]]
sapply(Format,function(x){
if(x == "pdf") temp = paste0("cairo_pdf(file = paste0(output,\'.\',x), bg = \'transparent\', width = width, height = height)")
else temp = paste0(x,"(file = paste0(output,\'.\',x), bg = \'transparent\', width = width*72, height = height*72, type = \'cairo\')")
eval(parse(text = temp))
# png(file = paste0(output_hc,".",x), bg = "transparent", width = width, height = height)
plot(p)
dev.off()
})
return(p)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.