1. Limma (Screening for Significantly Differentially Expressed Genes)
library(limma)
data<-read.table("expression-data.txt",sep = "\t",header=T,row.names=1)
data<- as.matrix(data)
mode(data)
group<-c(rep(1,85),rep(2,85)) 
design <- model.matrix(~ -1+factor(group)) 
colnames(design)<-c("Nonlesion","Lesion") 
design
contrast.matrix <- makeContrasts(Lesion-Nonlesion,levels=design)
fit <- lmFit(data, design)
fit1 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit1)
dif <- topTable(fit2, coef = 1,  n = nrow(fit2), lfc = 0) 
write.table(dif, "test.dif.txt", row.names = TRUE, sep = "\t")

2. clusterProfiler (Functional pathway analysis)
library(clusterProfiler)
library("org.Hs.eg.db")
A=read.table("input",header=T,sep="\t")
geneList = A[,1]
go <- enrichGO(geneList, OrgDb = org.Hs.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
kegg <- enrichKEGG(gene = genes,keyType = 'kegg', organism = 'hsa', pAdjustMethod = 'fdr', pvalueCutoff = 0.05,qvalueCutoff = 0.2)

3. one-way logistic regression
library(plyr)
library(rms)
library(epiDisplay)
A=read.table("input.txt",header=T,sep="\t")
glm1<- glm(A[,1]~A[,N],data=A,family = binomial)
glm2<- summary(glm1);glm2
OR<-round(exp(coef(glm1)),2)
SE<-glm2$coefficients[,2]
CI5<-round(exp(coef(glm1)-1.96*SE),2)
CI95<-round(exp(coef(glm1)+1.96*SE),2)
CI<-paste0(CI5,'-',CI95)
P<-round(glm2$coefficients[,4],3)
res1<-data.frame(OR,CI,P)[-1,];res1
summary(glm1)

4. Optimized Screening
library（lars）
data=read.table("input.txt",header=T)
x <- as.matrix(data[,2:N])
y <- data[,1]
alpha1_fit <- glmnet(x,y,alpha=1,family="gaussian")
plot(alpha1_fit,xvar="lambda",label=TRUE)
alpha1.fit <- cv.glmnet(x,y,type.measure = "mse",alpha=1,family="gaussian")
plot(alpha1.fit)

5. model construction
library(e1071)
library(pROC)
A=read.table("input.txt",header=T,sep="\t")
x <- A[,N]
y <- A[,1]
model <- svm(x, y) 
pred <- predict(model, x)
pred <- fitted(model)
model <- svm(x, y) 
modelroc <- roc(y,pred)
plot(modelroc, print.auc=TRUE, auc.polygon=TRUE, grid=c(0.1, 0.2),grid.col=c("grey", "grey"), max.auc.polygon=FALSE,auc.polygon.col="white", print.thres=TRUE,lwd=2)

6. Subtype analysis
A=read.table("input.txt",header=T)
d=data.matrix(A[,-1])
title=tempdir("Documents and Settings\Administrator")
results = ConsensusClusterPlus(d,maxK=6,reps=50,pItem=0.8,pFeature=1,title=title,clusterAlg="hc",distance="pearson",plot="pdf")

7.WGCNA
datExprA1p=read.table("input.txt",header=T,sep="\t")
datExprA1p=as.matrix(datExprA1p[,-1])
keepGenesExpr = rank(-rowMeans(datExprA1p))
datExprA1g = datExprA1p[keepGenesExpr,] 
adjacencyA1 = adjacency(t(datExprA1g),power=softPower,type="signed")
diag(adjacencyA1)=0
dissTOMA1 = 1-TOMsimilarity(adjacencyA1, TOMType="signed") 
geneTreeA1 = flashClust(as.dist(dissTOMA1), method="average") 
mColorh=NULL 
for (ds in 0:3){
tree = cutreeHybrid(dendro = geneTreeA1,pamStage=FALSE,minClusterSize = 30,cutHeight = 0.95,deepSplit = ds,distM = dissTOMA1)
mColorh=cbind(mColorh,labels2colors(tree$labels)); 
} 
plotDendroAndColors(geneTreeA1, mColorh, paste("dpSplt =",0:3), main = "",dendroLabels=FALSE)
modulesA1 =  mColorh[,1]
PCs1A = moduleEigengenes(t(datExprA1g), colors=modulesA1) 
ME_1A = PCs1A$eigengenes 
distPC1A = 1-abs(cor(ME_1A,use="p")) 
distPC1A = ifelse(is.na(distPC1A), 0, distPC1A) 
pcTree1A = hclust(as.dist(distPC1A),method="a") 
MDS_1A = cmdscale(as.dist(distPC1A),2) 
colorsA1 = names(table(modulesA1)) 
plotDendroAndColors(geneTreeA1, modulesA1, "Modules", dendroLabels=FALSE, hang=0.03, addGuide=TRUE, guideHang=0.05, main="RNAs dendrogram and module colors")