Hi James.

Very nice. Thank you. I was just looking at an email from that postdoc I consulted several weeks ago and this is what she said. Could you please explain as I am quite confused, working under an encroaching deadline, and without benefit the classmate or PI:

pData(gse) returns to the data frame of sample information only. A bit more idea of how gse ExpressionSet object looks like...you could imagine it is a multi-tab excel table. For example, when we extract the "sampledata <- pData(gse[[1]])", it is one of the tabs containing the pData data frame. For the downstream, you should look for "exprs(gse)" for the expression data, in which the rows are the genes and the columns are the sample IDs. Just do the differential expression analysis using the limma package instead of DESeq2. Now, you have the gse object which contains the assayData (expression data), phenoData (sample information), and featureData (gene annotation data).

1. Proceed to check if you need to normalize the expression data
2. Check the sample information, which column should be used as the "group"?
3. Proceed to the differential expression analysis using the limma package, you should get the "TopTable" results which include the gene IDs as the rows and the DE result as the column, we are interested in the columns of "log2FoldChange" and "adjusted p value". And then, you have to follow the instructions from clusterprofiler book. By using "log2FoldChange" and "adjusted p value", we can form a genelist object that is used in clusterprofiler package for the GO and KEGG enrichment analysis, as well as visualizing the results in the network figures.

I think the first step involves calling log2 on the ExpressionSet object? The second step is puzzling, since upon using the following code: sampleData <- pData(gse[[1]]), I don't know what column I should use as the "group" nor do I know what that means (GSE12657) although limma apparently depends upon the specification of one column name? Could you please explain 3? The postdoc appears to believe this step is obvious but it most certainly is not.

From the paper I am seeking to reproduce, the authors mentioned "hub genes" and the cytohubba package from Cytoscape. Would I need to transform the data at some point with the WGCNA package?

I spent the past several hours reading the following document and my code, thus far, is below.

https://github.com/Lindseynicer/How-to-analyze-GEO-microarray-data

Update:

library(GEOquery)
library(limma)
library(pheatmap)
library(dplyr)
library(ggplot2)
library(ggrepel)

gse <- getGEO("GSE12657", AnnotGPL = TRUE)
gse

# print the sample information
pData(gse[[1]])[1:2,]

# print the gene annotation
fData(gse[[1]])[1,]

# print the expression data
exprs(gse[[1]])[1,]

# Check the normalisation and scales used
pData(gse[[1]])$data_processing[1]

# Look at the expression value
summary(exprs(gse[[1]]))

# Perform a log2 transformation since the range (min and max) is not between 0 and 16
exprs(gse[[1]]) <- log2(exprs(gse[[1]]))

# Reassess the summary
summary(exprs(gse[[1]]))

# Plot the summary values
boxplot(exprs(gse),outline=F)

# Retrieve sampleInfo
sampleInfo <- pData(gse[[1]])

# Obtain table counts for Control, Glioblastoma, Oligodendroglioma, and 
#                         Pilocytic astrocytoma 
table(sampleInfo$characteristics_ch1)

# Select factors we need for the analysis
sampleInfo <- select(sampleInfo, characteristics_ch1)

# Rename the column to one that is more convenient
sampleInfo <- rename(sampleInfo, sample = characteristics_ch1)

library(pheatmap)
# The argument use="c" stops an error if there are any missing data points
# We can visualize the correlations between the samples by hierarchical clustering.
# The function ‘cor’ can calculate the correlation on the scale of 0 to 1, in a pairwise fashion between all samples, then 
# visualise on a heatmap. There are many ways to create heatmaps in R, here I use ‘pheatmap’, the only argument it requires # is a matrix of numeric values.

# We can add more sample info onto the plot to get a better pic of the group and clustering. Here, we make use of the # ‘sampleInfo’ file that was created earlier, to match with the columns of the correlation matrix.
corMatrix <- cor(exprs(gse[[1]]),use="c")
pheatmap(corMatrix)

# This this case, rownames(sampleInfo) == colnames(sampleInfo), so the correlation is 1.0

# Construct a design matrix will be created, this is a matrix of 0 and 1, one row for each 
# sample and one column for each sample group.
design <- model.matrix(~0 + sampleInfo$sample)

# Calculate the median expression level
threshold <- median(exprs(gse[[1]]))

#TRUE or FALSE for whether each gene is "expressed" in each sample
gene_is_expressed <- exprs(gse[[1]]) > threshold

# Identify genes that are expressed in more than 2 samples
retain <- rowSums(gene_is_expressed) > 3

# Check how many genes are removed / retained.
table(retain)

# Subset the data to keep only the expressed genes
**Error** 
gse <- gse[retain,]
Error in orig[[nm]][i, , ..., drop = drop] : 
  (subscript) logical subscript too long
gse <- gse[retain,]

Thank you James for your help. I greatly appreciate it.