Question: WGCNA code for rna-seq data SFT plot issue

0

krushnach80 •

**890**wrote:```
library(WGCNA)
library(flashClust)
options(stringsAsFactors = FALSE);
countdata <- read.csv('module500.csv', header=TRUE)
dim(countdata)
names(countdata)
femData <- countdata
names(femData)
datExpr = as.data.frame(t(femData[, -c(1)]))
View(datExpr)
dim(datExpr)
names(datExpr) = femData$Symbol
View(datExpr)
rownames(datExpr) = names(femData)[-c(1)]
dim(datExpr)
powers = c(c(1:10), seq(from = 12, to=20, by=2))
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,powerVector = powers,corFnc = cor,corOptions = list(use = 'p'),networkType = "unsigned")
sizeGrWindow(9, 5)
par(mfrow = c(1,2))
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")
# Red line corresponds to using an R^2 cut-off
abline(h=0.80,col="red")
# Mean connectivity as a function of the soft-thresholding power
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
adj= adjacency(datExpr,type = "unsigned", power = softPower)
#turn adjacency matrix into topological overlap to minimize the effects of noise and spurious associations
TOM=TOMsimilarityFromExpr(datExpr,networkType = "unsigned", TOMType = "unsigned", power = softPower)
colnames(TOM) =rownames(TOM) = datExpr
dissTOM=1-TOM
geneTree = flashClust(as.dist(dissTOM),method="average")
plot(geneTree, xlab="", sub="",cex=0.3)
minModuleSize = 20;
# Module identification using dynamic tree cut
dynamicMods = cutreeDynamic(dendro = geneTree, method="tree", minClusterSize = minModuleSize);
#dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, method="hybrid", deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize);
#the following command gives the module labels and the size of each module. Lable 0 is reserved for unassigned genes
table(dynamicMods)
dynamicColors = labels2colors(dynamicMods)
table(dynamicColors)
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
restGenes= (dynamicColors != "grey")
diss1=1-TOMsimilarityFromExpr(datExpr[,restGenes], power = softPower)
colnames(diss1) =rownames(diss1) = datExpr[restGenes]
hier1=flashClust(as.dist(diss1), method="average" )
plotDendroAndColors(hier1, dynamicColors[restGenes], "Dynamic Tree Cut", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05, main = "Gene dendrogram and module colors")
diag(diss1) = NA
sizeGrWindow(7,7)
TOMplot(diss1, hier1, as.character(dynamicColors[restGenes]))
module_colors= setdiff(unique(dynamicColors), "grey")
for (color in module_colors){
module=datExpr[which(dynamicColors==color)]
write.table(module, paste("module_",color, ".txt",sep=""), sep="\t", row.names=TRUE, col.names=TRUE,quote=FALSE)
}
module.order <- unlist(tapply(1:ncol(datExpr),as.factor(dynamicColors),I))
m<-t(t(datExpr[,module.order])/apply(datExpr[,module.order],2,max))
heatmap(t(m),zlim=c(0,1),col=gray.colors(100),Rowv=NA,Colv=NA,labRow=NA,scale="none",RowSideColors=dynamicColors[module.order])
MEList = moduleEigengenes(datExpr, colors = dynamicColors)
MEs = MEList$eigengenes
plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2))
```

So this is my code im using for WGCNA ,but at this point `abline(h=0.80,col="red")`

i not getting an abline in SFT plot ,am i doing something wrong or what is the issue.

Any help or suggestion would be highly appreciated

Did you check that your values are not all larger (or smaller) than 0.8?

Also, posting the plots would help

2.7kokay.i will post the plot .

890The y values that you are plotting in the graph in which you are trying to plot the red line

2.7kokay i got it my values starts from -0.8 to 0.0 the y axis , so when i change my value of h to -0.2 i get the line but how do i choose value of h ...?

890