### R code from vignette source 'vignettes/ProCoNA/inst/doc/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)