首先准备数据和分组文件
形式如下：
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
