This mornings post is an summery of the work I did on Friday and finished up over the weekend. There was a lot of coding and great visualizations! Lauren gave us a WGCNA tutorial. Below is the R code.

source(“http://bioconductor.org/biocLite.R”) biocLite(“impute”) install.packages(“WGCNA”) library(WGCNA); # The following setting is important, do not omit. options(stringsAsFactors = FALSE); #Read in the female liver data set SSData = read.csv(“Phel_rnaseq_normalized_expression.csv”); # Take a quick look at what is in the data set: dim(SSData); names(SSData); datExpr0 = as.data.frame(t(SSData )); names(datExpr0) = SSData$Contig; rownames(datExpr0) = names(SSData) ; datExpr0 datExpr=datExpr0 gsg = goodSamplesGenes(datExpr0, verbose = 3); gsg$allOK collectGarbage(); powers = c(c(1:10), seq(from = 12, to=20, by=2)) # Call the network topology analysis function sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5) # Plot the results: sizeGrWindow(9, 5) par(mfrow = c(1,2)); cex1 = 0.9; # Scale-free topology fit index as a function of the soft-thresholding power 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”); # this line corresponds to using an R^2 cut-off of h abline(h=0.90,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”) net = blockwiseModules(datExpr, power = 10, TOMType = “unsigned”, minModuleSize = 9, reassignThreshold = 0, mergeCutHeight = 0.25, numericLabels = TRUE, pamRespectsDendro = FALSE, saveTOMs = TRUE, saveTOMFileBase = “femaleMouseTOM”, verbose = 3) sizeGrWindow(12, 9) # Convert labels to colors for plotting mergedColors = labels2colors(net$colors) # Plot the dendrogram and the module colors underneath plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]], “Module colors”, dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)

After that we worked on finding DEG’s that were interesting and categorizing them based on function for those with a high fold number. Overall, just continuing to sort out what is happening in the massive immune response the exposed animals display.