R-Script To Analyze Rna-Seq Data Set Using Wgcna

Hi,

I am currently analyzing 90 RNA-seq samples. I am interested in using WGCNA for my analysis. Nonetheless, when I read the tutorials provided by Steve's group, I felt a little bit lost with some details like probe, and etc. Since most of the command lines are geared towards microarray data sets, I am having a little bit difficulty when I tried running the suggested command lines. I omitted the ones related to microarray in my analysis. And, I kept getting error or warning messages. I could not even continue past soft thresholding calculation step. If you have used WGCNA to analyze RNA-seq data sets using gene counts as input data, can you please share with me the R-script you used for your analysis just to give me an idea how to proceed with mine? I am using normalized and transformed (variance stabilizing transformation in DESeq2) gene counts as my input data.

Edit: Here's what I have so far:

#Normalizing and transforming gene counts
> WGCNATest = read.table("C:/Users/yfy/Desktop/NewMicroscopyDE.txt", row.names =1 , header = T, sep = "\t")
> colData=data.frame(row.names = colnames(WGCNATest),
    +temp=c("20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C","20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C","20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C"),
    +genotype=c("PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi"),
    +time=c("6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12"))
> dds = DESeqDataSetFromMatrix(countData = WGCNATest, colData = colData, design = ~genotype+time+temp)
> colData(dds)$time = relevel(colData(dds)$ time, "6")
> dds2 = DESeq(dds, betaPrior = FALSE)
> varianceStabilizingTransformation(dds3, blind=TRUE)

#Automatic construction of the gene network and identification of modules
> WGCNATree = flashClust(dist(colData), method = "average")
> sizeGrWindow(12,9)
NULL
> par(cex = 0.6)
> par(mar = c(0,4,2,0))
> plot(WGCNATree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2)
> powers = c(c(1:10), seq(from = 12, to=20, by=2))
> options(stringsAsFactors = FALSE)
> sft = pickSoftThreshold(colData, powerVector = powers, verbose = 5)
pickSoftThreshold: will use block size 3.
 pickSoftThreshold: calculating connectivity for given powers...
   ..working on genes 1 through 3 of  3
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
  0 (non-NA) cases
In addition: Warning messages:
1: executing %dopar% sequentially: no parallel backend registered 
2: In (function (x, y = NULL, use = "all.obs", method = c("pearson",  :
  NAs introduced by coercion
3: In (function (x, y = NULL, use = "all.obs", method = c("pearson",  :
  NAs introduced by coercion
4: In eval(expr, envir, enclos) :
  Some correlations are NA in block 1 : 3 .
5: In as.vector(log10(dk)) : NaNs produced

rnaseq

• 13,708 views • link

updated 12.4 years ago by Devon Ryan

• written 12.4 years ago by yfyian •

Why don't you post what you have so far up until the first error?

• 0 views • link

written 12.4 years ago by Devon Ryan

#Normalizing and transforming gene counts

WGCNATest = read.table("C:/Users/yfy/Desktop/NewMicroscopyDE.txt", row.names =1 , header = T, sep = "\t") colData=data.frame(row.names = colnames(WGCNATest), + temp=c("20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C","20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C","20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C"), + genotype=c("PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi"), + time=c("6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12")) dds = DESeqDataSetFromMatrix(countData = WGCNATest, colData = colData, design = ~genotype+time+temp) colData(dds)$time = relevel(colData(dds)$ time, "6") dds2 = DESeq(dds, betaPrior = FALSE) varianceStabilizingTransformation(dds3, blind=TRUE)

#Automatic construction of the gene network and identification of modules

WGCNATree = flashClust(dist(colData), method = "average") sizeGrWindow(12,9) NULL par(cex = 0.6) par(mar = c(0,4,2,0)) plot(WGCNATree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2) powers = c(c(1:10), seq(from = 12, to=20, by=2)) options(stringsAsFactors = FALSE) sft = pickSoftThreshold(colData, powerVector = powers, verbose = 5) pickSoftThreshold: will use block size 3. pickSoftThreshold: calculating connectivity for given powers... ..working on genes 1 through 3 of 3 Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases In addition: Warning messages: 1: executing %dopar% sequentially: no parallel backend registered 2: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 3: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 4: In eval(expr, envir, enclos) : Some correlations are NA in block 1 : 3 . 5: In as.vector(log10(dk)) : NaNs produced sizeGrWindow(9, 5) NULL par(mfrow = c(1,2)) 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")) Error in plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, : object 'sft' not found sft = pickSoftThreshold(colData, powerVector = powers, verbose = 5) pickSoftThreshold: will use block size 3. pickSoftThreshold: calculating connectivity for given powers... ..working on genes 1 through 3 of 3 Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases In addition: Warning messages: 1: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 2: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 3: In eval(expr, envir, enclos) : Some correlations are NA in block 1 : 3 . 4: In as.vector(log10(dk)) : NaNs produced text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])sft$fitIndices[,2],labels=powers,cex=cex1,col="red") Error in text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, : object 'sft' not found

• 0 views • link

written 12.4 years ago by yfyian •

I took the liberty of just updating your post, hope you're OK with that :)

• 0 views • link

written 12.4 years ago by Devon Ryan

Not at all! Thanks for doing that. Definitely easier to read now :-)

• 0 views • link

written 12.4 years ago by yfyian •

Also, please note that varianceStabilizingTransformation(dds3, blind=TRUE) should have been varianceStabilizingTransformation(dds2, blind=TRUE).

• 0 views • link

written 12.4 years ago by yfyian •

3 answers

Your problem appears to actually start here:

> varianceStabilizingTransformation(dds3, blind=TRUE)

which should be:

> vsd = getVarianceStabilizedData(dds2, blind=TRUE)

Then, remember that you normally want to work with the transformed counts rather than the information describing what group(s) each of the values in the count matrix is in (colData). So, this

> WGCNATree = flashClust(dist(colData), method = "average")

should likely be something like:

> vsd2 <- t(vsd) #Many functions expect the matrix to be transposed
> WGCNATree = flashClust(dist(vsd2), method = "average")

Similarly, instead of:

> sft = pickSoftThreshold(colData, powerVector = powers, verbose = 5)

you actually want something like:

> sft = pickSoftThreshold(vsd2, powerVector = powers, verbose = 5)

• 0 views • link

written 12.4 years ago by Devon Ryan

Hi dpryan!!

Thank you for your quick response. I tried the command line on vsd you suggested. Here's what I got:

vsd = getVarianceStabilizedData(dds2, blind=TRUE) Error in getVarianceStabilizedData(dds2, blind = TRUE) : unused argument (blind = TRUE)

It worked when I took out the command blind=TRUE. But, I'd like the transformation to be blind to the sample information specified by the design formula. Hence, I should use blind=TRUE. How do I get it to work together with blind=TRUE?

Thanks again, Yoong

• 0 views • link

written 12.4 years ago by yfyian •

Hi dpryan,

Here's what I did after sending the above message:

# to transform data without using information specified by design formula
> varianceStabilizingTransformation(dds2, blind=TRUE)
> vsd = getVarianceStabilizedData(dds2)

I think my data has been transformed without using any info specified by design formula. And, I assigned new object vsd into my R-analysis. I think that's what you wanted to achieve, right?

Yoong

• 0 views • link

updated 12.4 years ago by Devon Ryan

• written 12.4 years ago by yfyian •

I've generally assumed that getVarianceStabilizedData() always uses blind, since that's the default for varianceStabilizingTransformation(). BTW, just running varianceStabilizingTransformation(dds2, blind=TRUE) without assigning the resulting SummarizedExperiment object to something is pointless (i.e., you're not actually doing anything, since the function doesn't modify it's input). If you really want to be sure, then:

vsd <- varianceStabilizingTransformation(dds2, blind=TRUE)
vc <- assay(vsd)

or something like that

• 0 views • link

written 12.4 years ago by Devon Ryan

I managed to get my soft threshold (power) with your suggestions. Thanks a lot for your input!!!!!

• 0 views • link

written 12.4 years ago by yfyian •

Hi dpryan,

Here's what I have so far:

#Normalizing and transforming gene counts

WGCNATest = read.table("C:/Users/yfy/Desktop/NewMicroscopyDE.txt", row.names =1 , header = T, sep = "\t") colData=data.frame(row.names = colnames(WGCNATest), + temp=c("20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C","20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C","20C","20C","20C","20C","20C","20C","30C","30C","30C","30C","30C","30C"), + genotype=c("PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","PI","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","SAL","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi","RNAi"), + time=c("6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12","6","6","6","12","12","12")) dds = DESeqDataSetFromMatrix(countData = WGCNATest, colData = colData, design = ~genotype+time+temp) colData(dds)$time = relevel(colData(dds)$ time, "6") dds2 = DESeq(dds, betaPrior = FALSE) varianceStabilizingTransformation(dds3, blind=TRUE)

#Automatic construction of the gene network and identification of modules

WGCNATree = flashClust(dist(colData), method = "average") sizeGrWindow(12,9) NULL par(cex = 0.6) par(mar = c(0,4,2,0)) plot(WGCNATree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,cex.axis = 1.5, cex.main = 2) powers = c(c(1:10), seq(from = 12, to=20, by=2)) options(stringsAsFactors = FALSE) sft = pickSoftThreshold(colData, powerVector = powers, verbose = 5) pickSoftThreshold: will use block size 3. pickSoftThreshold: calculating connectivity for given powers... ..working on genes 1 through 3 of 3 Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases In addition: Warning messages: 1: executing %dopar% sequentially: no parallel backend registered 2: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 3: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 4: In eval(expr, envir, enclos) : Some correlations are NA in block 1 : 3 . 5: In as.vector(log10(dk)) : NaNs produced sizeGrWindow(9, 5) NULL par(mfrow = c(1,2)) 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")) Error in plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, : object 'sft' not found sft = pickSoftThreshold(colData, powerVector = powers, verbose = 5) pickSoftThreshold: will use block size 3. pickSoftThreshold: calculating connectivity for given powers... ..working on genes 1 through 3 of 3 Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 0 (non-NA) cases In addition: Warning messages: 1: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 2: In (function (x, y = NULL, use = "all.obs", method = c("pearson", : NAs introduced by coercion 3: In eval(expr, envir, enclos) : Some correlations are NA in block 1 : 3 . 4: In as.vector(log10(dk)) : NaNs produced text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])sft$fitIndices[,2],labels=powers,cex=cex1,col="red") Error in text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, : object 'sft' not found

Many thanks, Yoong

• 0 views • link

written 12.4 years ago by yfyian •

Sorry, the formatting is a little bit off. I will email it to you.

Yoong

• 0 views • link

written 12.4 years ago by yfyian •