Nothing
# TODO: Add comment
#
# Author: ahrnee-adm
###############################################################################
### load / source
library("limma")
library(gplots) # volcano plot
library(seqinr)
library(optparse)
library(data.table)
library(epiR)
library(corrplot)
library(Biobase)
### INIT
if(!grepl("SafeQuant\\.Rcheck",getwd())){ # DEV mode
wd <- dirname(sys.frame(1)$ofile)
setwd(dirname(sys.frame(1)$ofile))
sqRootDir <- dirname(getwd())
source(paste(sqRootDir,"/R/ExpressionAnalysis.R",collapse="",sep=""))
source(paste(sqRootDir,"/R/SafeQuantAnalysis.R",collapse="",sep=""))
source(paste(sqRootDir,"/R/Graphics.R",collapse="",sep=""))
source(paste(sqRootDir,"/R/IdentificationAnalysis.R",collapse="",sep=""))
source(paste(sqRootDir,"/R/Parser.R",collapse="",sep=""))
source(paste(sqRootDir,"/R/TMT.R",collapse="",sep=""))
source(paste(sqRootDir,"/R/UserOptions.R",collapse="",sep=""))
}else{ # CHECK mode
### wd already set to tests when running CHECK
library(SafeQuant)
}
### INIT
### VARIOUS TEST FILES
# progenesis
progenesisFeatureCsvFile1 <- "testData/progenesis_feature_export1.csv"
progenesisPeptideMeasurementCsvFile1 <- "testData/progenesis_pep_measurement1.csv"
#progenesisPeptideMeasurementCsvFile1 <- "testData/tmp.csv"
progenesisProteinCsvFile1 <- "testData//progenesis_protein_export1.csv"
progenesisPeptideMeasurementFractionatedCsvFile1 <- "testData/progenesis_pep_measurement_fractionated1.csv"
progenesisProteinCsvFile2 <- "testData/2014/proteins2.csv"
progenesisFeatureCsvFile2 <- "testData/2014/peptides2.csv"
# scaffold
scaffoldTmt6PlexRawTestFile <- "testData/scaffold_tmt6plex_raw.xls"
scaffoldTmt10PlexRawTestFile <- "testData/scaffold_tmt10plex_raw.xls"
scaffoldTmt10PlexCalibMixRawTestFile <- "testData/scaffold_tmt10plex_calibMix_raw.xls"
#scaffoldPtmTMTRawDataFile1 <- "testData/scaffoldPTM/Christoph-LE-Human-pH10fraction-TMT-20150630/Raw Data Report for Christoph-LE-Human-pH10fraction-TMT-20150630.xls"
#scaffoldPtmReportFile1 <- "testData/scaffoldPTM/Christoph-LE-Human-pH10fraction-TMT-20150630/Spectrum Report of Scaffold_PTM_P-TMT-pH10 Experiment.xls"
scaffoldPtmTMTRawDataFile1 <- "testData/scaffold_tmt10plex_raw_phospho.xls"
scaffoldPtmReportFile1 <- "testData/scaffoldPtm_spectrum_report.xls"
# maxquant
maxQuantProteinFileTxt <- "testData/maxquant_protein_groups.csv"
# db
fastaFile <- "testData/mouse_proteins.fasta"
### INIT END
## CREATE TEST DATA
set.seed(1234)
nbFeatures <- 900
peptide <- paste("pep",1:nbFeatures,sep="")
peptide[1] <- "VALGDGVQLPPGDYSTTPGGTLFSTTPGGTR"
peptide[2] <- "AQAGLTATDENEDDLGLPPSPGDSSYYQDQVDEFHEAR"
proteinName <- sort(rep(paste("prot",1:(nbFeatures/3),sep=""),3))
proteinName[1:200] <- paste("REV_",proteinName[1:200] ,sep="")
proteinName[1] <- "sp|Q60876|4EBP1_MOUSE"
proteinName[2] <- "sp|Q9JI13|SAS10_MOUSE"
idScore <- rep(0,length(proteinName))
idScore[1:200] <-rnorm(200,10,1)
idScore[c(1:2,201:900)] <- rnorm(702,15,1)
ptm <- rep("",900)
ptm[1] <- "[15] Phospho (ST)|[30] Phospho (ST)"
ptm[2] <- "[20] Phospho (ST)"
pMassError <- c(rnorm(200,0,1.5),rnorm(700,0,0.5))
charge <- round(runif(length(ptm),1.5,3.7))
peptideName <- paste(peptide,ptm)
proteinDescription <- sort(rep(paste("protDescription",1:(nbFeatures/3),sep=""),3))
isNormAnchor <- rep(T,nbFeatures)
isFiltered <- rep(F,nbFeatures)
m <- as.matrix( data.frame(rnorm(nbFeatures,1001),rnorm(nbFeatures,1001),rnorm(nbFeatures,1002),rnorm(nbFeatures,1002),rnorm(nbFeatures,1000),rnorm(nbFeatures,1000)) )
rownames(m) <- peptideName
colnames(m) <- c("A_rep_1","A_rep_2","B_rep_1","B_rep_2","C_rep_1","C_rep_2")
### phenoData: stores expDesign
#condition isControl
#A_rep_1 A FALSE
#A_rep_2 A FALSE
#B_rep_1 B FALSE
#B_rep_2 B FALSE
#C_rep_1 C TRUE
#C_rep_2 C TRUE
expDesign <- data.frame(condition=c("A","A","B","B","C","C"),isControl=c(F,F,F,F,T,T),row.names=colnames(m))
#expDesign <- data.frame(condition=c("A","A","B","B","C","C"),row.names=colnames(m))
featureAnnotations <- data.frame(
peptide
, charge
,proteinName
,proteinDescription
,idScore
,ptm
,pMassError
,isNormAnchor
,isFiltered
,row.names=peptideName)
eset <- createExpressionDataset(expressionMatrix=m,expDesign=expDesign,featureAnnotations=featureAnnotations)
sqa <- safeQuantAnalysis(eset)
# ABS. QUANT SIM. DATA
cpc <- rep(2^(1:5),10)
set.seed(1234)
signal <- rnorm(length(cpc),cpc,cpc/10)
absEstSimData <- data.frame(cpc = log10(cpc),signal = log10(signal))
absEstSimDataFit <- lm(cpc ~ signal, data=absEstSimData )
### CREATE PAIRED ESET
set.seed(1234)
esetPaired <- eset
exprs(esetPaired)[,1] <- rnorm(nrow(exprs(eset)),1500,1500/10)
exprs(esetPaired)[,2] <- rnorm(nrow(exprs(eset)),3000,3000/10)
exprs(esetPaired)[,3] <- rnorm(nrow(exprs(eset)),1200,1200/10)
exprs(esetPaired)[,4] <- rnorm(nrow(exprs(eset)),2400,2400/10)
exprs(esetPaired)[,5] <- rnorm(nrow(exprs(eset)),1000,1000/10)
exprs(esetPaired)[,6] <- rnorm(nrow(exprs(eset)),2000,2000/10)
esetPaired <- createPairedExpDesign(esetPaired)
##### TMT
### CREATE TEST DATA
tmtTestData6Plex <- matrix(rep(10,24),ncol=6)
tmtTestData6Plex[2,1:3] <- c(9,9,9)
tmtTestData6Plex[3,1:3] <- c(100,100,100)
tmtTestData6Plex[4,c(1,3,5)] <- c(100,100,100)
tmtTestData10Plex <- matrix(rep(10,100),ncol=10)
esetCalibMix <- parseScaffoldRawFile(file=scaffoldTmt10PlexCalibMixRawTestFile
,expDesign=data.frame(condition=paste("Condition",c(1,2,3,1,2,3,1,2,3,1),sep=""),isControl=c(T,F,F,T,F,F,T,F,F,T) ))
#esetCalibMixPair <- .getCalibMixPairedEset(.getCalibMixEset(esetCalibMix))
### CREATE TEST DATA END
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.