gistfile1.r

library(oligo)
library(genefilter)
library(RColorBrewer)
library(sva)
library(SpeCond)

basepath <- '/n/HSPH/projects/am_trap'

#
# Get array data
#
celFiles <- list.celfiles(file.path(basepath, 'data'), 
                          pattern='*urger*', full.names=T)
celFiles
affy <- read.celfiles(celFiles, verbose=T)

# Transcript (gene) level normalization using RMA
geneCore <- rma(affy, target='core')
geneCore

# Set up covariates based on sample information from the McMahon group
pDataFile <- file.path(basepath, 'data', 'noSurgery.txt')
pDataObj <- read.table(pDataFile, row.names=1, header=T, sep='\t')
pData(geneCore) <- pDataObj
pData(geneCore)

# Retrieving NetAffx Biological Annotation
featureData(geneCore) <- getNetAffx(geneCore, 'transcript')
varLabels(featureData(geneCore))

# Extract the 'gene assignment' annotation
annot <- pData(featureData(geneCore)[, c('geneassignment')])
head(annot[!is.na(annot), ], 2)

# Generate a list of gene symbols from the gene assignment
desc <- annot[, 1]
symbols <- unlist(lapply(desc, function(x) strsplit(x, ' // ')[[1]][2]))
length(featureData(geneCore)$probesetid) == length(symbols)
head(symbols[!is.na(symbols)])

# Merge probe identifier and symbols
desc <- paste(featureData(geneCore)$probesetid,
              symbols,
              sep='.')


#
# SpeCond testing
#
library(SpeCond)
factors <- pData(geneCore)$Celltype

# Default parameters are a bit too stringent
param.detection <- getDefaultParameter()
param.detection
param.detection2 <-createParameterMatrix(per.1=0.29)

# Start with the default parameters; get an expression matrix
Mexp <- exprs(geneCore)
generalResult <- SpeCond(Mexp, 
                         param.detection=param.detection2, 
                         multitest.correction.method="BY",
                         prefix.file="E", 
                         print.hist.pv=TRUE, 
                         fit1=NULL, 
                         fit2=NULL, 
                         specificOutlierStep1=NULL)
specificResult <- generalResult$specificResult

# Repeat with expression set
generalResultS <- SpeCond(geneCore, 
                         param.detection=param.detection2,
                         multitest.correction.method="BY",
                         prefix.file="E", 
                         print.hist.pv=TRUE, 
                         fit1=NULL,
                         fit2=NULL, 
                         specificOutlierStep1=NULL, 
                         condition.factor=factors,
                         condition.method="mean")
specificResultS <- generalResultS$specificResult

MexpS <- getMatrixFromExpressionSet(geneCore,
                                    condition.factor=factors,
                                    condition.method="mean")

# Take a look at the results
getFullHtmlSpeCondResult(SpeCondResult=generalResultS, 
                         param.detection=specificResultS$param.detection, 
                         page.name="Example_SpeCond_results",
                         page.title="Tissue specific results", 
                         sort.condition="all", 
                         heatmap.profile=TRUE,
                         heatmap.expression=TRUE, 
                         heatmap.unique.profile=TRUE, 
                         expressionMatrix=MexpS)
                                    
rownames(MexpS) <- desc
geneRows <- as.vector(row(MexpS)[specificResultS$L.specific.result$M.specific.sum.row == 1, ][, 1])
genePageInfo <- getGeneHtmlPage(MexpS, 
                                specificResultS, 
                                name.index.html="index_example_SpeCond_Results.html",
                                gene.html.ids=geneRows)

                         

# Export
L.specific.result.profile <- getProfile(specificResultS$L.specific.result$M.specific)

writeSpeCondResult(specificResultS$L.specific.result,
                   file.name.profile="Example_specific_profile.txt", 
                   file.specific.gene="Example_list_specific_gene.txt",
                   file.name.unique.profile="Example_specific_unique_profile.txt")

writeUniqueProfileSpecificResult(L.specific.result=specificResultS$L.specific.result,
                                 file.name.unique.profile="Example_specific_unique_profile.txt", 
                                 full.list.gene=FALSE)

writeGeneResult(specificResultS$L.specific.result, 
                file.name.result.gene="Example_gene_gummary_result.txt")