inst/doc/ProCoNA_Vignette.R

In ProCoNA: Protein co-expression network analysis (ProCoNA).

### R code from vignette source 'ProCoNA_Vignette.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: Removing Missing Data
###################################################
options(keep.source = TRUE, width = 70, stringsAsFactors=FALSE, digits=2)
library(ProCoNA)
data(ProCoNA_Data)
peptideData <- subsetPeptideData(peptideData, percentageNAsAllowed=0.2)
dim(peptideData)


###################################################
### code chunk number 2: Picking the scaling number.
###################################################
beta <- pickSoftThreshold(peptideData, networkType="signed", RsquaredCut=0.8)
beta$powerEstimate


###################################################
### code chunk number 3: Building the ProCoNA Network
###################################################
peptideNetwork <- buildProconaNetwork(networkName="my network", 
                                      pepdat=peptideData, 
                                      networkType="signed",
                                      pow=beta$powerEstmate, 
                                      pearson=FALSE, 
                                      toPermTestPermutes=1000)
peptideNetwork


###################################################
### code chunk number 4: Viewing proconaNet slots
###################################################
getSlots("proconaNet")


###################################################
### code chunk number 5: Network Summary
###################################################
printNet(peptideNetwork)


###################################################
### code chunk number 6: Accessor functions
###################################################
networkName(peptideNetwork)
samples(peptideNetwork)[1:5]
peptides(peptideNetwork)[1:5]
mergedColors(peptideNetwork)[1:5]


###################################################
### code chunk number 7: Module Significance
###################################################
peptideNetwork <- toPermTest(peptideNetwork, 100)
permtest(peptideNetwork)


###################################################
### code chunk number 8: Plot Network Dendrogram
###################################################
plotNet(peptideNetwork)


###################################################
### code chunk number 9: Operating on modules
###################################################
module1 <- which(mergedColors(peptideNetwork) == 1)
module1_TOM <- TOM(peptideNetwork)[module1, module1]
mean(utri(module1_TOM))


###################################################
### code chunk number 10: Correlation of module eigenvectors with phenotypes.
###################################################

phenotypeCor <- correlationWithPhenotypesHeatMap(net=peptideNetwork,
                                                 phenotypes=phenotypes[,1:5],
                                                 modules=1:5,                
                                                 plotName="my plot",
                                                 title="snazzy heatmap",
                                                 textSize=0.7)

pepcor <- moduleMemberCorrelations(pnet=peptideNetwork,
                                      pepdat=peptideData,
                                      phenotypes=phenotypes)

#########################################################################
# quick function to write out the tables for specific modules.
moduleData <- function(pepnet, pepcors, module, pepinfo, fileprefix) {
  moduleX <- peptides(pepnet)[which(mergedColors(pepnet)==module)]
  moduleInfo <- pepinfo[which(pepinfo$Mass_Tag_ID %in% moduleX),]
  moduleCors <- pepcors[which(pepcors$Module==module),]
  corname <- paste(fileprefix, "_correlations.csv", sep="")
  write.table(moduleCors, file=corname, sep=",", row.names=F)
  infoname <- paste(fileprefix, "_peptide_info.csv", sep="")
  write.table(moduleInfo, file=infoname, sep=",", row.names=F)
}
########################################################################

# WRITE OUT A TABLE WITH THE BELOW FUNCTION CALL#
# moduleData(peptideNetwork, pepcor, 1, masstagdb, "Module_1")


###################################################
### code chunk number 11: Working from raw data.
###################################################
library(MSnbase)
file <- dir(system.file(package = "MSnbase", dir = "extdata"),
            full.names = TRUE, pattern = "mzXML$")
rawdata <- readMSData(file, msLevel = 2, verbose = FALSE)
experiment <- removePeaks(itraqdata, t = 400, verbose = FALSE)
experiment <- trimMz(experiment, mzlim = c(112, 120))
qnt <- quantify(experiment,
                method = "trap",
                reporters = iTRAQ4,
                strict = FALSE,
                parallel = FALSE,
                verbose = FALSE)
qnt.quant <- normalise(qnt, "quantiles")
gb <- fData(qnt)$PeptideSequence
qnt2 <- combineFeatures(qnt.quant, groupBy = gb, fun = "median")


###################################################
### code chunk number 12: filterna
###################################################
sum(is.na(qnt2))
qnt2 <- filterNA(qnt2, pNA = 0)
sum(is.na(qnt2))


###################################################
### code chunk number 13: proconamsnset (eval = FALSE)
###################################################
## peptideNetwork <- buildProconaNetwork(networkName="my network", 
##                                       pepdat=qnt2, 
##                                       networkType="signed", 
##                                       pow=beta$powerEstmate, 
##                                       pearson=FALSE, 
##                                       toPermTestPermutes=1000)