library(affyPLM)
rm(list=ls())
rawdata <- ReadAffy() 
eset <- rma(rawdata)
exprSet<-exprs(eset)


source("http://bioconductor.org/biocLite.R")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")#
BiocManager::install("hgu133plus2.db")
library(hgu133plus2.db)
ls('package:hgu133plus2.db')
ids=toTable(hgu133plus2SYMBOL)
dim(ids)
length(unique(ids$symbol))
table(ids$symbol)
sort(table(ids$symbol))
tail(sort(table(ids$symbol)))
table(sort(table(ids$symbol)))
plot(table(sort(table(ids$symbol))))
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]     
dim(exprSet)    
ids=ids[match(rownames(exprSet), ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
tmp=by(exprSet, ids$symbol, function(x) rownames(x) [which.max(rowMeans(x))])
probes<- as.character(tmp)
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% probes, ]
dim(exprSet)
str(exprSet)
rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
exprSet[1:5,1:5]


exprSet["GAPDH",]
exprSet["ACTB",]
boxplot(exprSet)



exprSet1=exprSet[,1:10]
exprSet2=exprSet[,11:20]
exprSet3=exprSet[,21:30]
str(exprSet)
str(exprSet1)
str(exprSet2)
str(exprSet3)
colnames(exprSet1)
colnames(exprSet2)
colnames(exprSet3)
exprSet=cbind(exprSet3,exprSet1,exprSet2)

write.csv(exprSet,"exprSet.csv")


colnames(exprSet)
library(stringr)
colnames(exprSet)=str_split(colnames(exprSet),"\\.",simplify = T)[,1]
colnames(exprSet)


library(reshape2)
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('symbol','sample','value')


group_list=c(rep("NOR",10),rep("SSc",10),rep("SScPH",10))
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)



dim(exprSet)
20186*0.25
sd=apply(exprSet,1,sd)
length(sd)
head(sd)
sd=sort(sd,decreasing = T)
sd=names(sd[1:5046])



library(WGCNA)

datExpr<-exprSet[rownames(exprSet)%in%sd,]
colnames(datExpr)
colnames(exprSet)
rownames(datExpr)
datTraits<-as.data.frame(colnames(datExpr))

datTraits$group<-c(rep("NOR",10),rep("SSc",10),rep("SScPH",10))
colnames(datTraits)<-c("sample","group")
rownames(datTraits)<-datTraits[,1]

head(datExpr)
dim(datExpr)

datExpr0=as.data.frame(t(datExpr))
enableWGCNAThreads()

gsg = goodSamplesGenes(datExpr0, verbose = 3);
gsg$allOK


sampleTree = hclust(dist(datExpr0), method = "average")
#sizeGrWindow(12,9)
par(cex = 0.6)
par(mfrow=c(1,1))
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)
abline(h = 100, col = "red");
clust = cutreeStatic(sampleTree, cutHeight = 100, minSize = 10)
table(clust)
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)


powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
par(mfrow = c(1,1))
par(mar=c(5,12,4,12))
cex1 = 0.9;
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
abline(h=0.85,col="red")
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")


softPower=7
ADJ1_cor <- abs(WGCNA::cor( datExpr,use = "p" ))^softPower
k <- as.vector(apply(ADJ1_cor,2,sum,na.rm=T))
hist(k)
scaleFreePlot(k,main="Check Scale free topology\n")




net = blockwiseModules(datExpr, power = 7, maxBlockSize = 6000,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "AS-green-FPKM-TOM",
                       verbose = 3)
table(net$colors)
mergedColors = labels2colors(net$colors)
table(mergedColors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)



####
design=model.matrix(~0+ datTraits$group)
e=as.factor(datTraits$group)
str(e)
levels(e)
colnames(design)=levels(e)
moduleColors <- labels2colors(net$colors)
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, design , use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
sizeGrWindow(10,6)
textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) = dim(moduleTraitCor)
par(mar = c(6, 8.5, 3, 3));
str(design)
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = colnames(as.data.frame(design)),
               yLabels = names(MEs),
                ySymbols = names(MEs),
                colorLabels = FALSE,
                colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
              main = paste("Module-trait relationships"))




par(cex = 0.6)
par(mfrow=c(1,1))
par(mar = c(8,4,2,2))

GS1 <- as.numeric(WGCNA::cor(design[,2],datExpr,use="p",method="pearson"))
GeneSignificance <- abs(GS1)
ModuleSignificance <- tapply(GeneSignificance,mergedColors,mean,na.rm=T)
plotModuleSignificance(GeneSignificance, moduleColors, ylim=c(0,0.6), main="Module Significance",las=3)

GS2 <- as.numeric(WGCNA::cor(design[,3],datExpr,use="p",method="pearson"))
GeneSignificance <- abs(GS2)
ModuleSignificance <- tapply(GeneSignificance,mergedColors,mean,na.rm=T)
plotModuleSignificance(GeneSignificance, moduleColors, ylim=c(0,0.6), main="Module Significance",las=3)







modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
names(geneModuleMembership) = paste("MM", modNames, sep="")
head(MMPvalue)
names(MMPvalue) = paste("p.MM", modNames, sep="")









SScPH = as.data.frame(design[,3])
names(SScPH) = "SScPH"
geneTraitSignificance = as.data.frame(cor(datExpr, SScPH, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", names(SScPH), sep="")
names(GSPvalue) = paste("p.GS.", names(SScPH), sep="")


module = "midnightblue"
column = match(module, modNames)
moduleGenes = moduleColors==module
sizeGrWindow(7, 7)

par(cex = 0.6)
par(mfrow=c(1,1))
par(mar = c(6,6,6,2))

verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   
                   xlab = paste("Module Membership in", module, "module"),
                   
                   ylab = "Gene significance for SScPH",
                   
                   main = paste("Module membership vs. gene significance\n"),
                   
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)





SSc = as.data.frame(design[,2])
names(SSc) = "SSc"
geneTraitSignificance = as.data.frame(cor(datExpr, SSc, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
names(geneTraitSignificance) = paste("GS.", names(SSc), sep="")
names(GSPvalue) = paste("p.GS.", names(SSc), sep="")


module = "turquoise"
column = match(module, modNames)
moduleGenes = moduleColors==module

par(cex = 0.6)
par(mfrow=c(1,1))
par(mar = c(6,6,6,2))

verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                   
                   abs(geneTraitSignificance[moduleGenes, 1]),
                   
                   xlab = paste("Module Membership in", module, "module"),
                   
                   ylab = "Gene significance for SSc",
                   
                   main = paste("Module membership vs. gene significance\n"),
                   
                   cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)


module = "turquoise"
probes.wgcna = colnames(datExpr) 
inModule = (moduleColors==module)
modProbes = probes.wgcna[inModule]
length(modProbes)
write.table(modProbes,"turquoise.txt")

module = "midnightblue"
probes.wgcna = colnames(datExpr) 
inModule = (moduleColors==module)
modProbes = probes.wgcna[inModule]
length(modProbes)
write.table(modProbes,"midnightblue.txt")



###

SSc= as.data.frame(design[,2])
names(SSc) = "SSc"
GS.SSc = as.numeric(cor(SSc,datExpr, use="p"))

datKME=signedKME(datExpr, MEs, outputColumnName="MM.")
head(datKME)

hightraitssc=abs(GS.SSc)> .5
hightraitssc=datExpr[hightraitssc==T]
table(hightraitssc)
hightraitssc=colnames(hightraitssc)
write(hightraitssc,"hightraitssc.txt")


highconnectssc=abs(datKME$MM.turquoise)>.9
highconnectssc=datExpr[highconnectssc==T]
highconnectssc=colnames(highconnectssc)
write(highconnectssc,"highconnectssc.txt")

###
SScPH= as.data.frame(design[,3])
names(SScPH) = "SScPH"
GS.SScPH = as.numeric(cor(SScPH,datExpr, use="p"))

datKME=signedKME(datExpr, MEs, outputColumnName="MM.")
head(datKME)


hightraitsscph=abs(GS.SScPH)> .5
hightraitsscph=datExpr[hightraitsscph==T]
hightraitsscph=colnames(hightraitsscph)
write(hightraitsscph,"hightraitsscph.txt")


highconnectsscph=abs(datKME$MM.midnightblue)>.9
highconnectsscph=datExpr[highconnectsscph==T]
highconnectsscph=colnames(highconnectsscph)
write(highconnectsscph,"highconnectsscph.txt")


########

library(ggplot2)

library(ggpubr)

library(dplyr)


WDR36=exprSet_L[exprSet_L$symbol=="WDR36",]
SRSF1=exprSet_L[exprSet_L$symbol=="SRSF1",]
POLR1B=exprSet_L[exprSet_L$symbol=="POLR1B",]
TLR2=exprSet_L[exprSet_L$symbol=="TLR2",]
TNFAIP6=exprSet_L[exprSet_L$symbol=="TNFAIP6",]

testdata1=rbind(WDR36,POLR1B,TNFAIP6)
com_anova<- compare_means( value ~ group, data = testdata1, method = "anova" , group.by = "symbol")
com_anova

com_t<- compare_means( value ~ group,, data = testdata1, method = "t.test" , group.by = "symbol", paired = F)
com_t


ggboxplot(testdata1, x="group", y="value", color = "group", 
          palette = "jco",add = "jitter",facet.by = "symbol", short.panel.labs = FALSE)+
  stat_compare_means(method = "anova", label.y = 10.5)+ # Add global p-value
  stat_compare_means(label = "p.format", method = "t.test", comparisons = list(c("NOR","SSc")),label.y = 8.6)+
  stat_compare_means(label = "p.format", method = "t.test", comparisons = list(c("SSc","SScPH")))


testdata2=rbind(SRSF1,TLR2)
com_anova<- compare_means( value ~ group, data = testdata2, method = "anova" , group.by = "symbol")
com_anova

com_t<- compare_means( value ~ group,, data = testdata2, method = "t.test" , group.by = "symbol", paired = F)
com_t


ggboxplot(testdata2, x="group", y="value", color = "group", 
          palette = "jco",add = "jitter",facet.by = "symbol", short.panel.labs = FALSE)+
  stat_compare_means(method = "anova", label.y = 12)+ # Add global p-value
  stat_compare_means(label = "p.format", method = "t.test", comparisons = list(c("NOR","SSc")),label.y = 11.5)+
  stat_compare_means(label = "p.format", method = "t.test", comparisons = list(c("SSc","SScPH")))