Skip to content

Commit d827041

Browse files
TankMermaidTankMermaid
committed
test for mock community
0 parents  commit d827041

14 files changed

+566
-0
lines changed

‎application/20180809 干湿稻及小麦根腐病实验结果及讨论.docx‎

42.1 KB
Binary file not shown.

‎application/20180809 干湿稻根腐病 Figs.pptx‎

5.07 MB
Binary file not shown.

‎application/20180809 绝对定量-gxx4(文献引用版).docx‎

152 KB
Binary file not shown.

‎mock/20180817 MOCK -统计补充.pptx‎

789 KB
Binary file not shown.

‎mock/20180817 真菌 mock data.xlsx‎

19.6 KB
Binary file not shown.

‎mock/20180817 细菌mock data.xlsx‎

24.3 KB
Binary file not shown.

‎mock/FIG2b_FIGS1_test_bac.R‎

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
# originally by Tank
2+
3+
# Email:[email protected]
4+
5+
rm(list=ls())
6+
# install.packages("DescTools")
7+
library(DescTools)
8+
library(dplyr)
9+
library(vegan)
10+
library(reshape2)
11+
########### setting the working directory and print it ###################
12+
tem = "final_test"
13+
setwd("~/xiaoxuan/final/mock/")
14+
print(paste("Your working directory is in",getwd()))
15+
16+
17+
figures.dir <- paste("~/xiaoxuan/final/mock/bac/",tem,'/',sep = '')
18+
table.dir <- paste("~/xiaoxuan/final/mock/table/",tem,'/',sep = '')
19+
20+
fig_flag <- dir.exists(figures.dir)
21+
if( isTRUE(!fig_flag)){
22+
dir.create(figures.dir)
23+
}
24+
25+
tab_flag <- dir.exists(table.dir)
26+
if( isTRUE(!tab_flag)){
27+
dir.create(table.dir)
28+
}
29+
####################################
30+
31+
32+
#####spike-in design in this batch #####################
33+
spike <- c("BI-OS-11-3","BI-OS-12-4","BI-OS-10-2")
34+
35+
36+
#1 design Mapping file and Sample id
37+
design = read.table("bac/design.txt", header=T, row.names= 1, sep="\t")
38+
# design$SampleID <- row.names(design)
39+
# sample_list <- as.matrix(row.names(design))
40+
41+
otu_table = read.table("bac_RA_withoutSpike.txt", row.names= 1, header=1, sep="\t")
42+
otu_table_num=as.data.frame(lapply(otu_table, function(x) as.numeric(gsub("\\%", "", x))/100))
43+
row.names(otu_table_num)=row.names(otu_table)
44+
otu_table_t= as.data.frame(t(otu_table_num))
45+
# otu_table_t<-data.frame(lapply(otu_table_t, function(x) as.numeric(gsub("\\%", "", x))/100))
46+
otu_table_meta=merge(design,otu_table_t,by="row.names")
47+
colnames(otu_table_meta)[1]= "SampleID"
48+
49+
dat=melt(otu_table_meta,id.vars = c("SampleID", "Other","Description"),measure.vars = row.names(otu_table))
50+
# dat_s=dat[,c(3,5)]
51+
52+
# DunnettTest(value ~ Description, data = dat[,c(3:5)],control ="E00", conf.level=0.9)
53+
# boxplot(Ozone ~ Month, data = airquality)
54+
55+
#################补充FIG 2b 和FIG S1中，3种比例条件下，去除spike-in的菌株RA 和E00的菌株RA 有无显著差别
56+
57+
vec1=row.names(otu_table)
58+
vec2=unique(as.vector(dat$Other))
59+
library("multcomp")
60+
61+
for (i in 1:length(vec1)){
62+
for (j in 1:length(vec2)){
63+
64+
# data("recovery", package = "multcomp")
65+
RA.aov <- aov(value ~ Description, data = dat[dat$variable %in% row.names(otu_table)[i]& dat$Other %in% vec2[j],])
66+
RA.mc <- glht(RA.aov,
67+
linfct = mcp(Description = "Dunnett"),
68+
alternative = "less")
69+
70+
71+
print(summary(RA.mc))
72+
}
73+
}
74+
75+
#### to continue to rewrite with Rmarkdown
76+
77+
78+
79+
#### 补充不同spike-in梯度下9株菌在1:1:1 VS 2:2:2时RA值有无统计学差异
80+
dat_21=dat[dat$Other %in% c("1:1:1","2:2:2"), ]
81+
vec3=row.names(otu_table)
82+
vec4=unique(as.vector(dat_21$Description))
83+
84+
for (i in 1:length(vec3)){
85+
for (j in 1:length(vec4)){
86+
87+
# data("recovery", package = "multcomp")
88+
RA.t <- t.test(value ~ Other, data = dat_21[dat_21$variable %in% vec3[i]& dat_21$Description%in% vec4[j],])
89+
# RA.mc <- glht(RA.t,
90+
# linfct = mcp(Description = "t.test"),
91+
# alternative = "less")
92+
93+
94+
print(paste0("t-test for ",vec3[i]," ",vec4[j]," ratio for 2a:2b:2c over 1a:1b:1c is constructing"))
95+
print(RA.t)
96+
97+
}
98+
}
99+
100+
101+
######### 不同spike-in梯度下proteo的AA 在1:1:1 VS 2:2:1 时有无统计差异，理想情况是AA 无统计学差异
102+
103+
otu_table = read.table("bac_AA_mock.txt", row.names= 1, header=1, sep="\t")
104+
# row.names(otu_table_num)=row.names(otu_table)
105+
otu_table_t= as.data.frame(t(otu_table))
106+
# otu_table_t<-data.frame(lapply(otu_table_t, function(x) as.numeric(gsub("\\%", "", x))/100))
107+
otu_table_meta=merge(design,otu_table_t,by="row.names")
108+
colnames(otu_table_meta)[1]= "SampleID"
109+
110+
dat=melt(otu_table_meta,id.vars = c("SampleID", "Other","Description"),measure.vars = row.names(otu_table))
111+
112+
113+
dat_21=dat[dat$Other %in% c("1:1:1","2:2:1"), ]
114+
vec3=row.names(otu_table)
115+
vec4=unique(as.vector(dat_21$Description))
116+
117+
for (i in 6:length(vec3)){
118+
for (j in 1:length(vec4)){
119+
120+
# data("recovery", package = "multcomp")
121+
RA.t <- t.test(value ~ Other, data = dat_21[dat_21$variable %in% vec3[i]& dat_21$Description%in% vec4[j],])
122+
# RA.mc <- glht(RA.t,
123+
# linfct = mcp(Description = "t.test"),
124+
# alternative = "less")
125+
126+
127+
print(paste0("t-test for ",vec3[i]," ",vec4[j]," AAratio for 2a:2b:1c over 1a:1b:1c is constructing"))
128+
print(RA.t)
129+
130+
}
131+
}
132+
133+

‎mock/FIG3b_FIGS4_test_fungi.R‎

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
# originally by Tank
2+
3+
# Email:[email protected]
4+
5+
rm(list=ls())
6+
# install.packages("DescTools")
7+
library(DescTools)
8+
library(dplyr)
9+
library(vegan)
10+
library(reshape2)
11+
########### setting the working directory and print it ###################
12+
tem = "final_test"
13+
setwd("~/xiaoxuan/final/mock/")
14+
print(paste("Your working directory is in",getwd()))
15+
16+
17+
figures.dir <- paste("~/xiaoxuan/final/mock/fun/",tem,'/',sep = '')
18+
table.dir <- paste("~/xiaoxuan/final/mock/table/",tem,'/',sep = '')
19+
20+
fig_flag <- dir.exists(figures.dir)
21+
if( isTRUE(!fig_flag)){
22+
dir.create(figures.dir)
23+
}
24+
25+
tab_flag <- dir.exists(table.dir)
26+
if( isTRUE(!tab_flag)){
27+
dir.create(table.dir)
28+
}
29+
####################################
30+
31+
32+
#####spike-in design in this batch #####################
33+
spike <- c("BI-OS-11-3","BI-OS-12-4","BI-OS-10-2")
34+
35+
36+
#1 design Mapping file and Sample id
37+
design = read.table("fun/design.txt", header=T, row.names= 1, sep="\t")
38+
# design$SampleID <- row.names(design)
39+
# sample_list <- as.matrix(row.names(design))
40+
41+
otu_table = read.table("fungi_RA_withoutSpike.txt", row.names= 1, header=1, sep="\t")
42+
# otu_table_num=as.data.frame(lapply(otu_table, function(x) as.numeric(gsub("\\%", "", x))/100))
43+
# row.names(otu_table_num)=row.names(otu_table)
44+
otu_table_t= as.data.frame(t(otu_table))
45+
# otu_table_t<-data.frame(lapply(otu_table_t, function(x) as.numeric(gsub("\\%", "", x))/100))
46+
otu_table_meta=merge(design,otu_table_t,by="row.names")
47+
colnames(otu_table_meta)[1]= "SampleID"
48+
49+
dat=melt(otu_table_meta,id.vars = c("SampleID", "Other","Description"),measure.vars = row.names(otu_table))
50+
# dat_s=dat[,c(3,5)]
51+
52+
# DunnettTest(value ~ Description, data = dat[,c(3:5)],control ="E00", conf.level=0.9)
53+
# boxplot(Ozone ~ Month, data = airquality)
54+
55+
#################补充FIG 2b 和FIG S1中，3种比例条件下，去除spike-in的菌株RA 和E00的菌株RA 有无显著差别
56+
57+
vec1=row.names(otu_table)
58+
vec2=unique(as.vector(dat$Other))
59+
library("multcomp")
60+
61+
for (i in 1:length(vec1)){
62+
for (j in 1:length(vec2)){
63+
64+
# data("recovery", package = "multcomp")
65+
RA.aov <- aov(value ~ Description, data = dat[dat$variable %in% row.names(otu_table)[i]& dat$Other %in% vec2[j],])
66+
RA.mc <- glht(RA.aov,
67+
linfct = mcp(Description = "Dunnett"),
68+
alternative = "less")
69+
70+
print(paste0("dunnet-test for ",vec1[i]," ",vec2[j]," ratio is constructing"))
71+
print(summary(RA.mc))
72+
}
73+
}
74+
75+
#### to continue to rewrite with Rmarkdown
76+
77+
78+
79+
#### 补充不同spike-in梯度下9株菌在1:1:1 VS 2:2:2时RA值有无统计学差异
80+
dat_21=dat[dat$Other %in% c("1:1:1","2:2:2"), ]
81+
vec3=row.names(otu_table)
82+
vec4=unique(as.vector(dat_21$Description))
83+
84+
for (i in 1:length(vec3)){
85+
for (j in 1:length(vec4)){
86+
87+
# data("recovery", package = "multcomp")
88+
RA.t <- t.test(value ~ Other, data = dat_21[dat_21$variable %in% vec3[i]& dat_21$Description%in% vec4[j],])
89+
# RA.mc <- glht(RA.t,
90+
# linfct = mcp(Description = "t.test"),
91+
# alternative = "less")
92+
93+
94+
print(paste0("t-test for ",vec3[i]," ",vec4[j]," ratio for 2a:2b:2c over 1a:1b:1c is constructing"))
95+
print(RA.t)
96+
97+
}
98+
}
99+
100+
101+
######### 不同spike-in梯度下Asco的AA 在1:1:1 VS 2:2:1 时有无统计差异，理想情况是AA 无统计学差异
102+
103+
otu_table = read.table("fungi_AA_mock.txt", row.names= 1, header=1, sep="\t")
104+
# row.names(otu_table_num)=row.names(otu_table)
105+
otu_table_t= as.data.frame(t(otu_table))
106+
# otu_table_t<-data.frame(lapply(otu_table_t, function(x) as.numeric(gsub("\\%", "", x))/100))
107+
otu_table_meta=merge(design,otu_table_t,by="row.names")
108+
colnames(otu_table_meta)[1]= "SampleID"
109+
110+
dat=melt(otu_table_meta,id.vars = c("SampleID", "Other","Description"),measure.vars = row.names(otu_table))
111+
112+
113+
dat_21=dat[dat$Other %in% c("1:1:1","2:2:1")&!dat$Description %in% "E00", ]
114+
vec3=row.names(otu_table)
115+
vec4=unique(as.vector(dat_21$Description))
116+
117+
for (i in 1:2){
118+
for (j in 1:length(vec4)){
119+
120+
# data("recovery", package = "multcomp")
121+
AA.t <- t.test(value ~ Other, data = dat_21[dat_21$variable %in% vec3[i]& dat_21$Description%in% vec4[j],])
122+
# RA.mc <- glht(RA.t,
123+
# linfct = mcp(Description = "t.test"),
124+
# alternative = "less")
125+
126+
127+
print(paste0("t-test for ",vec3[i]," ",vec4[j]," AA ratio for 2a:2b:1c over 1a:1b:1c is constructing"))
128+
print(AA.t)
129+
130+
}
131+
}
132+
133+

0 commit comments

Comments
 (0)