{"id":310,"date":"2022-03-01T20:10:07","date_gmt":"2022-03-01T12:10:07","guid":{"rendered":"https:\/\/linguopeng.top\/?p=310"},"modified":"2022-03-06T23:13:26","modified_gmt":"2022-03-06T15:13:26","slug":"%e4%bb%a3%e8%b0%a2%e5%9b%be","status":"publish","type":"post","link":"https:\/\/linguopeng.top\/?p=310","title":{"rendered":"\u4ee3\u8c22\u56fe"},"content":{"rendered":"\n
######################################\u5b89\u88c5\u5305#########################\n#~\/MetaboAnalyst\uff0cMetaboAnalyst\/cq\n# # install.packages('devtools')\n# # options(BioC_mirror=\"https:\/\/mirrors.tuna.tsinghua.edu.cn\/bioconductor\/\")\n# # options(\"repos\" = c(CRAN=\"http:\/\/mirrors.cloud.tencent.com\/CRAN\/\"))\n# # options(download.file.method = 'libcurl')\n# # options(url.method='libcurl')\n# package_list <- c(\"tidyverse\",\"readxl\",\n#                   \"MetaboAnalystR\",\"FELLA\",\n#                   \"ropls\",\"ggrepel\",\n#                   \"ggforce\",\"httr\",\n#                   \"Hmisc\",\"tidytext\",\n#                   \"psych\",\"corrplot\")\n# for(p in package_list){\n#   if(!suppressWarnings(suppressMessages(require(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))){\n#     install.packages(p,  warn.conflicts = FALSE)\n#     suppressWarnings(suppressMessages(library(p, character.only = TRUE, quietly = TRUE, warn.conflicts = FALSE)))\n#   }\n# }\n# # library(devtools)\n# # install.packages(\"digest\")\n# # devtools::install_github(\"xia-lab\/MetaboAnalystR\", build = TRUE, build_vignettes = FALSE)\n# # BiocManager::install(\"FELLA\")\n# # BiocManager::install(\"ropls\")\n#######################################################################\nlibrary(readxl)\nlibrary(tidyverse)\nlibrary(MetaboAnalystR)\n# setwd(\"~\/Desktop\/2021b\/article_metirials\/metabolomics\/serum\/\")\n#PreparePDFReport(mSetObj = mSet_neg_serum, usrName = \"wangluyao\")\npos_serum <- read_xlsx(\".\/p_n_2.xlsx\",sheet = 1)\nneg_serum <- read_xlsx(\".\/p_n_2.xlsx\",sheet = 2)\nserum_metadata <- read_xlsx(\".\/metadata.xlsx\")\nneg_serum <- neg_serum %>% select(Name, starts_with(\"s\"), starts_with(\"QC\"))\npos_serum <- pos_serum %>% select(Name, starts_with(\"s\"), starts_with(\"QC\"))\nGroup <- c(\"Group\",serum_metadata$Group)\nserum <- rbind(select(neg_serum, -starts_with(\"QC\")),\n               select(pos_serum, -starts_with(\"QC\"))) %>%  \n  group_by(Name) %>% summarise_all(mean) %>% rbind(Group,.)\n\nGroup <- c(\"Group\",serum_metadata$Group,rep(\"QC\",9))\nneg_serum_QC <- neg_serum %>% group_by(Name) %>% summarise_all(mean) %>% rbind(Group,.)\nGroup <- c(\"Group\",serum_metadata$Group,rep(\"QC\",9))\npos_serum_QC <- pos_serum %>% group_by(Name) %>% summarise_all(mean) %>% rbind(Group,.)\nGroup <- c(\"Group\",serum_metadata$Group,rep(\"QC\",9))\nserum_QC <- rbind(select(neg_serum, starts_with(\"Name\"),starts_with(\"s\"),starts_with(\"QC\")),\n               select(pos_serum, starts_with(\"Name\"),starts_with(\"s\"),starts_with(\"QC\"))) %>%  \n  group_by(Name) %>% summarise_all(mean) %>% rbind(Group,.)\n\nwrite_csv(neg_serum_QC,\"neg_serum_QC.csv\")\nwrite_csv(pos_serum_QC,\"pos_serum_QC.csv\")\nwrite_csv(serum_QC,\"serum_QC.csv\")\n\n\n# Standardization and PCA analysis for ESI- and ESI+ respectively\nmSet_neg_serum <- InitDataObjects(data.type = \"pktable\", anal.type = \"stat\", paired = FALSE)\nmSet_neg_serum <- Read.TextData(mSet_neg_serum, \"neg_serum_QC.csv\", format = \"colu\", lbl.type = \"disc\")\nmSet_neg_serum <- SanityCheckData(mSet_neg_serum)\nmSet_neg_serum <- ReplaceMin(mSet_neg_serum);\nmSet_neg_serum <- SanityCheckData(mSet_neg_serum)\nmSet_neg_serum <- FilterVariable(mSet_neg_serum, \"none\", \"F\", 25)\nmSet_neg_serum <- PreparePrenormData(mSet_neg_serum)\nmSet_neg_serum <- Normalization(mSet_neg_serum, \"MedianNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\nmSet_neg_serum <- PlotNormSummary(mSet_neg_serum, \"norm_neg_serum\", \"png\", 300, width=NA)\nmSet_neg_serum <- PlotSampleNormSummary(mSet_neg_serum, \"snorm_neg_serum\", \"png\", 300, width=NA)\nneg_serum_norm <- mSet_neg_serum$dataSet$norm\nmSet_neg_serum <- PCA.Anal(mSet_neg_serum)\nmSet_neg_serum <- PlotPCA2DScore(mSet_neg_serum, \"pca_neg_serum\", \"png\", 300,\n                                 width=NA, pcx = 1, pcy = 2, reg = 0.95, show = 0, grey.scale = 0)\n\nmSet_pos_serum <- InitDataObjects(data.type = \"pktable\", anal.type = \"stat\", paired = FALSE)\nmSet_pos_serum <- Read.TextData(mSet_pos_serum, \"pos_serum_QC.csv\", format = \"colu\", lbl.type = \"disc\")\nmSet_pos_serum <- SanityCheckData(mSet_pos_serum)\nmSet_pos_serum <- ReplaceMin(mSet_pos_serum);\nmSet_pos_serum <- SanityCheckData(mSet_pos_serum)\nmSet_pos_serum <- FilterVariable(mSet_pos_serum, \"none\", \"F\", 25)\nmSet_pos_serum <- PreparePrenormData(mSet_pos_serum)\nmSet_pos_serum <- Normalization(mSet_pos_serum, \"MedianNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\nmSet_pos_serum <- PlotNormSummary(mSet_pos_serum, \"norm_pos_serum\", \"png\", 300, width=NA)\nmSet_pos_serum <- PlotSampleNormSummary(mSet_pos_serum, \"snorm_pos_serum\", \"png\", 300, width=NA)\npos_serum_norm <- mSet_pos_serum$dataSet$norm\nmSet_pos_serum <- PCA.Anal(mSet_pos_serum)\nmSet_pos_serum <- PlotPCA2DScore(mSet_pos_serum, \"pca_pos_serum\", \"png\", 300,\n                                 width=NA, pcx = 1, pcy = 2, reg = 0.95, show = 0, grey.scale = 0)\n\n# Standardization and differential analysis for concatenated table\nmSet_serum <- InitDataObjects(data.type = \"pktable\", anal.type = \"stat\", paired = FALSE)\nmSet_serum <- Read.TextData(mSet_serum, \"serum_QC.csv\", format = \"colu\", lbl.type = \"disc\")\nmSet_serum <- SanityCheckData(mSet_serum)\nmSet_serum <- ReplaceMin(mSet_serum);\nmSet_serum <- SanityCheckData(mSet_serum)\nmSet_serum <- FilterVariable(mSet_serum, \"none\", \"F\", 25)\nmSet_serum <- PreparePrenormData(mSet_serum)\nmSet_serum <- Normalization(mSet_serum, \"MedianNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\nmSet_serum <- PlotNormSummary(mSet_serum, \"norm_serum\", \"png\", 300, width=NA)\nmSet_serum <- PlotSampleNormSummary(mSet_serum, \"snorm_serum\", \"png\", 300, width=NA)\nmSet_serum <- PCA.Anal(mSet_serum)\nmSet_serum <- PlotPCA2DScore(mSet_serum, \"pca_serum\", \"png\", 300,\n                             width=NA, pcx = 1, pcy = 2, reg = 0.95, show = 0, grey.scale = 0)\n\n# comparisons <- unique(serum_metadata$Group)\n# for(comp in 1:length(comparisons)){\n#   tmp_tab <- serum[,-which(str_detect(serum[1,],comparisons[comp]))]\n#   tmp_grp <- unique(as.character(tmp_tab[1,-1]))\n#   tmp_comp <- paste(tmp_grp[1],tmp_grp[2],sep = \"vs\")\n#   write_csv(tmp_tab,paste(tmp_comp,\"csv\",sep = \".\"))\n# } \n###########mSet_CH,CONvsCCFM161########\n# mSet_CH <- PlotOPLS.Permutation(mSet_CH, \"opls_perm_1_\", \"png\", 72, width=NA)\n# mSet<-PlotCmpdSummary(mSet, \"L-Threonic acid\",\"NA\", 1, \"png\", 72, width=NA)\n#######################CONvsI##################\nCONvsI <- serum[,c(1,5:8,13:16)]\nwrite_csv(CONvsI,\"CONvsI.csv\")\nmSet_CI <- InitDataObjects(data.type = \"pktable\", anal.type = \"stat\", paired = FALSE)\nmSet_CI <- Read.TextData(mSet_CI, \"CONvsI.csv\", format = \"colu\", lbl.type = \"disc\")\nmSet_CI <- SanityCheckData(mSet_CI)\nmSet_CI <- ReplaceMin(mSet_CI)\nmSet_CI <- SanityCheckData(mSet_CI)\nmSet_CI <- FilterVariable(mSet_CI, \"none\", \"F\", 25)\nmSet_CI <- PreparePrenormData(mSet_CI)\nmSet_CI <- Normalization(mSet_CI, \"MedianNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\n\nmSet_CI <- Ttests.Anal(mSet_CI, nonpar = FALSE, paired = FALSE, equal.var = TRUE,\n                       pvalType = \"fdr\", all_results = TRUE, threshp = 0.05)\nmSet_CI <- FC.Anal(mSet_CI, fc.thresh = 2, cmp.type = 0, paired = FALSE)\nmSet_CI <- Volcano.Anal(mSet_CI, paired = FALSE, fcthresh = 2.0, cmpType = 0, \n                        nonpar = FALSE, threshp = 0.05, equal.var = TRUE, pval.type = \"raw\")\nmSet_CI <- PlotVolcano(mSet_CI , \"volcano_CI\", 1, 0, \"png\", 300, width=NA)\n\nmSet_CI <- OPLSR.Anal(mSetObj = mSet_CI, reg = TRUE)\nmSet_CI <- PlotOPLS2DScore(mSet_CI, \"opls_CI\", \"png\", 300, width=NA, 1, 2, 0.95, 0, 0)\n#######################CONvsT##################\nCONvsT <- serum[,c(1,5:8,22:24)]\nwrite_csv(CONvsT,\"CONvsT.csv\")\nmSet_CT <- InitDataObjects(data.type = \"pktable\", anal.type = \"stat\", paired = FALSE)\nmSet_CT <- Read.TextData(mSet_CT, \"CONvsT.csv\", format = \"colu\", lbl.type = \"disc\")\nmSet_CT <- SanityCheckData(mSet_CT)\nmSet_CT <- ReplaceMin(mSet_CT)\nmSet_CT <- SanityCheckData(mSet_CT)\nmSet_CT <- FilterVariable(mSet_CT, \"none\", \"F\", 25)\nmSet_CT <- PreparePrenormData(mSet_CT)\nmSet_CT <- Normalization(mSet_CT, \"MedianNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\n\nmSet_CT <- Ttests.Anal(mSet_CT, nonpar = FALSE, paired = FALSE, equal.var = TRUE,\n                       pvalType = \"fdr\", all_results = TRUE, threshp = 0.05)\nmSet_CT <- FC.Anal(mSet_CT, fc.thresh = 2, cmp.type = 0, paired = FALSE)\nmSet_CT <- Volcano.Anal(mSet_CT, paired = FALSE, fcthresh = 2.0, cmpType = 0, \n                        nonpar = FALSE, threshp = 0.05, equal.var = TRUE, pval.type = \"raw\")\nmSet_CT <- PlotVolcano(mSet_CT , \"volcano_CT\", 1, 0, \"png\", 300, width=NA)\n\nmSet_CT <- OPLSR.Anal(mSetObj = mSet_CT, reg = TRUE)\nmSet_CT <- PlotOPLS2DScore(mSet_CT, \"opls_T\", \"png\", 300, width=NA, 1, 2, 0.95, 0, 0)\n\n\nsig_ci <- as.data.frame(mSet_CI$analSet$volcano$sig.mat)\nsig_ct <- as.data.frame(mSet_CT$analSet$volcano$sig.mat)\n\nvip_ci <- as.data.frame(mSet_CI$analSet$oplsda$vipVn)\nvip_ci <- vip_ci[rownames(sig_ci),]\nsig_ci <- cbind(sig_ci,vip_ci)\n\nvip_ct <- as.data.frame(mSet_CT$analSet$oplsda$vipVn)\nvip_ct <- vip_ct[rownames(sig_ct),]\nsig_ct <- cbind(sig_ct,vip_ct)\n\n# metadata <- read_xlsx(\"~\/Desktop\/2021b\/article_metirials\/pheno_data.xlsx\")\n# metadata$Age <- as.numeric(metadata$Age)\n# metadata$BMI <- as.numeric(metadata$BMI)\n# impute <- function(x,x.impute){\n#   ifelse(is.na(x),x.impute,x)\n# }\n# metadata$Age <- impute(metadata$Age,mean(na.omit(metadata$Age)))\n# metadata$BMI <- impute(metadata$BMI,mean(na.omit(metadata$BMI)))\n# write_csv(x = metadata %>% filter(ID %in% colnames(serum)),file = \"metadata.csv\")\n# \n# mSet_cov <- InitDataObjects(\"pktable\",\"ts\",FALSE)\n# mSet_cov <- SetDesignType(mSet_cov,\"multi\")\n# mSet_cov <- Read.TextDataTs(mSet_cov,\"serum.csv\",\"colts\")\n# mSet_cov <- ReadMetaData(mSet_cov,metafilename = \".\/metadata.csv\")\n# mSet_cov <- SanityCheckData(mSet_cov)\n# mSet_cov <- ReplaceMin(mSet_cov)\n# mSet_cov <- SanityCheckMeta(mSet_cov, 1)\n# mSet_cov <- SetDataTypeOfMeta(mSet_cov)\n# mSet_cov <- SanityCheckData(mSet_cov)\n# mSet_cov <- FilterVariable(mSet_cov, \"none\", \"F\", 25)\n# mSet_cov <- PreparePrenormData(mSet_cov)\n# mSet_cov <- Normalization(mSet_cov, \"SumNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\n# adj.vec <- c(\"Gender\",\"BMI\",\"Age\")\n# mSet_cov <- CovariateScatter.Anal(mSet_cov, \"covariate_plot_0_dpi72.png\", \"png\", 72, \"default\", \"Group\", \"Health\", \"NA\" , 0.05, FALSE)\n# maaslin <- as.data.frame(mSet_cov$analSet$cov$sig.mat)\n\ndiff_ct <- sig_ct %>% rownames_to_column(var = \"Name\") %>% \n  filter(., abs(`log2(FC)`)>1, raw.pval<0.05, vip_ct>1) %>% \n  pull(Name)\ndiff_ci <- sig_ci %>% rownames_to_column(var = \"Name\") %>% \n  filter(., abs(`log2(FC)`)>1, raw.pval<0.05, vip_ci>1) %>% \n  pull(Name)\n\ndiff_cp_data <- serum %>% filter(Name %in% c(\"Group\",unique(c(diff_ct,diff_ci))))\n\n# Enrichment analysis\nlibrary(httr)\ncall <- \"http:\/\/api.xialab.ca\/mapcompounds\"\ntoSend = list(queryList = diff_cp_data$Name, inputType = \"name\")\nquery_results <- httr::POST(call, body = toSend, encode = \"json\")\nquery_results$status_code==200\nquery_results_text <- content(query_results, \"text\", encoding = \"UTF-8\")\nquery_results_json <- rjson::fromJSON(query_results_text,simplify = TRUE)\nquery_results_table <- cbind.data.frame(query_results_json$Query,query_results_json$Match,\n                                        as.character(query_results_json$KEGG)) %>% \n  rename(query=colnames(.)[1],match=colnames(.)[2],kegg=colnames(.)[3])\n\n# diff_cp_data$Name[which(diff_cp_data$Name %in% c(\"2-Methyl-4,6-dinitrophenol\",\"D-(-)-Quinic acid\",\n#                                                  \"Dodecyl sulfate\",\"L-Threonic acid\"))] <-\n#   c(\"4,6-Dinitro-o-cresol\",\"Quinate\",\"Erythromycin estolate\",\"Threonate\")\n\ntmp_query <- query_results_table %>% \n  filter(kegg!=\"NA\",kegg!=\"NULL\")\n# diff_cp_data <- diff_cp_data %>% filter(Name %in% c(\"Group\",\"4,6-Dinitro-o-cresol\",\"Quinate\",\n#                                                     \"Erythromycin estolate\",\"Threonate\",tmp_query$query))\ndiff_cp_data <- diff_cp_data %>% filter(Name %in% c(\"Group\",tmp_query$query))\nwrite_csv(diff_cp_data,file = \"diff_cp_data.csv\")\n\nmSet_en <- InitDataObjects(\"conc\",\"msetqea\",FALSE)\nmSet_en <- Read.TextData(mSet_en,\"diff_cp_data.csv\",\"colu\",\"disc\")\nmSet_en <- SanityCheckData(mSet_en)\nmSet_en <- ReplaceMin(mSet_en)\nmSet_en <- CrossReferencing(mSet_en, \"name\")\nmSet_en <- CreateMappingResultTable(mSet_en)\nmSet_en <- PreparePrenormData(mSet_en)\nmSet_en <- Normalization(mSet_en, \"MedianNorm\", \"LogNorm\", \"NULL\", ratio=FALSE, ratioNum=20)\nmSet_en <- SetMetabolomeFilter(mSet_en, F)\nmSet_en <- SetCurrentMsetLib(mSet_en, \"kegg_pathway\", 2)\nmSet_en <- CalculateGlobalTestScore(mSet_en)\nmSet_en <- PlotQEA.Overview(mSet_en, \"qea_bar\", \"net\", \"png\", 300, width=NA)\nmSet_en <- PlotEnrichDotPlot(mSet_en, \"qea\", \"qea_dot\", \"png\", 300, width=NA)\nmSet_en <- PlotQEA.MetSet(mSet_en, \"Biosynthesis of unsaturated fatty acids\", \"png\", 300, width=NA)\n\nmSet_pw <- InitDataObjects(\"conc\", \"pathora\", FALSE)\ncmpd_vec <- diff_cp_data$Name[-1]\nmSet_pw <- Setup.MapData(mSet_pw, cmpd_vec)\nmSet_pw <- CrossReferencing(mSet_pw, \"name\")\nmSet_pw <- CreateMappingResultTable(mSet_pw)\nmSet_pw <- SetKEGG.PathLib(mSet_pw, \"hsa\", \"current\")\nmSet_pw <- SetMetabolomeFilter(mSet_pw, F)\nmSet_pw <- CalculateOraScore(mSet_pw, \"dgr\", \"hyperg\")\nmSet_pw <- PlotPathSummary(mSet_pw, F, \"pathway\", \"png\", 300, width=NA,NA,NA)\n\n# mSet_nea <- InitDataObjects(\"conc\", \"msetora\", FALSE)\n# mSet_nea <- Setup.MapData(mSet_nea, diff_cp_data$Name[-1])\n# mSet_nea <- CrossReferencing(mSet_nea, \"name\")\n# mSet_nea <- CreateMappingResultTable(mSet_nea)\n# mSet_nea <- SetMetabolomeFilter(mSet_nea, F)\n# mSet_nea <- SetCurrentMsetLib(mSet_nea, \"kegg_pathway\", 2)\n# mSet_nea <- CalculateHyperScore(mSet_nea)\n# mSet_nea <- PlotORA(mSet_nea, \"bar_nea\", \"net\", \"png\", 300, width=NA)\n# mSet_nea <- PlotEnrichDotPlot(mSet_nea, \"ora\", \"dot_nea\", \"png\", 300, width=NA)\n###### FC vs. CT ######\n# diff_ch[which(diff_ch==\"2-Methyl-4,6-dinitrophenol\")] <- \"4,6-Dinitro-o-cresol\"\n# diff_ch[which(diff_ch==\"D-(-)-Quinic acid\")] <- \"Quinate\"\n# diff_ch[which(diff_ch==\"Dodecyl sulfate\")] <- \"Erythromycin estolate\"\n# diff_ch[which(diff_ch==\"L-Threonic acid\")] <- \"Threonate\"\n# mSet_en_ch <- InitDataObjects(\"conc\", \"msetora\", FALSE)\n# mSet_en_ch <- Setup.MapData(mSet_en_ch, diff_ch)\n# mSet_en_ch <- CrossReferencing(mSet_en_ch, \"name\")\n# mSet_en_ch <- CreateMappingResultTable(mSet_en_ch)\n# mSet_en_ch <- SetMetabolomeFilter(mSet_en_ch, F)\n# mSet_en_ch <- SetCurrentMsetLib(mSet_en_ch, \"kegg_pathway\", 2)\n# mSet_en_ch <- CalculateHyperScore(mSet_en_ch)\n# mSet_en_ch <- PlotORA(mSet_en_ch, \"bar_en_ch\", \"net\", \"png\", 300, width=NA)\n# mSet_en_ch <- PlotEnrichDotPlot(mSet_en_ch, \"ora\", \"dot_en_ch\", \"png\", 300, width=NA)\n\n###### FC vs. IBS-C ######\n# diff_ci[which(diff_ci==\"2-Methyl-4,6-dinitrophenol\")] <- \"4,6-Dinitro-o-cresol\"\n# diff_ci[which(diff_ci==\"D-(-)-Quinic acid\")] <- \"Quinate\"\n# diff_ci[which(diff_ci==\"Dodecyl sulfate\")] <- \"Erythromycin estolate\"\n# diff_ci[which(diff_ci==\"L-Threonic acid\")] <- \"Threonate\"\nmSet_en_ci <- InitDataObjects(\"conc\", \"msetora\", FALSE)\nmSet_en_ci <- Setup.MapData(mSet_en_ci, diff_ci)\nmSet_en_ci <- CrossReferencing(mSet_en_ci, \"name\")\nmSet_en_ci <- CreateMappingResultTable(mSet_en_ci)\nmSet_en_ci <- SetMetabolomeFilter(mSet_en_ci, F)\nmSet_en_ci <- SetCurrentMsetLib(mSet_en_ci, \"kegg_pathway\", 2)\nmSet_en_ci <- CalculateHyperScore(mSet_en_ci)\nmSet_en_ci <- PlotORA(mSet_en_ci, \"bar_en_ci\", \"net\", \"png\", 300, width=NA)\nmSet_en_ci <- PlotEnrichDotPlot(mSet_en_ci, \"ora\", \"dot_en_ci\", \"png\", 300, width=NA)\n\n###### IBS-C vs. HC ######\n# diff_ih[which(diff_ih==\"2-Methyl-4,6-dinitrophenol\")] <- \"4,6-Dinitro-o-cresol\"\n# diff_ih[which(diff_ih==\"D-(-)-Quinic acid\")] <- \"Quinate\"\n# diff_ih[which(diff_ih==\"Dodecyl sulfate\")] <- \"Erythromycin estolate\"\n# diff_ih[which(diff_ih==\"L-Threonic acid\")] <- \"Threonate\"\n\nmSet_en_ct <- InitDataObjects(\"conc\", \"msetora\", FALSE)\nmSet_en_ct <- Setup.MapData(mSet_en_ct, diff_ct)\nmSet_en_ct <- CrossReferencing(mSet_en_ct, \"name\")\nmSet_en_ct <- CreateMappingResultTable(mSet_en_ct)\nmSet_en_ct <- SetMetabolomeFilter(mSet_en_ct, F)\nmSet_en_ct <- SetCurrentMsetLib(mSet_en_ct, \"kegg_pathway\", 2)\nmSet_en_ct <- CalculateHyperScore(mSet_en_ct)\nmSet_en_ct <- PlotORA(mSet_en_ct, \"bar_en_ct\", \"net\", \"png\", 300, width=NA)\nmSet_en_ct <- PlotEnrictDotPlot(mSet_en_ct, \"ora\", \"dot_en_ct\", \"png\", 300, width=NA)\n\n### volcano plot visulization ###\nlibrary(ggrepel)\nlibrary(ggforce)\nlibrary(Hmisc)\n\nvolcano_ct <- cbind(mSet_CT$analSet$volcano$fc.all,mSet_CT$analSet$volcano$fc.log,\n                    mSet_CT$analSet$volcano$p.log,mSet_CT$analSet$oplsda$vipVn) %>% \n  as.data.frame() %>% \n  rownames_to_column(\"name\") %>% \n  rename(fc=V1,log_fc=V2,log_p=V3,vip=V4) %>% \n  mutate(cp_tp=if_else(\n    log_p < 1.30103,\"NoDiff\",if_else(\n      abs(log_fc) < 1,\"NoDiff\",\n      if_else(log_fc >=1,\"Upward\",\"Downward\"))),\n    comparison=\"Control vs. T\")\nvol_num_ct <- volcano_ct %>% group_by(cp_tp) %>% summarise(count=n())\n\nvolcano_ci <- cbind(mSet_CI$analSet$volcano$fc.all,mSet_CI$analSet$volcano$fc.log,\n                    mSet_CI$analSet$volcano$p.log,mSet_CI$analSet$oplsda$vipVn) %>% \n  as.data.frame() %>%\n  rownames_to_column(\"name\") %>% \n  rename(fc=V1,log_fc=V2,log_p=V3,vip=V4) %>% \n  mutate(cp_tp=if_else(\n    log_p < 1.30103,\"NoDiff\",if_else(\n      abs(log_fc) < 1,\"NoDiff\",\n      if_else(log_fc >=1,\"Upward\",\"Downward\"))),\n    comparison=\"Control vs. I\")\nvol_num_ci <- volcano_ci %>% group_by(cp_tp) %>% summarise(count=n())\n\nvolcano <- rbind(volcano_ci,volcano_ct)\nvolcano$cp_tp <- factor(volcano$cp_tp,levels = c(\"Upward\",\"NoDiff\",\"Downward\"))\nvolcano$name <- capitalize(volcano$name)\nvolcanolabeld<-volcano[volcano$cp_tp==\"Downward\",]\nvolcanolabelu<-volcano[volcano$cp_tp==\"Upward\",]\nvolcanolabel<-rbind(volcanolabeld,volcanolabelu)\np1 <- ggplot(volcano,aes(log_fc,log_p)) + \n  geom_point(aes(color=cp_tp,size=abs(vip)),alpha=0.8)+\n  geom_hline(yintercept = -log10(0.05),linetype =\"dashed\") + \n  geom_vline(xintercept = c(-1,1),linetype = \"dashed\")+\n  # geom_mark_circle(data=volcano %>% \n  #                    filter(name %in% mSet_en$analSet$qea.hits$`Biosynthesis of unsaturated fatty acids`),\n  #                  aes(fill = name,\n  #                      label = name),\n  #                  expand = unit(1.5, \"mm\"),\n  #                  alpha = 0,\n  #                  con.cap = 0,\n  #                  label.lineheight = 0.1,\n  #                  label.fontsize = c(5, 4),\n  #                  show.legend = FALSE)+\n  geom_label_repel(data=volcano %>%\n                     filter(name %in% mSet_en$analSet$qea.hits$`Biosynthesis of unsaturated fatty acids`),\n                   size=2,\n                   box.padding = 1, #\u5b57\u5230\u70b9\u7684\u8ddd\u79bb\n                   point.padding = 0.1, #\u5b57\u5230\u70b9\u7684\u8ddd\u79bb\uff0c\u70b9\u5468\u56f4\u7684\u7a7a\u767d\u5bbd\u5ea6\n                   min.segment.length = 0, #\u77ed\u7ebf\u6bb5\u53ef\u4ee5\u7701\u7565\n                   segment.color = \"black\", #segment.colour = NA, \u4e0d\u663e\u793a\u7ebf\u6bb5\n                   aes(label=name,color=cp_tp),\n                   arrow = arrow(length=unit(0.01, \"npc\")),\n                   show.legend = FALSE)+\n  facet_wrap(.~comparison,nrow=1,scales = \"free_y\")+\n  scale_color_manual(values = c(\"Upward\"=\"#DC0000FF\",\"NoDiff\"=\"grey\",\"Downward\"=\"#00A087FF\"))+\n  scale_size(range = c(0.5,3),breaks = c(0.5,1,2,3))+\n  guides(color=guide_legend(title = \"Regulation\",order = 1),\n         size=guide_legend(title = \"VIP\",order = 2))+\n  labs(x=\"log2(Fold Change)\",\n       y=\"-log10(pvalue)\")+\n  theme_bw()\np1\nggsave(\"diff_serum_volcano.pdf\",p1,path = \".\/\",units = \"mm\",width = 210,height = 99)\n\n### Enrichment analysis ###\nlibrary(tidytext)\nendot_ct <- as.data.frame(mSet_en_ct$analSet$ora.mat) %>% \n  rownames_to_column(\"name\") %>% \n  mutate(ratio=hits\/expected,\n         comparison=\"C vs. T\",\n         log_p=-log10(`Raw p`))\n# endot_ih <- as.data.frame(mSet_en_ih$analSet$ora.mat) %>% \n#   rownames_to_column(\"name\") %>% \n#   mutate(ratio=hits\/expected,\n#          comparison=\"IBS-C vs. HC\",\n#          log_p=-log10(`Raw p`))\n# endot_ci <- as.data.frame(mSet_en_ci$analSet$ora.mat) %>% \n#   rownames_to_column(\"name\") %>% \n#   mutate(ratio=hits\/expected,\n#          comparison=\"FC vs. IBS-C\",\n#          log_p=-log10(`Raw p`))\n# endot <- rbind(endot_ch[1:10,],endot_ci[1:10,],endot_ih[1:10,])\nendot <- endot_ct\nendot$name <- capitalize(endot$name)\n\np2 <- ggplot(endot,aes(x=log_p,y=reorder_within(name,log_p,comparison)))+\n  geom_point(aes(color=`Raw p`,size=ratio))+\n  geom_vline(xintercept = -log10(0.05),linetype = \"dashed\")+\n  facet_wrap(.~comparison,nrow = 3,strip.position = \"right\",scales = \"free_y\")+\n  labs(x=\"-log10(pvalue)\")+\n  scale_y_reordered()+\n  #scale_color_gradientn(colors = color_vector_use,trans=\"reverse\")+\n  guides(color=guide_colorbar(title = \"P value\"),\n         size=guide_legend(title = \"Enrichmenr ratio\"))+\n  theme_bw()+\n  theme(panel.grid.minor = element_blank(),\n        axis.line.x = element_line(size = 0.3, colour = \"black\"),#x\u8f74\u8fde\u7ebf\n        axis.ticks.length.x = unit(-0.20, \"cm\"),#\u4fee\u6539x\u8f74\u523b\u5ea6\u7684\u9ad8\u5ea6\uff0c\u8d1f\u53f7\u8868\u793a\u5411\u4e0a\n        axis.text.x = element_text(margin = margin(t = 0.3, unit = \"cm\")),#\u7ebf\u4e0e\u6570\u5b57\u4e0d\u8981\u91cd\u53e0\n        axis.ticks.x = element_line(colour = \"black\",size = 0.3),#\u4fee\u6539x\u8f74\u523b\u5ea6\u7684\u7ebf                         \n        axis.ticks.y = element_blank(),\n        axis.title.y  = element_blank())\np2\nggsave(\"serum_enrichment.pdf\",p2,path = \".\/\",units = \"mm\",width = 210,height = 139)\n\n### Unsaturated fatty acids ###\nlibrary(ggpubr)\n\nmypal <- c(\"#0072B5FF\",\"#BC3C29FF\",\"#E18727FF\",\"#20854EFF\")\ncomparisons <- list(c(\"CON\", \"I\"),\n                    c(\"CON\",\"T\"),\n                    c(\"I\",\"T\"))\n\ndiff_usfa <- cbind(mSet_en$dataSet$url.smp.nms,mSet_en$dataSet$orig.cls,mSet_en$analSet$msea.data) %>%\n  rename(ID=colnames(.)[1],Group=colnames(.)[2]) %>% \n  select(ID,Group,mSet_en$analSet$qea.hits$`Biosynthesis of unsaturated fatty acids`) %>% \n  gather(key = \"USFAs\",value = \"Abundance\",3)\ndiff_usfa$Group <- factor(diff_usfa$Group,levels = c(\"CON\", \"I\",\"T\"))\ndiff_usfa <- diff_usfa %>% \n  filter(Group!=\"NA\")\np3 <- ggboxplot(diff_usfa,x=\"Group\",y=\"Abundance\",color = \"Group\",\n                add = \"jitter\",width = .3,size = .7,palette = mypal)+\n  labs(x=\"Group\",y=\"Normalized intensity\")+\n  geom_signif(comparisons = comparisons,\n              map_signif_level=F,vjust=0.5,color=\"black\",\n              textsize=5,test=wilcox.test,step_increase=0.1,tip_length = 0.015)+\n  facet_wrap(.~USFAs,scales = \"free_y\")+\n  theme_bw()+\n  theme(axis.text.x=element_blank(), # remove the x axis scale\n        axis.ticks.x=element_blank(),\n        axis.title.x = element_blank(),\n        legend.title = element_blank(),\n        legend.key=element_blank(),\n        legend.text = element_text(color=\"black\",size=8),\n        legend.spacing.x=unit(0,'cm'),\n        legend.spacing.y = unit(0,'cm'),\n        legend.background=element_blank(),\n        legend.box.background=element_rect(colour = \"black\"),\n        legend.box.margin = margin(1,1,1,1),\n        legend.position = c(0.95,0.55),legend.justification = c(0.95,0.1))\np3\nggsave(\"diff_usfas.pdf\",p3,path = \".\/\",units = \"in\",width = 4.34,height = 3.22)\ndev.off()\n\n###### exercise for FELLA ######\nlibrary(FELLA)\n# A database built from KEGG has to be regenerated for once\nfella_graph <- buildGraphFromKEGGREST(organism = \"hsa\")\nbuildDataFromGraph(keggdata.graph = fella_graph,\n                   databaseDir = \".\/fella_database\",\n                   internalDir = FALSE,\n                   matrices = \"diffusion\",\n                   normality = \"diffusion\",\n                   niter = 50)\n\nfella_data <- loadKEGGdata(databaseDir = \".\/fella_database\",\n                           internalDir = FALSE,\n                           loadMatrix = \"diffusion\")\ncat(getInfo(fella_data))\n# A database built from KEGG has to be regenerated for once\nfella_df <- tmp_query$kegg\nfella_df\n\nfella_anal <- defineCompounds(compounds = fella_df,\n                              data = fella_data)\ngetExcluded(fella_anal)\n\nfella_anal <- runDiffusion(object = fella_anal,\n                           data = fella_data,\n                           approx = \"normality\")\nFELLA::plot(fella_anal,\n            method = \"diffusion\",\n            data = fella_data,\n            nlimit = 150,\n            vertex.label.cex = 0.5)\nexportResults(format = \"csv\",\n              file = \".\/ella_res.csv\",\n              method = \"diffusion\",\n              object = fella_anal,\n              data = fella_data)\n########################################################################################\n###### exercise for FELLA ######\n\n# ## Correlation analysis\n# # Pheno_data \n# library(psych)\n# library(corrplot)\n# cor_cp <- mSet_en$dataSet$norm\n# colnames(cor_cp) <- capitalize(colnames(cor_cp))\n# colnames(cor_cp)[which(str_detect(colnames(cor_cp),\"icosapentaenoic acid\"))] <- \"Eicosapentaenoic acid\"\n# cor_cp <- cor_cp %>% \n#   select(mSet_en$analSet$qea.hits$`Biosynthesis of unsaturated fatty acids`,everything())\n# pheno_data <- metadata %>% \n#   filter(ID %in% rownames(cor_cp)) %>% \n#   column_to_rownames(\"ID\") %>% select(SBM,BSFS,API)\n# \n# cor_list <- corr.test(pheno_data, cor_cp, method = \"spearman\", adjust = \"BH\", alpha = 0.05)\n# r_mat <- cor_list$r\n# p_mat <- cor_list$p.adj\n# \n# corrplot(corr = r_mat,p.mat = p_mat, insig = \"label_sig\",\n#          sig.level = c(.001, .01, .05), pch.cex = .6, pch.col = \"Black\", col = color_vector_use,\n#          tl.col = \"black\",tl.cex = 0.5,tl.srt = 45,diag = TRUE, cl.cex = 0.5, cl.length = 5,\n#          outline = \"Black\")\n# \n# # differential taxa\n# diff_16S <- readRDS(\"~\/Desktop\/2021b\/article_metirials\/microbiome\/diff_16S.rds\")\n# diff_16S <- diff_16S %>% column_to_rownames(\"ID\") %>% select(-group)\n# diff_16S <- log1p(diff_16S)\n# d <- merge(cor_cp, diff_16S, by = \"row.names\", all = FALSE) %>% \n#   column_to_rownames(\"Row.names\")\n# dat1 <- d[,1:47]\n# dat2 <- d[,48:54]\n# \n# cor_list <- corr.test(dat1, dat2, method = \"spearman\", adjust = \"BH\", alpha = 0.05)\n# r_mat <- cor_list$r\n# p_mat <- cor_list$p.adj\n# \n# corrplot(corr = r_mat,p.mat = p_mat, insig = \"label_sig\",\n#          sig.level = c(.001, .01, .05), pch.cex = .6, pch.col = \"Black\", col = color_vector_use,\n#          tl.col = \"black\",tl.cex = 0.6,tl.srt = 45,diag = TRUE, cl.cex = 0.5, outline = \"Black\")\n# \n# d2 <- pheno_data[rownames(d),]\n# \n# cor_list <- corr.test(d2, d, method = \"spearman\", adjust = \"BH\", alpha = 0.05)\n# r_mat <- cor_list$r\n# p_mat <- cor_list$p.adj\n# \n# corrplot(corr = r_mat,p.mat = p_mat, insig = \"label_sig\",\n#          sig.level = c(.001, .01, .05), pch.cex = .6, pch.col = \"Black\", col = color_vector_use,\n#          tl.col = \"black\",tl.cex = 0.6,tl.srt = 45,diag = TRUE, cl.cex = 0.5, cl.length = 5,\n#          outline = \"Black\")\n# \n# save(diff_cp_data,file = \"~\/Desktop\/2021b\/article_metirials\/diff_serum_cmp.rds\")\n# serum_norm <- mSet_serum$dataSet$norm\n# save(serum_norm,file = \"~\/Desktop\/2021b\/article_metirials\/serum_norm.rds\")\n# \n# # ggClusterNet\n# library(phyloseq)\n# library(igraph)\n# library(network)\n# library(sna)\n# library(ggClusterNet)\n# \n# ps_16S <- readRDS(\"~\/Desktop\/2021b\/article_metirials\/ps_16S.rds\")\n# ps_16S <- prune_samples(rownames(d), ps_16S)\n# ps_16S <- filter_OTU_ps(ps_16S, Top = 200)\n# net_grp <- as.data.frame(rep(\"serum_metabolites\",24))\n# rownames(net_grp) <- colnames(dat1)\n# colnames(net_grp)[1] <- \"group\"\n# dat1 <- as.data.frame(dat1) %>% rownames_to_column(\"ID\")\n# cor_net <- corBiostripe(data = dat1, group = net_grp, ps = ps_16S, method = \"spearman\")\n# \n# tmp_otu = as.data.frame(vegan_otu(ps_16S))\n# tmp_cor <- (dat1[-1])\n# row.names(tmp_cor) = dat1[[1]]\n# tmp_cor <- t(tmp_cor)\n# finaldata <- rbind(tmp_otu,tmp_cor)\n# occor = corr.test(dat1,tmp_otu,use=\"pairwise\",method= \"spearman\",adjust=\"fdr\",alpha=.05,)\n# occor.r = occor$r\n# occor.p = t(occor$p.adj)\n# occor.r[occor.p > 0.05|abs(occor.r)<0.6] = 0\n# \n# ######################### Abolished case #######################################\n# # # PCA\n# # library(FactoMineR)\n# # library(ggsci)\n# # library(scales)\n# # \n# # show_col(pal_nejm(palette = \"default\")(8))\n# # mypal <- c(\"#0072B5FF\",\"#BC3C29FF\",\"#E18727FF\",\"#20854EFF\")\n# # \n# # pca_neg_serum <- PCA(neg_serum_norm,graph = FALSE)\n# # pca_neg_serum$eig #\u663e\u793a\u7279\u5f81\u503c\u3001\u65b9\u5dee\u8d21\u732e\u7387\u548c\u7d2f\u8ba1\u65b9\u5dee\u8d21\u732e\u7387\n# # dat_1 <- as.data.frame(pca_neg_serum$ind$coord) #dat_1\u4e3a\u4e3b\u6210\u5206\u5f97\u5206\uff08\u6837\u672c\u5750\u6807\uff09\n# # name <- rownames(dat_1)\n# # Group_neg_serum <- c(serum_metadata$Group,rep(\"QC\",20))\n# # dat_1 <- cbind(dat_1,name,Group_neg_serum) %>% select(-Dim.3,-Dim.4,-Dim.5)\n# # colnames(dat_1) <- c(\"x\",\"y\",\"ID\",\"Group\")\n# # #dat_2 <- as.data.frame(pca$var$coord) #dat_2\u4e3a\u4e3b\u6210\u5206\u8f7d\u8377\uff08\u53d8\u91cf\u5750\u6807\uff09\n# # dat_1$Group <- factor(dat_1$Group,levels = c(\"Health\",\"Constipation\",\"IBS-C\",\"QC\"))\n# # dat_1$x<-as.numeric(dat_1$x)\n# # dat_1$y<-as.numeric(dat_1$y)\n# # \n# # p1 <- ggplot(dat_1, aes(x, y,color = Group)) + #\u6620\u5c04PC1,PC2\n# #   geom_point(size = 3, alpha=.7) + #\u7ed8\u5236\u5f97\u5206\u6563\u70b9\u56fe\n# #   labs(title = \"PCA biplot\") + #\u56fe\u8868\u6807\u9898\n# #   xlab(paste(\"PC1\", paste(round(pca_neg_serum$eig[1,2], 2), \"%\", sep = \"\"), sep = \"  \")) + #x\u8f74\u6807\u9898\n# #   ylab(paste(\"PC2\", paste(round(pca_neg_serum$eig[2,2], 2), \"%\", sep = \"\"), sep = \"  \")) + #y\u8f74\u6807\u9898\n# #   stat_ellipse( geom = \"polygon\", alpha = 0.1, level = 0.95, show.legend = F,\n# #                 col=\"black\",size=1)+\n# #   stat_ellipse(aes(x,y,color=Group),level = .68,inherit.aes = FALSE,size=.8)+#\u52a0\u7f6e\u4fe1\u692d\u5706\n# #   geom_vline(xintercept=mean(dat_1$x),lty=1,col=\"black\",lwd=0.5)+#\u6dfb\u52a0\u6a2a\u7ebf\n# #   geom_hline(yintercept =mean(dat_1$y),lty=1,col=\"black\",lwd=0.5)+#\u6dfb\u52a0\u7ad6\u7ebf\n# #   theme_bw()+ scale_color_manual(values = mypal)+\n# #   theme( plot.background = element_rect(fill=\"white\"),\n# #          panel.background = element_rect(fill='white', colour='gray'),          \n# #          text=element_text(family=\"sans\", size=12),\n# #          plot.title = element_text(hjust = .5))\n# # p1\n# #  \n# # pca_pos_serum <- PCA(pos_serum_norm,graph = FALSE)\n# # dat_2 <- as.data.frame(pca_pos_serum$ind$coord) #dat_2\u4e3a\u4e3b\u6210\u5206\u5f97\u5206\uff08\u6837\u672c\u5750\u6807\uff09\n# # name <- rownames(dat_2)\n# # Group_pos_serum <- c(serum_metadata$Group,rep(\"QC\",13))\n# # dat_2 <- cbind(dat_2,name,Group_pos_serum) %>% select(-Dim.3,-Dim.4,-Dim.5)\n# # colnames(dat_2) <- c(\"x\",\"y\",\"ID\",\"Group\")\n# # #dat_2 <- as.data.frame(pca$var$coord) #dat_2\u4e3a\u4e3b\u6210\u5206\u8f7d\u8377\uff08\u53d8\u91cf\u5750\u6807\uff09\n# # dat_2$Group <- factor(dat_2$Group,levels = c(\"Health\",\"Constipation\",\"IBS-C\",\"QC\"))\n# # dat_2$x<-as.numeric(dat_2$x)\n# # dat_2$y<-as.numeric(dat_2$y)\n# # \n# # p2 <- ggplot(dat_2, aes(x, y,color = Group)) + #\u6620\u5c04PC1,PC2\n# #   geom_point(size = 3, alpha=.7) + #\u7ed8\u5236\u5f97\u5206\u6563\u70b9\u56fe\n# #   labs(title = \"PCA biplot\") + #\u56fe\u8868\u6807\u9898\n# #   xlab(paste(\"PC1\", paste(round(pca_pos_serum$eig[1,2], 2), \"%\", sep = \"\"), sep = \"  \")) + #x\u8f74\u6807\u9898\n# #   ylab(paste(\"PC2\", paste(round(pca_pos_serum$eig[2,2], 2), \"%\", sep = \"\"), sep = \"  \")) + #y\u8f74\u6807\u9898\n# #   stat_ellipse( geom = \"polygon\", alpha = 0.1, level = 0.95, show.legend = F,\n# #                 col=\"black\",size=1)+\n# #   stat_ellipse(aes(x,y,color=Group),level = .68,inherit.aes = FALSE,size=.8)+#\u52a0\u7f6e\u4fe1\u692d\u5706\n# #   geom_vline(xintercept=mean(dat_2$x),lty=1,col=\"black\",lwd=0.5)+#\u6dfb\u52a0\u6a2a\u7ebf\n# #   geom_hline(yintercept =mean(dat_2$y),lty=1,col=\"black\",lwd=0.5)+#\u6dfb\u52a0\u7ad6\u7ebf\n# #   theme_bw()+ scale_color_manual(values = mypal)+\n# #   theme(plot.background = element_rect(fill=\"white\"),\n# #         panel.background = element_rect(fill='white', colour='gray'),          \n# #         text=element_text(family=\"sans\", size=12),\n# #         plot.title = element_text(hjust = .5))\n# # p2\n# # \n# # # OPLS-DA\n# # library(ropls)\n# # # library(o2plsda)\n# # # o2plsda makes no sense, not recommended.\n# # # detach(\"package:o2plsda\",unload = TRUE)\n# # ###### exercise for o2pls-da ######\n# # # cv <- o2cv(X = opl_neg_serum, Y = opl_pos_serum, \n# # #            nc = 1:5, nx = 1:3, ny = 1:3, group = Group_neg_serum, nr_folds = 10)\n# # # fit_serum <- o2pls(X = opl_neg_serum, Y = opl_pos_serum, nc = 5, nx = 2, ny = 3)\n# # # summary(fit_serum)\n# # # Xl <- loadings(fit_serum,loading=\"Xjoint\")\n# # # Xs <- scores(fit_serum,score=\"Xjoint\")\n# # # plot(fit_serum,type=\"score\",var=\"Xjoint\", group=Group_serum,repel = FALSE)\n# # # plot(fit_serum,type=\"loading\",var=\"Xjoint\", group=Group_serum,repel=F,rotation=TRUE)\n# # # res_serum <- oplsda(fit_serum, Group_serum, nc=5)\n# # # plot(res_serum, type=\"score\", group=Group_serum)\n# # # vip <- vip(res_serum)\n# # # plot(res_serum,type=\"vip\", group = Group_serum, repel = FALSE, order=TRUE)\n# # ###### exercise for o2pls-da ######\n# # \n# # \n# # \n# # Group_serum <- c(serum_metadata$Group)\n# # comparisons <- data.frame(c(\"Health\",\"Constipation\",\"IBS-C\"), \n# #                           c(\"Constipation\",\"IBS-C\",\"Health\"))\n# # \n# # for(comp in 1:nrow(comparisons)){\n# #     comp_name <- paste(comparisons[comp, 1], \n# #                      comparisons[comp, 2], \n# #                      sep = \"vs\")\n# #     tmp_tab <- serum\n# #     tmp_tab$Group <- Group_serum\n# #     tmp_tab <- tmp_tab %>% filter(Group %in% c(comparisons[comp,1],comparisons[comp,2]))\n# #     tmp_grp <- tmp_tab$Group\n# #     tmp_tab <- tmp_tab %>% select(-Group)\n# #   \n# #     comp_name <- opls(tmp_tab, tmp_grp, predI = 1, orthoI = 1, crossvalI =10,\n# #                       scaleC=\"none\")\n# #     \n# #     VIP<-comp_name@vipVn\n# #     t1<-comp_name@orthoScoreMN\n# #     to1<-comp_name@scoreMN\n# #     \n# #     j<-data.frame(t1,to1,colnames(scaledata[,-which(str_detect(colnames(scaledata), \"QC\"))]),rep(c(\"A\",\"B\"),each=3))\n# #     colnames(j)<-c(\"o1\",\"p1\",\"name\",\"group\")\n# #     ggplot(j, aes(p1, o1,group = group,col=group)) + #\u6620\u5c04PC1,PC2\n# #       geom_point(size = 7) + #\u7ed8\u5236\u5f97\u5206\u6563\u70b9\u56fe\n# #       labs(title = \"PCA - biplot\") + #\u56fe\u8868\u6807\u9898\n# #       xlab(paste(\"to1\")) + #x\u8f74\u6807\u9898\n# #       ylab(paste(\"t1\")) + #y\u8f74\u6807\u9898\n# #       stat_ellipse( geom = \"polygon\", alpha = 0.1, level = 0.95, show.legend = F,\n# #                     size=1)+#\u52a0\u7f6e\u4fe1\u692d\u5706\n# #       geom_vline(xintercept=mean(j$p1),lty=1,col=\"black\",lwd=0.5) +#\u6dfb\u52a0\u6a2a\u7ebf\n# #       geom_hline(yintercept =mean(j$o1),lty=1,col=\"black\",lwd=0.5)+#\u6dfb\u52a0\u7ad6\u7ebf\n# #       theme( plot.background = element_rect(fill=\"white\"),\n# #              panel.background = element_rect(fill='white', colour='gray'),          \n# #              strip.text.x=element_text(size=rel(1.2), family=\"serif\", angle=-90),\n# #              strip.text.y=element_text(size=rel(1.2),  family=\"serif\") ,\n# #              axis.text.x = element_text(size = 14,color=\"black\"),\n# #              axis.text.y = element_text(size = 20,color=\"black\"),axis.ticks.x=element_blank()\n# #       )\n# # }<\/code><\/pre>\n\n\n\n
\"\"<\/div><\/figure><\/div>\n\n\n\n
\"\"<\/div><\/figure><\/div>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n\n\n\n
\"\"<\/div><\/figure>\n","protected":false},"excerpt":{"rendered":"","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[1],"tags":[],"class_list":["post-310","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"_links":{"self":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/310","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=310"}],"version-history":[{"count":9,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/310\/revisions"}],"predecessor-version":[{"id":347,"href":"https:\/\/linguopeng.top\/index.php?rest_route=\/wp\/v2\/posts\/310\/revisions\/347"}],"wp:attachment":[{"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=310"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=310"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/linguopeng.top\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=310"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}