2019-04-22

rm(list=ls())
require("ggplot2")
require("reshape")
library(wesanderson)
require("plyr")

setwd("~/Documents/SD1/SD1_analyses/")

ts=read.csv("input_files/rep_set_tax_assignments_silva_clean2.csv")
colnames(ts)=c("OTUID", "kingdom_s", "phylum_s", "class_s", "order_s", "family_s", "genus_s", "species_s", "evalue", "ref")

tgg=read.csv("input_files/rep_set_tax_assignments_split2.csv")

s1=read.csv("input_files/SD1_map_clean.csv", header=TRUE, row.names=1)
s1=as.data.frame(s1)
c2norm=read.csv("input_files/c2norm.csv", header=TRUE, row.names=1)
#c2norm=as.matrix(sapply(c2norm, as.numeric)) 
check=cbind(rownames(s1), colnames(c2norm))

##Import all sig otu files 
Genotype_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Genotype_sigOTUs_alltax.csv")
T3_disease_state_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/T3_disease_state_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Dose_disease_state_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Dose_disease_state_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Timepoint_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Timepoint_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Dose_disease_state_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Dose_disease_state_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Timepoint_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Timepoint_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Dose_disease_state_Timepoint_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Dose_disease_state_Timepoint_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Dose_disease_state_Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Dose_disease_state_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Timepoint_Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Timepoint_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Dose_disease_state_Timepoint_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Dose_disease_state_Timepoint_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Dose_disease_state_Dose_site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Dose_disease_state_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Timepoint_Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Timepoint_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Dose_disease_state_Timepoint_Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Dose_disease_state_Timepoint_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)
Site_Dose_disease_state_Timepoint_Dose_Site_sigOTUs_alltaxr=read.csv("Output_files/GLMM_results/lists/Site_Dose_disease_state_Timepoint_Dose_Site_sigOTUs_alltax.csv", header=TRUE, row.names=1)


#####looking for otus that are more abundant in T3 disease state DD time 2 and time 3. 
T3ex=T3_disease_state_sigOTUs_alltaxr
T3ex2=T3ex[!T3ex$OTUID %in% Dose_Site_sigOTUs_alltaxr$OTUID,]
T3ex3=T3ex2[!T3ex2$OTUID %in% Site_Dose_Site_sigOTUs_alltaxr$OTUID,]

##Create sig otu otu table
c2norm_T3ex3=c2norm[row.names(c2norm) %in% T3ex3$OTUID, ]
b_T3ex3=t(c2norm_T3ex3)
d_T3ex3=cbind.data.frame(s1$T3_disease_state, s1$Dose_disease_state, s1$Timepoint, b_T3ex3)
colnames(d_T3ex3)[1]="T3_disease_state"
colnames(d_T3ex3)[2]="Dose_disease_state"
colnames(d_T3ex3)[3]="Timepoint"

#check order of samples
check=cbind(row.names(d_T3ex3), s1[1:1])

#average diseased and healthy
##T3 disease state
se=function(x) sd(x)/sqrt(length(x))
d_T3ex3_dh_mean=ddply(d_T3ex3, ~Dose_disease_state:T3_disease_state:Timepoint, colwise(mean))
d_T3ex3_dh_se=ddply(d_T3ex3, ~Dose_disease_state:T3_disease_state:Timepoint, colwise(se))
row.names(d_T3ex3_dh_mean)=gsub(":", "_", d_T3ex3_dh_mean[,1])
row.names(d_T3ex3_dh_se)=gsub(":", "_", d_T3ex3_dh_se[,1])
d_T3ex3_dh_mean=d_T3ex3_dh_mean[,-c(1)]
d_T3ex3_dh_se=d_T3ex3_dh_se[,-c(1)]
d_T3ex3_dh_mean=t(d_T3ex3_dh_mean)
d_T3ex3_dh_se=t(d_T3ex3_dh_se)
d_T3ex3_dh_mean=as.data.frame(d_T3ex3_dh_mean)
d_T3ex3_dh_se=as.data.frame(d_T3ex3_dh_se)
d_T3ex3_dh=d_T3ex3_dh_mean
#d_T3ex3_dh$DvH=d_T3ex3_dh$Diseased_Diseased_two-d_T3ex3_dh$Healthy_Healthy_two
UpT2D=subset(d_T3ex3_dh, d_T3ex3_dh$Diseased_Diseased_two>0&d_T3ex3_dh$Diseased_Diseased_three>0&d_T3ex3_dh$Diseased_Diseased_dose==0&d_T3ex3_dh$Healthy_Healthy_one>0)
UpT2Dtax=merge(UpT2D, ts, by.x="row.names", by.y="OTUID")
UpT2Dtax$DvHdose=UpT2Dtax$Diseased_Diseased_dose-UpT2Dtax$Healthy_Healthy_dose
UpT2Dtax$DvHT2=UpT2Dtax$Diseased_Diseased_two-UpT2Dtax$Healthy_Healthy_two
UpT2Dtax$DvHT3=UpT2Dtax$Diseased_Diseased_three-UpT2Dtax$Healthy_Healthy_three
UpT2DtaxUp=subset(UpT2Dtax, UpT2Dtax$DvHT2>0&UpT2Dtax$DvHT3>0)
UpT2DtaxUp$Healthy_Diseased_two=NULL
d_T3ex3_dh_se$Healthy_Diseased_two=NULL
T2D_se=d_T3ex3_dh_se[row.names(d_T3ex3_dh_se) %in% UpT2DtaxUp$Row.names,]
T2D_se=as.matrix(T2D_se)
T2D_se_filt=T2D_se[apply(T2D_se, 2, function(x) all(x<500)),]
UpT2DtaxUpse=UpT2DtaxUp[UpT2DtaxUp$Row.names %in% row.names(T2D_se_filt),]
UpT2DtaxUpse$DvH=NULL

colnames(UpT2DtaxUpse)= c("OTUID", "D_D_dose", "D_D_T3", "D_D_T2", 
                          "D_H_T3", "D_H_T2", "H_H_dose", "H_H_T1", 
                          "H_H_T3", "H_H_T2", "kingdom_s", "phylum_s", "class_s",
                          "order_s", "family_s", "genus_s", "species_s",               
                          "evalue", "ref", "DvHdose", "DvHT2", "DvHT3")  
colnames(T2D_se_filt)= c("D_D_dose", "D_D_T3", "D_D_T2", 
                         "D_H_T3", "D_H_T2", "H_H_dose", "H_H_T1", 
                         "H_H_T3", "H_H_T2") 


Primary=cbind(UpT2DtaxUpse[,1], UpT2DtaxUpse[,11:22])
colnames(Primary)[1]="OTUID"
write.csv(Primary, file="Output_files/GLMM_results/Primary_responders2.csv")
countso=as.data.frame(table(UpT2DtaxUpse$order_s))
countso=countso[order(-countso$Freq),]
counts=as.data.frame(table(UpT2DtaxUpse$family_s))
counts=counts[order(-counts$Freq),]

###pull out families and graph
#Rhodobacteraceae
UpT2DtaxUpr=subset(UpT2DtaxUpse, UpT2DtaxUpse$family_s=="Rhodobacteraceae")
ser=T2D_se_filt[row.names(T2D_se_filt) %in% UpT2DtaxUpr$OTUID,]
ser=as.data.frame(ser)
ser$OTUID=row.names(ser)

otu.long_d = melt(UpT2DtaxUpr, id = c("OTUID","family_s"), 
                  measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                              "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                  ))
otu.long_se = melt(ser, id = c("OTUID"), 
                   measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                               "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                   ))
colnames(otu.long_se)[3]="SE"
otu.long_se$OTU_group=paste(otu.long_se$OTUID, otu.long_d$variable, sep="_")


otu.long_d$OTU_group=paste(otu.long_d$OTUID, otu.long_d$variable, sep="_")

otu.long_full=merge(otu.long_d, otu.long_se, by="OTU_group")
otu.long_full=otu.long_full[,3:8]
colnames(otu.long_full)=c("family_s", "variable", "value", "OTUID", "variable.y", "SE")

otu.long_full$OTUID = factor(otu.long_full$OTUID, levels = otu.long_full$OTUID[order(otu.long_full$family_s)])
otu.long_full$variable = factor(otu.long_full$variable, levels = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                                                                   "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"))

pdf(file="Output_files/GLMM_results/plots/Primary_responders_r.pdf", width=6, height=10)
dodge=position_dodge(width=0.9)
plotgb=ggplot(data=otu.long_full, aes(x=OTUID, y=value ))+
  geom_bar(aes(fill=family_s),position=dodge, stat="identity") + theme_bw() +
  theme(legend.position = "left") +
  ggtitle("Primary pathogens") + 
  facet_grid(variable~., scale="free_y") + ylim(-1,15)
plotgb=plotgb + theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank()) +
  scale_x_discrete("OTUID")+scale_fill_manual(values="#E58601")
plotgb=plotgb+geom_errorbar(aes(ymin=value-SE, ymax=value+SE), position=dodge, width=0.001)
plotgb

dev.off()

####Plot cryomorphaceae##################
#separate Cryomorphaceae
UpT2DtaxUpp=subset(UpT2DtaxUpse, UpT2DtaxUpse$family_s=="Cryomorphaceae")
seps=T2D_se_filt[row.names(T2D_se_filt) %in% UpT2DtaxUpp$OTUID,]
seps=as.data.frame(seps)
seps$OTUID=row.names(seps)

#melt dataframes
otu.long_d = melt(UpT2DtaxUpp, id = c("OTUID","family_s"), 
                  measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                              "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                  ))
otu.long_se = melt(seps, id = c("OTUID"), 
                   measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                               "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                   ))
#combine means and ses
colnames(otu.long_se)[3]="SE"
otu.long_se$OTU_group=paste(otu.long_se$OTUID, otu.long_d$variable, sep="_")
otu.long_d$OTU_group=paste(otu.long_d$OTUID, otu.long_d$variable, sep="_")
otu.long_full=merge(otu.long_d, otu.long_se, by="OTU_group")
otu.long_full=otu.long_full[,3:8]
colnames(otu.long_full)=c("family_s", "variable", "value", "OTUID", "variable.y", "SE")

otu.long_full$OTUID = factor(otu.long_full$OTUID, levels = otu.long_full$OTUID[order(otu.long_full$family_s)])
otu.long_full$variable = factor(otu.long_full$variable, levels = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                                                                   "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"))

pdf(file="Output_files/GLMM_results/plots/Primary_responders_c.pdf", width=6, height=10)
dodge=position_dodge(width=0.9)
plotgb=ggplot(data=otu.long_full, aes(x=OTUID, y=value ))+
  geom_bar(aes(fill=family_s),position=dodge, stat="identity") + theme_bw() +
  theme(legend.position = "left") +
  ggtitle("Primary pathogens") + 
  facet_grid(variable~., scale="free_y") + ylim(-1,150)
plotgb=plotgb + theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank()) +
  scale_x_discrete("OTUID")+scale_fill_manual(values="#5B1A18")
plotgb=plotgb+geom_errorbar(aes(ymin=value-SE, ymax=value+SE), position=dodge, width=0.001)
plotgb

dev.off()


####Alteromonadaceae
UpT2DtaxUpa=subset(UpT2DtaxUpse, UpT2DtaxUpse$family_s=="Alteromonadaceae")
sea=T2D_se_filt[row.names(T2D_se_filt) %in% UpT2DtaxUpa$OTUID,]
sea=as.data.frame(sea)
sea$OTUID=row.names(sea)

otu.long_d = melt(UpT2DtaxUpa, id = c("OTUID","family_s"), 
                  measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                              "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                  ))
otu.long_se = melt(sea, id = c("OTUID"), 
                   measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                               "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                   ))
colnames(otu.long_se)[3]="SE"
otu.long_se$OTU_group=paste(otu.long_se$OTUID, otu.long_d$variable, sep="_")

otu.long_d$OTU_group=paste(otu.long_d$OTUID, otu.long_d$variable, sep="_")

otu.long_full=merge(otu.long_d, otu.long_se, by="OTU_group")
otu.long_full=otu.long_full[,3:8]
colnames(otu.long_full)=c("family_s", "variable", "value", "OTUID", "variable.y", "SE")

otu.long_full$OTUID = factor(otu.long_full$OTUID, levels = otu.long_full$OTUID[order(otu.long_full$family_s)])
otu.long_full$variable = factor(otu.long_full$variable, levels = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                                                                   "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"))

pdf(file="Output_files/GLMM_results/plots/Primary_responders_a.pdf", width=6, height=10)
dodge=position_dodge(width=0.9)
plotgb=ggplot(data=otu.long_full, aes(x=OTUID, y=value ))+
  geom_bar(aes(fill=family_s),position=dodge, stat="identity") + theme_bw() +
  theme(legend.position = "left") +
  ggtitle("Primary pathogens") + 
  facet_grid(variable~., scale="free_y") + ylim(-1,75)
plotgb=plotgb + theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank()) +
  scale_x_discrete("OTUID")+scale_fill_manual(values="#E6A0C4")
plotgb=plotgb+geom_errorbar(aes(ymin=value-SE, ymax=value+SE), position=dodge, width=0.001)
plotgb

dev.off()



####Flavobacteriaceae
UpT2DtaxUpf=subset(UpT2DtaxUpse, UpT2DtaxUpse$family_s=="Flavobacteriaceae")
sef=T2D_se_filt[row.names(T2D_se_filt) %in% UpT2DtaxUpf$OTUID,]
sef=as.data.frame(sef)
sef$OTUID=row.names(sef)

otu.long_d = melt(UpT2DtaxUpf, id = c("OTUID","family_s"), 
                  measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                              "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                  ))
otu.long_se = melt(sef, id = c("OTUID"), 
                   measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                               "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                   ))
colnames(otu.long_se)[3]="SE"
otu.long_se$OTU_group=paste(otu.long_se$OTUID, otu.long_d$variable, sep="_")

otu.long_d$OTU_group=paste(otu.long_d$OTUID, otu.long_d$variable, sep="_")

otu.long_full=merge(otu.long_d, otu.long_se, by="OTU_group")
otu.long_full=otu.long_full[,3:8]
colnames(otu.long_full)=c("family_s", "variable", "value", "OTUID", "variable.y", "SE")

otu.long_full$OTUID = factor(otu.long_full$OTUID, levels = otu.long_full$OTUID[order(otu.long_full$family_s)])
otu.long_full$variable = factor(otu.long_full$variable, levels = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                                                                   "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"))

pdf(file="Output_files/GLMM_results/Primary_responders_fl.pdf", width=6, height=10)
dodge=position_dodge(width=0.9)
plotgb=ggplot(data=otu.long_full, aes(x=OTUID, y=value ))+
  geom_bar(aes(fill=family_s),position=dodge, stat="identity") + theme_bw() +
  theme(legend.position = "left") +
  ggtitle("Primary pathogens") + 
  facet_grid(variable~., scale="free_y") + ylim(-1,70)
plotgb=plotgb + theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank()) +
  scale_x_discrete("OTUID")+scale_fill_manual(values="#78B7C5")
plotgb=plotgb+geom_errorbar(aes(ymin=value-SE, ymax=value+SE), position=dodge, width=0.001)
plotgb

dev.off()


###Primary colonizers not in T3
UpT2Donly=subset(d_T3ex3_dh, d_T3ex3_dh$Diseased_Diseased_two>0&d_T3ex3_dh$Diseased_Diseased_dose>0&d_T3ex3_dh$Healthy_Diseased_two==0)
UpT2Donlytax=merge(UpT2Donly, ts, by.x="row.names", by.y="OTUID")
UpT2Donlytax$DvHdose=UpT2Donlytax$Diseased_Diseased_dose-UpT2Donlytax$Healthy_Healthy_dose
UpT2Donlytax$DvHT2=UpT2Donlytax$Diseased_Diseased_two-UpT2Donlytax$Healthy_Healthy_two
UpT2Donlytax$DvHT3=UpT2Donlytax$Diseased_Diseased_three-UpT2Donlytax$Healthy_Healthy_three
UpT2DonlytaxUp=subset(UpT2Donlytax, UpT2Donlytax$DvHdose>0&UpT2Donlytax$DvHT2>0&UpT2Donlytax$DvHT3<=0&
                        UpT2Donlytax$Healthy_Healthy_dose==0)
UpT2DonlytaxUp$Healthy_Diseased_two=NULL
d_T3ex3_dh_se$Healthy_Diseased_two=NULL
T2D_only_se=d_T3ex3_dh_se[row.names(d_T3ex3_dh_se) %in% UpT2DonlytaxUp$Row.names,]
T2D_only_se=as.matrix(T2D_only_se)
T2D_only_se_filt=T2D_only_se[apply(T2D_only_se, 2, function(x) all(x<500)),]
UpT2DonlytaxUpse=UpT2DonlytaxUp[UpT2DonlytaxUp$Row.names %in% row.names(T2D_only_se_filt),]
UpT2DonlytaxUpse$DvH=NULL

Primary2=cbind(UpT2DonlytaxUpse[,1], UpT2DonlytaxUpse[,11:22])
colnames(Primary2)[1]="OTUID"
write.csv(Primary2, file="Output_files/GLMM_results/Primary_pathogens_2.csv")
countso=as.data.frame(table(UpT2DonlytaxUpse$order_s))
countso=countso[order(-countso$Freq),]
counts=as.data.frame(table(UpT2DonlytaxUpse$family_s))
counts=counts[order(-counts$Freq),]

colnames(UpT2DonlytaxUpse)= c("OTUID", "D_D_dose", "D_D_T3", "D_D_T2", 
                              "D_H_T3", "D_H_T2", "H_H_dose", "H_H_T1", 
                              "H_H_T3", "H_H_T2", "kingdom_s", "phylum_s", "class_s",
                              "order_s", "family_s", "genus_s", "species_s",               
                              "evalue", "ref", "DvHdose", "DvHT2", "DvHT3")  
colnames(T2D_only_se_filt)= c("D_D_dose", "D_D_T3", "D_D_T2", 
                              "D_H_T3", "D_H_T2", "H_H_dose", "H_H_T1", 
                              "H_H_T3", "H_H_T2")
T2D_only_se_filt=as.data.frame(T2D_only_se_filt)
T2D_only_se_filt$OTUID=row.names(T2D_only_se_filt)
###graph
otu.long_d = melt(UpT2DonlytaxUpse, id = c("OTUID","order_s"), 
                  measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                              "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                  ))
otu.long_se = melt(T2D_only_se_filt, id = c("OTUID"), 
                   measure = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                               "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"
                   ))
colnames(otu.long_se)[3]="SE"
otu.long_se$OTU_group=paste(otu.long_se$OTUID, otu.long_d$variable, sep="_")
otu.long_d$OTU_group=paste(otu.long_d$OTUID, otu.long_d$variable, sep="_")

otu.long_full=merge(otu.long_d, otu.long_se, by="OTU_group")
otu.long_full=otu.long_full[,3:8]
colnames(otu.long_full)=c("order", "variable", "value", "OTUID", "variable.y", "SE")

otu.long_full$OTUID = factor(otu.long_full$OTUID, levels = otu.long_full$OTUID[order(otu.long_full$order)])
otu.long_full$variable = factor(otu.long_full$variable, levels = c("H_H_dose", "D_D_dose", "H_H_T1", "H_H_T2", 
                                                                   "D_H_T2", "D_D_T2", "H_H_T3", "D_H_T3", "D_D_T3"))

pdf(file="Output_files/GLMM_results/Temporary_primary_colonizers.pdf", width=6, height=7)

dodge=position_dodge(width=0.9)
plotgb=ggplot(data=otu.long_full, aes(x=OTUID, y=value ))+
  geom_bar(aes(fill=order),position=dodge, stat="identity") + theme_bw() +
  theme(legend.position = "left") +
  ggtitle("Temporary primary colonizers") + 
  facet_grid(variable~., scale="free_y") + ylim(-1,40)
plotgb=plotgb + theme(axis.text.x = element_blank(), panel.grid.major.x = element_blank()) +
  scale_x_discrete("OTUID")
plotgb=plotgb+geom_errorbar(aes(ymin=value-SE, ymax=value+SE), position=dodge, width=0.001)
plotgb

dev.off()

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 225,764评论 6赞 524
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 96,848评论 3赞 409
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 173,181评论 0赞 370
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 61,430评论 1赞 303
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 70,451评论 6赞 403
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 53,879评论 1赞 316
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 42,163评论 3赞 432
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 41,189评论 0赞 281
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 47,727评论 1赞 328
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 39,741评论 3赞 350
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 41,852评论 1赞 358
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 37,444评论 5赞 352
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 43,162评论 3赞 341
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 33,569评论 0赞 25
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 34,735评论 1赞 278
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 50,443评论 3赞 383
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 46,931评论 2赞 368

2019-04-22

推荐阅读更多精彩内容