R语言进行最小显著性差异分析 （LSD）

首先准备数据和分组文件

形式如下：

image.png
image.png

#读取文件并处理
library(agricolae)
library("plyr")
feature <- read.csv("npgene.csv",row.names = 1)
group <- read.csv("group.CSV",row.names = 1)
a <- as.data.frame(t(feature))
#将分组合并
feature_group <- cbind(a,group[match(row.names(a),row.names(group)),])
#编写LSD两组分析函数
LSD_two <- function(data=data1,group = "site"){
#判断数据框内是否有和为零的列，有则删除（避免LSD（）出错）
  number1 = 1
  for (v in 1:(ncol(aaa)-2)) {
    if (sum(data[v] == 0)) {
      number1 <- c(number1,v)
    }else{
      abcd <- 0
    }
  }
  if (length(number1) == 1  ) {
    abcd <- 0
  }else{
    number1 <- number1[-1]
    data <- data[-number1] 
  }
#判断按照深度还是位置进行LSD
  if (group == "depth") {
    for (i in 1:(ncol(data)-2)) {
      aa <- data.frame()
      temps <- list()
      aa <- data[,c(i,ncol(data))]
      temps <- aov(data[,i] ~ groupd, aa)
      s[[i]] <- summary(temps)
      listgroups[[i]] <- LSD.test(temps, 'groupd', p.adj = 'none')
    }
  }
  if (group == "site") {
    for (i in 1:(ncol(data)-2)) {
      aa <- data.frame()
      temps <- list()
      aa <- data[,c(i,(ncol(data)-1))]
      temps <- aov(data[,i] ~ groups, aa)
      s[[i]] <- summary(temps)
      listgroups[[i]] <- LSD.test(temps, 'groups', p.adj = 'none')
    }
  }
#将没列LSD数据整合成数据框
  end <- listgroups[[1]]$groups
  zzz <- listgroups[[1]]$groups
  for (v in 1:(ncol(data)-2)) {
    vv <- listgroups[[v]]$groups
    vv <- vv[2]
    end <- cbind(end,vv[match(row.names(end),row.names(vv)),])
    end1 <- end
  }
  colnames(end1)[3:length(colnames(end))] <- colnames(data)[1:(ncol(data)-2)]
  end1 <- end1[,-c(1,2)]
  AEND <- end1
#通过差异筛选数据，只显示有差异的数据
  for (g in 1:length(colnames(AEND))) {
    if (length(levels(as.factor(AEND[,g]))) > 1) {
      zzf <- AEND[g]
      zzz <- cbind(zzz,zzf)
    }
  }
  zzz <- zzz[-c(1,2)]
  zzz$site <- row.names(zzz)
  diyi <- zzz[ncol(zzz)]
  dier <- zzz[-ncol(zzz)]
  zzz <- cbind(diyi,dier)
  return(zzz)
}
#运用
LSD_two(feature_group,"site")

输出文件如下

image.png

其他相关代码（对各个位置不同深度或各深度不同位置进行LSD）

#数据框同上，LSD_two()也同上，group = "depth" 代表计算各个位置不同深度的LSD；"site"代表各个深度不同位置的LSD。
LSD_onestep <- function(data=data1,group = "site"){
  LSD_two <- function(data=data1,group = "site"){
    number1 = 1
    for (v in 1:(ncol(data)-2)) {
      if (sum(data[v] == 0)) {
        number1 <- c(number1,v)
      }else{
        abcd <- 0
      }
    }
    if (length(number1) == 1  ) {
      abcd <- 0
    }else{
      number1 <- number1[-1]
      data <- data[-number1] 
    }
    if (group == "depth") {
      for (i in 1:(ncol(data)-2)) {
        aa <- data.frame()
        temps <- list()
        aa <- data[,c(i,ncol(data))]
        temps <- aov(data[,i] ~ groupd, aa)
        s[[i]] <- summary(temps)
        listgroups[[i]] <- LSD.test(temps, 'groupd', p.adj = 'none')
      }
    }
    if (group == "site") {
      for (i in 1:(ncol(data)-2)) {
        aa <- data.frame()
        temps <- list()
        aa <- data[,c(i,(ncol(data)-1))]
        temps <- aov(data[,i] ~ groups, aa)
        s[[i]] <- summary(temps)
        listgroups[[i]] <- LSD.test(temps, 'groups', p.adj = 'none')
      }
    }
    end <- listgroups[[1]]$groups
    zzz <- listgroups[[1]]$groups
    for (v in 1:(ncol(data)-2)) {
      vv <- listgroups[[v]]$groups
      vv <- vv[2]
      end <- cbind(end,vv[match(row.names(end),row.names(vv)),])
      end1 <- end
    }
    colnames(end1)[3:length(colnames(end))] <- colnames(data)[1:(ncol(data)-2)]
    end1 <- end1[,-c(1,2)]
    AEND <- end1
    for (g in 1:length(colnames(AEND))) {
      if (length(levels(as.factor(AEND[,g]))) > 1) {
        zzf <- AEND[g]
        zzz <- cbind(zzz,zzf)
      }
    }
    zzz <- zzz[-c(1,2)]
    zzz$site <- row.names(zzz)
    diyi <- zzz[ncol(zzz)]
    dier <- zzz[-ncol(zzz)]
    zzz <- cbind(diyi,dier)
    return(zzz)
  }
  #如果是深度，则按照位置切分数据框并计算LSD
  if(group=="depth"){
    list22 <- list()
    for ( i in levels(data[,"groups"])) {
      site11 <- LSD_two(data[data["groups"] == i,],group ="depth")
      colnames(site11)[1] <- i
      a <- as.data.frame(t(site11))
      a$ababa <- rownames(a)
      diyi <- a[ncol(a)]
      dier <- a[-ncol(a)]
      a <- cbind(diyi,dier)
      row.names(a)[1] <- "site"
      list22[[i]] <- as.data.frame(t(a))
    }
#利用rbind.fill()函数将不同列数的数据框按列合并，将NA值转换为""空值
    rbind.zzf <- function(data=a,data1=b){
      zzong <- rbind.fill(data,data1)
      kkk <- as.character(zzong[,1])
      for (i in 2:ncol(zzong)){
        zzf <- as.character(zzong[,i])
        zzf[is.na(zzf)] <- ""
        kkk <- cbind(kkk,zzf)
      }
      colnames(kkk) <- colnames(zzong)
      kkk <- as.data.frame(kkk)
      return(kkk)
    }
#通过递归方式将切分并计算的LSD合并
    kkka <- data.frame()
    for (i in 1:length(list22)) {
      kkka <- rbind.zzf(kkka, list22[[i]])
    }
    write.table(kkka,"LSDdepth.csv",sep = ",",row.names = F,col.names = F,quote=F)
  }else if(group == "site"){
    list22 <- list()
    for ( i in levels(data[,"groupd"])) {
      site11 <- LSD_two(data[data["groupd"] == i,],group ="site")
      colnames(site11)[1] <- i
      a <- as.data.frame(t(site11))
      a$ababa <- rownames(a)
      diyi <- a[ncol(a)]
      dier <- a[-ncol(a)]
      a <- cbind(diyi,dier)
      row.names(a)[1] <- "site"
      list22[[i]] <- as.data.frame(t(a))
    }
    rbind.zzf <- function(data=a,data1=b){
      zzong <- rbind.fill(data,data1)
      kkk <- as.character(zzong[,1])
      for (i in 2:ncol(zzong)){
        zzf <- as.character(zzong[,i])
        zzf[is.na(zzf)] <- ""
        kkk <- cbind(kkk,zzf)
      }
      colnames(kkk) <- colnames(zzong)
      kkk <- as.data.frame(kkk)
      return(kkk)
    }
    kkka <- data.frame()
    for (i in 1:length(list22)) {
      kkka <- rbind.zzf(kkka, list22[[i]])
    }
    write.table(kkka,"LSDsite.csv",sep = ",",row.names = F,col.names = F,quote=F)
  }
}

输出文件如下

image.png