Question: WGCNA code for rna-seq data SFT plot issue
0
gravatar for krushnach80
2.9 years ago by
krushnach80830
krushnach80830 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

rna-seq R • 1.8k views
ADD COMMENTlink modified 2.9 years ago • written 2.9 years ago by krushnach80830
1

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

Also, posting the plots would help

ADD REPLYlink written 2.9 years ago by Fabio Marroni2.6k

okay.i will post the plot .

Did you check that your values are not all larger (or smaller) than 0.8? what kind of value ?
ADD REPLYlink written 2.9 years ago by krushnach80830
1

what kind of value ?

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

ADD REPLYlink modified 2.9 years ago • written 2.9 years ago by Fabio Marroni2.6k

okay 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 ...?

ADD REPLYlink written 2.9 years ago by krushnach80830
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 633 users visited in the last hour