This repository contains the R code and package for the mi4p methodology (Multiple Imputation for Proteomics), proposed by Marie Chion, Christine Carapito and Frédéric Bertrand (2021) in Accounting for multiple imputation-induced variability for differential analysis in mass spectrometry-based label-free quantitative proteomics, https://arxiv.org/abs/2108.07086.
The following material is available on the Github repository of the package https://github.com/mariechion/mi4p/.
The Functions folder contains all the functions used
for the workflow.
The Simulation-1, Simulation-2 and
Simulation-3 folders contain all the R scripts and data
used to conduct simulated experiments and evaluate our
methodology.
The Arabidopsis_UPS and Yeast_UPS
folders contain all the R scripts and data used to challenge our
methodology on real proteomics datasets. Raw experimental data were
deposited with the ProteomeXchange Consortium via the PRIDE partner
repository with the dataset identifiers PXD003841 and
PXD027800.
This website and these examples were created by M. Chion, C. Carapito and F. Bertrand.
You can install the released version of mi4p from CRAN with:
install.packages("mi4p")You can install the development version of mi4p from github with:
devtools::install_github("mariechion/mi4p")library(mi4p)set.seed(4619)
datasim <- protdatasim()
str(datasim)
#> 'data.frame': 200 obs. of 11 variables:
#> $ id.obs: int 1 2 3 4 5 6 7 8 9 10 ...
#> $ X1 : num 99.6 99.9 100.2 99.8 100.4 ...
#> $ X2 : num 97.4 101.3 100.3 100.2 101.7 ...
#> $ X3 : num 100.3 100.9 99.1 101.2 100.6 ...
#> $ X4 : num 99.4 99.2 98.5 99.1 99.5 ...
#> $ X5 : num 98.5 99.7 100 100.2 100.7 ...
#> $ X6 : num 200 199 199 200 199 ...
#> $ X7 : num 200 200 202 199 199 ...
#> $ X8 : num 202 199 200 199 201 ...
#> $ X9 : num 200 200 199 201 200 ...
#> $ X10 : num 200 198 200 201 199 ...
#> - attr(*, "metadata")='data.frame': 10 obs. of 3 variables:
#> ..$ Sample.name: chr [1:10] "X1" "X2" "X3" "X4" ...
#> ..$ Condition : Factor w/ 2 levels "A","B": 1 1 1 1 1 2 2 2 2 2
#> ..$ Bio.Rep : int [1:10] 1 2 3 4 5 6 7 8 9 10It is the dataset shipped with package.
save(datasim, file="datasim.RData", compress = "xz")attr(datasim, "metadata")
#> Sample.name Condition Bio.Rep
#> 1 X1 A 1
#> 2 X2 A 2
#> 3 X3 A 3
#> 4 X4 A 4
#> 5 X5 A 5
#> 6 X6 B 6
#> 7 X7 B 7
#> 8 X8 B 8
#> 9 X9 B 9
#> 10 X10 B 10MV1pct.NA.data <- MVgen(dataset = datasim[,-1], prop_NA = 0.01)
MV1pct.NA.data[1:6,]
#> X1 X2 X3 X4 X5 X6 X7
#> 1 99.62136 97.3689 100.28075 99.37889 98.48006 200.2803 199.8042
#> 2 99.86342 101.2731 100.86058 99.19046 99.71866 199.2618 200.3934
#> 3 100.21425 100.3024 99.07865 98.51259 100.02675 199.2197 202.2855
#> 4 99.84381 100.2009 101.18119 99.05420 100.18467 200.0615 199.3424
#> 5 100.36412 101.7019 100.59543 99.49536 100.65887 199.2520 199.4715
#> 6 100.70214 NA 101.65048 99.11372 100.94608 201.7471 198.7024
#> X8 X9 X10
#> 1 201.6095 199.7451 199.7543
#> 2 199.4827 199.7297 198.3984
#> 3 199.7161 198.8883 200.4094
#> 4 199.1865 200.9596 201.2789
#> 5 200.9682 199.9838 198.7303
#> 6 201.8030 199.4030 199.9033MV1pct.impMLE <- multi.impute(data = MV1pct.NA.data, conditions = attr(datasim,"metadata")$Condition, method = "MLE", parallel = FALSE)print(paste(Sys.time(), "Dataset", 1, "out of", 1))
#> [1] "2024-09-30 23:18:29.427335 Dataset 1 out of 1"
MV1pct.impMLE.VarRubin.Mat <- rubin2.all(data = MV1pct.impMLE, metacond = attr(datasim, "metadata")$Condition) print(paste("Dataset", 1, "out of",1, Sys.time()))
#> [1] "Dataset 1 out of 1 2024-09-30 23:19:46.728503"
MV1pct.impMLE.VarRubin.S2 <- as.numeric(lapply(MV1pct.impMLE.VarRubin.Mat, function(aaa){
DesMat = mi4p::make.design(attr(datasim, "metadata"))
return(max(diag(aaa)%*%t(DesMat)%*%DesMat))
}))MV1pct.impMLE.mi4limma.res <- mi4limma(qData = apply(MV1pct.impMLE,1:2,mean),
sTab = attr(datasim, "metadata"),
VarRubin = sqrt(MV1pct.impMLE.VarRubin.S2))
rapply(MV1pct.impMLE.mi4limma.res,head)
#> logFC.A_vs_B_logFC1 logFC.A_vs_B_logFC2 logFC.A_vs_B_logFC3
#> -101.21269 -99.27197 -100.47686
#> logFC.A_vs_B_logFC4 logFC.A_vs_B_logFC5 logFC.A_vs_B_logFC6
#> -100.07284 -99.11801 -99.82359
#> P_Value.A_vs_B_pval1 P_Value.A_vs_B_pval2 P_Value.A_vs_B_pval3
#> 0.00000 0.00000 0.00000
#> P_Value.A_vs_B_pval4 P_Value.A_vs_B_pval5 P_Value.A_vs_B_pval6
#> 0.00000 0.00000 0.00000
(simplify2array(MV1pct.impMLE.mi4limma.res)$P_Value.A_vs_B_pval)[1:10]
#> [1] 0 0 0 0 0 0 0 0 0 0
(simplify2array(MV1pct.impMLE.mi4limma.res)$P_Value.A_vs_B_pval)[11:200]<=0.05
#> [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE
#> [12] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [23] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [34] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [45] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE
#> [56] FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
#> [67] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
#> [78] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [89] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [100] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [111] TRUE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
#> [122] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [133] FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE
#> [144] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [155] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [166] FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [177] FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
#> [188] FALSE FALSE FALSETrue positive rate
sum((simplify2array(MV1pct.impMLE.mi4limma.res)$P_Value.A_vs_B_pval)[1:10]<=0.05)/10
#> [1] 1False positive rate
sum((simplify2array(MV1pct.impMLE.mi4limma.res)$P_Value.A_vs_B_pval)[11:200]<=0.05)/190
#> [1] 0.05789474MV1pct.impMLE.dapar.res <-limmaCompleteTest.mod(qData = apply(MV1pct.impMLE,1:2,mean), sTab = attr(datasim, "metadata"))
rapply(MV1pct.impMLE.dapar.res,head)
#> res.l.logFC.A_vs_B_logFC1 res.l.logFC.A_vs_B_logFC2
#> -101.212692 -99.271966
#> res.l.logFC.A_vs_B_logFC3 res.l.logFC.A_vs_B_logFC4
#> -100.476863 -100.072843
#> res.l.logFC.A_vs_B_logFC5 res.l.logFC.A_vs_B_logFC6
#> -99.118007 -99.823590
#> res.l.P_Value.A_vs_B_pval1 res.l.P_Value.A_vs_B_pval2
#> 0.000000 0.000000
#> res.l.P_Value.A_vs_B_pval3 res.l.P_Value.A_vs_B_pval4
#> 0.000000 0.000000
#> res.l.P_Value.A_vs_B_pval5 res.l.P_Value.A_vs_B_pval6
#> 0.000000 0.000000
#> fit.s21 fit.s22
#> 1.017841 1.017841
#> fit.s23 fit.s24
#> 1.017841 1.017841
#> fit.s25 fit.s26
#> 1.017841 1.017841Simulate a list of 100 datasets.
set.seed(4619)
norm.200.m100.sd1.vs.m200.sd1.list <- lapply(1:100, protdatasim)
metadata <- attr(norm.200.m100.sd1.vs.m200.sd1.list[[1]],"metadata")It is the list of dataset shipped with package.
save(norm.200.m100.sd1.vs.m200.sd1.list, file="norm.200.m100.sd1.vs.m200.sd1.list.RData", compress = "xz")100 datasets with parallel comuting support. Quite long to run even with parallel computing support.
library(foreach)
doParallel::registerDoParallel(cores=NULL)
requireNamespace("foreach",quietly = TRUE)MV1pct.NA.data <- foreach::foreach(iforeach = norm.200.m100.sd1.vs.m200.sd1.list,
.errorhandling = 'stop', .verbose = T) %dopar%
MVgen(dataset = iforeach[,-1], prop_NA = 0.01)MV1pct.impMLE <- foreach::foreach(iforeach = MV1pct.NA.data,
.errorhandling = 'stop', .verbose = F) %dopar%
multi.impute(data = iforeach, conditions = metadata$Condition,
method = "MLE", parallel = F)MV1pct.impMLE.VarRubin.Mat <- lapply(1:length(MV1pct.impMLE), function(index){
print(paste(Sys.time(), "Dataset", index, "out of", length(MV1pct.impMLE)))
rubin2.all(data = MV1pct.impMLE[[index]], metacond = metadata$Condition)
})MV1pct.impMLE.VarRubin.S2 <- lapply(1:length(MV1pct.impMLE.VarRubin.Mat), function(id.dataset){
print(paste("Dataset", id.dataset, "out of",length(MV1pct.impMLE.VarRubin.Mat), Sys.time()))
as.numeric(lapply(MV1pct.impMLE.VarRubin.Mat[[id.dataset]], function(aaa){
DesMat = mi4p::make.design(metadata)
return(max(diag(aaa)%*%t(DesMat)%*%DesMat))
}))
})MV1pct.impMLE.mi4limma.res <- foreach(iforeach = 1:100, .errorhandling = 'stop', .verbose = T) %dopar%
mi4limma(qData = apply(MV1pct.impMLE[[iforeach]],1:2,mean),
sTab = metadata,
VarRubin = sqrt(MV1pct.impMLE.VarRubin.S2[[iforeach]]))
MV1pct.impMLE.dapar.res <- foreach(iforeach = 1:100, .errorhandling = 'stop', .verbose = T) %dopar%
limmaCompleteTest.mod(qData = apply(MV1pct.impMLE[[iforeach]],1:2,mean),
sTab = metadata)Complimentary useful tests
The g.test function of theProteoMM
Bioconductor package, implements the G-Test described in “A statistical
framework for protein quantitation in bottom-up MS based proteomics``
(Karpievitch et al. Bioinformatics 2009). For some experimental designs
of experiments, this test may be used to look for significant peptides
based on their absence/presence. For some designs, it will decrease the
precision of out methodology, see the arabidopsis example on github.
library(ProteoMM)
ProteoMM::g.test(c(TRUE, TRUE, FALSE, FALSE), as.factor(c('grp1', 'grp1', 'grp2', 'grp2')))
#>
#> Log likelihood ratio/G test of independence without correction
#>
#> data: c(TRUE, TRUE, FALSE, FALSE) and as.factor(c("grp1", "grp1", "grp2", "grp2"))
#> Log likelihood ratio statistic (G) = 5.5452, X-squared df = 1,
#> p-value = 0.01853
data("qData")
data("sTab")
tableNA.qData <- apply(is.na(qData),1,table,sTab$Condition)
id.mix <- unlist(lapply(tableNA.qData,function(res) nrow(res)>1))
# apply(is.na(qData[id.mix,]),1,g.test,sTab$Condition)
res.g.test <- cbind(rownames=as.data.frame(rownames(qData)[id.mix]),
p.val=apply(is.na(qData[id.mix,]),1,
function(tab) return(ProteoMM::g.test(x=tab,y=sTab$Condition)$p.value)))
res.g.test[res.g.test[,2]<0.05,]
#> rownames(qData)[id.mix] p.val
#> 29 29 0.003925917
#> 44 44 0.003925917
#> 105 105 0.003925917
#> 389 389 0.003925917
qData[rownames(res.g.test[res.g.test[,2]<0.05,]),]
#> Intensity_C_R1 Intensity_C_R2 Intensity_C_R3 Intensity_D_R1
#> 29 23.18826 22.74154 22.94633 NA
#> 44 NA NA NA 25.83019
#> 105 26.52709 26.37171 26.63201 NA
#> 389 NA NA NA 22.20932
#> Intensity_D_R2 Intensity_D_R3
#> 29 NA NA
#> 44 25.98796 26.03095
#> 105 NA NA
#> 389 21.93693 22.03093The eigen_pi function of the ProteoMM
Bioconductor package computes the proportion of observations missing
completely at random. It is used by the g.test function if
such an estimate is to be computed using the data .
library(ProteoMM)
data(mm_peptides)
intsCols = 8:13
metaCols = 1:7
m_Ints = mm_peptides[, intsCols]
m_prot.info = mm_peptides[, metaCols]
m_logInts = m_Ints
m_logInts[m_Ints==0] = NA
m_logInts = log2(m_logInts)
my.pi = ProteoMM::eigen_pi(m_logInts, toplot=TRUE)
plot of chunk unnamed-chunk-24