-
Notifications
You must be signed in to change notification settings - Fork 11
/
tukey_plots.R
72 lines (72 loc) · 2.77 KB
/
tukey_plots.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
#install.packages("klaR",dependencies=TRUE)
#install.packages("spdep",dependencies=TRUE)
#install.packages("AlgDesign",dependencies=TRUE)
#install.packages("agricolae",dependencies=TRUE)
#install.packages("multcomp",dependencies=TRUE)
#install.packages("repr")
library(agricolae)
library(multcomp)
library(repr)
library(RColorBrewer)
setwd("~/Documents/UCDavis/dib/MMETSP")
special_flowers = c("MMETSP0693","MMETSP1019","MMETSP0923","MMETSP0008","MMETSP1002","MMETSP1325","MMETSP1018",
"MMETSP1346","MMETSP0088","MMETSP0092","MMETSP0717","MMETSP0223","MMETSP0115","MMETSP0196",
"MMETSP0197","MMETSP0398","MMETSP0399","MMETSP0922")
tax_raw <- read.csv("assembly_evaluation_data/MMETSP_all_evaluation_matrix.csv")
#head(tax_raw)
#colors=palette(brewer.pal(n=7,name="Dark2"))
tax_raw <- tax_raw[!tax_raw$SampleName %in% special_flowers,]
phylum_data <-tax_raw[,c(2,34,45,64,61,49,114,41,45,54)]
phylum <- phylum_data$Phylum
# Restrict data to top 7 most common Phyla
sub_phy<-c("Bacillariophyta","Dinophyta","Ochrophyta","Haptophyta","Ciliophora","Chlorophyta","Cryptophyta")
sub<-phylum_data[phylum_data$Phylum %in% sub_phy,]
colors=rainbow(length(sub_phy))
# Mean %ORF
fit <- aov(Unique_kmers ~ Phylum,data=sub)
a<-HSD.test(fit,"Phylum",group=TRUE)
tuk<-glht(fit,linfct=mcp(Phylum="Tukey"))
summary(tuk)
tuk.cld<-cld(tuk)
opar<-par(mai=c(1,1,1.5,1))
options(repr.plot.width=12, repr.plot.height=10)
pdf("paper/Figure8_tukey_taxa_unique_kmers.pdf",width=11,height=8.5)
plot(tuk.cld,col=colors,ylab="Mean %ORF")
dev.off()
plot(tuk.cld,col=colors,ylab="Mean %ORF")
# Unique k-mers
fit2 <- aov(mean_orf_percent.x ~ Phylum,data=sub)
b<-TukeyHSD(fit2,"Phylum",conf.level=0.95)
tuk<-glht(fit2,linfct=mcp(Phylum="Tukey"))
summary(tuk)
tuk.cld<-cld(tuk)
opar<-par(mai=c(1,1,1.5,1))
options(repr.plot.width=12, repr.plot.height=10)
pdf("paper/Figure8_tukey_taxa_mean_orf.pdf",width=11,height=8.5)
plot(tuk.cld,col=colors)
dev.off()
plot(tuk.cld,col=colors, ylab="Unique k-mers")
# Number of Input reads
fit3 <- aov(Input.Reads ~ Phylum,data=sub)
b<-TukeyHSD(fit3,"Phylum",conf.level=0.95)
tuk<-glht(fit3,linfct=mcp(Phylum="Tukey"))
summary(tuk)
tuk.cld<-cld(tuk)
opar<-par(mai=c(1,1,1.5,1))
options(repr.plot.width=12, repr.plot.height=10)
pdf("paper/Figure8_tukey_taxa_input_Reads.pdf",width=11,height=8.5)
plot(tuk.cld,col=colors)
dev.off()
plot(tuk.cld,col=colors, ylab="Number of Input Reads")
# Number of contigs
fit4 <- aov(n_seqs.x ~ Phylum,data=sub)
b<-TukeyHSD(fit4,"Phylum",conf.level=0.95)
tuk<-glht(fit4,linfct=mcp(Phylum="Tukey"))
summary(tuk)
tuk.cld<-cld(tuk)
opar<-par(mai=c(1,1,1.5,1))
options(repr.plot.width=12, repr.plot.height=10)
pdf("paper/Figure8_tukey_taxa_n_contigs.pdf",width=11,height=8.5)
plot(tuk.cld,col=colors)
dev.off()
plot(tuk.cld,col=colors, ylab="Number of Contigs")