第二课题原始参考代码

机器学习部分
#Description
#The core R codes for statistical analyses of manuscript gut jnl-2019-318860
#If you have any question of our source code, please contact me at yuchen_li@foxmail.com
#>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#The tsne method for Dimension reduction and visualization 
library(Rtsne)
tsne <- Rtsne(matrix, dims=3, pca = F, initial_dims = 10, perplexity = 100, max_iter = 500)

#the lasso and random forest algorithm for exLR-seq marker selection 
library(glmnet)
lasso.results <- c()
for (j in 1:1000) {
  n1 <- n[sample(1:nrow(n),ceiling(dim(n)[1]*0.5),replace=F),]
  h1 <- h[sample(1:nrow(h),ceiling(dim(h)[1]*0.5),replace=F),]
  s <- rbind.data.frame(n1,h1)
  ts <- apply(s,2,as.numeric)
  y <- as.matrix(ts[,1])
  x <- as.matrix(ts[,2:dim(ts)[1]])
  cv.fit <- cv.glmnet(x,y,family="binomial", type.measure = "auc", nfolds=5)
  co<-coef.cv.glmnet(cv.fit,s="lambda.1se")
  name <- rownames(co)[co[,1]!=0]
  lasso.results <- c(lasso.results, name)
}
freq.lasso.results <- table(lasso.results)

library(varSelRF)
facy <- factor(training.data$status)
x <- training.data[,-1]
step=varSelRF(x, facy, c.sd = 1, mtryFactor = 1, ntree = 500,
              ntreeIterat = 500, vars.drop.num = NULL, vars.drop.frac = 0.1, #
              whole.range = TRUE, recompute.var.imp = FALSE, verbose = FALSE,
              returnFirstForest = TRUE, fitted.rf = NULL, keep.forest = FALSE)
select.history <- step$selec.history


#The training diganostic model construction and SVM algorithm operation
set.seed(1)
library(pROC)
library(caret)
library(e1071)
library(kernlab)
fitControl <- trainControl(method = "repeatedcv",
                           number = 5,
                           repeats = 10,
                           ## Estimate class probabilities
                           classProbs = TRUE,
                           ## Evaluate performance using 
                           ## the following function
                           summaryFunction = twoClassSummary)

svmFit <- train(status ~ ., data = training.data, 
                method = "svmRadial", 
                trControl = fitControl, 
                preProc = c("center", "scale"),
                tuneLength = 8,
                metric = "ROC")

#The establishment of d-signature from three cohort
train.pred <- predict(svmFit,training.data)
validation.pred <- predict(svmFit,validation.data)
testing.pred <- predict(svmFit,testing.data)

#The caculation of diagnostic efficacy of d-signature
train.con <- confusionMatrix(train.pred, training.data$status)
validation.con <- confusionMatrix(validation.pred, validation.data$status)
testing.con <- confusionMatrix(testing.pred, testing.data$status)
tr.acc <- train.con$overall[1]
tr.se <- train.con$byClass[1]
tr.sp <- train.con$byClass[2]
va.acc <- validation.con$overall[1]
va.se <- validation.con$byClass[1]
va.sp <- validation.con$byClass[2]
te.acc <- testing.con$overall[1]
te.se <- testing.con$byClass[1]
te.sp <- testing.con$byClass[2]

#The Normfinder method for RT-qPCR
source("r.NormOldStab5.txt")
Result=Normfinder("Datafile.txt",ctVal=FALSE)
selected.exLRs.normfinder.results <- Result$Ordered

性别和年龄调整部分

####### Differential gene expression analysis using a linear mixed-effects model

load("./datExpr.rda") # Load expression data (normalized for library size and technical variables, log2 transformed)
load("./datTraits.rda") # Load metadata

library(nlme)

# 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$Brain_ID)

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));
}

# Run the linear mixed-effects model for all miRNAs in the data set
for (i in 1:nrow(datExpr)){
  miR = as.numeric(datExpr[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(datExpr)[i]; DE_Diagnosis[i,"Beta"] = Beta_Diagnosis; DE_Diagnosis[i,"StdErr"] = StdErr_Diagnosis; DE_Diagnosis[i,"Pval"] = Pval_Diagnosis
  DE_Age[i,"miRID"] = rownames(datExpr)[i]; DE_Age[i,"Beta"] = Beta_Age; DE_Age[i,"StdErr"] = StdErr_Age; DE_Age[i,"Pval"] = Pval_Age
  DE_Sex[i,"miRID"] = rownames(datExpr)[i]; DE_Sex[i,"Beta"] = Beta_Sex; DE_Sex[i,"StdErr"] = StdErr_Sex; DE_Sex[i,"Pval"] = Pval_Sex
  DE_Region[i,"miRID"] = rownames(datExpr)[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))

write.csv(DE_Diagnosis,file = "./DGE_Diagnosis.csv")
write.csv(DE_Age,file = "./DGE_Age.csv")
write.csv(DE_Sex,file = "./DGE_Sex.csv")
write.csv(DE_Region,file = "./DGE_BrainRegion.csv")

WGCNA部分（未使用）

####### WGCNA analysis using bootstrapping

library(WGCNA)
library(flashClust)
options(stringsAsFactors  =  FALSE)

dir.create("./Network_resampled")
dir.create("./TOM_resampled")


### 1. Construct network using the original sample set
load("./datExpr.rda") # Load expression data (normalized for library size and technical variables, log2 transformed)
load("./datTraits.rda") # Load metadata

# Pick soft threshold
powers  =  c(c(1:10), seq(from = 12, to = 20, by = 2))
sft = pickSoftThreshold(datExpr,networkType = "signed",corFnc = "bicor",verbose = 5,powerVector = powers)

softPower = 8 # Choose based on fit to scale-free topology
adjacency = adjacency(datExpr, corFnc = "bicor", type = "signed", power = softPower)
TOM = TOMsimilarity(adjacency,TOMType = "signed", verbose = 0)
dissTOM = 1-TOM
geneTree = flashClust(as.dist(dissTOM), method = "average")

minModuleSize = 10
ds = 3
cutHeight = 0.99999
dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, method = "hybrid",
                            deepSplit = ds, pamRespectsDendro = T,pamStage = T,
                            minClusterSize = minModuleSize, cutHeight = cutHeight)
dynamicColors = labels2colors(dynamicMods)

# Calculate eigengenes
MEList = moduleEigengenes(datExpr, colors = dynamicColors,softPower = softPower)
MEs = MEList$eigengenes
# Merge modules whose module eigengenes are highly correlated
MEDissThres = 0.15
merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)
moduleColors = merge$colors
MEs = merge$newMEs

save(MEs, TOM, moduleColors, geneTree, file = "./Network.RData")

### 2. Perform bootstrapping 
datTraits_ASD = datTraits[datTraits$Diagnosis == "ASD",]
datTraits_Ctrl = datTraits[datTraits$Diagnosis == "CTL",]
datExpr_ASD = datExpr[match(rownames(datTraits_ASD),rownames(datExpr)),]
datExpr_Ctrl = datExpr[match(rownames(datTraits_Ctrl),rownames(datExpr)),]

softPower = 8 
ds = 3
cutHeight = 0.99999
minModuleSize = 10

nSets = 200

for (i in 1:nSets){
  if (i%%10 == 0) { print(paste(i,"th resampling",sep = "")) }
  # Resample within ASD and CTL groups
  ASD_resample = sample(1:nrow(datExpr_ASD),size = nrow(datExpr_ASD),replace = T) 
  datExpr_ASD_resample = datExpr_ASD[ASD_resample,]
  Ctrl_resample = sample(1:nrow(datExpr_Ctrl),size = nrow(datExpr_Ctrl),replace = T)
  datExpr_Ctrl_resample = datExpr_Ctrl[Ctrl_resample,]
  datExpr_resample = rbind(datExpr_ASD_resample,datExpr_Ctrl_resample)
  
  adjacency = adjacency(datExpr_resample, corFnc = "bicor", type = "signed", power = softPower)
  TOM = TOMsimilarity(adjacency,TOMType = "signed")
  dissTOM = 1-TOM
  geneTree = flashClust(as.dist(dissTOM), method = "average")
  dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, method="hybrid",
                              deepSplit = ds, pamRespectsDendro = T, pamStage=T,
                              minClusterSize = minModuleSize, cutHeight = cutHeight, verbose=0)
  dynamicColors = labels2colors(dynamicMods)
  
  MEList = moduleEigengenes(datExpr_resample, colors = dynamicColors,softPower = softPower)
  MEs = MEList$eigengenes
  
  # Merge modules whose module eigengenes are highly correlated
  MEDissThres = 0.15
  merge = mergeCloseModules(datExpr_resample, dynamicColors, cutHeight = MEDissThres, verbose = 0)
  moduleColors = merge$colors
  MEs = merge$newMEs
  
  save(MEs, TOM, moduleColors, geneTree, file = paste("./Network_resampled/Network_resampled",i,".RData",sep = ""))
}

### 3. Calculate consensus TO matrix 
nSets = 200
prob = 0.95 
nSubset = 20000

load("./Network.RData")
TOM.mat = as.matrix(TOM)
nGenes = nrow(TOM.mat)

subsetEntries = sample(nGenes*(nGenes-1)/2, size = nSubset)
TOMsubset = list()

TOMsubset.main = vectorizeMatrix(TOM.mat)[subsetEntries] 
quantile.TOM.main = quantile(TOMsubset.main,probs = prob,type = 8) 

quantile.TOM = rep(1, nSets) ## vector of quantiles of the individual TO matrices
beta.prob = rep(1, nSets) ## Scaling powers to calibrate TOM values                                        

for (set in 1:nSets) { 
  if (!file.exists(paste("./TOM_resampled/TOM_resampled",set,".RData",sep = ""))) {
    print(paste("On resampled TOM ",set,"...",sep = ""))
    load(paste("./Network_resampled/Network_resampled",set,".RData",sep = ""))
    tmpTOM = as.matrix(TOM)
    TOMsubset[[set]] = vectorizeMatrix(tmpTOM)[subsetEntries]  
    quantile.TOM[set] = quantile(TOMsubset[[set]],probs = prob,type = 8) 
    beta.prob[set] = log(quantile.TOM.main)/log(quantile.TOM[set]) ## calculate the power of the adjacency function
    print(paste("Scaling power for this TOM is ",signif(beta.prob[set],3),sep = ""))
    TOM = tmpTOM^(beta.prob[set]) ## use the power adjacency function for calibrating the TOM
    save(TOM,file = paste("./TOM_resampled/TOM_resampled",set,".RData",sep = ""))
  } else {
    print(paste("./Network_resampled/Network_resampled",set,".RData",sep = ""))
    print("Already processed!")
  }
}

rm(TOM.mat,tmpTOM)

consTOM = vector("list",nSets)
for (set in 1:nSets) {
  print(paste("loading TOM",set,sep = " "))
  load(paste("./TOM_resampled/TOM_resampled",set,".RData",sep = ""))
  consTOM[[set]]$TOM = TOM
}
print("Calculting the median quantile across the TOMs")
consensusTOM = matrix(NA,nrow = nrow(TOM),ncol = ncol(TOM))
for (i in 1:nrow(TOM) ) {
  if (i%%100 == 0) { print(paste("On gene",i,sep = " ")) }          
  for (j in 1:ncol(TOM)) {
    thisvec = vector(mode = "numeric",length = nSets)
    for (k in 1:nSets){
      thisvec[k]=consTOM[[k]]$TOM[i,j]
    }
    consensusTOM[i,j] = median(thisvec)
  }
}
save(consensusTOM,file = "./TOM_resampled/ConsensusTOM.RData")

# Define modules based on consensus TO matrix
cons_dissTOM = 1-consensusTOM
cons_geneTree = flashClust(as.dist(cons_dissTOM), method = "average")

softPower = 8
ds = 3
cutHeight = 0.99999
minModuleSize = 10
cons_dynamicMods = cutreeDynamic(dendro = cons_geneTree, distM = cons_dissTOM, method = "hybrid",
                                 deepSplit = ds, pamRespectsDendro = T, pamStage = T,
                                 minClusterSize = minModuleSize, cutHeight = cutHeight, verbose = 0)
cons_dynamicColors = labels2colors(cons_dynamicMods)

MEDissThres = 0.15
cons_merge = mergeCloseModules(datExpr, cons_dynamicColors, cutHeight = MEDissThres, verbose = 0)
cons_moduleColors = cons_merge$colors
cons_moduleColors = matchLabels(cons_moduleColors, moduleColors, pThreshold = 5e-2) # Match consensus module colors with original module colors

cons_MEList = moduleEigengenes(datExpr, colors = cons_moduleColors, softPower = softPower)
cons_MEs = cons_MEList$eigengenes

save(cons_MEs, consensusTOM, cons_moduleColors, cons_geneTree, file = "ConsensusNetwork.RData")

# Calculate module membership
modNames=substring(names(cons_MEs),3)
geneModuleBicor=bicorAndPvalue(datExpr, cons_MEs, use = "p")
geneModuleMembership = as.data.frame(geneModuleBicor$bicor)
MMPvalue = as.data.frame(geneModuleBicor$p)
names(geneModuleMembership) = paste("MM", modNames, sep = "")
names(MMPvalue) = paste("p.MM", modNames, sep = "")

geneInfo = data.frame(miRID = colnames(datExpr), moduleColor = cons_moduleColors, geneModuleMembership, MMPvalue)
write.csv(geneInfo, file = "./ConsensusNetwork_ModuleMembership.csv")

参考文献：

  • Plasma extracellular vesicle long RNA profiling identifies a diagnostic signature for the detection of pancreatic ductal adenocarcinoma.
  • Extracellular Vesicles Long RNA Sequencing Reveals Abundant mRNA, circRNA, and lncRNA in Human Blood as Potential Biomarkers for Cancer Diagnosis.
  • Genome-wide, integrative analysis implicates microRNA dysregulation in autism spectrum disorder
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,547评论 6赞 477
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,399评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,428评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,599评论 1赞 274
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,612评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,577评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,941评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,603评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,852评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,605评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,693评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,375评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,955评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,936评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,172评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 43,970评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,414评论 2赞 342