xinv xinv%*%I pbinom(7, size=12, prob=0.02) 1-pbinom(7, size=12, prob=0.02) -10*log10(1-pbinom(7, size=12, prob=0.02)) -10*log10(1-pbinom(6, size=12, prob=0.02)) plot q q() message(log(-1)) traceback() f = function(x) {r = x-g(x)} r library(isobar) data(ibspiked_set1) ib = ibspiked_set1; rm(ibspiked_set1) debug(normalize) normalize(ib) undebug(normalize) traceback() library(vsn) data("kidney") kidney justvsn vsn2 showMethods("vsn2") qt(0.95, 10) 0.95/108 qt(0.008796296, 52) source('C:/Research/Projects/QUANrevision/201508/generalizedESD.R') y source('C:/Research/Projects/QUANrevision/201508/generalizedESD.R') source('C:/Research/Projects/QUANrevision/201508/generalizedESD.R') source('C:/Research/Projects/QUANrevision/201508/generalizedESD.R') y x = rnorm(100) x dcauchy(x) y = dcauchy(x) plot(x,y) x = rcauchy(100, location = 0, scale = 1) y1 = dcauchy(x, location = 0, scale = 1) y2 = dcauchy(x, location = 1, scale = 1) plot(x, y1) plot(x, y2) plot(x, y1) sum(y1) sum(y2) y2 = dcauchy(x, location = 0, scale = 5) plot(x, y2) library(limma) removeBatchEffect() removeBatchEffect pnorm ??pnorm ?pnorm pnorm(0,0,1) getwd() return(est.err) # simulate N sets of observations N = 1000; # number of sets of observations data.bl = cbind(rep(1,N),rep(3,N),rep(10,N)) # baseline data err = 1; # noise level to estimate err1 = matrix(rnorm(3*N,sd=0.1),N,3) # random noise data.obs = data.bl+err+err1; # strategy #1 func.1 = function(data){ mn = apply(data,2,mean) X=mn[1]; Y=mn[2]; Z=mn[3]; a1=(2/3)^2; b1=-2*2/3*(X-1/3*Y); a2=(9/10)^2; b2=-2*9/10*(X-1/10*Z); a=a1+a2; b=b1+b2; est.err = -b/(2*a); return(est.err) } func.1(data.obs) # strategy #2 Actually, it is the same with strategy #1 func.2 = function(data){ a1=(2/3)^2; a2=(9/10)^2; X=data[,1];Y=data[,2];Z=data[,3]; b1= -2*2/3*(X-1/3*Y); b2=-2*9/10*(X-1/10*Z); a=a1+a2; b=mean(b1)+mean(b2) est.err = -b/(2*a); return(est.err) } func.2(data.obs) data data.bl data.obs mean(data.obs) apply(data.obs, 'mean', 2) apply(data.obs, 'mean', 1) apply(data.obs, 1, mean) apply(data.obs, 2, mean) # simulate N sets of observations N = 1000; # number of sets of observations data.bl = cbind(rep(1,N),rep(3,N),rep(10,N)) # baseline data err = 1; # noise level to estimate err1 = matrix(rnorm(3*N,sd=0.1),N,3) # random noise data.obs = data.bl+err+err1; # strategy #1 func.1 = function(data){ mn = apply(data,2,mean) X=mn[1]; Y=mn[2]; Z=mn[3]; a1=(2/3)^2; b1=-2*2/3*(X-1/3*Y); a2=(9/10)^2; b2=-2*9/10*(X-1/10*Z); a=a1+a2; b=b1+b2; est.err = -b/(2*a); return(est.err) } func.1(data.obs) # strategy #2 Actually, it is the same with strategy #1 func.2 = function(data){ a1=(2/3)^2; a2=(9/10)^2; X=data[,1];Y=data[,2];Z=data[,3]; b1= -2*2/3*(X-1/3*Y); b2=-2*9/10*(X-1/10*Z); a=a1+a2; b=mean(b1)+mean(b2) est.err = -b/(2*a); return(est.err) } func.2(data.obs) source("https://bioconductor.org/biocLite.R") biocLite("viper") biocLite('mixtools') biocLite('bcellViper') library(viper) biocLite('mixtools') source("https://bioconductor.org/biocLite.R") source("https://bioconductor.org/biocLite.R") source("https://bioconductor.org/biocLite.R") sum(c(1:7)^2) rank.mat = matrix(data = c(3, 2.5, 3, 5, 6, 4, 2, 4, 2.5, 5, 5, 5, 4, 1, 7, 6, 7, 5, 1, 4, 3, 1.5, 2.5, 2, 1.5, 2, 4, 4, 1.5, 6, 6, 1.5, 4, 4, 5, 5.5, 6, 4, 5, 7, 4, 7, 5.5, 2.5, 1, 5, 3,4,6), nrow = 7, ncol = 7) rank.mat rank.mat = t(rank.mat) rank.mat feature.weight = c(7:1) feature.weight rank.mat2 <- (t(apply(rank.mat, 1, function(x) feature.weight * x))) rank.mat2 rmat.n <- apply(rank.mat2, 2, function(x) x/length(which(!is.na(x)))) rmat.n pep.sumrank <- rowSums(rmat.n, na.rm = TRUE) pep.sumrank av.rank <- pep.sumrank/sum(c(1:7)^2) ## worst av.rank = 1 av.rank <- 1 - av.rank ## reverse: best = 1 "weight" = high = reliable peptide av.rank <- av.rank^2 av.rank weight <- av.rank weight <- weight / sum(weight) trend <- apply(Mat, 2, function(x) weighted.mean(x, w=weight)) weight feature.weight = c(7,6,5,4,3,2,1)^2 feature.weight library(msnbase) source("http://bioconductor.org/biocLite.R") biocLite("Msnbase") biocLite("MSnbase") library(MSnbase) data(msnset2) head(exprs(msnset2)) debug(combineFeatures) prot = combineFeatures(msnset2, groupBy = fData(msnset2)$accession, fun = "iPQF") debug(iPQF) debug(iPQF.method) pos mat ffeature.weight features mat Mat rank.mat rank.mat2 49*4 feature.weight rmat.n 108/21.6 pep.sumrank 19.6+ 21.6+9.6+9+2.4+0.6+25 av.rank av.rank av.rank sum(av.rank) install.packages("protiq") library(protiq) library("protiq") library("protiq") install.packages("protiq") library("protiq") install.packages("graph") library("protiq") install.packages("graph") source("https://bioconductor.org/biocLite.R") biocLite("graph") library("protiq") biocLite("RBGL") library("protiq") ls data("leptoSRM") View(leptoSRMedgespp) View(leptoSRMpeptides) View(leptoSRMproteins) View(leptoSRMpeptides) dataChecked <- checkInputData(scampiData=scampi(peptides=leptoSRMpeptides, proteins=leptoSRMproteins, edgespp=leptoSRMedgespp), rescaling=FALSE) tmpPrepro <- preprocessInputData(scampiData=dataChecked) dataPrepro <- tmpPrepro[["dataPrepro"]] myCCList <- tmpPrepro[["ccList"]] rm(tmpPrepro) myCCList <- lapply(myCCList, getCovU, beta=0.2, tau=0.5) pepAbundanceScore <- quantifyPeptide(pepInfo=dataPrepro@peptides[13, c("pepId", "ccInd")], ccList=myCCList, param=c(alphaH=0, betaH=0.2, muH=0.3, tauH=0.2)) View(pepAbundanceScore) scampiRes <- quantifyProteins(scampiData=dataPrepro, ccList=myCCList, paramList=list(LSE=c(alphaH=0, betaH=0.2, muH=0.3, tauH=0.2)), quantifyPeptides=FALSE) names(scampiRes) scampiRes scampiRes@proteins View(leptoSRMpeptides) View(leptoSRMedgespp) View(leptoSRMproteins) protAbundanceScore <- quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId", "ccInd")], ccList=myCCList, param=c(alphaH=0, betaH=0.2, muH=0.3, tauH=0.2))[1] protAbundanceScore dataPrepro@proteins[7, c("protId", "ccInd")] debug(quantifyProtein) protAbundanceScore <- + quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId", "ccInd")], + ccList=myCCList, param=c(alphaH=0, betaH=0.2, muH=0.3, + tauH=0.2))[1] quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId", "ccInd")], ccList=myCCList, param=c(alphaH=0, betaH=0.2, muH=0.3, tauH=0.2))[1] cc$pepObs param["betaH"] Gcu scampiParam <- estimateModelParameters(method="all", ccList=myCCList, peptides=dataPrepro@peptides, numIter=10) scampiParam a = quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId","ccInd")], ccList = myCCList, param = c(alphaH = 5.88, betaH = 0.6, muH = 0, tauH = 0.48) ) a undebug(quantifyProtein) quantifyProtein a = quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId","ccInd")], ccList = myCCList, param = c(alphaH = 5.88, betaH = 0.6, muH = 0, tauH = 0.48) a a = quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId","ccInd")], ccList = myCCList, param = c(alphaH = 5.88, betaH = 0.6, muH = 0, tauH = 0.48)) a scampiParam <- estimateModelParameters(method="all", ccList=myCCList, peptides=dataPrepro@peptides, numIter=10) rm(list = ls()) data("leptoSRM") scampiOut <- runScampi(peptides=leptoSRMpeptides, proteins=leptoSRMproteins, edgespp=leptoSRMedgespp, rescaling=FALSE, method="LSE") scampiOut show(scampiOut) summary(scampiOut) scampiOut@proteins scampiOut@edgespp library(protiq) data("leptoSRM") scampiOut <- runScampi(peptides=leptoSRMpeptides, proteins=leptoSRMproteins, edgespp=leptoSRMedgespp, rescaling=FALSE, method="LSE") summary(scampiOut) library(protiq) data("leptoSRM") dataChecked <- checkInputData(scampiData=scampi(peptides=leptoSRMpeptides, proteins=leptoSRMproteins, edgespp=leptoSRMedgespp), rescaling = FALSE) tmpPrepro <- preprocessInputData(scampiData=dataChecked) dataPrepro <- tmpPrepro[["dataPrepro"]] myCCList <- tmpPrepro[["ccList"]] rm(tmpPrepro) myCCList <- lapply(myCCList, getCovU, beta=0.2, tau=0.5) debug(quantifyProtein) protAbundanceScore <- quantifyProtein(protInfo=dataPrepro@proteins[7, c("protId", "ccInd")], ccList=myCCList, param=c(alphaH=0, betaH=0.2, muH=0.3, tauH=0.2))[1] myCCList$distance ccList shiny::runApp('C:/Research/Projects/R_Shiny/jump_q') runApp('C:/Research/Projects/R_Shiny/jump_q') runApp('C:/Research/Projects/R_Shiny/jump_q') runApp('C:/Research/Projects/R_Shiny/jump_q') source('C:/Research/Projects/R_Shiny/jump_q/R/statTest.R') shiny::runApp('C:/Research/Projects/R_Shiny/jump_q') runApp('C:/Research/Projects/R_Shiny/jump_q') runApp('C:/Research/Projects/R_Shiny/jump_q') res res[1,] res[1,] res[1,] ggplot(data = res, aes(logfc, significance)) + geom_point() + geom_points() ggplot(data = res, aes(logfc, significance)) + geom_point() xlab ylab runApp('C:/Research/Projects/R_Shiny/jump_q') source('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmaply.R') install("d3heatmap") install.packages("d3heatmap") install.packages("d3heatmap") source('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmaply.R') library(d3heatmap) d3heatmap() d3heatmap(mtcars) d3heatmap(mtcars, scale = "column") source('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmaply.R') source('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmap.R') mtcars source('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmap.R') library(shiny); runApp('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmap.R') runApp('C:/Research/Projects/R_Shiny/jump_q/test_d3heatmap.R') runApp('C:/Research/Projects/R_Shiny/jump_q') runApp('C:/Research/Projects/R_Shiny/jump_q') class(data) class(mat) mat[1,] class(mtcars) mtcars mtcars[1,] runApp('C:/Research/Projects/R_Shiny/jump_q') mat[1,] d3heatmap(mat) d3heatmap(mat, labRow = NA) d3heatmap(mat, labRow = F) d3heatmap(mat, labRow = NULL) graphics.off() runApp('C:/Research/Projects/R_Shiny/jump_q') library(plotly) library(plotly) set.seed(100) d <- diamonds[sample(nrow(diamonds), 1000), ] plot_ly(d, x = carat, y = price, text = paste("Clarity: ", clarity), mode = "markers", color = carat, size = carat) source('~/test_plotly.R') d source('~/test_plotly.R') plot_ly carat d d$carat source('~/test_plotly.R') source('~/test_plotly.R') runApp('C:/Research/Projects/R_Shiny/jump_q') debugSource('C:/Research/Projects/SuffixArray/sa.R') seq seq[1] class(seq) substr(seq, 462590, 1) substr(seq, 462590, 2) substr(seq, 462590, 3) nchar(seq) substr(seq, 1, 3) substr(seq, 462590, 3) substr(seq, 462610, 3) substr(seq, 46000, 3) substr(seq, 1000, 1) substr(seq, 100, 1) substr(seq, 1, 1) substr(seq, 3, 1) substr(seq, 4, 1) substr(seq, 462590, 462600) loess() matrix(rnorm, ncol = 3, nrow = 2) matrix(rnorm(), ncol = 3, nrow = 2) rnorm(3, 2) rnorm(0, 1) rnorm(0, 1) rnorm(1, 1) rnorm(1, 0) rnorm(1, 0) rnorm(1, 0) rnorm(1, 0) rnorm(1, 0) rnorm(1, 0) rnorm(10) debugSource('C:/Research/Projects/11Collaborations/KunduLab/20181129/DEPs.R') debugSource('C:/Research/Projects/11Collaborations/KunduLab/20181129/DEPs.R') dir() cos(0.2) cos(c(0.1,0.2)) a = 0.2*pi*seq(0, 9) a cos(a) sum(cos(a)) seq(0, pi) seq(0, pi, n = 2501) linspace(0, pi, n = 2501) seq(0, pi, length.out= 2501) a = seq(0, pi, length.out= 2501) a[1:10] a[2499:2501] shiny::runApp('U:/Research/Projects/9Visualization/shiny/jump_q') runApp('U:/Research/Projects/9Visualization/shiny/jump_q') runApp('U:/Research/Projects/9Visualization/shiny/jump_q') qnorm() qnorm qnorm(0.95, 0, 1) qnorm(0.975, 0, 1) a = NULL for (i in 1:10000) {a = c(a, i)} a = list() for (i in 1:10000) {a = c(a, i)} a a = NULL for (i in 1:10000) {a = c(a, i)} a source('~/test.R') length(aP) ppois ppois(10, 1787 * 0.01) pnorm pnorm(1200, mean = 1020, sd = 50) 1 - pnorm(1200, mean = 1020, sd = 50) ppois ppois(3, lambda = 2.5 * 4) ppois(10, lambda = 1787 * 0.01) pbinom pmultnom pmultinom pmultinomial phyper(x, 13, 39, 5, lower.tail=TRUE) phyper(1, 13, 39, 5, lower.tail=TRUE) phyper(1, 13, 39, 5, lower.tail=FALSE) qhyper m = matrix(c(17, 21, 34, 5263), nrow = 2) m fisher.test(m) phyper(17, 34, 21, 5263) phyper(17, 34, 5350-21, 5263) phyper showMethods(phyper) phyper(17, 34, 5350-21, 5263) qhyper(17, 34, 5350-21, 5263) phyper(17,38,5335-38,51) 1- phyper(17,38,5335-38,51) dhyper(17,38,5335-38,51) dhyper(0,38,5335-38,51) dhyper(0:16,38,5335-38,51) sum(dhyper(0:16,38,5335-38,51)) double(sum(dhyper(0:16,38,5335-38,51))) float(sum(dhyper(0:16,38,5335-38,51))) 1 - sum(dhyper(0:16,38,5335-38,51)) a = sum(dhyper(0:16,38,5335-38,51)) a a = dhyper(0:16,38,5335-38,51) a sum(a) 1 - sum(a) cumsum(a) 1 - cumsum(a) a sum(a[1:5]) sum(a[1:6]) sum(a[1:7]) 1 - sum(a[1:7]) 1 - sum(a[1:8]) sum(a[1:7]) - sum(a[1:8]) sum(a[1:8]) - sum(a[1:7]) sum(a[1:6]) sum(a[1:7]) 1-sum(a[1:6]) 1 - sum(a[1:5]) 1 - sum(a[1]) 1 - sum(a[1:2]) 1 - sum(a[1:3]) 1 - sum(a[1:4]) 1 - sum(a[1:5]) 1 - sum(a[1:6]) 1 - sum(a[1:7]) a[7] dhyper(7,28,5573-28,51) dhyper(0:6,28,5573-28,51) sum(dhyper(0:6,28,5573-28,51)) 1-sum(dhyper(0:6,28,5573-28,51)) library(devtools) install_github("c-ruttkies/MetFragR/metfRag") install.packages("rJava") install_github("c-ruttkies/MetFragR/metfRag") setwd("U:/Research/Projects/7Metabolomics/JUMPm/R") dir() source('U:/Research/Projects/7Metabolomics/JUMPm/R/alignment.R') source('U:/Research/Projects/7Metabolomics/JUMPm/R/alignment.R') source('U:/Research/Projects/7Metabolomics/JUMPm/R/alignment.R')