【统计】相关与回归：混合效应模型

• RNA-seq数据调整性别和年龄，及一个病例对应多个样本。

library(nlme)
datTraits<-read_csv("2-datTraits_95.csv", col_names = T)
head(names(datTraits))

# Fixed effects:
Diagnosis = as.factor(datTraits$Diagnosis)
Diagnosis = relevel(Diagnosis,ref = "CTL") # Using control as reference level
Age = as.numeric(datTraits$Age)
Sex = as.factor(datTraits$Sex)
Region = as.factor(datTraits$Region)
# Random effect:
BrainID = as.factor(datTraits$BrainID)
DE_Diagnosis = data.frame(miRID = NA,Beta = NA,StdErr = NA,Pval = NA)
DE_Age = data.frame(miRID = NA,Beta = NA,StdErr = NA,Pval = NA)
DE_Sex = data.frame(miRID = NA,Beta = NA,StdErr = NA,Pval = NA)
DE_Region = data.frame(miRID = NA,Beta = NA,StdErr = NA,Pval = NA)

# Define a function that runs the linear mixed-effects model and extract Beta values, standard errors, and P values for the fixed effects
runlme = function(miR) {
  lme = lme(miR~Diagnosis + Age + Sex + Region, rand = ~1|BrainID, na.action = na.exclude)
  slme = summary(lme);
  Beta_Diagnosis = slme$tTable[2,1]
  Beta_Age = slme$tTable[3,1]
  Beta_Sex = slme$tTable[4,1]
  Beta_Region = slme$tTable[5,1]
  StdErr_Diagnosis = slme$tTable[2,2]
  StdErr_Age = slme$tTable[3,2]
  StdErr_Sex = slme$tTable[4,2]
  StdErr_Region = slme$tTable[5,2]
  Pval_Diagnosis = slme$tTable[2,5]
  Pval_Age = slme$tTable[3,5]
  Pval_Sex = slme$tTable[4,5]
  Pval_Region = slme$tTable[5,5]
  LME_Diagnosis = list(Beta = Beta_Diagnosis,StdErr = StdErr_Diagnosis,Pval = Pval_Diagnosis)
  LME_Age = list(Beta = Beta_Age,StdErr = StdErr_Age,Pval = Pval_Age)
  LME_Sex = list(Beta = Beta_Sex,StdErr = StdErr_Sex,Pval = Pval_Sex)
  LME_Region = list(Beta = Beta_Region,StdErr = StdErr_Region,Pval = Pval_Region)
  return(list(LME_Diagnosis = LME_Diagnosis,LME_Age = LME_Age,LME_Sex = LME_Sex,LME_Region = LME_Region));
}

# 计算
experset<-read.csv('3-normdatExpr_DEseq_95_miRBAse_batch.csv',header = T, row.names = 1, sep = ",") 
names(experset)[c(1:10)]

for (i in 1:nrow(experset)){
  #i=1
  miR = as.numeric(experset[i,])
  result = try(runlme(miR),silent = F)
  if (length(result) == 4) {
    Beta_Diagnosis = result[[1]]$Beta
    StdErr_Diagnosis = result[[1]]$StdErr
    Pval_Diagnosis = result[[1]]$Pval
    Beta_Age = result[[2]]$Beta
    StdErr_Age = result[[2]]$StdErr
    Pval_Age = result[[2]]$Pval
    Beta_Sex = result[[3]]$Beta
    StdErr_Sex = result[[3]]$StdErr
    Pval_Sex = result[[3]]$Pval
    Beta_Region = result[[4]]$Beta
    StdErr_Region = result[[4]]$StdErr
    Pval_Region = result[[4]]$Pval
  } else {
    cat('Error in LME for microRNA',rownames(datExpr)[i],'\n')
    cat(' Setting Beta value=0, SE=Inf, and P-value=1\n')
    Beta_Diagnosis = Beta_Age = Beta_Sex = Beta_Region = 0
    StdErr_Diagnosis = StdErr_Age = StdErr_Sex = StdErr_Region = Inf
    Pval_Diagnosis = Pval_Age = Pval_Sex = Pval_Region = 1
  }
  DE_Diagnosis[i,"miRID"] = rownames(experset)[i]; DE_Diagnosis[i,"Beta"] = Beta_Diagnosis; DE_Diagnosis[i,"StdErr"] = StdErr_Diagnosis; DE_Diagnosis[i,"Pval"] = Pval_Diagnosis
  DE_Age[i,"miRID"] = rownames(experset)[i]; DE_Age[i,"Beta"] = Beta_Age; DE_Age[i,"StdErr"] = StdErr_Age; DE_Age[i,"Pval"] = Pval_Age
  DE_Sex[i,"miRID"] = rownames(experset)[i]; DE_Sex[i,"Beta"] = Beta_Sex; DE_Sex[i,"StdErr"] = StdErr_Sex; DE_Sex[i,"Pval"] = Pval_Sex
  DE_Region[i,"miRID"] = rownames(experset)[i]; DE_Region[i,"Beta"] = Beta_Region; DE_Region[i,"StdErr"] = StdErr_Region; DE_Region[i,"Pval"] = Pval_Region
}

# Perform FDR adjustment using Benjamini-Hochberg correction
DE_Diagnosis$adjPval = p.adjust(DE_Diagnosis$Pval,method = "BH",n = nrow(DE_Diagnosis))
DE_Age$adjPval = p.adjust(DE_Age$Pval,method = "BH",n = nrow(DE_Age))
DE_Sex$adjPval = p.adjust(DE_Sex$Pval,method = "BH",n = nrow(DE_Sex))
DE_Region$adjPval = p.adjust(DE_Region$Pval,method = "BH",n = nrow(DE_Region))

DE_Diagnosis<-DE_Diagnosis[order(DE_Diagnosis$Pval),] #排序
DE_Age<-DE_Diagnosis[order(DE_Age$Pval),]
DE_Sex<-DE_Diagnosis[order(DE_Sex$Pval),]
DE_Region<-DE_Diagnosis[order(DE_Region$Pval),]

write.csv(DE_Diagnosis,file = "./4-batch_DGE_Diagnosis.csv",row.names = F) 
write.csv(DE_Age,file = "./4-batch_DGE_Age_batch.csv",row.names = F)
write.csv(DE_Sex,file = "./4-batch_DGE_Sex.csv",row.names = F)
write.csv(DE_Region,file = "./4-batch_DGE_BrainRegion.csv",row.names = F)