Gut microbial metabolic disorder in depression: insights from computational modeling and mediation analysis.

1. rm(list = ls())
2. library(Maaslin2)
3. library(MicrobiotaProcess)
4. library(tidyverse)
5. library(compositions)
6. library(pheatmap)
7. library(phyloseq)
8. library(ggalluvial)
9. library(dplyr)
10. library(qvalue)
11. library(gridExtra)
12. library(ggsignif)
13. library(ggpubr)
14. library(data.table)
15. library(ggplot2)
16. library(ggExtra)
17. library(ggtext)
18. library(vegan)
19. library(caret)
20. library(data.table)
21. library(ggtext)
22. library(ggbreak)
23. meta_D<-read.table('D:/microbiome/micom/202404/paper/data/meta_Depression.txt',sep = '\t',header = T)
24. meta_H<-read.table('D:/microbiome/micom/202404/paper/data/meta_healthy.txt',sep = '\t',header = T)
25. colnames(meta_H)[8]<-'disease'
26. colnames(meta_D)[8]<-'disease'
27. meta_H$disease<-'Healthy'
28. meta_D$disease<-'Depression'
29. age<-rbind(meta_H[,c(6,8),drop=F],meta_D[,c(6,8),drop=F])
30. p1<-ggplot(age, aes(x = Host.age, fill = disease)) +
31. geom_density(alpha = 0.5) +
32. labs(x = "Age(Before matching)") +
33. theme_minimal() +
34. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
35. ggtitle('C')+
36. theme(legend.position = 'none')+
37. theme(plot.title = element_text(face = 'bold'))+
38. theme(axis.text.y = element_blank(),
39. axis.title = element_blank(),
40. panel.grid = element_blank())
41. p1##p.val:6.265e-15
42. BMI<-rbind(meta_H[,c(7,8),drop=F],meta_D[,c(7,8),drop=F])
43. p2<-ggplot(BMI, aes(x = BMI, fill = disease)) +
44. geom_density(alpha = 0.5) +
45. labs(x = "BMI(Before matching)") +
46. theme_minimal() +
47. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
48. xlim(c(0,55))+
49. ggtitle('E')+
50. theme(plot.title = element_text(face = 'bold'))+
51. theme(legend.position = 'none')+
52. theme(axis.text.y = element_blank(),
53. axis.title = element_blank(),
54. panel.grid = element_blank())
55. p2##p.val:0.4529
56. meta_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_H_used.csv',row.names = 1)
57. meta_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_D_used.csv',row.names = 1)
58. colnames(meta_H)[5]<-'disease'
59. colnames(meta_D)[5]<-'disease'
60. meta_H$disease<-'Healthy'
61. meta_D$disease<-'Depression'
62. age<-rbind(meta_H[,c(3,5),drop=F],meta_D[,c(3,5),drop=F])
63. p3<-ggplot(age, aes(x = Host.age, fill = disease)) +
64. geom_density(alpha = 0.5) +
65. labs(x = "Age(After matching)") +
66. theme_minimal() +
67. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
68. ggtitle('D')+
69. theme(plot.title = element_text(face = 'bold'))+
70. theme(legend.position = 'none')+
71. theme(axis.text.y = element_blank(),
72. axis.title = element_blank(),
73. panel.grid = element_blank())
74. p3##p.val:0.5574
75. BMI<-rbind(meta_H[,c(4,5),drop=F],meta_D[,c(4,5),drop=F])
76. p4<-ggplot(BMI, aes(x = BMI, fill = disease)) +
77. geom_density(alpha = 0.5) +
78. labs(x = "BMI(After matching)") +
79. theme_minimal() +
80. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
81. ggtitle('F')+
82. theme(plot.title = element_text(face = 'bold'))+
83. theme(axis.text.y = element_blank(),
84. axis.title = element_blank(),
85. panel.grid = element_blank())
86. p4##p.val:0.8053
87. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig1.pdf',width = 15,height = 5.82)
88. grid.arrange(p1, p3,p2,p4,nrow = 1)
89. dev.off()
91. abundance_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/depression_genus.csv',row.names = 1)
92. abundance_D<-abundance_D[match(rownames(meta_D),colnames(abundance_D))]
93. temp<-abundance_D
94. temp[temp!=0]<-1
95. abundance_D<-abundance_D[which(rowSums(temp)>ncol(abundance_D)*0.1)%>%as.numeric(),]
96. abundance_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/healthy_genus.csv',row.names = 1)
97. abundance_H<-abundance_H[match(rownames(meta_H),colnames(abundance_H))]
98. temp<-abundance_H
99. temp[temp!=0]<-1
100. #abundance_H<-abundance_H[which(rowSums(temp)>ncol(abundance_H)*0.1)%>%as.numeric(),]
101. abundance<-merge(abundance_D,abundance_H,by = "row.names",all = TRUE)
102. abundance[is.na(abundance)]<-0
103. rownames(abundance)<-abundance$Row.names
104. abundance$Row.names<-NULL
105. #abundance_clr<-t(abundance/100)%>%clr()%>%as.matrix()%>%t()
106. meta<-rbind(meta_D,meta_H)
107. files<-list.files('D:/microbiome/micom/genus_models/genus_models/')
108. for(i in 1:length(files)){
109. files[i]<-substr(files[i],1,nchar(files[i])-4)
110. }
111. abundance1<-abundance_H[files,]
112. abundance1<-na.omit(abundance1)
113. colSums(abundance1)%>%hist()
114. ##########alpha
115. OTU<-phyloseq::otu_table(abundance,taxa_are_rows = T)
116. meta1<-phyloseq::sample_data(meta)
117. physeq<-phyloseq(OTU,meta1)
118. mpse<-as.MPSE(physeq)
119. mpse@assays@data@listData[["Abundance"]]<-mpse@assays@data@listData[["Abundance"]]*1e10
120. mpse %<>%
121. mp_decostand(.abundance=Abundance)
122. cols<-c('Depression'='#E7B800', 'Healthy'='#00AFBB')
123. mpse %<>% mp_cal_alpha(.abundance = Abundance,force = T)
124. p5 <- mpse %>%
125. mp_plot_alpha(
126. .alpha = c(Observe, Shannon,Chao1,Simpson, Pielou),
127. .group = disease,
128. ) +
129. scale_fill_manual(values=cols) +
130. scale_color_manual(values=cols) +
131. theme(
132. legend.position="none",
133. strip.background = element_rect(colour=NA, fill="grey")
134. )+
135. ggtitle('A')+
136. theme(plot.title = element_text(face = 'bold'))
137. p5
138. ##########anosim
139. mpse %>%
140. mp_anosim(.abundance=Abundance, .group=disease, action="get")
141. mpse %<>% mp_cal_dist(.abundance=hellinger, distmethod="bray")
142. p6 <- mpse %>% mp_plot_dist(.distmethod = bray, .group = disease, group.test=TRUE, textsize=2)+
143. ggtitle('C')+
144. theme(plot.title = element_text(face = 'bold'))
145. p6
146. ##########pcoa
147. dist_bray <- phyloseq::distance(physeq, method = "bray")
148. pcoa <- ordinate(physeq, method = "PCoA", distance = dist_bray)
149. pcoa_df <- data.frame(pcoa$vectors)
150. pcoa_df$SampleID <- rownames(pcoa_df)
151. sample_data<-cbind(rownames(meta),meta$disease)%>%as.data.frame()
152. colnames(sample_data)<-c('SampleID','disease')
153. rownames(sample_data)<-sample_data$SampleID
154. pcoa_df <- merge(pcoa_df, sample_data, by.x = "SampleID", by.y = "SampleID")
155. p7<-ggplot(pcoa_df, aes(x = Axis.1, y = Axis.2, color = disease)) +
156. geom_point(alpha = 0.4) +
157. labs(
158. x = paste0("PCOA1 (", round(pcoa[["values"]][["Eigenvalues"]][1], 2), "%)"),
159. y = paste0("PCOA2 (", round(pcoa[["values"]][["Eigenvalues"]][2], 2), "%)"),
160. color = "label"
161. ) +
162. theme_minimal()+
163. theme(legend.position = 'none')+
164. stat_ellipse(aes(color = disease),type = "t", level = 0.95)+
165. scale_color_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
166. annotate("text", x = -0.3, y = -0.3, label = "ANOSIM R: 0.007962
167. P: 0.004")+
168. ggtitle('B')+
169. theme(plot.title = element_text(face = 'bold'))
170. p7
171. p7<-ggMarginal(p7, type = "boxplot", margins = "both", size = 5, groupColour = TRUE, groupFill = TRUE)
172. p7
173. #abundance1<-t(abundance)%>%as.data.frame()
174. #env<-cbind(meta$Host.age,meta$BMI)%>%as.data.frame()
175. #colnames(env)<-c('Age','BMI')
176. #rownames(env)<-rownames(abundance1)
177. #cca <- cca(abundance1 ~ ., data = env)
178. #cca_summary<-summary(cca)
179. #cca_df<-cca_summary$sites%>%as.data.frame()
180. #cca_df$disease<-meta$disease[match(rownames(cca_df),rownames(meta))]
181. #scores <- as.data.frame(scores(cca)$biplot)
182. #p7<-ggplot(cca_df, aes(x = CCA1, y = CCA2, color = disease)) +
183. # geom_point(alpha = 0.4) +
184. # geom_segment(data = scores, aes(x = 0, xend = CCA1*20, y = 0, yend = CCA2*20),
185. # arrow = arrow(length = unit(0.2, "cm")), color = "red")+
186. # labs(
187. # x = paste0("CCA1 (", round(cca$CCA$eig[1] * 100, 2), "%)"),
188. # y = paste0("CCA2 (", round(cca$CCA$eig[2] * 100, 2), "%)"),
189. # color = "label"
190. # ) +
191. # theme_minimal()+
192. # theme(legend.position = 'none')+
193. # stat_ellipse(aes(color = disease),type = "t", level = 0.95)+
194. # scale_color_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
195. # annotate("text", x = -10, y = -20, label = "ANOSIM R: 0.007962
196. # P: 0.004")+
197. # ggtitle('B')
198. #p7<-ggMarginal(p7, type = "boxplot", margins = "both", size = 5, groupColour = TRUE, groupFill = TRUE)
199. #p7
201. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig2_1.pdf',width = 15,height = 4)
202. grid.arrange(p5,p7,ncol=2)
203. dev.off()
204. #########
205. #top5_D<-rownames(abundance_D)[order(rowSums(abundance_D),decreasing = T)[1:5]]
206. #df_D<-rowSums(abundance_D[top5_D,])%>%as.data.frame()
207. #colnames(df_D)<-'Depression'
208. #df_D<-rbind(df_D,abundance_D[-order(rowSums(abundance_D),decreasing = T)[1:5],]%>%rowSums()%>%sum())
209. #rownames(df_D)[6]<-'Others'
210. #top5_H<-rownames(abundance_H)[order(rowSums(abundance_H),decreasing = T)[1:5]]
211. #df_H<-rowSums(abundance_H[top5_H,])%>%as.data.frame()
212. #colnames(df_H)<-'Healthy'
213. #df_H<-rbind(df_H,abundance_H[-order(rowSums(abundance_H),decreasing = T)[1:5],]%>%rowSums()%>%sum())
214. #rownames(df_H)[6]<-'Others'
215. #df<-merge(df_D,df_H,by = "row.names",all = TRUE)
216. #rownames(df)<-df$Row.names
217. #df$Row.names<-NULL
218. #df[is.na(df)]<-0
219. #df$genus<-rownames(df)
220. #temp<-colSums(df[,1:2])%>%as.numeric()
221. #df[,1]<-df[,1]/temp[1]
222. #df[,2]<-df[,2]/temp[2]
223. #data_long <- df %>%
224. # tidyr::gather(key = "group", value = "value", -genus)
225. #data_alluvial <- data_long %>%
226. # dplyr::mutate(flow = rep(1:length(unique(data_long$genus)),2))
227. #desired_order<-c('Healthy','Depression')
228. #colnames(data_alluvial)[1]<-'Genus'
229. #p8<-ggplot(data_alluvial, aes(x = group, y = value, fill = Genus, stratum = Genus, alluvium = flow)) +
230. # geom_stratum(width = 0.4) +
231. # geom_flow(stat = "alluvium", lode.guidance = "forward", aes.flow = "backward") +
232. # theme_minimal() +
233. # labs(x = "", y = "abundance") +
234. # scale_fill_brewer(type = "qual", palette = "Set3") +
235. # scale_x_discrete(limits=desired_order)+
236. # theme(legend.text=element_text(face="italic"))+
237. # ggtitle('D')+
238. # theme(plot.title = element_text(face = 'bold'))
239. #p8
240. ############Wilcoxon signed rank test
241. abundance_D<-abundance[,match(rownames(meta_D),colnames(abundance))]
242. abundance_H<-abundance[,match(rownames(meta_H),colnames(abundance))]
243. res<-matrix(0,nrow(abundance),2)%>%as.data.frame()
244. rownames(res)<-rownames(abundance)
245. colnames(res)<-c('p.val','fdr')
246. for(i in 1:nrow(abundance)){
247. res[i,1]<-wilcox.test(abundance_D[i,]%>%as.numeric(),abundance_H[i,]%>%as.numeric(),paired = T)$p.val
248. }
249. res[,2]<-p.adjust(res[,1],method = 'BH')
250. diff_top50<-rownames(res)[order(res[,2],decreasing = F)[1:50]]
251. diff_top50
252. diff<-rownames(res)[which(res[,2]<0.05)]
253. #diff1<-rownames(res)[which(res[,2]>1e-5&res[,2]<0.05)]
254. #diff1
255. #exp<-abundance_clr[diff_top50,]
256. #label<-meta[,5,drop=F]
257. #ann_colors = list(disease=c(Healthy="#00AFBB",Depression="#E7B800"))
258. #col <- c(colorRampPalette(c("#4DBBD5FF", "white"))(abs(min(exp))*100),
259. # colorRampPalette('white')(1),
260. # colorRampPalette(c("white", "#F39B7FFF"))(max(exp)*100))
261. #p9<-pheatmap(exp,show_colnames = F,show_rownames = T,
262. # annotation_col = label,
263. # annotation_colors = ann_colors,
264. # color = col,
265. # cluster_cols = F,
266. # fontfamily="serif")
267. #p9
268. #pdf(file = 'diff_genus.pdf',width = 8,height = 15)
269. #pheatmap(exp,show_colnames = F,show_rownames = T,
270. # annotation_col = label,
271. # annotation_colors = ann_colors,
272. # color = col,
273. # cluster_cols = F,
274. # fontfamily="serif")
275. #dev.off()
276. #diff<-rownames(res)[which(res[,2]<0.005)]
277. #a<-c('Bifidobacterium','Blautia','Ruminococcus','Oscillibacter','Faecalibacterium','Corprococcus','Alistipes','Dialister')
278. #intersect(a,diff)
279. #diff_abundance<-abundance[diff,]
280. #top5_abundance<-rowSums(diff_abundance)%>%as.data.frame()
281. #top5_genus<-diff_abundance[rownames(top5_abundance)[order(top5_abundance[,1],decreasing = T)[1:5]],]
282. #top5_genus<-t(top5_genus)%>%as.data.frame()
283. #top5_genus$id<-rownames(top5_genus)
284. #top5_genus<-pivot_longer(top5_genus,
285. # cols = -id,
286. # names_to = "Genus",
287. # values_to = "Relative Abundance"
288. #)
289. #top5_genus$label<-meta$disease[match(top5_genus$id,rownames(meta))]
290. #p10<-ggplot(top5_genus, aes(x = `Relative Abundance`, y = Genus, fill = label)) +
291. # geom_boxplot(outlier.shape = NA) +
292. # theme_minimal() +
293. # scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
294. # theme(axis.text.y = element_text(face = "italic"))+
295. # ggtitle('C')
296. #p10
297. #pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig2.pdf',width = 12,height = 10)
298. #grid.arrange(arrangeGrob(p5, p7,ncol=2),p10, nrow = 2)
299. #dev.off()
300. #exp1<-abundance_clr[diff1,]
301. #col <- c(colorRampPalette(c("#4DBBD5FF", "white"))(abs(min(exp1))*100),
302. # colorRampPalette('white')(1),
303. # colorRampPalette(c("white", "#F39B7FFF"))(max(exp1)*100))
304. #p5<-pheatmap(exp1,show_colnames = F,show_rownames = F,
305. # annotation_col = label,
306. # annotation_colors = ann_colors,
307. # color = col)
308. ##########Maaslin2
309. #meta_D1<-meta_D
310. #meta_D1$disease[1:100]<-'Healthy'
311. #fit_data = Maaslin2(input_data = t(abundance_D),
312. # input_metadata = meta_D1,
313. # min_prevalence = 0,
314. # normalization = "NONE",
315. # max_significance = 1,
316. # output = "D:/microbiome/micom/202404/paper/result/1-1/Maaslin2_result/Depression",
317. # fixed_effects = c("Host.age",'BMI'))
318. #
319. #meta_H1<-meta_H
320. #meta_H1$disease[1:100]<-'Depression'
321. #fit_data = Maaslin2(input_data = t(abundance_H),
322. # input_metadata = meta_H1,
323. # min_prevalence = 0,
324. # normalization = "NONE",
325. # max_significance = 1,
326. # output = "D:/microbiome/micom/202404/paper/result/1-1/Maaslin2_result/Healthy",
327. # fixed_effects = c("Host.age",'BMI'))
328. ############
329. #coef_D<-read.table('D:/microbiome/micom/202404/paper/result/1-1/Maaslin2_result/Depression/all_results.tsv',header = T)
330. #coef_H<-read.table('D:/microbiome/micom/202404/paper/result/1-1/Maaslin2_result/Healthy/all_results.tsv',header = T)
331. #coef_D<-coef_D[which(coef_D$metadata=='Host.age'),]
332. #coef_H<-coef_H[which(coef_H$metadata=='Host.age'),]
333. #coef_D<-coef_D[match(intersect(coef_D$feature,coef_H$feature),coef_D$feature),]
334. #coef_H<-coef_H[match(intersect(coef_D$feature,coef_H$feature),coef_H$feature),]
335. #coef<-cbind(coef_D$coef,coef_H$coef)%>%as.data.frame()
336. #colnames(coef)<-c('Depression','Healthy')
337. #rownames(coef)<-coef_D$feature
338. #coef$index<-abs(coef[,1]-coef[,2])
339. #genus<-c('Ruminococcus','Lactonifactor')
340. ######
341. #a<-cor(t(abundance_D),meta_D$Host.age)
342. #b<-cor(t(abundance_H),meta_H$Host.age)
343. ###########
344. Reactionabundance<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/Results/ReactionAbundance.csv',row.names = 1)
345. #Reactionpresence<-Reactionabundance
346. #Reactionpresence[Reactionpresence!=0]<-1
347. #Subsystemabundance<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/Results/SubsystemAbundance.csv',row.names = 1)
348. #Subsystempresence<-Subsystemabundance
349. #Subsystempresence[Subsystempresence!=0]<-1
350. #genuspresence<-abundance
351. #genuspresence[genuspresence!=0]<-1
352. #df<-rbind(genuspresence,Reactionpresence,Subsystempresence)
353. #df<-t(df)
354. #label_col<-c(rep('Genus',nrow(genuspresence)),
355. # rep('Reactions',nrow(Reactionpresence)),
356. # rep('Subsystem',nrow(Subsystempresence)))%>%as.data.frame()
357. #colnames(label_col)<-'label'
358. #rownames(label_col)<-colnames(df)
359. #label_row<-meta[,5,drop=F]
360. #col <-c("white", "#4DBBD5FF")
361. #p11<-pheatmap(df,show_colnames = F,show_rownames = F,
362. # annotation_row = label_row,
363. # annotation_col = label_col,
364. # annotation_colors = ann_colors,
365. # color = col,
366. # cluster_rows = F,
367. # cluster_cols = F,
368. # legend = F)
369. #p11
370. ########Wilcoxon signed rank test (Reaction)
371. res_reaction<-matrix(0,nrow(Reactionabundance),2)
372. colnames(res_reaction)<-c('p.val','fdr')
373. rownames(res_reaction)<-rownames(Reactionabundance)
374. Reactionabundance_D<-Reactionabundance[,match(rownames(meta_D),colnames(Reactionabundance))]
375. Reactionabundance_H<-Reactionabundance[,match(rownames(meta_H),colnames(Reactionabundance))]
376. for(i in 1:nrow(Reactionabundance)){
377. print(i)
378. res_reaction[i,1]<-wilcox.test(Reactionabundance_D[i,]%>%as.numeric(),Reactionabundance_H[i,]%>%as.numeric(),paired = T)$p.val
379. }
380. res_reaction[,2]<-p.adjust(res_reaction[,1],method = 'BH')
381. diff_reaction<-rownames(res_reaction)[which(res_reaction[,2]<0.05)]
382. diff_reaction
383. df<-t(Reactionabundance[rownames(res_reaction)[order(res_reaction[,1],decreasing = F)[1:20]],,drop=F])%>%as.data.frame()
384. df$label<-meta$disease[match(rownames(df),rownames(meta))]
385. df1<-pivot_longer(df,
386. cols = -label,
387. names_to = "reaction",
388. values_to = "value"
389. )
390. p12<-ggplot(df1, aes(x = value, y = reaction, fill = label)) +
391. geom_boxplot(outlier.shape = NA) +
392. theme_minimal() +
393. scale_x_break(c(0.1,0.99), scales = 0.5,ticklabels=c(0.99,0.995,1))+
394. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
395. scale_x_continuous(labels = function(x) sprintf("%.2f", x))+
396. ggtitle('D')+
397. theme(plot.title = element_text(face = 'bold'))
398. p12
401. #pdf(file = 'diff_reactions.pdf',width = 8,height = 15)
402. #pheatmap(exp,show_colnames = F,show_rownames = T,
403. # annotation_col = label,
404. # annotation_colors = ann_colors,
405. # color = col,
406. # cluster_cols = F)
407. #dev.off()
408. ########Wilcoxon signed rank test (Subsystem)
409. #res_Subsystem<-matrix(0,nrow(Subsystemabundance),2)
410. #colnames(res_Subsystem)<-c('p.val','fdr')
411. #rownames(res_Subsystem)<-rownames(Subsystemabundance)
412. #Subsystemabundance_D<-Subsystemabundance[,match(rownames(meta_D),colnames(Subsystemabundance))]
413. #Subsystemabundance_H<-Subsystemabundance[,match(rownames(meta_H),colnames(Subsystemabundance))]
414. #for(i in 1:nrow(Subsystemabundance)){
415. # print(i)
416. # res_Subsystem[i,1]<-wilcox.test(Subsystemabundance_D[i,]%>%as.numeric(),Subsystemabundance_H[i,]%>%as.numeric(),paired = T)$p.val
417. #}
418. #res_Subsystem[,2]<-p.adjust(res_Subsystem[,1],method = 'BH')
419. #diff_Subsystem<-rownames(res_Subsystem)[which(res_Subsystem[,1]<0.05)]
420. #diff_Subsystem
421. #exp<-Subsystemabundance[diff_Subsystem,]
422. #col <- c(colorRampPalette('white')(1),
423. # colorRampPalette(c("white", "#F39B7FFF"))(max(exp)*100))
424. #p13<-pheatmap(exp,show_colnames = F,show_rownames = T,
425. # annotation_col = label_row,
426. # annotation_colors = ann_colors,
427. # color = col,
428. # cluster_cols = F)
429. #p13
430. ############alpha/beta(reaction)
431. temp<-colSums(Reactionabundance)%>%as.numeric()
432. for(i in 1:ncol(Reactionabundance)){
433. Reactionabundance[,i]<-Reactionabundance[,i]/temp[i]
434. }
435. OTU<-phyloseq::otu_table(Reactionabundance,taxa_are_rows = T)
436. #rownames(meta)<-meta$sampleid
437. meta1<-phyloseq::sample_data(meta)
438. physeq<-phyloseq(OTU,meta1)
439. mpse<-as.MPSE(physeq)
440. mpse@assays@data@listData[["Abundance"]]<-mpse@assays@data@listData[["Abundance"]]*1e20
441. mpse %<>%
442. mp_decostand(.abundance=Abundance)
443. cols<-c('Depression'='#E7B800', 'Healthy'='#00AFBB')
444. mpse %<>% mp_cal_alpha(.abundance = Abundance,force = T)
445. p14 <- mpse %>%
446. mp_plot_alpha(
447. .alpha = c(Observe, Shannon,Chao1,Simpson, Pielou),
448. .group = disease,
449. ) +
450. scale_fill_manual(values=cols) +
451. scale_color_manual(values=cols) +
452. theme(
453. legend.position="none",
454. strip.background = element_rect(colour=NA, fill="grey")
455. )+
456. ggtitle('A')+
457. theme(plot.title = element_text(face = 'bold'))
459. mpse %>%
460. mp_anosim(.abundance=Abundance, .group=disease, action="get")
462. p14
464. dist_bray <- phyloseq::distance(physeq, method = "bray")
465. pcoa <- ordinate(physeq, method = "PCoA", distance = dist_bray)
466. pcoa_df <- data.frame(pcoa$vectors)
467. pcoa_df$SampleID <- rownames(pcoa_df)
468. sample_data<-cbind(rownames(meta),meta$disease)%>%as.data.frame()
469. colnames(sample_data)<-c('SampleID','disease')
470. rownames(sample_data)<-sample_data$SampleID
471. pcoa_df <- merge(pcoa_df, sample_data, by.x = "SampleID", by.y = "SampleID")
472. p15<-ggplot(pcoa_df, aes(x = Axis.1, y = Axis.2, color = disease)) +
473. geom_point(alpha = 0.4) +
474. labs(
475. x = paste0("PCOA1 (", round(pcoa[["values"]][["Eigenvalues"]][1], 2), "%)"),
476. y = paste0("PCOA2 (", round(pcoa[["values"]][["Eigenvalues"]][2], 2), "%)"),
477. color = "label"
478. ) +
479. theme_minimal()+
480. theme(legend.position = 'none')+
481. stat_ellipse(aes(color = disease),type = "t", level = 0.95)+
482. scale_color_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
483. annotate("text", x = -0.1, y = -0.2, label = "ANOSIM R: -0.0003651
484. P: 0.493")+
485. ggtitle('B')+
486. theme(plot.title = element_text(face = 'bold'))
487. p15
488. p15<-ggMarginal(p15, type = "boxplot", margins = "both", size = 5, groupColour = TRUE, groupFill = TRUE)
489. p15
490. #abundance1<-t(Reactionabundance)%>%as.data.frame()
491. #env<-cbind(meta$Host.age,meta$BMI)%>%as.data.frame()
492. #colnames(env)<-c('Age','BMI')
493. #rownames(env)<-rownames(abundance1)
494. #cca <- cca(abundance1 ~ ., data = env)
495. #cca_summary<-summary(cca)
496. #cca_df<-cca_summary$sites%>%as.data.frame()
497. #cca_df$disease<-meta$disease[match(rownames(cca_df),rownames(meta))]
498. #scores <- as.data.frame(scores(cca)$biplot)
499. #p15<-ggplot(cca_df, aes(x = CCA1, y = CCA2, color = disease)) +
500. # geom_point(alpha = 0.4) +
501. # geom_segment(data = scores, aes(x = 0, xend = CCA1*20, y = 0, yend = CCA2*20),
502. # arrow = arrow(length = unit(0.2, "cm")), color = "red")+
503. # labs(
504. # x = paste0("CCA1 (", round(cca$CCA$eig[1] * 100, 2), "%)"),
505. # y = paste0("CCA2 (", round(cca$CCA$eig[2] * 100, 2), "%)"),
506. # color = "label"
507. # ) +
508. # theme_minimal()+
509. # theme(legend.position = 'none')+
510. # stat_ellipse(aes(color = disease),type = "t", level = 0.95)+
511. # scale_color_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
512. # annotate("text", x = -10, y = -20, label = "ANOSIM R: 0.007962
513. # P: 0.004")+
514. # ggtitle('B')
515. #p15<-ggMarginal(p15, type = "boxplot", margins = "both", size = 5, groupColour = TRUE, groupFill = TRUE)
516. #p15
517. #########
518. fc_reaction<-rowMeans(Reactionabundance_D)/rowMeans(Reactionabundance_H)%>%as.data.frame()
519. res_reaction<-cbind(res_reaction,fc_reaction)
520. colnames(res_reaction)[3]<-'log2fc'
521. res_reaction[which(res_reaction$log2fc!=0),3]<-log2(res_reaction[which(res_reaction$log2fc!=0),3])
522. res_reaction$`-log10p.val`<-(-log10(res_reaction$p.val))
524. cut_off_log2FC<-1
525. cut_off_p<-0.05
526. res_reaction$sig<-ifelse(res_reaction$p.val < cut_off_p &
527. abs(res_reaction$log2fc) >= cut_off_log2FC,
528. ifelse(res_reaction$log2fc > cut_off_log2FC ,'Up','Down'),'No significance')
529. p16<-ggplot(res_reaction, aes(x =log2fc, y=`-log10p.val`,colour = sig)) +
530. geom_point(alpha=0.65, size=2) +
531. scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757")) +
532. geom_vline(xintercept=c(-cut_off_log2FC,cut_off_log2FC),lty=4,col="black",lwd=0.8) +
533. geom_hline(yintercept = -log10(cut_off_p), lty=4,col="black",lwd=0.8) +
534. labs(x="log2FC", y="-log10(p.val)") +
535. theme_bw() +
536. theme(plot.title = element_text(hjust = 0),
537. legend.position="right",
538. legend.title = element_blank()
539. )+
540. ggtitle('C')+
541. theme(plot.title = element_text(face = 'bold'))
542. p16
544. ############alpha/beta(Subsystemabundance)
545. #temp<-colSums(Subsystemabundance)%>%as.numeric()
546. #for(i in 1:ncol(Subsystemabundance)){
547. # Subsystemabundance[,i]<-Subsystemabundance[,i]/temp[i]
548. #}
549. #OTU<-phyloseq::otu_table(Subsystemabundance,taxa_are_rows = T)
550. #meta1<-phyloseq::sample_data(meta)
551. #physeq<-phyloseq(OTU,meta1)
552. #mpse<-as.MPSE(physeq)
553. #mpse@assays@data@listData[["Abundance"]]<-mpse@assays@data@listData[["Abundance"]]*1e20
554. #mpse %<>%
555. # mp_decostand(.abundance=Abundance)
556. #cols<-c('Depression'='#E7B800', 'Healthy'='#00AFBB')
557. #mpse %<>% mp_cal_alpha(.abundance = Abundance,force = T)
558. #p16 <- mpse %>%
559. # mp_plot_alpha(
560. # .alpha = c(Observe, Shannon,Chao1,Simpson, Pielou),
561. # .group = disease,
562. # ) +
563. # scale_fill_manual(values=cols) +
564. # scale_color_manual(values=cols) +
565. # theme(
566. # legend.position="none",
567. # strip.background = element_rect(colour=NA, fill="grey")
568. # )
569. #p16
570. #dist_bray <- phyloseq::distance(physeq, method = "bray")
571. #pcoa <- ordinate(physeq, method = "PCoA", distance = dist_bray)
572. #pcoa_df <- data.frame(pcoa$vectors)
573. #pcoa_df$SampleID <- rownames(pcoa_df)
574. #sample_data<-cbind(rownames(meta),meta$disease)%>%as.data.frame()
575. #colnames(sample_data)<-c('SampleID','disease')
576. #rownames(sample_data)<-sample_data$SampleID
577. #pcoa_df <- merge(pcoa_df, sample_data, by.x = "SampleID", by.y = "SampleID")
578. #p17<-ggplot(pcoa_df, aes(x = Axis.1, y = Axis.2, color = disease)) +
579. # geom_point(size = 2) +
580. # labs(x = paste0("PCoA1 (", round(pcoa$values$Relative_eig[1] * 100, 2), "%)"),
581. # y = paste0("PCoA2 (", round(pcoa$values$Relative_eig[2] * 100, 2), "%)")) +
582. # theme_minimal() +
583. # stat_ellipse(aes(color = disease),type = "t", level = 0.95) +
584. # scale_color_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))
585. #p17<-ggMarginal(p17, type = "boxplot", margins = "both", size = 5, groupColour = TRUE, groupFill = TRUE)
586. #p17
587. ###############
588. #abundance_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/healthy_genus.csv',row.names = 1)
589. #files<-list.files('D:/microbiome/micom/genus_models/genus_models')
590. #for(i in 1:length(files)){
591. # files[i]<-substr(files[i],1,nchar(files[i])-4)
592. #}
593. #abundance_H<-abundance_H[intersect(rownames(abundance_H),files),]
594. #sampleids<-list.files('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/healthy')
595. #for(i in 1:length(sampleids)){
596. # sampleids[i]<-substr(sampleids[i],1,nchar(sampleids[i])-4)
597. #}
598. #for(sampleid in sampleids){
599. # print(sampleid)
600. # flux_file<-paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/healthy/',sampleid,'.csv')
601. # flux<-read.csv(flux_file,row.names = 1)
602. # flux<-flux[which(flux[,1]!=0),,drop=F]
603. # temp<-abundance_H[,match(colnames(flux),colnames(abundance_H)),drop=F]
604. # temp<-temp[which(temp[,1]!=0),,drop=F]
605. # flux$genus<-0
606. # flux$metabolite<-0
607. # for(i in 1:nrow(temp)){
608. # temp1<-grep(rownames(temp)[i],rownames(flux))
609. # flux[temp1,2]<-rownames(temp)[i]
610. # }
611. # for(i in 1:nrow(flux)){
612. # if(flux$genus[i]!='0'){
613. # flux$metabolite[i]<-substr(rownames(flux)[i],nchar(flux$genus[i])+1,nchar(rownames(flux)[i]))
614. # }else{
615. # flux$metabolite[i]<-rownames(flux)[i]
616. # }
617. # }
618. # flux$genus[which(flux$genus=='0')]<-'thylakoid'
619. # metabolite<-unique(flux$metabolite)
620. # genus<-unique(flux$genus)
621. # temp3<-matrix(0,length(metabolite),1)%>%as.data.frame()
622. # rownames(temp3)<-metabolite
623. # temp_a<-abundance_H[,sampleid,drop=F]
624. # temp_a<-temp_a/colSums(temp_a)
625. # for(i in 1:length(genus)){
626. # if(genus[i]!='thylakoid'){
627. # flux[which(flux$genus==genus[i]),1]<-flux[which(flux$genus==genus[i]),1]*temp_a[genus[i],]
628. # }
629. # }
630. # for(i in 1:length(metabolite)){
631. # temp3[i,1]<-flux[which(flux$metabolite==metabolite[i]),1]%>%sum()
632. # }
633. # colnames(temp3)<-sampleid
634. # output_file<-paste0('D:/microbiome/micom/202404/paper/result/1-1//FBA_result_p/healthy/',sampleid,'.csv')
635. # write.csv(temp3,output_file)
636. #}
638. #temp4<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result_p/healthy/',sampleids[1],'.csv'),row.names = 1)
639. #temp5<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result_p/healthy/',sampleids[2],'.csv'),row.names = 1)
640. #merged_df <- merge(temp4, temp5, by = "row.names", all = TRUE)
641. #merged_df[is.na(merged_df)] <- 0
642. #rownames(merged_df) <- merged_df$Row.names
643. #merged_df$Row.names <- NULL
644. #for(i in 3:length(sampleids)){
645. # print(i)
646. # temp6<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result_p/healthy/',sampleids[i],'.csv'),row.names = 1)
647. # merged_df <- merge(merged_df, temp6, by = "row.names", all = TRUE)
648. # merged_df[is.na(merged_df)] <- 0
649. # rownames(merged_df) <- merged_df$Row.names
650. # merged_df$Row.names <- NULL
651. #}
652. #write.csv(merged_df,'D:/microbiome/micom/202404/paper/result/1-1/metabolite_H.csv')
653. ########
654. #abundance_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/depression_genus.csv',row.names = 1)
655. #files<-list.files('D:/microbiome/micom/genus_models/genus_models')
656. #for(i in 1:length(files)){
657. # files[i]<-substr(files[i],1,nchar(files[i])-4)
658. #}
659. #abundance_D<-abundance_D[intersect(rownames(abundance_D),files),]
660. #sampleids<-list.files('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression')
661. #for(i in 1:length(sampleids)){
662. # sampleids[i]<-substr(sampleids[i],1,nchar(sampleids[i])-4)
663. #}
664. #for(sampleid in sampleids){
665. # print(sampleid)
666. # flux_file<-paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleid,'.csv')
667. # flux<-read.csv(flux_file,row.names = 1)
668. # flux<-flux[which(flux[,1]!=0),,drop=F]
669. # temp<-abundance_D[,match(colnames(flux),colnames(abundance_D)),drop=F]
670. # temp<-temp[which(temp[,1]!=0),,drop=F]
671. # flux$genus<-0
672. # flux$metabolite<-0
673. # for(i in 1:nrow(temp)){
674. # temp1<-grep(rownames(temp)[i],rownames(flux))
675. # flux[temp1,2]<-rownames(temp)[i]
676. # }
677. # for(i in 1:nrow(flux)){
678. # if(flux$genus[i]!='0'){
679. # flux$metabolite[i]<-substr(rownames(flux)[i],nchar(flux$genus[i])+1,nchar(rownames(flux)[i]))
680. # }else{
681. # flux$metabolite[i]<-rownames(flux)[i]
682. # }
683. # }
684. # flux$genus[which(flux$genus=='0')]<-'thylakoid'
685. # metabolite<-unique(flux$metabolite)
686. # genus<-unique(flux$genus)
687. # temp3<-matrix(0,length(metabolite),1)%>%as.data.frame()
688. # rownames(temp3)<-metabolite
689. # temp_a<-abundance_D[,sampleid,drop=F]
690. # temp_a<-temp_a/colSums(temp_a)
691. # for(i in 1:length(genus)){
692. # if(genus[i]!='thylakoid'){
693. # flux[which(flux$genus==genus[i]),1]<-flux[which(flux$genus==genus[i]),1]*temp_a[genus[i],]
694. # }
695. # }
696. # for(i in 1:length(metabolite)){
697. # temp3[i,1]<-flux[which(flux$metabolite==metabolite[i]),1]%>%sum()
698. # }
699. # colnames(temp3)<-sampleid
700. # output_file<-paste0('D:/microbiome/micom/202404/paper/result/1-1//FBA_result_p/depression/',sampleid,'.csv')
701. # write.csv(temp3,output_file)
702. #}
703. #
704. #temp4<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result_p/depression/',sampleids[1],'.csv'),row.names = 1)
705. #temp5<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result_p/depression/',sampleids[2],'.csv'),row.names = 1)
706. #merged_df <- merge(temp4, temp5, by = "row.names", all = TRUE)
707. #merged_df[is.na(merged_df)] <- 0
708. #rownames(merged_df) <- merged_df$Row.names
709. #merged_df$Row.names <- NULL
710. #for(i in 3:length(sampleids)){
711. # print(i)
712. # temp6<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result_p/depression/',sampleids[i],'.csv'),row.names = 1)
713. # merged_df <- merge(merged_df, temp6, by = "row.names", all = TRUE)
714. # merged_df[is.na(merged_df)] <- 0
715. # rownames(merged_df) <- merged_df$Row.names
716. # merged_df$Row.names <- NULL
717. #}
718. #write.csv(merged_df,'D:/microbiome/micom/202404/paper/result/1-1/metabolite_D.csv')
721. ########
722. metabolite_D<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/metabolite_D.csv',row.names = 1)
723. metabolite_H<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/metabolite_H.csv',row.names = 1)
724. metabolite_D<-metabolite_D[intersect(rownames(metabolite_D),rownames(metabolite_H)),]
725. metabolite_H<-metabolite_H[intersect(rownames(metabolite_D),rownames(metabolite_H)),]
726. metabolite_D<-metabolite_D[,match(rownames(meta_D),colnames(metabolite_D))]
727. metabolite_H<-metabolite_H[,match(rownames(meta_H),colnames(metabolite_H))]
728. metabolite<-cbind(metabolite_D,metabolite_H)
729. meta_D<-read.table('D:/microbiome/micom/202404/paper/data/meta_Depression.txt',sep = '\t',header = T)
730. meta_H<-read.table('D:/microbiome/micom/202404/paper/data/meta_healthy.txt',sep = '\t',header = T)
731. colnames(meta_H)[8]<-'disease'
732. colnames(meta_D)[8]<-'disease'
733. meta_H$disease<-'Healthy'
734. meta_D$disease<-'Depression'
735. meta<-rbind(meta_D,meta_H)
737. res_metabolite<-matrix(0,nrow(metabolite_D),2)%>%as.data.frame()
738. rownames(res_metabolite)<-rownames(metabolite_D)
739. colnames(res_metabolite)<-c('p.val','fdr')
740. for(i in 1:nrow(metabolite_D)){
741. print(i)
742. res_metabolite[i,1]<-wilcox.test(metabolite_D[i,]%>%as.numeric(),metabolite_H[i,]%>%as.numeric(),paired = T)$p.val
743. }
744. res_metabolite[,2]<-p.adjust(res_metabolite[,1],method = 'BH')
745. diff_metabolite<-rownames(res_metabolite)[which(res_metabolite[,2]<0.05)]
746. diff_metabolite
747. df<-t(metabolite[rownames(res_metabolite)[which(res_metabolite[,2]<0.05)],,drop=F])%>%as.data.frame()
748. df<-scale(df,scale = T,center = T)%>%as.data.frame()
749. metabolite_all<-fread('D:/microbiome/micom/202404/paper/data/MetaboliteDatabase.txt')%>%as.data.frame()
750. compartment<-c()
751. for(i in 1:ncol(df)){
752. compartment<-c(compartment,substr(colnames(df)[i],nchar(colnames(df)[i])-2,nchar(colnames(df)[i])))
753. }
754. name<-c()
755. for(i in 1:ncol(df)){
756. name<-c(name,substr(colnames(df)[i],1,nchar(colnames(df)[i])-3))
757. }
758. name1<-metabolite_all[match(name,metabolite_all[,1]),2]
759. name1<-paste0(name1,compartment)
760. colnames(df)<-name1
761. df$label<-meta$disease[match(rownames(df),meta$Run.ID)]
764. df1<-pivot_longer(df,
765. cols = -label,
766. names_to = "metabolite",
767. values_to = "flux"
768. )
770. #p<-list()
771. #for(i in 1:(ncol(df)-1)){
772. # print(i)
773. # temp<-df1[which(df1$metabolite==colnames(df)[i]),]
774. # p[[i]]<-ggplot(temp, aes(x = label, y = flux, fill = label)) +
775. # geom_boxplot(outlier.shape = NA) +
776. # theme_minimal() +
777. # scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
778. # coord_cartesian(ylim=c(quantile(temp$flux,0.1),quantile(temp$flux,0.9)+0.2))+
779. #
780. # ggtitle(colnames(df)[i])
781. #}
782. p18<-ggplot(df1, aes(x = flux, y = metabolite, fill = label)) +
783. geom_boxplot(outlier.alpha=0.2) +
784. theme_minimal() +
785. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
786. #xlim(c(-5,5))+
787. labs(x='flux\n mmol/[gDW h](scaled)',y='')+
788. theme(axis.text.y = element_text(size = 12),
789. axis.title.y = element_text(size = 12),
790. axis.text.x = element_text(size = 12),
791. axis.title.x = element_text(size = 12),
792. legend.title = element_blank(),
793. legend.position = 'right',
794. plot.title = element_text(hjust = -2.5,face = 'bold'))+
795. ggtitle('A')
797. p18
798. ###############
799. allmetabolite<-fread('D:/microbiome/micom/202404/paper/data/MetaboliteDatabase.txt')%>%as.data.frame()
800. diff_metabolite<-strsplit(diff_metabolite,'\\[')%>%unlist
801. diff_metabolite<-diff_metabolite[-grep('\\]',diff_metabolite)]%>%unique
802. id<-allmetabolite$V6[match(diff_metabolite,allmetabolite$V1)]
803. id<-id[id!='']
804. id
805. write.table(id,'D:/microbiome/micom/202404/paper/result/1-1/keggid_diff.txt',quote = F,sep = '\t',row.names = F,col.names = F)
807. stitch<-fread('D:/microbiome/micom/202404/paper/result/1-1/stitch_interactions (1).tsv')%>%as.data.frame()
808. protein<-c(stitch$`#node1`,
809. stitch$node2)
810. names(protein)<-c(stitch$node1_external_id,stitch$node2_external_id)
811. protein<-protein[grep('ENSP',names(protein))]
812. protein<-unique(protein)
813. write.table(protein,'D:/microbiome/micom/202404/paper/result/1-1/protein.txt',quote = F,sep = '\t',row.names = F,col.names = F)
815. protein<-read.table('D:/microbiome/micom/202404/paper/result/1-1/protein.txt',sep='\t')
816. gene_list<-protein$V1
817. library(httr)
818. library(jsonlite)
819. library(stringr)
820. get_uniprot_id <- function(gene) {
821. url <- paste0("https://rest.uniprot.org/uniprotkb/search?query=gene_exact:",
822. gene, "+AND+organism_id:9606+AND+reviewed:true&format=json&size=1")
823. res <- GET(url)
824. if (res$status_code != 200) return(NA)
825. data <- fromJSON(content(res, "text", encoding = "UTF-8"))
826. if (length(data$results) == 0) return(NA)
827. return(data[["results"]][["primaryAccession"]])
828. }
829. get_fasta <- function(uniprot_id) {
830. url <- paste0("https://rest.uniprot.org/uniprotkb/", uniprot_id, ".fasta")
831. res <- GET(url)
832. if (res$status_code != 200) return(NULL)
833. return(content(res, "text", encoding = "UTF-8"))
834. }
835. fasta_all <- ""
836. for (gene in gene_list) {
837. cat("Processing:", gene, "\n")
838. uid <- get_uniprot_id(gene)
839. if (is.na(uid)) {
840. warning(paste("No UniProt ID found for", gene))
841. next
842. }
843. fasta_seq <- get_fasta(uid)
844. if (is.null(fasta_seq)) {
845. warning(paste("Failed to retrieve FASTA for", uid))
846. next
847. }
848. fasta_all <- paste0(fasta_all, fasta_seq, "\n")
849. }
851. # 5. ?????? FASTA ??????
852. writeLines(fasta_all, "output_sequences.fasta")
853. ##############
854. #metabolite_used_D<-metabolite_D[diff_metabolite,]%>%t()%>%as.data.frame()
855. #metabolite_used_H<-metabolite_H[diff_metabolite,]%>%t()%>%as.data.frame()
856. #metabolite_used_D$label<-'Depression'
857. #metabolite_used_H$label<-'Healthy'
858. #set.seed(123)
859. #temp<-sample(1:nrow(metabolite_used_D),0.7*nrow(metabolite_used_D))
860. #train_D<-metabolite_used_D[temp,]
861. #train_H<-metabolite_used_H[temp,]
862. #data_train<-rbind(train_D,train_H)
863. #test_D<-metabolite_used_D[-temp,]
864. #test_H<-metabolite_used_H[-temp,]
865. #data_test<-rbind(test_D,test_H)
866. #train_index <- createFolds(data_train$label, k = 10)
867. #randomForestFit <- data_train %>% train(label ~ .,
868. # method = "rf",
869. # data = .,
870. # tuneLength = 5,
871. # trControl = trainControl(method = "cv", indexOut = train_index))
872. #randomForestFit
873. #pr_rf<-predict(randomForestFit,data_test,type = 'prob')
874. #library(pROC)
875. #roc_curve<-roc(data_test$label,pr_rf[,1])
876. #plot(roc_curve)
877. #sp.obj <- ci.sp(roc_curve, sensitivities=seq(0, 1, .01), boot.n=100)
878. #plot(sp.obj, type="shape", col="#5bd1d7")
879. #text(0,0.5,paste0("AUC:",round(roc_curve[["auc"]],4)))
890. #pdf(file = 'diff_metabolite.pdf',width = 10,height = 5)
891. #ggplot(df1, aes(x = metabolite, y = value, fill = label)) +
892. # geom_boxplot() +
893. # theme_minimal() +
894. # scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))
895. #dev.off()
896. #######34dhphe[e]
897. #sampleids<-list.files('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression')
898. #for(i in 1:length(sampleids)){
899. # sampleids[i]<-substr(sampleids[i],1,nchar(sampleids[i])-4)
900. #}
901. #temp7<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleids[1],'.csv'),row.names = 1)
902. #temp7<-temp7[grep('34dhphe\\[e\\]',rownames(temp7)),,drop=F]
903. #temp8<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleids[2],'.csv'),row.names = 1)
904. #temp8<-temp8[grep('34dhphe\\[e\\]',rownames(temp8)),,drop=F]
905. #
906. #merged_df <- merge(temp7, temp8, by = "row.names", all = TRUE)
907. #merged_df[is.na(merged_df)] <- 0
908. #rownames(merged_df) <- merged_df$Row.names
909. #merged_df$Row.names <- NULL
910. #for(i in 3:length(sampleids)){
911. # print(i)
912. # temp9<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleids[i],'.csv'),row.names = 1)
913. # temp9<-temp9[grep('34dhphe\\[e\\]',rownames(temp9)),,drop=F]
914. # merged_df <- merge(merged_df, temp9, by = "row.names", all = TRUE)
915. # merged_df[is.na(merged_df)] <- 0
916. # rownames(merged_df) <- merged_df$Row.names
917. # merged_df$Row.names <- NULL
918. #}
919. #write.csv(merged_df,'D:/microbiome/micom/202404/paper/result/1-1/34dhphe[e]_D.csv')
920. #######
921. #L_Dopa_D<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/34dhphe[e]_D.csv',row.names = 1)
922. #L_Dopa_H<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/34dhphe[e]_H.csv',row.names = 1)
923. #for(i in 1:nrow(L_Dopa_D)){
924. # rownames(L_Dopa_D)[i]<-substr(rownames(L_Dopa_D)[i],1,nchar(rownames(L_Dopa_D)[i])-10)
925. #}
926. #for(i in 1:nrow(L_Dopa_H)){
927. # rownames(L_Dopa_H)[i]<-substr(rownames(L_Dopa_H)[i],1,nchar(rownames(L_Dopa_H)[i])-10)
928. #}
929. #L_Dopa_H<-abs(L_Dopa_H)
930. #L_Dopa_D<-abs(L_Dopa_D)
931. #top20_D<-rownames(L_Dopa_D)[order(rowSums(L_Dopa_D),decreasing = T)[1:20]]
932. #df_D<-rowSums(L_Dopa_D[top20_D,])%>%as.data.frame()
933. #colnames(df_D)<-'Depression'
934. #df_D<-rbind(df_D,L_Dopa_D[-order(rowSums(L_Dopa_D),decreasing = T)[1:20],]%>%rowSums()%>%sum())
935. #rownames(df_D)[21]<-'Others'
936. #top20_H<-rownames(L_Dopa_H)[order(rowSums(L_Dopa_H),decreasing = T)[1:20]]
937. #df_H<-rowSums(L_Dopa_H[top20_H,])%>%as.data.frame()
938. #colnames(df_H)<-'Healthy'
939. #df_H<-rbind(df_H,L_Dopa_H[-order(rowSums(L_Dopa_H),decreasing = T)[1:20],]%>%rowSums()%>%sum())
940. #rownames(df_H)[21]<-'Others'
941. #df<-merge(df_D,df_H,by = "row.names",all = TRUE)
942. #rownames(df)<-df$Row.names
943. #df$Row.names<-NULL
944. #df[is.na(df)]<-0
945. #df$genus<-rownames(df)
946. ##temp<-colSums(df[,1:2])%>%as.numeric()
947. ##df[,1]<-df[,1]/temp[1]
948. ##df[,2]<-df[,2]/temp[2]
949. #data_long <- df %>%
950. # tidyr::gather(key = "group", value = "value", -genus)
951. #data_alluvial <- data_long %>%
952. # dplyr::mutate(flow = rep(1:length(unique(data_long$genus)),2))
953. #desired_order<-c('Healthy','Depression')
954. #col<-rgb(runif(21), runif(21), runif(21))
955. #p19<-ggplot(data_alluvial, aes(x = group, y = value, fill = genus, stratum = genus, alluvium = flow)) +
956. # geom_stratum(width = 0.4) +
957. # geom_flow(stat = "alluvium", lode.guidance = "forward", aes.flow = "backward") +
958. # theme_minimal() +
959. # labs(x = "", y = "L-Dopa[e]") +
960. # scale_x_discrete(limits=desired_order)+
961. # ggtitle('B')
962. #p19
963. #
964. #
965. #tmao_D<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/tmao[e]_D.csv',row.names = 1)
966. #tmao_H<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/tmao[e]_H.csv',row.names = 1)
967. #for(i in 1:nrow(tmao_D)){
968. # rownames(tmao_D)[i]<-substr(rownames(tmao_D)[i],1,nchar(rownames(tmao_D)[i])-7)
969. #}
970. #for(i in 1:nrow(tmao_H)){
971. # rownames(tmao_H)[i]<-substr(rownames(tmao_H)[i],1,nchar(rownames(tmao_H)[i])-7)
972. #}
973. #tmao_H<-abs(tmao_H)
974. #tmao_D<-abs(tmao_D)
975. #top20_D<-rownames(tmao_D)[order(rowSums(tmao_D),decreasing = T)[1:20]]
976. #df_D<-rowSums(tmao_D[top20_D,])%>%as.data.frame()
977. #colnames(df_D)<-'Depression'
978. #df_D<-rbind(df_D,tmao_D[-order(rowSums(tmao_D),decreasing = T)[1:20],]%>%rowSums()%>%sum())
979. #rownames(df_D)[21]<-'Others'
980. #top20_H<-rownames(tmao_H)[order(rowSums(tmao_H),decreasing = T)[1:20]]
981. #df_H<-rowSums(tmao_H[top20_H,])%>%as.data.frame()
982. #colnames(df_H)<-'Healthy'
983. #df_H<-rbind(df_H,tmao_H[-order(rowSums(tmao_H),decreasing = T)[1:20],]%>%rowSums()%>%sum())
984. #rownames(df_H)[21]<-'Others'
985. #df<-merge(df_D,df_H,by = "row.names",all = TRUE)
986. #rownames(df)<-df$Row.names
987. #df$Row.names<-NULL
988. #df[is.na(df)]<-0
989. #df$genus<-rownames(df)
990. ##temp<-colSums(df[,1:2])%>%as.numeric()
991. ##df[,1]<-df[,1]/temp[1]
992. ##df[,2]<-df[,2]/temp[2]
993. #data_long <- df %>%
994. # tidyr::gather(key = "group", value = "value", -genus)
995. #data_alluvial <- data_long %>%
996. # dplyr::mutate(flow = rep(1:length(unique(data_long$genus)),2))
997. #desired_order<-c('Healthy','Depression')
998. #col<-rgb(runif(22), runif(22), runif(22))
999. #p20<-ggplot(data_alluvial, aes(x = group, y = value, fill = genus, stratum = genus, alluvium = flow)) +
1000. # geom_stratum(width = 0.4) +
1001. # geom_flow(stat = "alluvium", lode.guidance = "forward", aes.flow = "backward") +
1002. # theme_minimal() +
1003. # labs(x = "", y = "Trimethylamine N-oxide[e]") +
1004. # scale_x_discrete(limits=desired_order)+
1005. # ggtitle('C')
1006. #p20
1007. #
1008. #pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig4.pdf',width = 10,height = 7)
1009. #grid.arrange(p18,arrangeGrob(p19, p20,ncol=2), nrow = 2)
1010. #dev.off()
1011. #
1012. ########
1013. #L_Dopa_D<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/34dhphe[e]_D.csv',row.names = 1)
1014. #L_Dopa_H<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/34dhphe[e]_H.csv',row.names = 1)
1015. #L_Dopa <- merge(L_Dopa_D, L_Dopa_H, by = "row.names", all = TRUE)
1016. #L_Dopa[is.na(L_Dopa)] <- 0
1017. #rownames(L_Dopa) <- L_Dopa$Row.names
1018. #L_Dopa$Row.names <- NULL
1019. #for(i in 1:nrow(L_Dopa)){
1020. # rownames(L_Dopa)[i]<-substr(rownames(L_Dopa)[i],1,nchar(rownames(L_Dopa)[i])-10)
1021. #}
1022. #L_Dopa_D<-L_Dopa[,match(colnames(L_Dopa_D),colnames(L_Dopa))]
1023. #L_Dopa_H<-L_Dopa[,match(colnames(L_Dopa_H),colnames(L_Dopa))]
1024. #res_L_Dopa<-matrix(0,nrow(L_Dopa),2)%>%as.data.frame()
1025. #rownames(res_L_Dopa)<-rownames(L_Dopa)
1026. #colnames(res_L_Dopa)<-c('p.val','fdr')
1027. #for(i in 1:nrow(L_Dopa_D)){
1028. # res_L_Dopa[i,1]<-wilcox.test(L_Dopa_D[i,]%>%as.numeric(),L_Dopa_H[i,]%>%as.numeric(),paired = T)$p.val
1029. #}
1030. #res_L_Dopa[,2]<-p.adjust(res_L_Dopa[,1],method = 'BH')
1031. #diff_L_Dopa<-rownames(res_L_Dopa)[which(res_L_Dopa[,2]<0.1)]
1032. #setdiff(diff_L_Dopa,diff)
1033. ##########
1034. #tmao_D<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/tmao[e]_D.csv',row.names = 1)
1035. #tmao_H<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/tmao[e]_H.csv',row.names = 1)
1036. #tmao <- merge(tmao_D, tmao_H, by = "row.names", all = TRUE)
1037. #tmao[is.na(tmao)] <- 0
1038. #rownames(tmao) <- tmao$Row.names
1039. #tmao$Row.names <- NULL
1040. #for(i in 1:nrow(tmao)){
1041. # rownames(tmao)[i]<-substr(rownames(tmao)[i],1,nchar(rownames(tmao)[i])-7)
1042. #}
1043. #tmao_D<-tmao[,match(colnames(tmao_D),colnames(tmao))]
1044. #tmao_H<-tmao[,match(colnames(tmao_H),colnames(tmao))]
1045. #res_tmao<-matrix(0,nrow(tmao),2)%>%as.data.frame()
1046. #rownames(res_tmao)<-rownames(tmao)
1047. #colnames(res_tmao)<-c('p.val','fdr')
1048. #for(i in 1:nrow(tmao_D)){
1049. # res_tmao[i,1]<-wilcox.test(tmao_D[i,]%>%as.numeric(),tmao_H[i,]%>%as.numeric(),paired = T)$p.val
1050. #}
1051. #res_tmao[,2]<-p.adjust(res_tmao[,1],method = 'BH')
1052. #diff_tmao<-rownames(res_tmao)[which(res_tmao[,1]<0.05)]
1053. #setdiff(diff_tmao,diff)
1057. #######FEA
1058. reaction_all<-read.csv('D:/microbiome/micom/202404/paper/data/ReactionDatabase.csv',row.names = 1)
1059. reaction_used<-reaction_all[match(rownames(res_reaction),rownames(reaction_all)),,drop=F]
1060. diff_reaction_used<-reaction_all[match(diff_reaction,rownames(reaction_all)),,drop=F]
1061. reaction_used<-na.omit(reaction_used)
1062. reaction_used$Subsystem[which(reaction_used$Subsystem=='')]<-'Other'
1063. Subsystem_list <- split(rownames(reaction_used), reaction_used$Subsystem)
1064. all_reactions<-rownames(reaction_used)
1065. selected_reactions<-rownames(diff_reaction_used)
1066. res_FEA <- data.frame(
1067. Subsystem = character(),
1068. PValue = numeric(),
1069. stringsAsFactors = FALSE
1070. )
1071. N <- length(all_reactions)
1072. M <- length(selected_reactions)
1073. for (Subsystem in names(Subsystem_list)) {
1074. K <- length(Subsystem_list[[Subsystem]])
1075. x <- sum(selected_reactions %in% Subsystem_list[[Subsystem]])
1076. p_value <- phyper(x - 1, K, N - K, M, lower.tail = FALSE)
1077. res_FEA <- rbind(res_FEA, data.frame(Subsystem = Subsystem, PValue = p_value))
1078. }
1079. res_FEA <- res_FEA %>%
1080. mutate(AdjustedPValue = p.adjust(PValue, method = "BH")) %>%
1081. arrange(AdjustedPValue)
1082. res_FEA_sig<-res_FEA[which(res_FEA$AdjustedPValue<0.05),]
1083. index<-match(res_FEA_sig$Subsystem,names(Subsystem_list))
1084. ratio<-c()
1085. for(i in 1:length(index)){
1086. ratio<-c(ratio,length(which(diff_reaction_used$Subsystem==names(Subsystem_list)[[index[i]]]))/length(Subsystem_list[[index[i]]]))
1087. }
1088. res_FEA_sig$neg_log_padj <- -log10(res_FEA_sig$AdjustedPValue)
1089. res_FEA_sig$ratio<-ratio
1090. p21<-ggplot(res_FEA_sig, aes(x = ratio, y = Subsystem, size = neg_log_padj, color = AdjustedPValue)) +
1091. geom_point(alpha = 0.7) +
1092. scale_size_continuous(name = "-log10(p.adjust)") +
1093. scale_color_gradient(low = "#4DBBD5FF", high = "#F39B7FFF", name = "fdr") +
1094. labs(title = "",
1095. x = "Reaction Ratio",
1096. y = "Subsystem") +
1097. theme_minimal()+
1098. theme(axis.text.y = element_text(size = 12),
1099. axis.text.x = element_text(size = 12),
1100. plot.title = element_text(hjust = -1.05,face = 'bold'))+
1101. guides(color = guide_legend(ncol = 2),
1102. size = guide_legend(ncol = 2))+
1103. ggtitle('E')
1105. p21
1106. a<-ggplot()+theme_minimal()
1107. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig3_12.pdf',width = 16,height = 14)
1108. grid.arrange(arrangeGrob(p14, p15,nrow=2),arrangeGrob(p16,a,p21,nrow=3,heights = c(0.3,0.4,0.3)), ncol = 2)
1109. dev.off()
1110. pdf(file = 'D:/microbiome/micom/202404/paper/figures/p12.pdf',width = 9,height = 6)
1111. p12
1112. dev.off()
1114. ##########
1115. #metabolite_D<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/metabolite_D.csv',row.names = 1)
1116. #metabolite_H<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/metabolite_H.csv',row.names = 1)
1117. #merged_df <- merge(metabolite_D, metabolite_H, by = "row.names", all = TRUE)
1118. #merged_df[is.na(merged_df)] <- 0
1119. #rownames(merged_df) <- merged_df$Row.names
1120. #merged_df$Row.names <- NULL
1121. #merged_df<-abs(merged_df)
1122. #OTU<-phyloseq::otu_table(merged_df,taxa_are_rows = T)
1123. #meta1<-phyloseq::sample_data(meta)
1124. #physeq<-phyloseq(OTU,meta1)
1125. #mpse<-as.MPSE(physeq)
1126. #mpse@assays@data@listData[["Abundance"]]<-mpse@assays@data@listData[["Abundance"]]*1e10
1127. #mpse %<>%
1128. # mp_decostand(.abundance=Abundance)
1129. #cols<-c('Depression'='#E7B800', 'Healthy'='#00AFBB')
1130. #mpse %<>% mp_cal_alpha(.abundance = Abundance,force = T)
1131. #p22 <- mpse %>%
1132. # mp_plot_alpha(
1133. # .alpha = c(Observe, Shannon,Chao1,Simpson, Pielou),
1134. # .group = disease,
1135. # ) +
1136. # scale_fill_manual(values=cols) +
1137. # scale_color_manual(values=cols) +
1138. # theme(
1139. # legend.position="none",
1140. # strip.background = element_rect(colour=NA, fill="grey")
1141. # )+
1142. # ggtitle('A')
1143. #p22
1144. ###########anosim
1145. #
1146. ###########pcoa
1147. #dist_bray <- phyloseq::distance(physeq, method = 'bray')
1148. #pcoa <- ordinate(physeq, method = "PCoA", distance = dist_bray)
1149. #pcoa_df <- data.frame(pcoa$vectors)
1150. #pcoa_df$SampleID <- rownames(pcoa_df)
1151. #sample_data<-cbind(rownames(meta),meta$disease)%>%as.data.frame()
1152. #colnames(sample_data)<-c('SampleID','disease')
1153. #rownames(sample_data)<-sample_data$SampleID
1154. #pcoa_df <- merge(pcoa_df, sample_data, by.x = "SampleID", by.y = "SampleID")
1155. #pcoa_df$age<-meta$Host.age[match(pcoa_df$SampleID,rownames(meta))]
1156. #pcoa_df$BMI<-meta$BMI[match(pcoa_df$SampleID,rownames(meta))]
1157. #pcoa_df$age[which(pcoa_df$age>=0&pcoa_df$age<=20)]<-'0-20'
1158. #pcoa_df$age[which(pcoa_df$age>20&pcoa_df$age<=40)]<-'20-40'
1159. #pcoa_df$age[which(pcoa_df$age>40&pcoa_df$age<=60)]<-'40-60'
1160. #pcoa_df$age[which(pcoa_df$age>60&pcoa_df$age<=80)]<-'60-80'
1161. #p23<-ggplot(pcoa_df, aes(x = Axis.1, y = Axis.2, color = disease, shape = age, size = BMI)) +
1162. # geom_point(alpha = 0.4) +
1163. # labs(
1164. # x = paste0("PCoA1 (", round(pcoa$values$Relative_eig[1] * 100, 2), "%)"),
1165. # y = paste0("PCoA2 (", round(pcoa$values$Relative_eig[2] * 100, 2), "%)"),
1166. # color = "label",
1167. # shape = "Age",
1168. # size = "BMI"
1169. # ) +
1170. # theme_minimal()+
1171. # scale_color_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
1172. # annotate("text", x = 0.4, y = -0.4, label = "ANOSIM R: 0.007962
1173. # P: 0.004")+
1174. # ggtitle('B')
1175. #p23<-ggMarginal(p23, type = "boxplot", margins = "both", size = 5, groupColour = TRUE, groupFill = TRUE)
1176. #p23
1177. ########
1178. #a<-res[rownames(res_L_Dopa),]
1179. #a<-na.omit(a)
1180. #b<-res_L_Dopa[match(rownames(a),rownames(res_L_Dopa)),]
1181. #aa<-cbind(a[,1],b[1])
1182. #colnames(aa)<-c('p.val1','p.val2')
1183. #write.csv(aa,'D:/microbiome/micom/202404/paper/result/1-1/res_genus-ldopa.csv')
1184. #
1185. #a<-res[rownames(res_tmao),]
1186. #a<-na.omit(a)
1187. #b<-res_tmao[match(rownames(a),rownames(res_tmao)),]
1188. #aa<-cbind(a[,1],b[1])
1189. #colnames(aa)<-c('p.val1','p.val2')
1190. #write.csv(aa,'D:/microbiome/micom/202404/paper/result/1-1/res_genus-tmao.csv')
1191. ########
1193. metabolites<-c('15dap\\[e\\]','34dhphe\\[c\\]','3c3hmp\\[c\\]','5dglcn\\[c\\]','acACP\\[c\\]',
1194. 'C02528\\[c\\]','cinnm\\[c\\]','cu2\\[e\\]',
1195. 'gdpddman\\[c\\]','gdpfuc\\[c\\]','idon_L\\[c\\]')
1197. #for(i in 1:length(metabolites)){
1198. # print(i)
1199. # sampleids_D<-list.files('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression')
1200. # for(j in 1:length(sampleids_D)){
1201. # sampleids_D[j]<-substr(sampleids_D[j],1,nchar(sampleids_D[j])-4)
1202. # }
1203. # temp7<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleids_D[1],'.csv'),row.names = 1)
1204. # temp7<-temp7[grep(metabolites[i],rownames(temp7)),,drop=F]
1205. # temp8<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleids_D[2],'.csv'),row.names = 1)
1206. # temp8<-temp8[grep(metabolites[i],rownames(temp8)),,drop=F]
1207. #
1208. # merged_df <- merge(temp7, temp8, by = "row.names", all = TRUE)
1209. # merged_df[is.na(merged_df)] <- 0
1210. # rownames(merged_df) <- merged_df$Row.names
1211. # merged_df$Row.names <- NULL
1212. # for(k in 3:length(sampleids_D)){
1213. # temp9<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/',sampleids_D[k],'.csv'),row.names = 1)
1214. # temp9<-temp9[grep(metabolites[i],rownames(temp9)),,drop=F]
1215. # merged_df <- merge(merged_df, temp9, by = "row.names", all = TRUE)
1216. # merged_df[is.na(merged_df)] <- 0
1217. # rownames(merged_df) <- merged_df$Row.names
1218. # merged_df$Row.names <- NULL
1219. # }
1220. # write.csv(merged_df,paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'))
1221. #
1222. #
1223. # sampleids_H<-list.files('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/healthy')
1224. # for(j1 in 1:length(sampleids_H)){
1225. # sampleids_H[j1]<-substr(sampleids_H[j1],1,nchar(sampleids_H[j1])-4)
1226. # }
1227. # temp7<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/healthy/',sampleids_H[1],'.csv'),row.names = 1)
1228. # temp7<-temp7[grep(metabolites[i],rownames(temp7)),,drop=F]
1229. # temp8<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/healthy/',sampleids_H[2],'.csv'),row.names = 1)
1230. # temp8<-temp8[grep(metabolites[i],rownames(temp8)),,drop=F]
1231. #
1232. # merged_df <- merge(temp7, temp8, by = "row.names", all = TRUE)
1233. # merged_df[is.na(merged_df)] <- 0
1234. # rownames(merged_df) <- merged_df$Row.names
1235. # merged_df$Row.names <- NULL
1236. # for(k1 in 3:length(sampleids_H)){
1237. # temp9<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/healthy/',sampleids_H[k1],'.csv'),row.names = 1)
1238. # temp9<-temp9[grep(metabolites[i],rownames(temp9)),,drop=F]
1239. # merged_df <- merge(merged_df, temp9, by = "row.names", all = TRUE)
1240. # merged_df[is.na(merged_df)] <- 0
1241. # rownames(merged_df) <- merged_df$Row.names
1242. # merged_df$Row.names <- NULL
1243. # }
1244. # write.csv(merged_df,paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'))
1245. #}
1247. #########
1248. i<-1
1249. temp_D<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'),row.names = 1)
1250. temp_H<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'),row.names = 1)
1251. meta_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_H_used.csv',row.names = 1)
1252. meta_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_D_used.csv',row.names = 1)
1253. colnames(meta_H)[5]<-'disease'
1254. colnames(meta_D)[5]<-'disease'
1255. meta_H$disease<-'Healthy'
1256. meta_D$disease<-'Depression'
1257. meta<-rbind(meta_D,meta_H)
1258. for(j in 1:nrow(temp_D)){
1259. rownames(temp_D)[j]<-substr(rownames(temp_D)[j],1,nchar(rownames(temp_D)[j])-(nchar(metabolites[i])-2))
1260. }
1261. for(j in 1:nrow(temp_H)){
1262. rownames(temp_H)[j]<-substr(rownames(temp_H)[j],1,nchar(rownames(temp_H)[j])-(nchar(metabolites[i])-2))
1263. }
1265. temp <- merge(temp_D, temp_H, by = "row.names", all = TRUE)
1266. temp[is.na(temp)] <- 0
1267. rownames(temp) <- temp$Row.names
1268. temp$Row.names <- NULL
1269. temp<-t(temp)%>%as.data.frame()
1270. temp<-temp[,-which(colSums(temp)==0)]
1271. temp$sampleid<-rownames(temp)
1273. df1<-pivot_longer(temp,
1274. cols = -sampleid,
1275. names_to = "genus",
1276. values_to = "flux"
1277. )
1278. df1$metabolite<-substr(metabolites[i],1,nchar(metabolites[i])-5)
1279. df1$label<-meta$disease[match(df1$sampleid,rownames(meta))]
1280. for(i in 2:length(metabolites)){
1281. print(i)
1282. temp_D<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'),row.names = 1)
1283. temp_H<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'),row.names = 1)
1284. for(j in 1:nrow(temp_D)){
1285. rownames(temp_D)[j]<-substr(rownames(temp_D)[j],1,nchar(rownames(temp_D)[j])-(nchar(metabolites[i])-2))
1286. }
1287. for(j in 1:nrow(temp_H)){
1288. rownames(temp_H)[j]<-substr(rownames(temp_H)[j],1,nchar(rownames(temp_H)[j])-(nchar(metabolites[i])-2))
1289. }
1291. temp <- merge(temp_D, temp_H, by = "row.names", all = TRUE)
1292. temp[is.na(temp)] <- 0
1293. rownames(temp) <- temp$Row.names
1294. temp$Row.names <- NULL
1295. temp<-t(temp)%>%as.data.frame()
1296. temp<-temp[,-which(colSums(temp)==0)]
1297. temp$sampleid<-rownames(temp)
1299. df3<-pivot_longer(temp,
1300. cols = -sampleid,
1301. names_to = "genus",
1302. values_to = "flux"
1303. )
1304. df3$metabolite<-substr(metabolites[i],1,nchar(metabolites[i])-5)
1305. df3$label<-meta$disease[match(df3$sampleid,rownames(meta))]
1306. df1<-rbind(df1,df3)
1307. }
1308. #pdf(file = 'test.pdf',width = 10,height = 1200)
1309. #ggplot(df1, aes(x = label, y = flux, fill = label)) +
1310. # geom_boxplot(outlier.shape = NA) +
1311. # theme(axis.text.x = element_blank(), # ihxh=4ff,f g->
1312. # axis.ticks.x = element_blank())+
1313. # stat_compare_means(comparisons = list(c("Depression", "Healthy")),label.y = 1.5)+
1314. # scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
1315. # coord_cartesian(ylim = c(-1, 1))+
1316. # facet_grid(genus ~ metabolite)
1317. #dev.off()
1320. #write.csv(df1,'data.csv')
1321. #########
1322. #diff<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/diff_taxon.csv',header = F)%>%as.character()
1323. library(tidyverse)
1324. library(RColorBrewer)
1325. library(effsize)
1326. #df1$flux<-scale(df1$flux,scale = T,center = T)
1327. metabolite_used<-unique(df1$metabolite)
1328. metabolite_all<-fread('D:/microbiome/micom/202404/paper/data/MetaboliteDatabase.txt')%>%as.data.frame()
1329. metabolite_used1<-metabolite_all[match(metabolite_used,metabolite_all[,1]),2]
1330. df1$label<-meta$disease[match(df1$sampleid,rownames(meta))]
1331. meta<-df1[!duplicated(df1$sampleid),c(1,5)]
1332. res_all<-data.frame()
1334. for(i in 1:length(metabolite_used)){
1335. print(i)
1336. temp_m<-df1[df1$metabolite==metabolite_used[i],]
1337. con<-meta$sampleid[meta$label=="Healthy"]
1338. ds<-meta$sampleid[meta$label=="Depression"]
1339. m1<-df1[df1$metabolite==metabolite_used[i],]
1340. mat1<-pivot_wider(m1[,1:3],names_from = sampleid,values_from = flux)
1341. mat1<-column_to_rownames(mat1,var = "genus")
1342. mat1<-mat1[,meta$sampleid]
1343. fc1<-rowMeans(mat1[,ds])/rowMeans(mat1[,con])
1344. p1<-c()
1345. cliff_delta<-c()
1346. a<-c()
1347. for (j in 1:nrow(mat1)) {
1348. ptmp<-wilcox.test(as.numeric(mat1[j,ds]),as.numeric(mat1[j,con]),paired = T)$p.val
1349. cliff_delta_temp<-cliff.delta(as.numeric(mat1[j,ds]),as.numeric(mat1[j,con]))$estimate
1350. p1<-c(p1,ptmp)
1351. cliff_delta<-c(cliff_delta,cliff_delta_temp)
1352. }
1353. res1<-data.frame(genus=rownames(mat1),fc=fc1,p=p1,metabolite=metabolite_used1[i],cliff_delta=cliff_delta,a=rowMeans(mat1[,ds]))
1354. res_all<-rbind(res_all,res1)
1355. }
1356. res_all$sig<-'no'
1357. res_all$sig[which(res_all$p<0.05)]<-'yes'
1358. res_all1<-res_all[res_all$sig=='yes',]
1359. #setdiff(res_all1$genus%>%unique(),diff_top50)
1360. cu_genus<-res_all1$genus[which(res_all1$metabolite=='Cu2+')]
1362. res_all2<-res_all1[which(res_all1$metabolite!='Cu2+'),]
1363. res_cu<-res_all1[which(res_all1$metabolite=='Cu2+'),]
1364. res_all2$a[res_all2$a>0]<-'export'
1365. res_all2$a[res_all2$a<0]<-'import'
1366. colnames(res_all2)[6]<-'direction'
1367. res_all2 <- res_all2 %>%
1368. mutate(significance = case_when(
1369. p < 0.001 ~ "***",
1370. p < 0.01 ~ "**",
1371. p < 0.05 ~ "*",
1372. TRUE ~ ""
1373. ))
1375. p22<-ggplot2::ggplot()+
1376. geom_point(data = res_all2,
1377. aes(x = genus, y = metabolite,size = fc,color = significance,shape=direction))+
1378. theme_minimal()+
1379. #scale_fill_manual(c('*'='red','**'='green'))+
1380. theme(axis.text.x = element_text(angle = 90,size = 12,face = "italic"),
1381. axis.title.x = element_text(size = 12),
1382. axis.text.y = element_text(size = 12),
1383. legend.position = 'left',
1384. plot.title = element_text(hjust = -2,face = 'bold'))+
1385. scale_y_discrete(position = 'right')+
1386. labs(x='',y='')+
1387. guides(color = guide_legend(ncol = 1),
1388. size = guide_legend(ncol = 2)) +
1389. scale_x_discrete(labels = function(y) {
1390. ifelse(y %in% setdiff(res_all2$genus%>%unique(),diff),
1391. paste0("<span style='color:red;'>", y, "</span>"),
1392. y)
1393. })+
1394. scale_shape_manual(values = c("export" = 24, "import" = 25))+
1395. theme(axis.text.x=element_markdown())+
1396. ggtitle('B')
1397. p22
1398. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig4.pdf',width = 15,height = 12)
1399. grid.arrange(p18,p22, ncol = 2)
1400. dev.off()
1401. #ggsave("dotp.pdf",height = 15,width = 6,limitsize = FALSE)
1402. ######
1404. i<-8
1405. cu_D<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'),row.names = 1)
1406. cu_H<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'),row.names = 1)
1407. taxon<-read.csv('D:/microbiome/micom/202404/paper/data/taxa.csv')
1408. for(j in 1:nrow(cu_D)){
1409. rownames(cu_D)[j]<-substr(rownames(cu_D)[j],1,nchar(rownames(cu_D)[j])-(nchar(metabolites[i])-2))
1410. }
1411. for(j in 1:nrow(cu_H)){
1412. rownames(cu_H)[j]<-substr(rownames(cu_H)[j],1,nchar(rownames(cu_H)[j])-(nchar(metabolites[i])-2))
1413. }
1415. cu <- merge(cu_D, cu_H, by = "row.names", all = TRUE)
1416. cu[is.na(cu)] <- 0
1417. rownames(cu) <- cu$Row.names
1418. cu$Row.names <- NULL
1419. cu_data<-cu
1420. cu<-t(cu)%>%as.data.frame()
1421. cu<-cu[,-which(colSums(cu)==0)]
1422. cu<-abs(cu)
1423. for(i in 1:nrow(cu)){
1424. cu[i,]<-cu[i,]/rowSums(cu)[i]
1425. }
1426. cu$sampleid<-rownames(cu)
1428. df<-pivot_longer(cu,
1429. cols = -sampleid,
1430. names_to = "genus",
1431. values_to = "flux"
1432. )
1433. df$phyla<-taxon$p__[match(df$genus,taxon$g__)]
1434. df$flux<-abs(df$flux)
1435. df$flux<-log2(df$flux+1)
1436. df$metabolic<-'Cu2+'
1437. df$group<-meta$disease[match(df$sampleid,rownames(meta))]
1438. set.seed(1)
1439. col<-rgb(runif(length(unique(df$phyla))),runif(length(unique(df$phyla))),runif(length(unique(df$phyla))))
1440. names(col)<-unique(df$phyla)
1441. ggplot(df, aes(x = sampleid, y = flux, fill = phyla)) +
1442. geom_bar(stat = "identity") +
1443. scale_fill_manual(values = col)+
1444. theme(axis.text.x = element_blank(),
1445. axis.ticks.x = element_blank(),
1446. panel.grid = element_blank(),
1447. strip.text = element_text(size = 12))+
1448. facet_grid(group~.,scales = "free_x",as.table = F)
1453. #######
1454. p_phyla<-list()
1455. p_genus<-list()
1456. meta_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_H_used.csv',row.names = 1)
1457. meta_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_D_used.csv',row.names = 1)
1458. colnames(meta_H)[5]<-'disease'
1459. colnames(meta_D)[5]<-'disease'
1460. meta_H$disease<-'Healthy'
1461. meta_D$disease<-'Depression'
1462. meta<-rbind(meta_D,meta_H)
1463. i=8
1464. for(i in 1:length(metabolites)){
1465. cu_D<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'),row.names = 1)
1466. cu_H<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'),row.names = 1)
1467. cu_D<-cu_D[,-which(colnames(cu_H)=='ERR1090035')]
1468. cu_H<-cu_H[,-which(colnames(cu_H)=='ERR1090035')]
1469. taxon<-read.csv('D:/microbiome/micom/202404/paper/data/taxa.csv')
1470. for(j in 1:nrow(cu_D)){
1471. rownames(cu_D)[j]<-substr(rownames(cu_D)[j],1,nchar(rownames(cu_D)[j])-(nchar(metabolites[i])-2))
1472. }
1473. for(j in 1:nrow(cu_H)){
1474. rownames(cu_H)[j]<-substr(rownames(cu_H)[j],1,nchar(rownames(cu_H)[j])-(nchar(metabolites[i])-2))
1475. }
1477. cu <- merge(cu_D, cu_H, by = "row.names", all = TRUE)
1478. cu[is.na(cu)] <- 0
1479. rownames(cu) <- cu$Row.names
1480. cu$Row.names <- NULL
1481. cu<-t(cu)%>%as.data.frame()
1482. cu<-cu[,-which(colSums(cu)==0)]
1483. cu$sampleid<-rownames(cu)
1484. df<-pivot_longer(cu,
1485. cols = -sampleid,
1486. names_to = "genus",
1487. values_to = "flux"
1488. )
1489. df$phyla<-taxon$p__[match(df$genus,taxon$g__)]
1490. df$flux<-abs(df$flux)
1491. df$flux<-log2(df$flux+1)
1492. df$metabolic<-metabolites[i]
1493. set.seed(1)
1494. col<-rgb(runif(length(unique(df$phyla))),runif(length(unique(df$phyla))),runif(length(unique(df$phyla))))
1495. names(col)<-unique(df$phyla)
1496. col_genus<-rgb(runif(length(unique(df$genus))),runif(length(unique(df$genus))),runif(length(unique(df$genus))))
1497. names(col_genus)<-unique(df$genus)
1498. p_genus[[i]]<-ggplot(df, aes(x = sampleid, y = flux, fill = genus)) +
1499. geom_bar(stat = "identity") +
1500. scale_fill_manual(values = col_genus)+
1501. theme(axis.text.x = element_blank(),
1502. axis.ticks.x = element_blank(),
1503. panel.grid = element_blank(),
1504. strip.text = element_text(size = 12),
1505. legend.position = 'none')+
1506. facet_wrap(~ metabolic)
1507. }
1508. for(i in 1:length(metabolites)){
1509. cu_D<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'),row.names = 1)
1510. cu_H<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'),row.names = 1)
1511. cu_D<-cu_D[,-which(colnames(cu_H)=='ERR1090035')]
1512. cu_H<-cu_H[,-which(colnames(cu_H)=='ERR1090035')]
1513. taxon<-read.csv('D:/microbiome/micom/202404/paper/data/taxa.csv')
1514. for(j in 1:nrow(cu_D)){
1515. rownames(cu_D)[j]<-substr(rownames(cu_D)[j],1,nchar(rownames(cu_D)[j])-(nchar(metabolites[i])-2))
1516. }
1517. for(j in 1:nrow(cu_H)){
1518. rownames(cu_H)[j]<-substr(rownames(cu_H)[j],1,nchar(rownames(cu_H)[j])-(nchar(metabolites[i])-2))
1519. }
1521. cu <- merge(cu_D, cu_H, by = "row.names", all = TRUE)
1522. cu[is.na(cu)] <- 0
1523. rownames(cu) <- cu$Row.names
1524. cu$Row.names <- NULL
1525. cu<-t(cu)%>%as.data.frame()
1526. cu<-cu[,-which(colSums(cu)==0)]
1527. cu$sampleid<-rownames(cu)
1528. df<-pivot_longer(cu,
1529. cols = -sampleid,
1530. names_to = "genus",
1531. values_to = "flux"
1532. )
1533. df$phyla<-taxon$p__[match(df$genus,taxon$g__)]
1534. df$flux<-abs(df$flux)
1535. df$flux<-log2(df$flux+1)
1536. df$metabolic<-metabolites[i]
1537. Firm<-df[which(df$phyla=='Proteobacteria'),]
1538. Firm_summarized <- Firm %>%
1539. group_by(sampleid) %>%
1540. summarize(
1541. flux = sum(flux),
1542. )
1543. set.seed(1)
1544. col_phyla<-rgb(runif(length(unique(df$phyla))),runif(length(unique(df$phyla))),runif(length(unique(df$phyla))))
1545. names(col_phyla)<-unique(df$phyla)
1546. col_phyla<-c(col_phyla,'#E7B800','#00AFBB')
1547. names(col_phyla)[19:20]<-c('Depression','Healthy')
1548. col_genus<-rgb(runif(length(unique(df$genus))),runif(length(unique(df$genus))),runif(length(unique(df$genus))))
1549. names(col_genus)<-unique(df$genus)
1550. Firm_summarized_D<-Firm_summarized[match(rownames(meta)[which(meta$disease=='Depression')],Firm_summarized$sampleid),]
1551. Firm_summarized_H<-Firm_summarized[match(rownames(meta)[which(meta$disease=='Healthy')],Firm_summarized$sampleid),]
1552. levels<-c(Firm_summarized_D$sampleid[order(Firm_summarized_D$flux)],rownames(meta_H)[match(Firm_summarized_D$sampleid[order(Firm_summarized_D$flux)],rownames(meta_D))])
1553. ##df$sampleid<-factor(df$sampleid,levels = Firm_summarized$sampleid[order(Firm_summarized$flux)])
1554. df$sampleid<-factor(df$sampleid,levels = levels)
1555. df$label<-meta$disease[match(df$sampleid,rownames(meta))]
1556. df<-na.omit(df)
1557. #tile_data <- df %>%
1558. # mutate(xmin = as.numeric(factor(sampleid)) - 0.5,
1559. # xmax = as.numeric(factor(sampleid)) + 0.5,
1560. # ymin = -1,
1561. # ymax = -0.5)
1562. p_phyla[[i]]<-ggplot(df, aes(x = sampleid, y = flux, fill = phyla)) +
1563. geom_bar(stat = "identity") +
1564. theme(axis.text.x = element_blank(),
1565. axis.ticks.x = element_blank(),
1566. panel.grid = element_blank(),
1567. strip.text = element_text(size = 12)) +
1568. scale_fill_manual(values = col_phyla)+
1569. facet_wrap(~ label, scales = "free_x", ncol = 1)
1571. }
1572. p25<-ggplot(df, aes(x = sampleid, y = flux, fill = phyla)) +
1573. geom_bar(stat = "identity") +
1574. theme(axis.text.x = element_blank(),
1575. axis.ticks.x = element_blank(),
1576. panel.grid = element_blank(),
1577. strip.text = element_text(size = 12),
1578. plot.title = element_text(face = 'bold')) +
1579. scale_fill_manual(values = col_phyla)+
1580. labs(y='flux\n mmol/[gDW h](scaled)',x='Sample id')+
1581. facet_wrap(~ label, scales = "free_x", ncol = 1)+
1582. ggtitle('D')
1583. p25
1584. ########
1586. cu_data$genus<-rownames(cu_data)
1587. cu_data<-pivot_longer(cu_data,
1588. cols = -genus,
1589. names_to = "sampleid",
1590. values_to = "flux"
1591. )
1592. cu_data$label<-meta$label[match(cu_data$sampleid,meta$sampleid)]
1593. p<-list()
1594. for(i in 1:length(cu_genus)){
1595. temp<-cu_data[which(cu_data$genus==cu_genus[i]),]
1596. p[[i]]<-ggplot(temp, aes(x = label, y = flux, fill = label)) +
1597. geom_violin() +
1598. labs( x = "", y = '') +
1599. facet_wrap(~ genus)+
1600. stat_compare_means(aes(group = label),label = "p.signif",label.x = 1.5) +
1601. theme(axis.text.x = element_blank(),
1602. strip.text = element_text(size=14),
1603. panel.grid.major = element_blank(),
1604. panel.grid.minor = element_blank(),
1605. axis.line.x = element_blank(),
1606. axis.ticks.x = element_blank(),
1607. axis.text.y = element_blank(),
1608. axis.ticks.y = element_blank(),
1609. plot.margin = unit(c(0.01, 0.01, 0.01, 0.01), "cm"),#top, right, bottom, and left
1610. legend.position = 'none')+
1611. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))
1612. }
1613. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig6.pdf',width = 36,height = 20)
1614. grid.arrange(arrangeGrob(p[[1]],p[[2]],p[[3]],p[[4]],p[[5]],p[[6]],p[[7]],p[[8]],p[[9]],ncol=9),
1615. arrangeGrob(p[[10]],p[[11]],p[[12]],p[[13]],p[[14]],p[[15]],p[[16]],p[[17]],p[[18]],ncol=9),
1616. arrangeGrob(p[[19]],p[[20]],p[[21]],p[[22]],p[[23]],p[[24]],p[[25]],p[[26]],p[[27]],ncol=9),
1617. arrangeGrob(p[[28]],p[[29]],p[[30]],p[[31]],p[[32]],p[[33]],p[[34]],p[[35]],p[[36]],ncol=9),
1618. arrangeGrob(p[[37]],p[[38]],p[[39]],p[[40]],p[[41]],p[[42]],p[[43]],p[[44]],p[[45]],ncol=9),
1619. nrow = 5)
1620. dev.off()
1621. #########
1622. Firm<-df[which(df$phyla=='Firmicutes'),]
1623. Firm_summarized <- Firm %>%
1624. group_by(sampleid) %>%
1625. summarize(
1626. flux = sum(flux),
1627. )
1628. ########
1629. i=8
1630. cu_D<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_D.csv'),row.names = 1)
1631. cu_H<-read.csv(paste0('D:/microbiome/micom/202404/paper/result/1-1/',substr(metabolites[i],1,nchar(metabolites[i])-5),'_H.csv'),row.names = 1)
1632. cu <- merge(cu_D, cu_H, by = "row.names", all = TRUE)
1633. cu[is.na(cu)] <- 0
1634. rownames(cu) <- cu$Row.names
1635. cu$Row.names <- NULL
1636. cu<-t(cu)%>%as.data.frame()
1637. cu<-cu[,-which(colSums(cu)==0)]
1638. cu<-abs(cu)
1639. for(i in 1:ncol(cu)){
1640. colnames(cu)[i]<-substr(colnames(cu)[i],1,nchar(colnames(cu)[i])-6)
1641. }
1642. cu<-cu[,match(cu_genus,colnames(cu))]
1643. cu$group<-meta$disease[match(rownames(cu),rownames(meta))]
1644. cu<-cbind(cu[,ncol(cu)],cu[,-ncol(cu)])
1645. colnames(cu)[1]<-'label'
1646. cu_D<-cu[which(cu$label=='Depression'),]
1647. cu_H<-cu[which(cu$label=='Healthy'),]
1648. temp1<-colMeans(cu_D[,-1])
1649. temp2<-colMeans(cu_H[,-1])
1650. temp<-rbind(temp1,temp2)%>%as.data.frame()
1651. temp$label<-c('Depression','Healthy')
1652. temp<-cbind(temp[,ncol(temp)],temp[,-ncol(temp)])
1654. p23<-ggradar(temp,
1655. grid.line.width = 0,
1656. axis.label.size= 0,
1657. group.line.width = 0.2,
1658. group.point.size = 0.2,
1659. plot.extent.x.sf = 1.6,
1660. background.circle.colour = 'white',
1661. grid.label.size = 0,
1662. grid.min = 0,
1663. grid.mid = 0.2,
1664. grid.max = 0.3,
1665. group.colours = c('#E7B800','#00AFBB'),
1666. background.circle.transparency = 0,
1667. legend.position = 'none')
1668. p23
1669. p24<-ggplot()+theme_minimal()
1670. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig4_12.pdf',width = 15,height = 18)
1671. grid.arrange(arrangeGrob(p18,arrangeGrob(p22,p24,nrow = 2),ncol=2),p25,nrow=2)
1672. dev.off()
1674. pdf(file = 'D:/microbiome/micom/202404/paper/figures/cu_radar.pdf',width = 7,height = 4)
1675. p23
1676. dev.off()
1677. ###########
1678. library(mediation)
1679. data<-rbind(abundance,metabolite)
1680. data<-t(data)%>%as.data.frame()
1681. for(i in 1:ncol(data)){
1682. colnames(data)[i]<-paste0('f_',colnames(data)[i])
1683. }
1684. data$group<-meta$disease[match(rownames(data),rownames(meta))]
1685. data$gender<-meta$Sex[match(rownames(data),rownames(meta))]
1686. data$age<-meta$Host.age[match(rownames(data),rownames(meta))]
1687. data$BMI<-meta$BMI[match(rownames(data),rownames(meta))]
1688. genusnames<-diff
1689. metabolitenames<-diff_metabolite
1690. data$group<-(as.factor(data$group)%>%as.numeric()-1)
1691. colnames(data)<-gsub("\\[|\\]","_",colnames(data))
1692. f1<-paste0('group','~','f_',genusnames[i])
1693. f1<-as.formula(f1)
1694. Y1<-glm(f1,data=data,family = binomial("probit"))
1695. f2<-paste0('group','~','f_',metabolitenames[j],'+gender+BMI+age')
1696. f2<-as.formula(f2)
1697. Y2<-lm(f2,data=data)
1698. ############
1699. library(mediation)
1700. data_a<-abundance[diff,]
1701. data_m<-metabolite[diff_metabolite,]
1702. data<-rbind(data_a,data_m)
1703. data<-t(data)%>%as.data.frame()
1704. allfeatures<-colnames(data)
1705. for(i in 1:ncol(data)){
1706. colnames(data)[i]<-paste0('features',i)
1707. }
1708. data$group<-meta$disease[match(rownames(data),rownames(meta))]
1709. data$gender<-meta$Sex[match(rownames(data),rownames(meta))]
1710. data$age<-meta$Host.age[match(rownames(data),rownames(meta))]
1711. data$BMI<-meta$BMI[match(rownames(data),rownames(meta))]
1712. data$group<-(as.factor(data$group)%>%as.numeric()-1)
1713. result_mediate<-c()
1714. for(i in 1:nrow(data_a)){
1715. print(i)
1716. for(j in (nrow(data_a)+1):(ncol(data)-4)){
1717. print(j)
1718. f1<-paste0('features',j,'~','features',i,'+gender+BMI+age')
1719. f1<-as.formula(f1)
1720. Y1<-lm(f1,data=data)
1721. f2<-paste0('group','~','features',i,'+','features',j,'+gender+BMI+age')
1722. f2<-as.formula(f2)
1723. Y2<-glm(f2,data=data,family = binomial("probit"))
1724. r_causal2 <- mediate(Y1, Y2, treat=paste0('features',i), mediator=paste0('features',j), covariates=c('gender','BMI','age'), boot=TRUE, sims=1000)
1725. res_temp<-c(paste0('features',i),paste0('features',j),
1726. r_causal2$d1.p,r_causal2$d0.p,
1727. r_causal2$d1,r_causal2$d0)
1728. result_mediate <- append(result_mediate, list(res_temp))
1729. }
1730. }
1731. features1<-c()
1732. features2<-c()
1733. res_m<-c()
1734. for(i in 1:length(result_mediate)){
1735. if(result_mediate[[i]][3]%>%as.numeric()<0.001&
1736. result_mediate[[i]][4]%>%as.numeric()<0.001){
1737. features1<-append(features1,allfeatures[match(result_mediate[[i]][1],colnames(data))])
1738. features2<-append(features2,allfeatures[match(result_mediate[[i]][2],colnames(data))])
1739. temp<-c(allfeatures[match(result_mediate[[i]][1],colnames(data))],
1740. allfeatures[match(result_mediate[[i]][2],colnames(data))],
1741. result_mediate[[i]][3]%>%as.numeric())%>%t()%>%as.data.frame()
1742. res_m<-rbind(res_m,temp)
1743. }
1744. }
1746. features1<-c()
1747. features2<-c()
1749. for(i in 1:length(result_mediate)){
1750. if(result_mediate[[i]][3]%>%as.numeric()<0.001&
1751. result_mediate[[i]][4]%>%as.numeric()<0.001){
1752. features1<-append(features1,result_mediate[[i]][1])
1753. features2<-append(features2,result_mediate[[i]][2])
1754. }
1755. }
1757. colnames(res_m)<-c('genus','metabolites','p.val')
1758. res_m[,3]<-as.numeric(res_m[,3])
1759. res_m1<-res_m
1760. #########
1761. res_m1<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/res_m2.csv',row.names = 1)
1762. features1<-res_m1$genus
1763. features2<-res_m1$metabolites
1764. features1<-match(features1,allfeatures)
1765. features1<-paste0('features',features1)
1766. features2<-match(features2,allfeatures)
1767. features2<-paste0('features',features2)
1769. result_mediate1<-c()
1771. for(i in 1:nrow(res_m1)){
1772. print(i)
1773. f0<-paste0('group','~',features1[i],'+gender+BMI+age')#d8-d;ei~h*ei+f77fe g4
1774. f0<-as.formula(f0)
1775. Y0<-glm(f0,data=data,family = binomial("probit"))
1776. f1<-paste0(features2[i],'~',features1[i],'+gender+BMI+age')#d8-d;ei~h*ei+f77fe g4
1777. f1<-as.formula(f1)
1778. Y1<-lm(f1,data=data)
1779. print(summary(Y1))
1780. f2<-paste0('group','~',features1[i],'+',features2[i],'+gender+BMI+age')#e ei~h*ei+d8-d;ei+f77fe g4
1781. f2<-as.formula(f2)
1782. Y2<-glm(f2,data=data,family = binomial("probit"))
1783. r_causal2 <- mediate(Y1, Y2, treat=features1[i], mediator=features2[i], covariates=c('gender','BMI','age'), boot=TRUE, sims=1000)
1784. model_summary<-summary(r_causal2)
1785. res_temp<-c(features1[i],
1786. features2[i],
1787. r_causal2$d.avg,
1788. r_causal2$d.avg.p,
1789. r_causal2$z0.p,
1790. r_causal2$z1.p,
1791. r_causal2$z.avg,
1792. r_causal2$z.avg.p,
1793. Y1$coefficients[2]%>%as.numeric(),
1794. Y2$coefficients[2]%>%as.numeric(),
1795. Y2$coefficients[3]%>%as.numeric(),
1796. r_causal2$n.avg,
1797. Y0$coefficients[2]
1798. )
1799. result_mediate1 <- append(result_mediate1, list(res_temp))
1800. }
1801. features1<-c()
1802. features2<-c()
1803. d.avg<-c()
1804. z.avg<-c()
1805. coefficientsy1<-c()
1806. coefficientsy2<-c()
1807. coefficientsy3<-c()
1808. d.avg.p<-c()
1809. z.avg.p<-c()
1810. n.avg<-c()
1812. res_m2<-c()
1813. for(i in 1:length(result_mediate1)){
1814. if(result_mediate1[[i]][8]%>%as.numeric()<0.01&result_mediate1[[i]][4]%>%as.numeric()<0.01){
1815. print(i)
1816. features1<-result_mediate1[[i]][1]
1817. features2<-result_mediate1[[i]][2]
1818. d.avg<-result_mediate1[[i]][3]
1819. z.avg<-result_mediate1[[i]][7]
1820. d.avg.p<-result_mediate1[[i]][4]
1821. z.avg.p<-result_mediate1[[i]][8]
1822. coefficientsy1<-result_mediate1[[i]][9]
1823. coefficientsy2<-result_mediate1[[i]][10]
1824. coefficientsy3<-result_mediate1[[i]][13]
1825. n.avg<-result_mediate1[[i]][12]
1826. temp<-c(allfeatures[match(features1,colnames(data))],
1827. allfeatures[match(features2,colnames(data))],
1828. coefficientsy1,
1829. coefficientsy2,
1830. coefficientsy3,
1831. d.avg.p,
1832. z.avg.p,
1833. d.avg,
1834. z.avg,
1835. n.avg)%>%t()%>%as.data.frame()
1836. colnames(temp)<-c('genus','metabolites','coefficients1','coefficients2','coefficients3','d.avg.p','z.avg.p','d.avg','z.avg','n.avg')
1837. res_m2<-rbind(res_m2,temp)
1838. }
1839. }
1840. res_m2<-res_m2[res_m2$n.avg%>%as.numeric()>0.1,]
1841. result_mediate2<-c()
1842. for(i in 1:nrow(res_m1)){
1843. print(i)
1844. f1<-paste0('group','~',features1[i],'+gender+BMI+age')
1845. f1<-as.formula(f1)
1846. Y1<-glm(f1,data=data,family = binomial("probit"))
1847. f2<-paste0(features2[i],'~',features1[i],'+','group','+gender+BMI+age')
1848. f2<-as.formula(f2)
1849. Y2<-lm(f2,data=data)
1850. r_causal2 <- mediate(Y1, Y2, treat=features1[i], mediator=features2[i], covariates=c('gender','BMI','age'), boot=TRUE, sims=1000)
1851. res_temp<-c(features1[i],
1852. features2[i],
1853. r_causal2$d.avg,
1854. r_causal2$d.avg.p,
1855. r_causal2$z0.p,
1856. r_causal2$z1.p,
1857. r_causal2$z.avg,
1858. r_causal2$z.avg.p,
1859. #r_causal2$n.avg.p
1860. Y1$coefficients[2]%>%as.numeric(),
1861. Y2$coefficients[3]%>%as.numeric())
1862. result_mediate2 <- append(result_mediate2, list(res_temp))
1863. }
1865. res_m3<-c()
1866. for(i in 1:length(result_mediate2)){
1867. print(i)
1868. features1<-result_mediate2[[i]][1]
1869. features2<-result_mediate2[[i]][2]
1870. #d.avg<-result_mediate2[[i]][3]
1871. #z.avg<-result_mediate2[[i]][7]
1872. d.avg.p<-result_mediate2[[i]][4]
1873. #z.avg.p<-result_mediate2[[i]][8]
1874. #coefficientsy1<-result_mediate2[[i]][9]
1875. #coefficientsy2<-result_mediate2[[i]][10]
1876. #coefficientsy3<-result_mediate2[[i]][13]
1877. #n.avg<-result_mediate2[[i]][12]
1878. temp<-c(allfeatures[match(features1,colnames(data))],
1879. allfeatures[match(features2,colnames(data))],
1880. #coefficientsy1,
1881. #coefficientsy2,
1882. #coefficientsy3,
1883. d.avg.p
1884. #z.avg.p,
1885. #d.avg,
1886. #z.avg,
1887. #n.avg)
1888. )%>%t()%>%as.data.frame()
1889. colnames(temp)<-c('genus','metabolites','d.avg.p')
1890. res_m3<-rbind(res_m3,temp)
1891. }
1897. write.csv(res_m2,'D:/microbiome/micom/202404/paper/result/1-1/res_m2.csv')
1899. #########
1900. library(igraph)
1901. library(ggraph)
1902. metabolite<-unique(res_m1$metabolites)
1903. nodetype<-data.frame(nodes=c(unique(res_m1$genus),
1904. unique(res_m1$metabolites),
1905. "Depression"),
1906. type=c(rep(c("Independent","Mediator","Dependent"),
1907. c(length(unique(res_m1$genus)),length(unique(res_m1$metabolites)),1))))
1908. res_m1<-rbind(res_m1,data.frame(genus=metabolite,metabolites="Depression",p.val=0))
1910. graph<-graph_from_data_frame(res_m1,vertices = nodetype)
1911. pos<-nodetype
1912. pos$x<-c(seq(0.1,0.1*length(unique(res_m1$genus)),0.1),seq(0.2,0.2*length(unique(res_m1$metabolites)),0.2),1.2)
1913. pos$y<-c(rep(0.8,23),rep(0.5,11),0.2)
1915. p26<-ggraph(graph,
1916. layout = "manual",
1917. x = pos$x,
1918. y = pos$y) +
1919. geom_edge_link0(width=1,arrow = arrow(length = unit(3, 'mm'), type = "closed"),alpha=0.5) +
1920. geom_node_point(aes(fill = type), shape = 21, size = 6,alpha=0.5) +
1921. #scale_color_manual(values = c("Independent"="orangered4","Mediator"="orchid4","Dependent"="red")) +
1922. scale_fill_manual(values = c("Independent"="orangered4","Mediator"="orchid4","Dependent"="red")) +
1923. #scale_edge_colour_gradientn(
1924. # colors = c("aquamarine3","grey90"),
1925. # #limits = c(0, 0.28),
1926. # space = "Lab",
1927. # na.value = "grey50")+
1928. guides(fill = guide_legend(ncol = 1))+
1929. theme_graph(base_family = "sans")+#theme(legend.box = 'horizontal',
1930. # legend.box.just = 'top')+
1931. theme(plot.title = element_text(hjust = -0.065,vjust = 0.25),title = element_text(size = 12))+
1932. geom_node_text(aes(label = pos$nodes), size = 5, vjust = ifelse(1:35 <= 23, -0.5, 1),
1933. angle = ifelse(1:35 <= 23, 45, 0),hjust = ifelse(1:35 <= 23, 0, 0.5))+
1934. expand_limits(y=c(0.2,1.2),x=c(0,2.6))
1935. #p26<-p26+ggtitle('E')
1936. p26
1937. pdf(file = 'D:/microbiome/micom/202404/paper/figures/fig4_9.pdf',width = 15,height = 25)
1938. grid.arrange(arrangeGrob(p18,arrangeGrob(p22,p24,nrow = 2),ncol=2),p25,nrow=2)
1939. dev.off()
1940. ################
1941. meta_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_H_used.csv',row.names = 1)
1942. meta_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_D_used.csv',row.names = 1)
1943. colnames(meta_H)[5]<-'disease'
1944. colnames(meta_D)[5]<-'disease'
1945. meta_H$disease<-'Healthy'
1946. meta_D$disease<-'Depression'
1947. meta<-rbind(meta_D,meta_H)
1948. abundance_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/depression_genus.csv',row.names = 1)
1949. abundance_D<-abundance_D[match(rownames(meta_D),colnames(abundance_D))]
1950. temp<-abundance_D
1951. temp[temp!=0]<-1
1952. abundance_D<-abundance_D[which(rowSums(temp)>ncol(abundance_D)*0.1)%>%as.numeric(),]
1953. abundance_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/healthy_genus.csv',row.names = 1)
1954. abundance_H<-abundance_H[match(rownames(meta_H),colnames(abundance_H))]
1955. temp<-abundance_H
1956. temp[temp!=0]<-1
1957. #abundance_H<-abundance_H[which(rowSums(temp)>ncol(abundance_H)*0.1)%>%as.numeric(),]
1958. abundance<-merge(abundance_D,abundance_H,by = "row.names",all = TRUE)
1959. abundance[is.na(abundance)]<-0
1960. rownames(abundance)<-abundance$Row.names
1961. abundance$Row.names<-NULL
1962. abundance_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/depression_genus.csv',row.names = 1)
1963. abundance_D<-abundance_D[match(rownames(meta_D),colnames(abundance_D))]
1964. temp<-abundance_D
1965. temp[temp!=0]<-1
1966. abundance_D<-abundance_D[which(rowSums(temp)>ncol(abundance_D)*0.1)%>%as.numeric(),]
1967. abundance_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/healthy_genus.csv',row.names = 1)
1968. abundance_H<-abundance_H[match(rownames(meta_H),colnames(abundance_H))]
1969. temp<-abundance_H
1970. temp[temp!=0]<-1
1971. #abundance_H<-abundance_H[which(rowSums(temp)>ncol(abundance_H)*0.1)%>%as.numeric(),]
1972. abundance<-merge(abundance_D,abundance_H,by = "row.names",all = TRUE)
1973. abundance[is.na(abundance)]<-0
1974. rownames(abundance)<-abundance$Row.names
1975. abundance$Row.names<-NULL
1976. files <- list.files('D:/microbiome/micom/202404/paper/result/1-1/FBA_result/depression/', full.names = TRUE)
1977. df_list <- vector("list", length(files))
1978. for(i in seq_along(files)){
1979. print(i)
1980. sampleid <- tools::file_path_sans_ext(basename(files[i]))
1981. temp_m <- read.csv(files[i], row.names = 1)
1982. colnames(temp_m) <- 'flux'
1983. temp_m$sampleid <- sampleid
1984. temp_m$genus <- NA
1985. temp_m$metabolite <- NA
1986. temp_a <- abundance[which(abundance[, sampleid, drop = FALSE] > 0), , drop = FALSE]
1987. genus <- rownames(temp_a)
1988. for(gen in genus){
1989. indices <- grep(gen, rownames(temp_m))
1990. if(length(indices) > 0){
1991. temp_m$genus[indices] <- gen
1992. temp_m$metabolite[indices] <- substr(rownames(temp_m)[indices], nchar(gen) + 1, nchar(rownames(temp_m)[indices]))
1993. }
1994. }
1995. df_list[[i]] <- temp_m
1996. }
1997. df_depression <- do.call(rbind, df_list)
1998. ##########
1999. taxon<-read.csv('D:/microbiome/micom/202404/paper/data/taxa.csv')
2000. tsnedata<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/tsne_data.csv',row.names = 1)
2001. tsnedata$group<-meta$disease[match(tsnedata$sampleid,rownames(meta))]
2002. #tsnedata<-tsnedata[tsnedata$genus%in%diff,]
2003. tsnedata$phyla<-taxon$p__[match(tsnedata$genus,taxon$g__)]
2004. ggplot(tsnedata, aes(x = TSNE.1, y = TSNE.2,color=genus, shape = group)) +
2005. geom_point(size = 1,alpha=0.6) +
2006. theme_minimal() +
2007. theme(legend.position = 'none')
2008. ##########
2009. mes<-read.csv('D:/microbiome/micom/202404/paper/result/1-1/mes.csv',row.names = 1)
2010. mes$group<-meta$disease[match(mes$sampleid,rownames(meta))]
2011. mes_diff<-mes[mes$metabolite%in%diff_metabolite,]
2012. mes_diff<-mes_diff[mes_diff$MES>0,]
2013. mes<-mes[mes$MES>0,]
2014. ggplot(mes_diff,aes(x=MES,fill=group))+
2015. geom_density(alpha=0.6)+
2016. theme_minimal()+
2017. theme(panel.grid = element_blank(),
2018. axis.line = element_line(),
2019. axis.ticks = element_line())+
2020. scale_fill_manual(values=c('Depression'='#E7B800', 'Healthy'='#00AFBB'))+
2021. xlim(c(0,6))
2024. metabolite_D_diff<-metabolite_D[diff_metabolite,]
2025. metabolite_H_diff<-metabolite_H[diff_metabolite,]
2026. rowMeans(metabolite_D_diff)/rowMeans(metabolite_H_diff)
2028. ###########
2029. all_genus<-unique(tsnedata$genus)
2030. res_tsne<-c()
2031. for(i in 1:length(all_genus)){
2032. print(i)
2033. tsne1_d<-tsnedata_d$TSNE.1[which(tsnedata_d$genus==all_genus[i])]
2034. tsne1_h<-tsnedata_h$TSNE.1[which(tsnedata_h$genus==all_genus[i])]
2035. freq_d<-length(tsne1_d)
2036. freq_h<-length(tsne1_h)
2037. if(freq_d>0&freq_h>0){
2038. print(shapiro.test(tsne1_d)$p.val)
2039. print(shapiro.test(tsne1_h)$p.val)
2040. temp<-c(t.test(tsne1_d,tsne1_h)$p.val,freq_d,freq_h)%>%as.data.frame()
2041. }else{
2042. temp<-c(NA,freq_d,freq_h)%>%as.data.frame()
2043. }
2044. temp<-t(temp)%>%as.data.frame()
2045. colnames(temp)<-c('p.val','freq_d','freq_h')
2046. rownames(temp)<-all_genus[i]
2047. res_tsne<-rbind(res_tsne,temp)
2048. }
2049. res_tsne$p.adj<-NA
2050. res_tsne$p.adj<-p.adjust(res_tsne$p.val,method = 'BH')
2051. ######
2052. taxon<-read.csv('D:/microbiome/micom/202404/paper/data/taxa.csv')
2053. meta_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_H_used.csv',row.names = 1)
2054. meta_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/meta_D_used.csv',row.names = 1)
2055. colnames(meta_H)[5]<-'disease'
2056. colnames(meta_D)[5]<-'disease'
2057. meta_H$disease<-'Healthy'
2058. meta_D$disease<-'Depression'
2059. meta<-rbind(meta_D,meta_H)
2060. abundance_D<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/depression_genus.csv',row.names = 1)
2061. abundance_D<-abundance_D[match(rownames(meta_D),colnames(abundance_D))]
2062. temp<-abundance_D
2063. temp[temp!=0]<-1
2064. #abundance_D<-abundance_D[which(rowSums(temp)>ncol(abundance_D)*0.1)%>%as.numeric(),]
2065. abundance_H<-read.csv('D:/microbiome/micom/202404/paper/data/1-1/healthy_genus.csv',row.names = 1)
2066. abundance_H<-abundance_H[match(rownames(meta_H),colnames(abundance_H))]
2067. temp<-abundance_H
2068. temp[temp!=0]<-1
2069. #abundance_H<-abundance_H[which(rowSums(temp)>ncol(abundance_H)*0.1)%>%as.numeric(),]
2070. abundance<-merge(abundance_D,abundance_H,by = "row.names",all = TRUE)
2071. abundance[is.na(abundance)]<-0
2072. rownames(abundance)<-abundance$Row.names
2073. abundance$Row.names<-NULL
2074. abundance$phyla<-taxon$p__[match(rownames(abundance),taxon$g__)]
2075. abundance<-na.omit(abundance)
2076. Firm<-colSums(abundance[abundance$phyla=='Firmicutes',1:ncol(abundance)-1])
2077. Bact<-colSums(abundance[abundance$phyla=='Bacteroidetes',1:ncol(abundance)-1])
2078. F_B<-(Firm/Bact)%>%as.data.frame()
2079. colnames(F_B)<-'Firmicutes/Bacteroidetes'
2080. F_B$group<-meta$disease[match(rownames(F_B),rownames(meta))]
2081. ggplot(F_B,aes(x=group,y=`Firmicutes/Bacteroidetes`,))+
2082. geom_boxplot()+
2083. ylim(c(0,10))
2084. abundance$phyla<-NULL
2085. shannon<-diversity(t(abundance))%>%as.data.frame()
2086. simpson<-diversity(t(abundance),index = 'simpson')%>%as.data.frame()
2087. abundance[abundance<0.01]<-0
2088. temp<-abundance
2089. temp$phyla<-NULL
2090. temp[temp!=0]<-1
2091. temp1<-rowSums(temp)
2092. abundance<-abundance[rownames(temp)[temp1>53],]
2093. abundance1<-rbind(abundance,shannon,simpson)
2094. a<-colSums(abundance)
2096. b<-colSums(abundance[,1:ncol(abundance)-1])
2097. (a-b)%>%sort()%>%plot()
2098. (a-b)%>%sort()%>%hist()
2100. #t.test(F_B$`Firmicutes/Bacteroidetes`[F_B$group=='Depression'],
2101. # F_B$`Firmicutes/Bacteroidetes`[F_B$group=='Healthy'],paired = T)
2102. #####
2103. abundance1<-t(abundance1)%>%as.data.frame()
2104. abundance1$`F/B`<-F_B$`Firmicutes/Bacteroidetes`
2105. exchange_mat<-read.csv('D:/xy/20250108/exchanges_mat.csv',row.names = 1)
2106. topo<-read.csv('D:/microbiome/Depression/topolog.csv',row.names = 1)
2107. id<-intersect(rownames(abundance1),rownames(topo))
2108. abundance1<-abundance1[id,]
2109. exchange_mat<-exchange_mat[id,]
2110. topo<-topo[id,]
2111. shannon<-shannon[id,]
2112. simpson<-simpson[id,]
2113. data<-cbind(#abundance1,
2114. exchange_mat,
2115. topo
2116. )
2117. #data<-cbind(shannon,simpson,exchange_mat,topo)
2119. data$group<-meta$disease[match(rownames(data),rownames(meta))]
2120. data$Disease.name<-NULL
2121. data_split <- createDataPartition(data$group, p = .70, list = FALSE)
2122. training_data <- data[ data_split,]
2123. testing_data <- data[-data_split,]
2124. train_index <- createFolds(training_data$group, k = 10)
2126. randomForestFit <- training_data %>% train(group ~ .,
2127. method = "rf",
2128. data = .,
2129. tuneLength = 5,
2130. trControl = trainControl(method = "cv", indexOut = train_index))
2131. randomForestFit
2132. pr_rf<-predict(randomForestFit,testing_data)
2133. confusionMatrix(pr_rf,reference = testing_data$group%>%as.factor())
2134. pr_rf<-predict(randomForestFit,testing_data,type = 'prob')
2135. library(pROC)
2136. roc_curve<-roc(testing_data$group,pr_rf[,1])
2137. plot(roc_curve)
2138. sp.obj <- ci.sp(roc_curve, sensitivities=seq(0, 1, .01), boot.n=100)
2139. plot(sp.obj, type="shape", col="#5bd1d7")
2140. text(0,0.5,paste0("AUC:",round(roc_curve[["auc"]],4)))
2141. varImp(randomForestFit)

code1.R at commit ffa78a5, no license · at the source