微生物联合挖掘教程连载-PC-6-差异分析多种展示方式

diffpath = paste(otupath,"/diff_Manhattan/",sep = "")dir.create(diffpath)
  ps = ps  pvalue = 0.05  lfc = 0  diffpath = diffpath 
# otu化count = ps %>%     ggClusterNet::vegan_otu() %>%    t()# create DGE list# 制作DGE序列# DGEList是一个可以包含多种内容和统计的列表。DGEList至少需要的元素：counts、samples（包含group分组信息和lib.size文库大小），counts用来存放表达矩阵，samples用来标记样本信息和库的大小，group声明组别  d = edgeR::DGEList(counts=count, group= phyloseq::sample_data(ps)$Group)  d = edgeR::calcNormFactors(d)  design.mat = model.matrix(~ 0 + d$samples$group)  colnames(design.mat)=levels(as.factor(phyloseq::sample_data(ps)$Group))  d2 = edgeR::estimateGLMCommonDisp(d, design.mat)  d2 = edgeR::estimateGLMTagwiseDisp(d2, design.mat)  fit = edgeR::glmFit(d2, design.mat)# fit里有系数、拟合值等数值  
  Desep_group <- as.character(levels(as.factor(phyloseq::sample_data(ps)$Group)))  Desep_group  # 显示分类的组名
  aaa = combn(Desep_group,2)# 两两间的组合# 开始做分析  for (i in 1:dim(aaa)[2]) {    # i = 1    Desep_group = aaa[,i]    print( Desep_group)    group = paste(Desep_group[1],Desep_group[2],sep = "-")    group    BvsA <- limma::makeContrasts(contrasts =  group,levels=design.mat)#注意是以GF1为对照做的比较    # 组间比较,统计Fold change, Pvalue    lrt = edgeR::glmLRT(fit,contrast=BvsA)        # FDR检验，控制假阳性率小于5%    de_lrt = edgeR::decideTestsDGE(lrt, adjust.method="fdr", p.value=pvalue,lfc=lfc)#lfc=0这个是默认值    summary(de_lrt)    # 导出计算结果    x=lrt$table    x$sig=de_lrt        x$sig=de_lrt    head(x)    #------差异结果符合otu表格的顺序    row.names(count)[1:6]        x <- cbind(x, padj = p.adjust(x$PValue, method = "fdr"))    enriched = row.names(subset(x,sig==1))    depleted = row.names(subset(x,sig==-1))    x$level = as.factor(ifelse(as.vector(x$sig) ==1, "enriched",ifelse(as.vector(x$sig)==-1, "depleted","nosig")))    x$otu = rownames(x)    x$neglogp = -log(x$PValue)
    tax = ps %>%       ggClusterNet::vegan_tax() %>%      as.data.frame()    head(tax)    x = cbind(x,tax)    head(x)        top_phylum=c("Bacteroidetes","Firmicutes","Planctomycetes","Proteobacteria","Verrucomicrobia")    # 在这里改想要呈现的门    x[!(x$Phylum %in% top_phylum),]$Phylum = "Low Abundance" # no level can get value    x$otu = factor(x$otu, levels=x$otu)   # set x order    x$level = factor(x$level, levels=c("enriched","depleted","nosig"))    levels(x$Phylum)=c(top_phylum,"Low Abundance")    x = x %>% arrange(Phylum)    head(x)    x$otu = factor(x$otu, levels=x$otu)           # if (tem != 0) {    #   tem = x[x$neglogp>15,] %>% nrow()    # }        FDR = min(x$neglogp[x$level=="depleted"])    p = ggplot(x, aes(x=otu, y=neglogp, color=Phylum, size=logCPM, shape=level)) +      geom_point(alpha=.7) +       geom_hline(yintercept=FDR, linetype=2, color="lightgrey") +      scale_shape_manual(values=c(17, 25, 20))+      scale_size(breaks=c(5, 10, 15)) +      labs(x="OTU", y="-loge(P)") +      theme(axis.ticks.x=element_blank(),axis.text.x=element_blank(),legend.position="top") +      scale_color_manual(values = c(RColorBrewer::brewer.pal(9,"Set1")))    p        filename = paste(diffpath,"/",paste(Desep_group[1],Desep_group[2],sep = "_"),"Manhattan_plot.pdf",sep = "")    ggsave(filename,p,width = 16,height = 6)    filename = paste(diffpath,"/",paste(Desep_group[1],Desep_group[2],sep = "_"),"Manhattan_plot.png",sep = "")    ggsave(filename,p,width = 16,height = 6,dpi = 72)      }  

ps = readRDS("./data/dataNEW/ps_16s.rds")#==这里有点小问题，就是我的电脑没有ps_16s.rds文件，我用的是另一个rds文件，可能出图有点区别==heatpath = paste(otupath,"/heapmap_boplot/",sep = "")dir.create(heatpath)
map = phyloseq::sample_data(ps)map$ID = row.names(map)phyloseq::sample_data(ps) = mapj = 2 #做门ps_tem = ps %>%     ggClusterNet::scale_micro(method = "TMM") %>%    ggClusterNet::tax_glom_wt(ranks = j) 
rowSD = function(x){    apply(x,1, sd)  }  rowCV = function(x){    rowSD(x)/rowMeans(x)  }  id <- ps %>%     ggClusterNet::scale_micro(method = "TMM") %>%    ggClusterNet::tax_glom_wt(ranks = j) %>%    ggClusterNet::filter_OTU_ps(100) %>%    ggClusterNet::vegan_otu() %>%    t() %>% as.data.frame() %>%rowCV %>%    sort(decreasing = TRUE) %>%    head(20) %>%    names()# 得到门类的名字  
ps_rela= ps_temheatnum = 20label=TRUEcol_cluster=TRUErow_cluster=TRUE

  map = phyloseq::sample_data(ps_rela)  map$ID = row.names(map)  phyloseq::sample_data(ps_rela) = map  otu = as.data.frame(t(ggClusterNet::vegan_otu(ps_rela)))  otu = as.matrix(otu[id,])  ps_heatm = phyloseq::phyloseq(    phyloseq::otu_table(otu,taxa_are_rows = TRUE),    phyloseq::tax_table(ps_rela),    phyloseq::sample_data(ps_rela)      )    # print(ps_heatm)  # 得到otu表    datah <- as.data.frame(t(ggClusterNet::vegan_otu(ps_heatm)))  head(datah)   tax = as.data.frame(ggClusterNet::vegan_tax(ps_heatm))    otutaxh = cbind(datah,tax)  head(otutaxh)    otutaxh$id = paste(row.names(otutaxh),otutaxh$Genus,sep = "_")  # otutaxh$id  = row.names(otutaxh)  row.names(otutaxh) = otutaxh$id   # 转化成dataframe    data <- otutaxh[,c("id",phyloseq::sample_names(ps))]
  rig <- data[,phyloseq::sample_names(ps)] %>% rowMeans() %>% as.data.frame()  head(rig)  # 算出均值并排序  colnames(rig) = "MeanAbundance"  rig$id = row.names(rig)  rig = rig %>% dplyr::arrange(MeanAbundance)  rig$id = factor(rig$id,levels = rig$id)  p_rig = ggplot(rig) + geom_bar(aes(y = id,x = MeanAbundance),                                 fill = "#A54657",                                 stat = "identity") + theme_void()    tem = data[,phyloseq::sample_names(ps)] %>% as.matrix()     tem = scale(t(tem)) %>% t() %>%      as.data.frame()   data[,phyloseq::sample_names(ps)] = tem     # 计算距离矩阵，出分类树     # data[data > 0.3]<-0.3  mat <- data[,-1] #drop gene column as now in rows    if (col_cluster ==  TRUE) {    clust <- hclust(dist(mat %>% as.matrix())) # hclust with distance matrix    ggtree_plot <- ggtree::ggtree(clust)  }  if (row_cluster ==  TRUE) {    v_clust <- hclust(dist(mat %>% as.matrix() %>% t()))    ggtree_plot_col <- ggtree::ggtree(v_clust) + ggtree::layout_dendrogram()  }    if (label ==  TRUE) {    map = as.data.frame(phyloseq::sample_data(ps))    map$ID = row.names(map)    labels= ggplot(map, aes(x = ID, y=1, fill=Group)) + geom_tile() +      scale_fill_brewer(palette = 'Set1',name="Cell Type") +      theme_void()  }    map = phyloseq::sample_data(ps) %>% as.tibble() %>%    dplyr::arrange(Group) %>% as.data.frame()  map$ID  pcm = reshape2::melt(data, id = c("id"))  pcm$variable = factor(pcm$variable,levels = map$ID)  pcm$id = factor(pcm$id,levels = rig$id)      p1 = ggplot(pcm, aes(y = id, x = variable)) +     # geom_point(aes(size = value,fill = value), alpha = 0.75, shape = 21) +     geom_tile(aes(size = value,fill = value))+    scale_size_continuous(limits = c(0.000001, 100), range = c(2,25), breaks = c(0.1,0.5,1)) +     labs( y= "", x = "", size = "Relative Abundance (%)", fill = "")  +     # scale_fill_manual(values = colours, guide = FALSE) +     scale_x_discrete(limits = rev(levels(pcm$variable)))  +     scale_y_discrete(position = "right") +    scale_fill_gradientn(colours =colorRampPalette(RColorBrewer::brewer.pal(11,"Spectral")[11:1])(60))+    theme(      panel.background=element_blank(),      panel.grid=element_blank(),      axis.text.x = element_text(colour = "black",angle = 90)          )    colours = c( "#A54657",  "#582630", "#F7EE7F", "#4DAA57","#F1A66A","#F26157", "#F9ECCC", "#679289", "#33658A",               "#F6AE2D","#86BBD8")  #----样本在y轴上  p2 = ggplot(pcm, aes(y = id, x = variable)) +     geom_point(aes(size = value,fill = value), alpha = 0.75, shape = 21) +     scale_size_continuous(limits = c(0.000001, 100), range = c(2,25), breaks = c(0.1,0.5,1)) +     labs( y= "", x = "", size = "Relative Abundance (%)", fill = "")  +     # scale_fill_manual(values = colours, guide = FALSE) +     scale_x_discrete(limits = rev(levels(pcm$variable)))  +     scale_y_discrete(position = "right")  +    scale_fill_gradientn(colours =colorRampPalette(RColorBrewer::brewer.pal(11,"Spectral")[11:1])(60)) +     theme(      panel.background=element_blank(),      panel.grid=element_blank(),      axis.text.x = element_text(colour = "black",angle = 90)            )
  p1 <- p1  %>%    aplot::insert_right(p_rig, width=.2)   p2 <- p2  %>%    aplot::insert_right(p_rig, width=.2)   if (col_cluster ==  TRUE) {    p1 <- p1  %>%      aplot::insert_left(ggtree_plot, width=.2)    p2 <- p2  %>%      aplot::insert_left(ggtree_plot, width=.2)  }
  if (row_cluster ==  TRUE) {    p1 <- p1  %>%      aplot::insert_top(labels, height=.02)     p2 <- p2  %>%      aplot::insert_top(labels, height=.02)   }    if (label ==  TRUE) {    p1 <- p1  %>%      aplot::insert_top(ggtree_plot_col, height=.1)    p2 <- p2  %>%      aplot::insert_top(ggtree_plot_col, height=.1)  }

  filename = paste(heatpath,"/",phyloseq::rank.names(ps)[j],"Topggheatmap.pdf",sep = "")  ggsave(filename,p1,width = 14,height = (6 + heatnum/10))    filename = paste(heatpath,phyloseq::rank.names(ps)[j],"Topggbubble.pdf",sep = "")  ggsave(filename,p2,width = 14,height = (6 + heatnum/10))    filename = paste(heatpath,"/",phyloseq::rank.names(ps)[j],"Topggheatmap.png",sep = "")  ggsave(filename,p1,width = 14,height = (6 + heatnum/10))    filename = paste(heatpath,phyloseq::rank.names(ps)[j],"Topggbubble.png",sep = "")  ggsave(filename,p2,width = 14,height = (6 + heatnum/10))

# BiocManager::install("MetBrewer")

diffpath.1 = paste(otupath,"/Mui.Group.v/",sep = "")dir.create(diffpath.1)
source("./function/EdgerSuper2.R")res = EdgerSuper2 (ps = ps,group  = "Group",artGroup =NULL,                   j = "OTU",                   path = diffpath.1)
head(res)# res里有logFC值，level(enrich/nosing...),p值等相关数据
  res$ID = row.names(res)  datv = res   # for循环挑选每个cluster的top前5 gene symbol  tm.g <- function(data){    id = data$group %>% unique()        for (i in 1:length(id)) {      tem = filter(data,group==id[i],level != "nosig") %>%         distinct(ID,.keep_all = TRUE) %>%         top_n(5,abs(logFC))#想改在这里改      if (i == 1) {        tem2 = tem      } else {        tem2 = rbind(tem2,tem)      }    }    return(tem2)  }    top <- tm.g(datv)  # 先画背景柱，根据数据log2FC的max值,min值来确定  #根据数据中log2FC区间确定背景柱长度：  # 长度越长，说明差异越大  head(datv)    tem = datv %>% group_by(group) %>% summarise(max = max(logFC),min = min(logFC)) %>% as.data.frame()    col1<-data.frame(x=tem$group,                   y=tem$max)  col2<-data.frame(x=tem$group,                   y=tem$min)  # 绘制背景柱  p1 <- ggplot()+    geom_col(data = col1,             mapping = aes(x = x,y = y),             fill = "#dcdcdc",alpha = 0.6)+    geom_col(data = col2,             mapping = aes(x = x,y = y),             fill = "#dcdcdc",alpha = 0.6)  p1        #把散点火山图叠加到背景柱上：  head(datv)    p2 <- ggplot()+    geom_col(data = col1,             mapping = aes(x = x,y = y),             fill = "#dcdcdc",alpha = 0.6)+    geom_col(data = col2,             mapping = aes(x = x,y = y),             fill = "#dcdcdc",alpha = 0.6)+    geom_jitter(data = datv,                aes(x =group , y = logFC, color =level ),                size = 1,                width =0.4)+    scale_color_manual(name=NULL,                       values = c("#4393C3","#FC4E2A","grey40"))+    labs(x="",y="log2(FoldChange)")    #改颜色，第一个是deplete,第二个是enrich，第三个是nosing  p2    # 添加X轴的分组色块标签：  dfcol<-data.frame(x=tem$group,                    y=0,                    label=tem$group)  # 添加分组色块标签  dfcol$group <- tem$group  # 加载包  library(RColorBrewer)  library(MetBrewer)  # BiocManager::install("MetBrewer")  # 自定义分组色块的颜色  tile_color <- met.brewer("Thomas",length(tem$group))    # 在图中镶嵌色块  p3 <- p2 + geom_tile(data = dfcol,                       aes(x=x,y=y),                       height=1.75,                       color = "black",                       fill = tile_color,                       alpha = 0.6,                       show.legend = F)+    geom_text(data=dfcol,              aes(x=x,y=y,label=group),              size =3.5,              color ="white") + theme_classic()  p3    library(ggrepel)  p4<-p3+geom_text_repel(    data=top,    aes(x=group,y=logFC,label=ID),    force = 1.2,    arrow = arrow(length = unit(0.008, "npc"),                  type = "open", ends = "last"))  p4  # 去除背景，美化图片  p5 <- p4+    theme_minimal()+    theme(      axis.title = element_text(size = 18,                                color = "black",                                face = "bold"),      axis.line.y = element_line(color = "black",                                 size = 1.2),      axis.line.x = element_blank(),      axis.text.x = element_blank(),      panel.grid = element_blank(),      legend.position = "top",      legend.direction = "vertical",      legend.justification = c(1,0),      legend.text = element_text(size = 12)    )  p5
# return(list(p5,p3,datv,top))p = p3pfilename = paste(diffpath.1,"/","Mui.group.volcano.pdf",sep = "")ggsave(filename,p,width = 12,height = 6,limitsize = FALSE)

> source("./function/EdgerSuper2.R")
Error in file(filename, "r", encoding = encoding) :cannot open the connectionIn addition: Warning message:In file(filename, "r", encoding = encoding) :cannot open file './function/EdgerSuper2.R': No such file or directory