annotate quantitation.r @ 0:bb199421f731 draft default tip

planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
author galaxyp
date Tue, 02 Oct 2018 16:30:33 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
1 #################################
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
2 library(tidyverse)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
3 library(furrr)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
4 library(lme4)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
5 library(MSnbase)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
6 library(MSqRob)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
7
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
8 ##Import and preprocess data
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
9 ############################
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
10 MSnSet2df = function(msnset){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
11 ## Converts Msnset to a tidy dataframe
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
12 ## Always creates feature and vector column so these shouldn't be defined by user.
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
13 ## convenient for downstream analysis steps.
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
14 if(any(c("sample", "feature", "expression") %in% c(colnames(fData(msnset)),colnames(pData(msnset))))){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
15 stop("Column names in the \"fData\" or \"pData\" slot of the \"msnset\" object cannot be named
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
16 \"sample\", \"feature\" or \"expression\". Please rename these columns before running the analysis.")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
17 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
18
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
19 dt <- as.data.frame(Biobase::exprs(msnset)) %>% mutate(feature = rownames(.)) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
20 gather(sample, expression, - feature, na.rm=TRUE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
21 dt <- fData(msnset) %>% mutate(feature = rownames(.)) %>% left_join(dt,. , by = 'feature')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
22 dt <- pData(msnset) %>% mutate(sample = rownames(.)) %>% left_join(dt,. , by = 'sample')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
23 as_data_frame(dt)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
24 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
25
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
26 ## robust summarisation
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
27 do_robust_summaristion = function(msnset, group_var = Proteins, keep_fData_cols = NULL, nIter = 20,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
28 sum_fun = summarizeRobust){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
29 ## TODO use funture_map instead of mutate to speed up
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
30 ## Uses assumption that featureNames and sampleNames exist in every msnset
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
31 ## Can also be used for multiple rounds of normalization, e.g. first from PSMs to peptides, then from peptides to proteins
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
32 system.time({## Time how long it takes
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
33 group_var <- enquo(group_var) ;#group_var = quo(Proteins)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
34 ## Make tidy dataframe from Msnset
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
35 df <- MSnSet2df(msnset)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
36 ## Do summarisision according defined groups
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
37 dt <- filter(df, !is.na(expression)) %>% group_by(!!group_var) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
38 mutate(expression = sum_fun(expression, feature, sample, nIter = nIter)) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
39 dplyr::select(!!group_var, sample, expression) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
40 ## collapse to one value per group
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
41 distinct
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
42 ## Construct an Msnset object from dataframe
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
43 dt_exprs <- spread(dt, sample, expression) %>% ungroup
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
44 exprs_data <- dplyr::select(dt_exprs, - !!group_var) %>% as.matrix
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
45 rownames(exprs_data) <- as.character(pull(dt_exprs, !!group_var))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
46
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
47 fd <- dplyr::select(dt_exprs,!!group_var)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
48
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
49 ##Select the group variable and all variables you want to keep
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
50 if (!is.null(keep_fData_cols)){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
51 fd_ext <- dplyr::select(df, !!group_var, one_of(keep_fData_cols)) %>% distinct
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
52 if(nrow(fd)!=nrow(fd_ext)){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
53 stop("Values in the \"group_var\" column can only correspond to a single value in the \"keep_fData_cols\" column.")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
54 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
55 fd <- left_join(fd,fd_ext)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
56 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
57
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
58 fd <- as.data.frame(fd)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
59 rownames(fd) <- as.character(pull(fd, !!group_var))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
60 out <- MSnSet(exprs_data, fData = AnnotatedDataFrame(fd) , pData = pData(msnset)[colnames(exprs_data),,drop = FALSE])
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
61 }) %>% print
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
62 out
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
63 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
64
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
65 summarizeRobust <- function(expression, feature, sample, nIter=100,...) {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
66 ## Assumes that intensities mx are already log-transformed
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
67 ## characters are faster to construct and work with then factors
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
68 feature <- as.character(feature)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
69 ##If there is only one 1 peptide for all samples return expression of that peptide
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
70 if (length(unique(feature)) == 1L) return(expression)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
71 sample <- as.character(sample)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
72 ## modelmatrix breaks on factors with 1 level so make vector of ones (will be intercept)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
73 if (length(unique(sample)) == 1L) sample = rep(1,length(sample))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
74
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
75 ## sum contrast on peptide level so sample effect will be mean over all peptides instead of reference level
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
76 X = model.matrix(~ -1 + sample + feature,contrasts.arg = list(feature = 'contr.sum'))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
77 ## MasS::rlm breaks on singulare values.
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
78 ## check with base lm if singular values are present.
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
79 ## if so, these coefficients will be zero, remove this collumn from model matrix
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
80 ## rinse and repeat on reduced modelmatrx till no singular values are present
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
81 repeat {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
82 fit = .lm.fit(X,expression)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
83 id = fit$coefficients != 0
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
84 X = X[ , id, drop =FALSE]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
85 if (!any(!id)) break
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
86 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
87 ## Last step is always rlm
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
88 fit = MASS::rlm(X,expression,maxit = nIter,...)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
89 ## Only return the estimated effects effects as summarised values
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
90 sampleid = seq_along(unique(sample))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
91 return(X[,sampleid,drop = FALSE] %*% fit$coefficients[sampleid])
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
92 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
93
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
94
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
95
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
96
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
97
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
98
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
99
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
100
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
101 ## mixed models
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
102 ###############
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
103 setGeneric (
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
104 name= "getBetaB",
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
105 def=function(model,...){standardGeneric("getBetaB")}
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
106 )
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
107
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
108 .getBetaBMermod = function(model) {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
109 betaB <- c(as.vector(getME(model,"beta")),as.vector(getME(model,"b")))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
110 names(betaB) <- c(colnames(getME(model,"X")),rownames(getME(model,"Zt")))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
111 betaB
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
112 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
113 setMethod("getBetaB", "lmerMod", .getBetaBMermod)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
114
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
115 .getBetaBGlm = function(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
116 model$coefficients
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
117
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
118 setGeneric (
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
119 name= "getVcovBetaBUnscaled",
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
120 def=function(model,...){standardGeneric("getVcovBetaBUnscaled")}
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
121 )
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
122
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
123 setMethod("getBetaB", "glm", .getBetaBGlm)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
124
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
125 .getVcovBetaBUnscaledMermod = function(model){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
126 ## TODO speed up (see code GAM4)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
127 p <- ncol(getME(model,"X"))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
128 q <- nrow(getME(model,"Zt"))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
129 Ct <- rbind2(t(getME(model,"X")),getME(model,"Zt"))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
130 Ginv <- solve(tcrossprod(getME(model,"Lambda"))+Diagonal(q,1e-18))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
131 vcovInv <- tcrossprod(Ct)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
132 vcovInv[((p+1):(q+p)),((p+1):(q+p))] <- vcovInv[((p+1):(q+p)),((p+1):(q+p))]+Ginv
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
133
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
134 solve(vcovInv)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
135 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
136
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
137 setMethod("getVcovBetaBUnscaled", "lmerMod", .getVcovBetaBUnscaledMermod)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
138
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
139 .getVcovBetaBUnscaledGlm = function(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
140 ## cov.scaled is scaled with the dispersion, "cov.scaled" is without the dispersion!
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
141 ## MSqRob::getSigma is needed because regular "sigma" function can return "NaN" when sigma is very small!
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
142 ## This might cause contrasts that can be estimated using summary() to be NA with our approach!
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
143 summary(model)$cov.scaled/MSqRob::getSigma(model)^2
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
144
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
145 setMethod("getVcovBetaBUnscaled", "glm", .getVcovBetaBUnscaledGlm)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
146
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
147 ## Estimate pvalues contrasts
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
148 contrast_helper = function(formula, msnset, contrast = NULL){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
149 ## Gives back the coefficients you can use to make contrasts with given the formula and dataset
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
150 ## If a factor variable is specified (that is present in the formula) all the possible contrasts
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
151 ## within this variable are returned
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
152 contrast <- enquo(contrast) ;#contrast = quo(condition)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
153 df <- MSnSet2df(msnset)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
154 all_vars <- formula %>% terms %>% delete.response %>% all.vars
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
155 names(all_vars) <- all_vars
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
156 df[,all_vars] <- map2_dfr(all_vars,df[,all_vars],paste0)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
157 coefficients <- c("(Intercept)", df %>% dplyr::select(all_vars) %>% unlist %>% unique %>% as.character)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
158 if (contrast != ~NULL) {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
159 c <- pull(df,!! contrast) %>% unique %>% sort %>% as.factor
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
160 comp <- combn(c,2,simplify = FALSE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
161 ## condIds = map(comp,~which( coefficients %in% .x))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
162 ## L = rep(0,length(coefficients))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
163 ## L = sapply(condIds,function(x){L[x]=c(-1,1);L})
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
164 ## rownames(L) = coefficients
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
165 ## colnames(L) = map_chr(comp, ~paste(.x,collapse = '-'))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
166 condIds <- map(comp, ~which(coefficients %in% .x))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
167 L <- rep(0,nlevels(c))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
168 L <- sapply(comp,function(x){L[x]=c(-1,1);L})
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
169 rownames(L) <- levels(c)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
170 colnames(L) <- map_chr(comp, ~paste(rev(.x),collapse = '-'))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
171 L
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
172 } else coefficients
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
173 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
174
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
175 setGeneric (
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
176 name= "getXLevels",
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
177 def=function(model,...){standardGeneric("getXLevels")}
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
178 )
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
179
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
180 .getXLevelsGlm = function(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
181 map2(names(model$xlevels), model$xlevels, paste0) %>% unlist
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
182
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
183 setMethod("getXLevels", "glm", .getXLevelsGlm)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
184
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
185 .getXLevelsMermod = function(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
186 getME(model,"flist") %>% map(levels) %>% unlist %>% unname
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
187
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
188 setMethod("getXLevels", "lmerMod", .getXLevelsMermod)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
189
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
190 contEst <- function(model, contrasts, var, df, lfc = 0){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
191 #TODO only contrast of random effect possible and not between fixed regression terms
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
192 betaB <- getBetaB(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
193 vcov <- getVcovBetaBUnscaled(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
194 coefficients <- names(betaB)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
195 id <- coefficients %in% rownames(contrasts)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
196
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
197 coefficients <- coefficients[id]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
198 vcov <- vcov[id,id]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
199 betaB <- betaB[id]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
200
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
201 xlevels <- getXLevels(model)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
202 id <- !apply(contrasts,2,function(x){any(x[!(rownames(contrasts) %in% xlevels)] !=0)})
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
203 contrasts <- contrasts[coefficients, id, drop = FALSE]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
204 ## If no contrasts could be found, terminate
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
205 if (is.null(colnames(contrasts))) return(new_tibble(list()))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
206
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
207 se <- sqrt(diag(t(contrasts)%*%vcov%*%contrasts)*var)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
208 logFC <- (t(contrasts)%*%betaB)[,1]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
209
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
210 ### Addition to allow testing against another log FC (lfc)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
211 ### See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2654802/
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
212
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
213 lfc <- abs(lfc)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
214 aest <- abs(logFC)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
215
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
216 Tval <- setNames(rep(0, length(logFC)),names(se))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
217 tstat.right <- (aest - lfc)/se
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
218 tstat.left <- (aest + lfc)/se
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
219 pval <- pt(tstat.right, df = df, lower.tail = FALSE) +
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
220 pt(tstat.left, df = df, lower.tail = FALSE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
221 tstat.right <- pmax(tstat.right, 0)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
222 fc.up <- (logFC >= lfc)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
223 fc.up[is.na(fc.up)] <- FALSE
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
224 fc.down <- (logFC < -lfc)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
225 fc.down[is.na(fc.down)] <- FALSE
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
226 Tval[fc.up] <- tstat.right[fc.up]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
227 Tval[fc.down] <- -tstat.right[fc.down]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
228 Tval[is.na(logFC)] <- NA
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
229
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
230 new_tibble(list(contrast = colnames(contrasts),
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
231 logFC = logFC,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
232 se = se,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
233 t = Tval,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
234 df = rep(df, length(se)),
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
235 pvalue = pval))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
236 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
237
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
238 do_lmerfit = function(df, form, nIter = 10, tol = 1e-6, control = lmerControl(calc.derivs = FALSE)){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
239 fit <- lmer(form, data = df, control = control)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
240 ##Initialize SSE
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
241 res <- resid(fit)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
242 ## sseOld=sum(res^2)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
243 sseOld <- fit@devcomp$cmp['pwrss']
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
244 while (nIter > 0){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
245 nIter = nIter-1
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
246 fit@frame$`(weights)` <- MASS::psi.huber(res/(mad(res)))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
247 fit <- refit(fit)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
248 res <- resid(fit)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
249 ## sse=sum(res^2)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
250 sse <- fit@devcomp$cmp['pwrss']
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
251 if(abs(sseOld-sse)/sseOld <= tol) break
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
252 sseOld <- sse
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
253 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
254 return(fit)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
255 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
256
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
257 calculate_df = function(df, model, vars){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
258 ## Get all the variables in the formula that are not defined in vars
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
259 form <- attributes(model@frame)$formula
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
260 vars_formula <- all.vars(form)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
261 vars_drop <- vars_formula[!vars_formula %in% vars]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
262 ## Sum of number of columns -1 of Zt mtrix of each random effect that does not involve a variable in vars_drop
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
263 mq <- getME(model,'q_i')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
264 id <- !map_lgl(names(mq),~{any(stringr::str_detect(.x,vars_drop))})
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
265 p <- sum(mq[id]) - sum(id)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
266 ## Sum of fixed effect parameters that do not involve a variable in vars_drop
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
267 mx <- getME(model,'X')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
268 id <- !map_lgl(colnames(mx),~{any(stringr::str_detect(.x,vars_drop))})
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
269 p <- p + sum(id)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
270
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
271 ## n is number of sample because 1 protein defined per sample
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
272 n <- n_distinct(df$sample)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
273 n-p
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
274 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
275
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
276 do_mm = function(formulas, msnset, group_vars = feature,type_df = 'traceHat'
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
277 , contrasts = NULL, lfc = 0, p.adjust.method = "BH", max_iter = 20L
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
278 , squeeze_variance = TRUE
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
279 , control = lmerControl(calc.derivs = FALSE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
280 ## choose parallel = plan(sequential) if you don't want parallelisation
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
281 ## , parallel_plan = plan(cluster, workers = makeClusterPSOCK(availableCores()))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
282 , parallel = TRUE, cores = availableCores()
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
283 ){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
284 if(!(type_df %in% c("conservative", "traceHat")))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
285 stop("Invalid input `type_df`.")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
286
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
287 system.time({## can take a while
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
288 if (parallel){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
289 cl <- makeClusterPSOCK(cores)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
290 plan(cluster, workers = cl)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
291 } else {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
292 plan(sequential)}
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
293
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
294 ## future::plan(parallel_plan,gc = TRUE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
295 formulas <- map(c(formulas), ~update(.,expression ~ . ))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
296 group_vars <- enquo(group_vars) # group_var = quo(protein)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
297 df <- MSnSet2df(msnset)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
298
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
299 ## Glm adds variable name to levels in catogorical (eg for contrast)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
300 ## lme4 doesnt do this for random effect, so add beforehand
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
301 ## Ludger needs this for Hurdle
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
302 df = formulas %>% map(lme4:::findbars) %>% unlist %>% map_chr(all.vars) %>% unique %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
303 purrr::reduce(~{mutate_at(.x,.y,funs(paste0(.y,.)))}, .init=df)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
304
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
305 cat("Fitting mixed models\n")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
306 ## select only columns needed for fitting
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
307 df_prot <- df %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
308 group_by_at(vars(!!group_vars)) %>% nest %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
309 mutate(model = furrr::future_map(data,~{
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
310 for (form in formulas){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
311 fit = try(do_lmerfit(.x, form, nIter = max_iter,control = control))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
312 if (class(fit) == "lmerMod") return(fit)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
313 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
314 fit
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
315 }))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
316
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
317 ## Return also failed ones afterward
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
318 df_prot_failed <- filter(df_prot, map_lgl(model,~{class(.x) != "lmerMod"}))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
319 df_prot <- filter(df_prot, map_lgl(model, ~{class(.x)=="lmerMod"}))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
320
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
321 if(nrow(df_prot) == 0) {print("No models could be fitted"); return(df_prot_failed)}
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
322
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
323 df_prot <-
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
324 mutate(df_prot
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
325 , formula = map(model,~{attributes(.@frame)$formula})
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
326 ## get trace hat df for squeezeVar
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
327 , df = map_dbl(model, ~getDf(.x))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
328 , sigma = map_dbl(model,~{MSqRob::getSigma(.x)})) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
329 ## Squeeze variance
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
330 bind_cols(as_data_frame(MSqRob::squeezeVarRob(.$sigma^2, .$df, robust = TRUE))) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
331 ## mutate(var_protein = ifelse(squeeze_variance,var.post,sigma^2),
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
332 mutate(var_protein = if (squeeze_variance) var.post else sigma^2,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
333 df_post = df + df.prior
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
334 , df_protein =
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
335 if (type_df == "conservative")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
336 ## Calculate df on protein level, assumption is that there is only one protein value/run,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
337 map2_dbl(data, model,~calculate_df(.x,.y, vars = colnames(pData(msnset))))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
338 else if (type_df == "traceHat")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
339 ## Alternative: MSqRob implementation with trace(Hat):
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
340 if(squeeze_variance) df_post else df
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
341 )
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
342
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
343 ## Calculate fold changes and p values for contrast
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
344 cat("Estimating p-values contrasts\n")
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
345 df_prot <- df_prot %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
346 mutate(contrasts = furrr::future_pmap(list(model = model, contrasts = list(contrasts),
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
347 var = var_protein,
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
348 df = df_protein, lfc = lfc), contEst)) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
349 ## Calculate qvalues BH
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
350 select_at(vars(!!group_vars, contrasts)) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
351 unnest %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
352 group_by(contrast) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
353 mutate(qvalue = p.adjust(pvalue, method = p.adjust.method)) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
354 group_by_at(vars(!!group_vars)) %>% nest(.key = contrasts) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
355 left_join(df_prot,.)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
356 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
357 ) %>% print
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
358 if (parallel) stopCluster(cl)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
359 bind_rows(df_prot,df_prot_failed)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
360 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
361
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
362 read_moff = function(moff,meta){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
363 print('START READING MOFF DATA')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
364 set = readMSnSet2(moff, ecol = -c(1,2),fnames = 'peptide',
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
365 sep = '\t',stringsAsFactors = FALSE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
366 colnames(fData(set)) = c('peptide','protein')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
367 pd = read_tsv(meta) %>% column_to_rownames('sample') %>% as.data.frame
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
368
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
369 ## fix msnbase bug 1
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
370 ## if there is only 1 sample. Msnbase doesn't name it
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
371 if (length(sampleNames(set) ==1))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
372 sampleNames(set) = rownames(pd)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
373
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
374 pData(set) = pd
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
375 ## fix msnbase bug 2
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
376 ## bug in msnbase in summarisation (samplenames should be alphabetically)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
377 sample_order = order(sampleNames(set))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
378 set = MSnSet(exprs(set)[,sample_order,drop = FALSE]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
379 , fData = AnnotatedDataFrame(fData(set))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
380 , pData = AnnotatedDataFrame(pData(set)[sample_order,,drop = FALSE]))
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
381 print('END READING MOFF DATA')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
382 set
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
383 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
384
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
385 preprocess = function(set){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
386 print('START PREPROCESSING')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
387 if (ncol(set) == 1){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
388 exprs(set)[0 == (exprs(set))] <- NA
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
389 set = log(set, base = 2)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
390 ## keep smallest unique groups
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
391 groups2 <- smallestUniqueGroups(fData(set)$protein,split = ',')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
392 sel <- fData(set)$protein %in% groups2
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
393 set <- set[sel,]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
394 } else {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
395 ## normalisation
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
396 exprs(set)[0 == (exprs(set))] <- NA
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
397 set <- normalize(set, 'vsn')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
398 ## keep smallest unique groups
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
399 groups2 <- smallestUniqueGroups(fData(set)$protein,split = ',')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
400 sel <- fData(set)$protein %in% groups2
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
401 set <- set[sel,]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
402 ## remove peptides with less then 2 observations
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
403 sel <- rowSums(!is.na(exprs(set))) >= 2
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
404 set <- set[sel]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
405 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
406 print('END PREPROCESSING')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
407 set
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
408 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
409
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
410 summarise = function(set){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
411 print('START SUMMARISATION')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
412 ## Summarisation
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
413 if (ncol(set) == 1){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
414 set = combineFeatures(set,fun="median", groupBy = fData(set)$protein,cv = FALSE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
415 } else {
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
416 ## set = combineFeatures(set,fun="robust", groupBy = fData(set)$protein,cv = FALSE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
417 set = do_robust_summaristion(set,protein)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
418 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
419 exprs(set) %>% as.data.frame %>% rownames_to_column('protein') %>% write_tsv('summarised_proteins.tsv')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
420 print('END SUMMARISATION')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
421 set
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
422 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
423
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
424 quantify = function(set, cpu = 0){
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
425 print('START QUANTITATION')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
426 if ((cpu == 0) | is.na(cpu)) cpu = availableCores()
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
427 print(cpu)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
428 form = colnames(pData(set)) %>% paste0('(1|',.,')',collapse='+') %>% paste('~',.) %>% as.formula
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
429 contrasts <- contrast_helper(form, set, condition)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
430 res = do_mm(formulas = form, set, group_vars = c(protein)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
431 , contrasts = contrasts,type_df = 'traceHat', parallel = TRUE,cores = cpu) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
432 filter(!map_lgl(contrasts, is.null)) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
433 transmute(protein, contrasts) %>% unnest %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
434 transmute(protein
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
435 , comparison = str_replace_all(contrast, 'condition', '')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
436 , logFC,pvalue,qvalue) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
437 write_tsv('quantitation.tsv')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
438 print('END QUANTITATION')
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
439 }
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
440
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
441 args = commandArgs(trailingOnly=TRUE)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
442 moff = args[1]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
443 meta = args[2]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
444 summarise_only = args[3]
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
445 cpu = strtoi(args[4])
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
446
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
447 res = read_moff(moff, meta) %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
448 preprocess %>%
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
449 summarise
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
450 if (summarise_only != 1)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
451 quantify(res, cpu)
bb199421f731 planemo upload for repository https://github.com/galaxyproteomics/tools-galaxyp/tree/master/tools/lfq_protein_quant commit 26ff08776f90f96646598a19cfcf57d42aa4a43b
galaxyp
parents:
diff changeset
452