R/evalfq.R
In JianboFu0406/EVALFQ111: R Package for Evaluating Label-Free Proteome Quantification
Defines functions
lfqevalueall
Documented in
lfqevalueall
#' @title Conduct LFQ and assess performance of all possible LFQ workflows.
#' @description The ProteoLFQ enables the label-free quantification of proteomic
#' data and the performance assessment of each LFQ workflow from multiple perspectives. Moreover, it
#' provides the unique function of ranking all possible LFQ workflows (>3,000 random combinations of
#' transformation, normalization and imputation methods) based on their performances.
#' All in all, this tool makes the performance assessment of whole LFQ workflow possible
#' (collectively assessed by five well-established criteria with distinct underlying theories) and
#' gives the ranking results of all possible workflows based on the criteria preferred and selected by the users.
#' For function definitions and descriptions please use "??ProteoLFQ" command in R.
#' @param data_q This input file should be numeric type except the first and second column containing the names and label (control or case) of the studied samples, respectively. The intensity data should be provided in this input file with the following order: samples in row and proteins/peptides in column. Missing value (NA) of protein intensity are allowed.
#' @param assum_a all proteins were assumed to be equally important.The authors will be asked to input a letter “Y” to indicate the corresponding assumption is held for the studied dataset and a letter “N” to denote the opposite.
#' @param assum_b The level of protein abundance was assumed to be constant among all samples. The authors will be asked to input a letter “Y” to indicate the corresponding assumption is held for the studied dataset and a letter “N” to denote the opposite.
#' @param assum_c The intensities of the vast majority of the proteins were assumed to be unchanged under the studied conditions. The authors will be asked to input a letter “Y” to indicate the corresponding assumption is held for the studied dataset and a letter “N” to denote the opposite.
#' @param Ca Criterion (a): precision of LFQ based on the proteomes among replicates (Proteomics. 15:3140-51, 2015). If set 1, the user chooses to assess LFQ workflows using Criterion (a). If set 0, the user excludes Criterion (a) from performance assessment. The default setting of this value is “1”.
#' @param Cb Criterion (b): classification ability of LFQ between distinct sample groups (Nat Biotechnol. 28:83-9, 2010). If set 1, the user chooses to assess LFQ workflows using Criterion (b). If set 0, the user excludes Criterion (b) from performance assessment. The default setting of this value is “1”.
#' @param Cc Criterion (c): differential expression analysis by reproducibility-optimization (Nat Biotechnol. 32:896-902, 2014). If set 1, the user chooses to assess LFQ workflows using Criterion (c). If set 0, the user excludes Criterion (c) from performance assessment. The default setting of this value is “1”.
#' @param Cd Criterion (d): reproducibility of the identified protein markers among different datasets (Mol Biosyst. 11:1235-40, 2015). If set 1, the user chooses to assess LFQ workflows using Criterion (d). If set 0, the user excludes Criterion (d) from performance assessment. The default setting of this value is “1”.
#' @return preprocessed matrix
#' @import utils stats
#' @import metabolomics
#' @import affy vsn
#' @import MASS limma
#' @import ProteoMM ROTS
#' @importFrom grDevices colorRampPalette dev.off pdf
#' @useDynLib EVALFQ
#' @importFrom Rcpp sourceCpp
#' @rawNamespace import(dplyr, except=c(filter,lag,select,combine))
#' @rawNamespace import(gplots, except=lowess)
#' @importFrom pcaMethods pca
#' @importFrom pcaMethods completeObs
#' @import impute
#' @usage lfqevalueall(data_q,
#' assum_a="Y", assum_b="Y", assum_c="Y", Ca="1", Cb="1", Cc="1", Cd="1")
#' @export lfqevalueall
lfqevalueall <- function(data_q, assum_a="Y", assum_b="Y", assum_c="Y", Ca="1", Cb="1", Cc="1", Cd="1"){
#### 样本在行，特征在列
trans <- function(data,n){
matrix <- switch(
n,
Box_Cox(data),
log2(data),
data
)
return(matrix)
}
cents <- function(data,n){
matrix <- switch(
n,
MEC(data),
MDC(data),
data
)
return(matrix)
}
scals<- function(data,n){
matrix <- switch(
n,
1,
AUTO1(data),
PARETO1(data),
VAST1(data),
RANGE1(data)
)
return(matrix)
}
norm <- function(data,n){
matrix <- switch(
n,
t(fastlo(as.matrix(data))),
t(EIGENMS(data, label)),
t(LOWESS(data)),
t(SMAD(data)),
t(MEAN(data)),
t(MEDIAN(data)),
t(data),
t(PQN(data)),
t(QUANTILE(as.matrix(data))),
t(RLR1(data)),
t(MSTUS(data)),
t(TMM(data)),
t(VSN(as.matrix(data)))
)
return(matrix)
}
impute <- function(data,n){
matrix <- switch(
n,
filter_train_data,
t(back(filter_train_data)),
t(bpca(filter_train_data,nPcs=3)),
t(censor(filter_train_data)),
t(knn(filter_train_data,k=10)),
t(svdm(filter_train_data,nPcs=3)),
t(zero(filter_train_data))
)
return(matrix)
}
consistency <- function(fold = 5, top = 20) {
folds <- fold
control.label <- control.y # variable-1
test.fold1 <- split(sample(1:length(control.label)), 1:folds) #ignore warning
case.label <- case.y # variable-2
test.fold2 <- split(sample(1:length(case.label)), 1:folds) #ignore warning
DEG <- list()
for (i in 1:folds) {
com.x <- cbind(control.x[, test.fold1[[i]]], case.x[, test.fold2[[i]]]) # variable-3 & 4.
lab.ct <- test.fold1[[i]]
lab.ca <- test.fold2[[i]]
design <- cbind(Grp1 = 1, Grp2vs1 = c(rep(0, length(lab.ct)), rep(1, length(lab.ca))))
fit <- limma::lmFit(com.x, design)
fit <- limma::eBayes(fit)
DEG[[i]] <- rownames(limma::topTable(fit, coef = 2, number = nrow(com.x)))
}
names(DEG) <- LETTERS[1:folds]
top.n <-top # Extracting the top n genes.
DEG.list <- DEG
for (g in 1:length(DEG.list)) {
DEG.list[[g]] <- DEG.list[[g]][1:top.n]
}
# Calculating consistency score:
setlist <- DEG.list
OLlist <- overLapper(setlist=setlist, sep="", type="vennsets")
con.score <- 0
VennList <- OLlist$Venn_List
for (i in 1:length(VennList)) {
insect.n <- nchar(names(VennList[i]))
if (insect.n < 2) next
num.i <- 2^(insect.n - 2) * length(VennList[[i]])
con.score <- con.score + num.i
}
return(con.score) # consistense score
}
# Stable consistense score, with 20 repeats
stabel.score <- function(repeats = 20, fold = 5, top = 10) {
score <- 0
for (r in 1:repeats) {
score <- score + consistency(fold, top)
}
return(score/repeats)
}
#####################################################################
iName<-c("BOX","LOG","NON")
oName<-c("MEC","MDC","NON")
pName<-c("NON","ATO","PAR","RAN","VAS")
jName<-c("CYC","EIG","LOW","MAD","MEA","MED","NON","PQN","QUA","RLR","TIC","TMM","VSN")
gName<-c("NON","BAK","BPC","CEN","KNN","SVD","ZER")
dataa<-data_q
rownames(dataa)<-dataa[,1]
label <- dataa[,2]
frame <- dataa[, -(1:2)]
frame <- t(frame)
frame <- data.frame(frame)
train_data <- data.matrix(frame, rownames.force = NA)
train_data_t <- train_data
Fpcv <- list()
Fscore <- list()
Faccuracy <- list()
Fpmad <- list()
Fbar1num <- list()
Fbaronummean <- list()
Fbarosd <- list()
Fbarorsd <- list()
Fbaro_rsdr_to_bar1num <- list()
spike <- list()
backgound <- list()
time <- 0
center <- c(1,2,3)
scaling <- c(1,2,3,4,5)
normalization <- c(1,2,3,4,5,6,7,8,9,10,11,12)
cat("Assumption A: all proteins were equally important (Y/N): ", assum_a, "\n")
cat("Assumption B: the level of protein abundance was constant among all samples (Y/N): ", assum_b, "\n")
cat("Assumption C: the intensity of the majority of proteins were unchanged (Y/N): ", assum_c, "\n")
if(is.na(match("Y",assum_a))){
center<-3
scaling<-1
}
if(is.na(match("Y",assum_b)) & is.na(match("Y",assum_c))) {
normalization <- c(2,7)
}
if(is.na(match("Y",assum_b)) & !is.na(match("Y",assum_c))){
normalization <- c(1,2,7,8,9,10,11)
}
if(!is.na(match("Y",assum_b)) & is.na(match("Y",assum_c))){
normalization <- c(1,4,5,6,7,11)
}
for(i in 1:3){
train_data_tran <- try(trans(train_data_t,i))
if (inherits(train_data_tran, "try-error"))
{ next }
if(i!=3){
for(o in center){
tran_train_data <- train_data_tran
tran_train_data[is.infinite(data.matrix(tran_train_data))] <- NA
cen_train_data <- try(cents(tran_train_data,o))
if (inherits(cen_train_data, "try-error"))
{ next }
for(p in scaling){
scal_factor <- try(scals(tran_train_data,p))
if (inherits(scal_factor, "try-error"))
#if(class(scal_factor) == "try-error")
{ next }
scal_train_data <- cen_train_data/scal_factor
if (inherits(scal_train_data, "try-error"))
#if(class(scal_train_data)=="try-error")
{ next }
for(j in normalization){
scal_train_data[is.nan(scal_train_data)] <- NA
scal_train_data[is.infinite(scal_train_data)] <- NA
normalized_data <- try(norm(scal_train_data,j))
if (inherits(normalized_data, "try-error"))
#if(class(normalized_data)=="try-error")
{ next }
label_c <- as.factor(label)
g1 <- table(label_c)[levels(label_c)[1]]*0.8
g2 <- table(label_c)[levels(label_c)[2]]*0.8
train_data_filtering <- try(Basicfilter(normalized_data,label,g1=2,g2=2))
if (inherits(train_data_filtering, "try-error"))
#if(class(train_data_filtering) == "try-error")
{ next }
filter_train_data <- train_data_filtering
for(g in 1:7){
if(g==1){
imputed_data <- filter_train_data
}
if(g!=1){
imputed_data <- try(impute(filter_train_data,g))
if (inherits(imputed_data, "try-error"))
#if(class(imputed_data)=="try-error")
{ next }
}
##### Feature Selection
time = time+1
dataa <- data_q
frame <- imputed_data
label <- dataa[,2]
rots.out <-try(ROTS(data = t(frame), groups = as.character(label), B = 200, K = 500 , seed = 1234,log = FALSE))
if (inherits(rots.out, "try-error"))
#if(class(rots.out)=="try-error")
{ next }
frame <- imputed_data
label <- as.factor(as.character(label))
im.data <- data.frame(label=label,frame)
#(a) Precision of LFQ Based on the Proteomes among Replicates
if( Ca == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria A: Precision"),"\n")
data <- im.data
result <- PCV1(data)
pcv <- sapply(1:3, function(i){round(1000*mean(as.numeric(result[,i])))/1000})[3]
pmad <- try(mean(PMAD(data)))
if (inherits(pmad, "try-error"))
#if(class(pmad) == "try-error")
{ next }
Fpmad[paste(i,o,p,j,g,sep="")] <- pmad
}else{
message("'Criteria A: Precision' cannot be evaluated, Please Check!")
}
#(b) Classification Ability of LFQ between Distinct Sample Groups
if( Cb == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria B: Classification.Ability"),"\n")
data <- im.data
rots.out <- rots.out
col_pos <- which(rots.out$FDR<0.05)
if(length(col_pos) <= 10) {
markerid <-order(rots.out$pvalue)[1:20]
}else {
markerid <- which(rots.out$FDR<0.05)+1
}
clusters <- hclust(dist(data[,markerid]))
clusterCut <- cutree(clusters, 2)
dataa <- data_q
label <- dataa[,2]
tmatrix <- table(clusterCut, label)
tru <- as.numeric(data[,1])
accuracy <- (tmatrix[1,1]+ tmatrix[2,2])/length(label)
Faccuracy[paste(i,o,p,j,g,sep="")] <- accuracy
}else{
message("'Criteria B: Classification.Ability' cannot be evaluated, Please Check!")
}
#(c) Differential Expression Analysis Based on Reproducibility-optimization
if( Cc == 1 ){
cat(paste("Assessing" ,paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria C: Differential.Expression"),"\n")
rots.out <- rots.out
breaks <- seq(0,1,0.05)
sdres <- affy::hist(rots.out$pvalue,breaks=breaks)
bar1num <- sdres$counts[1]
baronummean <- mean(sdres$counts[-1])
barosd <- sd(sdres$counts[-1])
barorsd <- barosd/baronummean
Fbar1num[paste(i,o,p,j,g,sep="")] <- bar1num
Fbaronummean[paste(i,o,p,j,g,sep="")] <- baronummean
Fbarosd[paste(i,o,p,j,g,sep="")] <- barosd
Fbarorsd[paste(i,o,p,j,g,sep="")] <- barorsd
Fbaro_rsdr_to_bar1num[paste(i,o,p,j,g,sep="")] <- barorsd/bar1num
}else{
message("'Criteria C: Differential.Expression' cannot be evaluated, Please Check!")
}
#(d) Reproducibility of the Identified Protein Markers among Different Datasets
if( Cd == 1 && length(label) >= 20){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria D: Reproducibility"),"\n")
test_data <- imputed_data
label.vector <- names(table(label))
control.x <- as.data.frame(t(test_data[label == label.vector[1], -1]))
case.x <- as.data.frame(t(test_data[label == label.vector[2], -1]))
control.y <- rep(0, table(label)[1])
case.y <- rep(1, table(label)[2])
score <- try(stabel.score(repeats = 200, fold = 5, top = 20))
if (inherits(score, "try-error"))
print(score)
Fscore[paste(i,o,p,j,g,sep="")] <- score
}else{
message("'Criteria.D-Reproducibility' cannot be evaluated, Please Check!")
}
}
}
}
}
}else{
for(o in 3){
tran_train_data <- train_data_tran
tran_train_data[is.infinite(data.matrix(tran_train_data))] <- NA
cen_train_data <- try(cents(tran_train_data,o))
if (inherits(cen_train_data, "try-error"))
#if(class(cen_train_data) == "try-error")
{ next }
for(p in 1){
scal_factor <- try(scals(tran_train_data,p))
if (inherits(scal_factor, "try-error"))
#if(class(scal_factor)=="try-error")
{ next }
scal_train_data <- cen_train_data/scal_factor
if (inherits(scal_train_data, "try-error"))
#if(class(scal_train_data) == "try-error")
{ next }
for(j in 13){
scal_train_data[is.nan(scal_train_data)] <- NA
scal_train_data[is.infinite(scal_train_data)] <- NA
normalized_data <- try(norm(scal_train_data,j))
if (inherits(normalized_data, "try-error"))
#if(class(normalized_data) == "try-error")
{ next }
label_c <- as.factor(label)
g1 <- table(label_c)[levels(label_c)[1]]*0.8
g2 <- table(label_c)[levels(label_c)[2]]*0.8
train_data_filtering <- try(Basicfilter(normalized_data,label,g1=2,g2=2))
if (inherits(train_data_filtering, "try-error"))
#if(class(train_data_filtering)=="try-error")
{ next }
filter_train_data <- train_data_filtering
for(g in 1:7){
if(g==1){
imputed_data <- filter_train_data
}
if(g!=1){
imputed_data <- try(impute(filter_train_data,g))
if (inherits(imputed_data, "try-error"))
#if(class(imputed_data) == "try-error")
{ next }
}
##### Feature Selection
time = time+1
dataa <- data_q
frame <- imputed_data
label <- dataa[,2]
rots.out <- try(ROTS(data = t(frame), groups = as.character(label), B = 200, K = 500 , seed = 1234,log = FALSE))
if (inherits(rots.out, "try-error"))
#if(class(rots.out) == "try-error")
{ next }
frame <- imputed_data
label <- as.factor(as.character(label))
im.data <- data.frame(label=label,frame)
#(a) Precision of LFQ Based on the Proteomes among Replicates
if( Ca == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria A: Precision"),"\n")
data <- im.data
result <- PCV1(data)
pcv <- sapply(1:3, function(i){round(1000*mean(as.numeric(result[,i])))/1000})[3]
pmad <- try(mean(PMAD(data)))
if (inherits(pmad, "try-error"))
#if(class(pmad) == "try-error")
{ next }
Fpmad[paste(i,o,p,j,g,sep="")] <- pmad
}else{
message("'Criteria A: Precision' cannot be evaluated, Please Check!")
}
#(b) Classification Ability of LFQ between Distinct Sample Groups
if( Cb == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria B: Classification.Ability"),"\n")
data <- im.data
rots.out <- rots.out
col_pos <- which(rots.out$FDR<0.05)
if(length(col_pos) <= 10) {
markerid <- order(rots.out$pvalue)[1:20]
}else {
markerid <- which(rots.out$FDR<0.05)+1
}
clusters <- hclust(dist(data[,markerid]))
clusterCut <- cutree(clusters, 2)
dataa <- data_q
label <- dataa[,2]
tmatrix <- table(clusterCut, label)
tru <- as.numeric(data[,1])
accuracy <- (tmatrix[1,1]+ tmatrix[2,2])/length(label)
Faccuracy[paste(i,o,p,j,g,sep="")] <- accuracy
}else{
message("'Criteria B: Classification.Ability' cannot be evaluated, Please Check!")
}
#(c) Differential Expression Analysis Based on Reproducibility-optimization
if( Cc == 1 ){
cat(paste("Assessing" ,paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria C: Differential.Expression"),"\n")
rots.out <- rots.out
breaks <- seq(0,1,0.05)
sdres <- affy::hist(rots.out$pvalue,breaks=breaks)
bar1num <- sdres$counts[1]
baronummean <- mean(sdres$counts[-1])
barosd <- sd(sdres$counts[-1])
barorsd <- barosd/baronummean
Fbar1num[paste(i,o,p,j,g,sep="")] <- bar1num
Fbaronummean[paste(i,o,p,j,g,sep="")] <- baronummean
Fbarosd[paste(i,o,p,j,g,sep="")] <- barosd
Fbarorsd[paste(i,o,p,j,g,sep="")] <- barorsd
Fbaro_rsdr_to_bar1num[paste(i,o,p,j,g,sep="")] <- barorsd/bar1num
}else{
message("'Criteria C: Differential.Expression' cannot be evaluated, Please Check!")
}
#(d) Reproducibility of the Identified Protein Markers among Different Datasets
if( Cd == 1 && length(label) >= 20){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria D: Reproducibility"),"\n")
test_data <- imputed_data
label.vector <- names(table(label))
control.x <- as.data.frame(t(test_data[label == label.vector[1], -1]))
case.x <- as.data.frame(t(test_data[label == label.vector[2], -1]))
control.y <- rep(0, table(label)[1])
case.y <- rep(1, table(label)[2])
score <- try(stabel.score(repeats = 200, fold = 5, top = 20))
Fscore[paste(i,o,p,j,g,sep="")] <- score
}else{
message("'Criteria.D-Reproducibility' cannot be evaluated, Please Check!")
}
}
}
}
}
}
}
###################################################Step-2
if(!is.null(unlist(Faccuracy))){
Acc_revison <- unlist(Faccuracy)
Acc_revison_ID <- which(Acc_revison<0.5)
Acc_revison[Acc_revison_ID] <- 1-Acc_revison[Acc_revison_ID]
}
result <- dplyr::bind_rows("Precision"=unlist(Fpmad),"Classcification.Ability"=Acc_revison,"Differential.Expression"=unlist(Fbaro_rsdr_to_bar1num),"Reproducibility"=unlist(Fscore),.id = "id")
result1 <- t(result)
colnames(result1) <- result1["id",]
result2 <- result1[-1,]
result2 <- data.frame(result2, check.names=FALSE)
rownames(result2) <- names[match(rownames(result2), names[,1]), 2]
Rank <- apply(result2, 2, function(x){rank(as.numeric(as.character(x)), ties.method="min", na.last = "keep")})
rownames(Rank) <- rownames(result2)
if(length(grep("Reproducibility", colnames(Rank))) == 1){
Rank[,"Reproducibility"] <- rank(-as.numeric(as.character(result2[,"Reproducibility"])), ties.method="min", na.last = "keep")
}else{
Rank <- Rank
}
if(length(grep("Classcification.Ability", colnames(Rank))) == 1){
Rank[,"Classcification.Ability"] <- rank(-as.numeric(as.character(result2[,"Classcification.Ability"])),ties.method="min",na.last = "keep")
}else{
Rank <- Rank
}
Rank_revision <- apply(Rank, 2, function(x){x[is.na(x)] <- nrow(Rank);return(x)})
Ranksum0 <- apply(Rank_revision, 1, sum)
Rankres0 <- cbind("OverallRank" = Ranksum0,Rank_revision)
zuihou0 <- cbind("Rank"=Rankres0,"Value"=result2)
if(length(grep("Precision",colnames(result2))) == 1){
ID <- which(as.numeric(as.character(result2[,"Precision"]))>0.7)
AA <- zuihou0[-ID,]
AA0 <- AA[order(AA[,1], decreasing = FALSE),]
AA0[,1] <- rank(AA0[,1], ties.method = "min")
BB <- zuihou0[ID,]
BB0 <- BB[order(BB[,1], decreasing = FALSE),]
BB0[,1] <- rank(BB0[,1], ties.method = "min") + max(AA0[,1])
zuihou2 <- rbind(AA0,BB0)
}else{
zuihou1 <- zuihou0[order(Rankres0[,"OverallRank"], decreasing = FALSE),]
zuihou1[,1] <- rank(zuihou1[,1], ties.method = "min")
zuihou2 <- zuihou1
}
return(zuihou2)
###############################################End##########################################################################
}
JianboFu0406/EVALFQ111 documentation built on Aug. 10, 2020, 1:49 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions lfqevalueall
Documented in lfqevalueall
#' @title Conduct LFQ and assess performance of all possible LFQ workflows.
#' @description The ProteoLFQ enables the label-free quantification of proteomic
#' data and the performance assessment of each LFQ workflow from multiple perspectives. Moreover, it
#' provides the unique function of ranking all possible LFQ workflows (>3,000 random combinations of
#' transformation, normalization and imputation methods) based on their performances.
#' All in all, this tool makes the performance assessment of whole LFQ workflow possible
#' (collectively assessed by five well-established criteria with distinct underlying theories) and
#' gives the ranking results of all possible workflows based on the criteria preferred and selected by the users.
#' For function definitions and descriptions please use "??ProteoLFQ" command in R.
#' @param data_q This input file should be numeric type except the first and second column containing the names and label (control or case) of the studied samples, respectively. The intensity data should be provided in this input file with the following order: samples in row and proteins/peptides in column. Missing value (NA) of protein intensity are allowed.
#' @param assum_a all proteins were assumed to be equally important.The authors will be asked to input a letter “Y” to indicate the corresponding assumption is held for the studied dataset and a letter “N” to denote the opposite.
#' @param assum_b The level of protein abundance was assumed to be constant among all samples. The authors will be asked to input a letter “Y” to indicate the corresponding assumption is held for the studied dataset and a letter “N” to denote the opposite.
#' @param assum_c The intensities of the vast majority of the proteins were assumed to be unchanged under the studied conditions. The authors will be asked to input a letter “Y” to indicate the corresponding assumption is held for the studied dataset and a letter “N” to denote the opposite.
#' @param Ca Criterion (a): precision of LFQ based on the proteomes among replicates (Proteomics. 15:3140-51, 2015). If set 1, the user chooses to assess LFQ workflows using Criterion (a). If set 0, the user excludes Criterion (a) from performance assessment. The default setting of this value is “1”.
#' @param Cb Criterion (b): classification ability of LFQ between distinct sample groups (Nat Biotechnol. 28:83-9, 2010). If set 1, the user chooses to assess LFQ workflows using Criterion (b). If set 0, the user excludes Criterion (b) from performance assessment. The default setting of this value is “1”.
#' @param Cc Criterion (c): differential expression analysis by reproducibility-optimization (Nat Biotechnol. 32:896-902, 2014). If set 1, the user chooses to assess LFQ workflows using Criterion (c). If set 0, the user excludes Criterion (c) from performance assessment. The default setting of this value is “1”.
#' @param Cd Criterion (d): reproducibility of the identified protein markers among different datasets (Mol Biosyst. 11:1235-40, 2015). If set 1, the user chooses to assess LFQ workflows using Criterion (d). If set 0, the user excludes Criterion (d) from performance assessment. The default setting of this value is “1”.
#' @return preprocessed matrix
#' @import utils stats
#' @import metabolomics
#' @import affy vsn
#' @import MASS limma
#' @import ProteoMM ROTS
#' @importFrom grDevices colorRampPalette dev.off pdf
#' @useDynLib EVALFQ
#' @importFrom Rcpp sourceCpp
#' @rawNamespace import(dplyr, except=c(filter,lag,select,combine))
#' @rawNamespace import(gplots, except=lowess)
#' @importFrom pcaMethods pca
#' @importFrom pcaMethods completeObs
#' @import impute
#' @usage lfqevalueall(data_q,
#' assum_a="Y", assum_b="Y", assum_c="Y", Ca="1", Cb="1", Cc="1", Cd="1")
#' @export lfqevalueall
lfqevalueall <- function(data_q, assum_a="Y", assum_b="Y", assum_c="Y", Ca="1", Cb="1", Cc="1", Cd="1"){
#### 样本在行，特征在列
trans <- function(data,n){
matrix <- switch(
n,
Box_Cox(data),
log2(data),
data
)
return(matrix)
}
cents <- function(data,n){
matrix <- switch(
n,
MEC(data),
MDC(data),
data
)
return(matrix)
}
scals<- function(data,n){
matrix <- switch(
n,
1,
AUTO1(data),
PARETO1(data),
VAST1(data),
RANGE1(data)
)
return(matrix)
}
norm <- function(data,n){
matrix <- switch(
n,
t(fastlo(as.matrix(data))),
t(EIGENMS(data, label)),
t(LOWESS(data)),
t(SMAD(data)),
t(MEAN(data)),
t(MEDIAN(data)),
t(data),
t(PQN(data)),
t(QUANTILE(as.matrix(data))),
t(RLR1(data)),
t(MSTUS(data)),
t(TMM(data)),
t(VSN(as.matrix(data)))
)
return(matrix)
}
impute <- function(data,n){
matrix <- switch(
n,
filter_train_data,
t(back(filter_train_data)),
t(bpca(filter_train_data,nPcs=3)),
t(censor(filter_train_data)),
t(knn(filter_train_data,k=10)),
t(svdm(filter_train_data,nPcs=3)),
t(zero(filter_train_data))
)
return(matrix)
}
consistency <- function(fold = 5, top = 20) {
folds <- fold
control.label <- control.y # variable-1
test.fold1 <- split(sample(1:length(control.label)), 1:folds) #ignore warning
case.label <- case.y # variable-2
test.fold2 <- split(sample(1:length(case.label)), 1:folds) #ignore warning
DEG <- list()
for (i in 1:folds) {
com.x <- cbind(control.x[, test.fold1[[i]]], case.x[, test.fold2[[i]]]) # variable-3 & 4.
lab.ct <- test.fold1[[i]]
lab.ca <- test.fold2[[i]]
design <- cbind(Grp1 = 1, Grp2vs1 = c(rep(0, length(lab.ct)), rep(1, length(lab.ca))))
fit <- limma::lmFit(com.x, design)
fit <- limma::eBayes(fit)
DEG[[i]] <- rownames(limma::topTable(fit, coef = 2, number = nrow(com.x)))
}
names(DEG) <- LETTERS[1:folds]
top.n <-top # Extracting the top n genes.
DEG.list <- DEG
for (g in 1:length(DEG.list)) {
DEG.list[[g]] <- DEG.list[[g]][1:top.n]
}
# Calculating consistency score:
setlist <- DEG.list
OLlist <- overLapper(setlist=setlist, sep="", type="vennsets")
con.score <- 0
VennList <- OLlist$Venn_List
for (i in 1:length(VennList)) {
insect.n <- nchar(names(VennList[i]))
if (insect.n < 2) next
num.i <- 2^(insect.n - 2) * length(VennList[[i]])
con.score <- con.score + num.i
}
return(con.score) # consistense score
}
# Stable consistense score, with 20 repeats
stabel.score <- function(repeats = 20, fold = 5, top = 10) {
score <- 0
for (r in 1:repeats) {
score <- score + consistency(fold, top)
}
return(score/repeats)
}
#####################################################################
iName<-c("BOX","LOG","NON")
oName<-c("MEC","MDC","NON")
pName<-c("NON","ATO","PAR","RAN","VAS")
jName<-c("CYC","EIG","LOW","MAD","MEA","MED","NON","PQN","QUA","RLR","TIC","TMM","VSN")
gName<-c("NON","BAK","BPC","CEN","KNN","SVD","ZER")
dataa<-data_q
rownames(dataa)<-dataa[,1]
label <- dataa[,2]
frame <- dataa[, -(1:2)]
frame <- t(frame)
frame <- data.frame(frame)
train_data <- data.matrix(frame, rownames.force = NA)
train_data_t <- train_data
Fpcv <- list()
Fscore <- list()
Faccuracy <- list()
Fpmad <- list()
Fbar1num <- list()
Fbaronummean <- list()
Fbarosd <- list()
Fbarorsd <- list()
Fbaro_rsdr_to_bar1num <- list()
spike <- list()
backgound <- list()
time <- 0
center <- c(1,2,3)
scaling <- c(1,2,3,4,5)
normalization <- c(1,2,3,4,5,6,7,8,9,10,11,12)
cat("Assumption A: all proteins were equally important (Y/N): ", assum_a, "\n")
cat("Assumption B: the level of protein abundance was constant among all samples (Y/N): ", assum_b, "\n")
cat("Assumption C: the intensity of the majority of proteins were unchanged (Y/N): ", assum_c, "\n")
if(is.na(match("Y",assum_a))){
center<-3
scaling<-1
}
if(is.na(match("Y",assum_b)) & is.na(match("Y",assum_c))) {
normalization <- c(2,7)
}
if(is.na(match("Y",assum_b)) & !is.na(match("Y",assum_c))){
normalization <- c(1,2,7,8,9,10,11)
}
if(!is.na(match("Y",assum_b)) & is.na(match("Y",assum_c))){
normalization <- c(1,4,5,6,7,11)
}
for(i in 1:3){
train_data_tran <- try(trans(train_data_t,i))
if (inherits(train_data_tran, "try-error"))
{ next }
if(i!=3){
for(o in center){
tran_train_data <- train_data_tran
tran_train_data[is.infinite(data.matrix(tran_train_data))] <- NA
cen_train_data <- try(cents(tran_train_data,o))
if (inherits(cen_train_data, "try-error"))
{ next }
for(p in scaling){
scal_factor <- try(scals(tran_train_data,p))
if (inherits(scal_factor, "try-error"))
#if(class(scal_factor) == "try-error")
{ next }
scal_train_data <- cen_train_data/scal_factor
if (inherits(scal_train_data, "try-error"))
#if(class(scal_train_data)=="try-error")
{ next }
for(j in normalization){
scal_train_data[is.nan(scal_train_data)] <- NA
scal_train_data[is.infinite(scal_train_data)] <- NA
normalized_data <- try(norm(scal_train_data,j))
if (inherits(normalized_data, "try-error"))
#if(class(normalized_data)=="try-error")
{ next }
label_c <- as.factor(label)
g1 <- table(label_c)[levels(label_c)[1]]*0.8
g2 <- table(label_c)[levels(label_c)[2]]*0.8
train_data_filtering <- try(Basicfilter(normalized_data,label,g1=2,g2=2))
if (inherits(train_data_filtering, "try-error"))
#if(class(train_data_filtering) == "try-error")
{ next }
filter_train_data <- train_data_filtering
for(g in 1:7){
if(g==1){
imputed_data <- filter_train_data
}
if(g!=1){
imputed_data <- try(impute(filter_train_data,g))
if (inherits(imputed_data, "try-error"))
#if(class(imputed_data)=="try-error")
{ next }
}
##### Feature Selection
time = time+1
dataa <- data_q
frame <- imputed_data
label <- dataa[,2]
rots.out <-try(ROTS(data = t(frame), groups = as.character(label), B = 200, K = 500 , seed = 1234,log = FALSE))
if (inherits(rots.out, "try-error"))
#if(class(rots.out)=="try-error")
{ next }
frame <- imputed_data
label <- as.factor(as.character(label))
im.data <- data.frame(label=label,frame)
#(a) Precision of LFQ Based on the Proteomes among Replicates
if( Ca == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria A: Precision"),"\n")
data <- im.data
result <- PCV1(data)
pcv <- sapply(1:3, function(i){round(1000*mean(as.numeric(result[,i])))/1000})[3]
pmad <- try(mean(PMAD(data)))
if (inherits(pmad, "try-error"))
#if(class(pmad) == "try-error")
{ next }
Fpmad[paste(i,o,p,j,g,sep="")] <- pmad
}else{
message("'Criteria A: Precision' cannot be evaluated, Please Check!")
}
#(b) Classification Ability of LFQ between Distinct Sample Groups
if( Cb == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria B: Classification.Ability"),"\n")
data <- im.data
rots.out <- rots.out
col_pos <- which(rots.out$FDR<0.05)
if(length(col_pos) <= 10) {
markerid <-order(rots.out$pvalue)[1:20]
}else {
markerid <- which(rots.out$FDR<0.05)+1
}
clusters <- hclust(dist(data[,markerid]))
clusterCut <- cutree(clusters, 2)
dataa <- data_q
label <- dataa[,2]
tmatrix <- table(clusterCut, label)
tru <- as.numeric(data[,1])
accuracy <- (tmatrix[1,1]+ tmatrix[2,2])/length(label)
Faccuracy[paste(i,o,p,j,g,sep="")] <- accuracy
}else{
message("'Criteria B: Classification.Ability' cannot be evaluated, Please Check!")
}
#(c) Differential Expression Analysis Based on Reproducibility-optimization
if( Cc == 1 ){
cat(paste("Assessing" ,paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria C: Differential.Expression"),"\n")
rots.out <- rots.out
breaks <- seq(0,1,0.05)
sdres <- affy::hist(rots.out$pvalue,breaks=breaks)
bar1num <- sdres$counts[1]
baronummean <- mean(sdres$counts[-1])
barosd <- sd(sdres$counts[-1])
barorsd <- barosd/baronummean
Fbar1num[paste(i,o,p,j,g,sep="")] <- bar1num
Fbaronummean[paste(i,o,p,j,g,sep="")] <- baronummean
Fbarosd[paste(i,o,p,j,g,sep="")] <- barosd
Fbarorsd[paste(i,o,p,j,g,sep="")] <- barorsd
Fbaro_rsdr_to_bar1num[paste(i,o,p,j,g,sep="")] <- barorsd/bar1num
}else{
message("'Criteria C: Differential.Expression' cannot be evaluated, Please Check!")
}
#(d) Reproducibility of the Identified Protein Markers among Different Datasets
if( Cd == 1 && length(label) >= 20){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria D: Reproducibility"),"\n")
test_data <- imputed_data
label.vector <- names(table(label))
control.x <- as.data.frame(t(test_data[label == label.vector[1], -1]))
case.x <- as.data.frame(t(test_data[label == label.vector[2], -1]))
control.y <- rep(0, table(label)[1])
case.y <- rep(1, table(label)[2])
score <- try(stabel.score(repeats = 200, fold = 5, top = 20))
if (inherits(score, "try-error"))
print(score)
Fscore[paste(i,o,p,j,g,sep="")] <- score
}else{
message("'Criteria.D-Reproducibility' cannot be evaluated, Please Check!")
}
}
}
}
}
}else{
for(o in 3){
tran_train_data <- train_data_tran
tran_train_data[is.infinite(data.matrix(tran_train_data))] <- NA
cen_train_data <- try(cents(tran_train_data,o))
if (inherits(cen_train_data, "try-error"))
#if(class(cen_train_data) == "try-error")
{ next }
for(p in 1){
scal_factor <- try(scals(tran_train_data,p))
if (inherits(scal_factor, "try-error"))
#if(class(scal_factor)=="try-error")
{ next }
scal_train_data <- cen_train_data/scal_factor
if (inherits(scal_train_data, "try-error"))
#if(class(scal_train_data) == "try-error")
{ next }
for(j in 13){
scal_train_data[is.nan(scal_train_data)] <- NA
scal_train_data[is.infinite(scal_train_data)] <- NA
normalized_data <- try(norm(scal_train_data,j))
if (inherits(normalized_data, "try-error"))
#if(class(normalized_data) == "try-error")
{ next }
label_c <- as.factor(label)
g1 <- table(label_c)[levels(label_c)[1]]*0.8
g2 <- table(label_c)[levels(label_c)[2]]*0.8
train_data_filtering <- try(Basicfilter(normalized_data,label,g1=2,g2=2))
if (inherits(train_data_filtering, "try-error"))
#if(class(train_data_filtering)=="try-error")
{ next }
filter_train_data <- train_data_filtering
for(g in 1:7){
if(g==1){
imputed_data <- filter_train_data
}
if(g!=1){
imputed_data <- try(impute(filter_train_data,g))
if (inherits(imputed_data, "try-error"))
#if(class(imputed_data) == "try-error")
{ next }
}
##### Feature Selection
time = time+1
dataa <- data_q
frame <- imputed_data
label <- dataa[,2]
rots.out <- try(ROTS(data = t(frame), groups = as.character(label), B = 200, K = 500 , seed = 1234,log = FALSE))
if (inherits(rots.out, "try-error"))
#if(class(rots.out) == "try-error")
{ next }
frame <- imputed_data
label <- as.factor(as.character(label))
im.data <- data.frame(label=label,frame)
#(a) Precision of LFQ Based on the Proteomes among Replicates
if( Ca == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria A: Precision"),"\n")
data <- im.data
result <- PCV1(data)
pcv <- sapply(1:3, function(i){round(1000*mean(as.numeric(result[,i])))/1000})[3]
pmad <- try(mean(PMAD(data)))
if (inherits(pmad, "try-error"))
#if(class(pmad) == "try-error")
{ next }
Fpmad[paste(i,o,p,j,g,sep="")] <- pmad
}else{
message("'Criteria A: Precision' cannot be evaluated, Please Check!")
}
#(b) Classification Ability of LFQ between Distinct Sample Groups
if( Cb == 1 ){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-",oName[o],"-",pName[p],"-", jName[j], "-",gName[g],sep=""),"Under Criteria B: Classification.Ability"),"\n")
data <- im.data
rots.out <- rots.out
col_pos <- which(rots.out$FDR<0.05)
if(length(col_pos) <= 10) {
markerid <- order(rots.out$pvalue)[1:20]
}else {
markerid <- which(rots.out$FDR<0.05)+1
}
clusters <- hclust(dist(data[,markerid]))
clusterCut <- cutree(clusters, 2)
dataa <- data_q
label <- dataa[,2]
tmatrix <- table(clusterCut, label)
tru <- as.numeric(data[,1])
accuracy <- (tmatrix[1,1]+ tmatrix[2,2])/length(label)
Faccuracy[paste(i,o,p,j,g,sep="")] <- accuracy
}else{
message("'Criteria B: Classification.Ability' cannot be evaluated, Please Check!")
}
#(c) Differential Expression Analysis Based on Reproducibility-optimization
if( Cc == 1 ){
cat(paste("Assessing" ,paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria C: Differential.Expression"),"\n")
rots.out <- rots.out
breaks <- seq(0,1,0.05)
sdres <- affy::hist(rots.out$pvalue,breaks=breaks)
bar1num <- sdres$counts[1]
baronummean <- mean(sdres$counts[-1])
barosd <- sd(sdres$counts[-1])
barorsd <- barosd/baronummean
Fbar1num[paste(i,o,p,j,g,sep="")] <- bar1num
Fbaronummean[paste(i,o,p,j,g,sep="")] <- baronummean
Fbarosd[paste(i,o,p,j,g,sep="")] <- barosd
Fbarorsd[paste(i,o,p,j,g,sep="")] <- barorsd
Fbaro_rsdr_to_bar1num[paste(i,o,p,j,g,sep="")] <- barorsd/bar1num
}else{
message("'Criteria C: Differential.Expression' cannot be evaluated, Please Check!")
}
#(d) Reproducibility of the Identified Protein Markers among Different Datasets
if( Cd == 1 && length(label) >= 20){
cat(paste("Assessing" , paste(time,".",sep=""), paste(iName[i],"-", oName[o],"-",pName[p],"-",jName[j], "-",gName[g],sep=""),"Under Criteria D: Reproducibility"),"\n")
test_data <- imputed_data
label.vector <- names(table(label))
control.x <- as.data.frame(t(test_data[label == label.vector[1], -1]))
case.x <- as.data.frame(t(test_data[label == label.vector[2], -1]))
control.y <- rep(0, table(label)[1])
case.y <- rep(1, table(label)[2])
score <- try(stabel.score(repeats = 200, fold = 5, top = 20))
Fscore[paste(i,o,p,j,g,sep="")] <- score
}else{
message("'Criteria.D-Reproducibility' cannot be evaluated, Please Check!")
}
}
}
}
}
}
}
###################################################Step-2
if(!is.null(unlist(Faccuracy))){
Acc_revison <- unlist(Faccuracy)
Acc_revison_ID <- which(Acc_revison<0.5)
Acc_revison[Acc_revison_ID] <- 1-Acc_revison[Acc_revison_ID]
}
result <- dplyr::bind_rows("Precision"=unlist(Fpmad),"Classcification.Ability"=Acc_revison,"Differential.Expression"=unlist(Fbaro_rsdr_to_bar1num),"Reproducibility"=unlist(Fscore),.id = "id")
result1 <- t(result)
colnames(result1) <- result1["id",]
result2 <- result1[-1,]
result2 <- data.frame(result2, check.names=FALSE)
rownames(result2) <- names[match(rownames(result2), names[,1]), 2]
Rank <- apply(result2, 2, function(x){rank(as.numeric(as.character(x)), ties.method="min", na.last = "keep")})
rownames(Rank) <- rownames(result2)
if(length(grep("Reproducibility", colnames(Rank))) == 1){
Rank[,"Reproducibility"] <- rank(-as.numeric(as.character(result2[,"Reproducibility"])), ties.method="min", na.last = "keep")
}else{
Rank <- Rank
}
if(length(grep("Classcification.Ability", colnames(Rank))) == 1){
Rank[,"Classcification.Ability"] <- rank(-as.numeric(as.character(result2[,"Classcification.Ability"])),ties.method="min",na.last = "keep")
}else{
Rank <- Rank
}
Rank_revision <- apply(Rank, 2, function(x){x[is.na(x)] <- nrow(Rank);return(x)})
Ranksum0 <- apply(Rank_revision, 1, sum)
Rankres0 <- cbind("OverallRank" = Ranksum0,Rank_revision)
zuihou0 <- cbind("Rank"=Rankres0,"Value"=result2)
if(length(grep("Precision",colnames(result2))) == 1){
ID <- which(as.numeric(as.character(result2[,"Precision"]))>0.7)
AA <- zuihou0[-ID,]
AA0 <- AA[order(AA[,1], decreasing = FALSE),]
AA0[,1] <- rank(AA0[,1], ties.method = "min")
BB <- zuihou0[ID,]
BB0 <- BB[order(BB[,1], decreasing = FALSE),]
BB0[,1] <- rank(BB0[,1], ties.method = "min") + max(AA0[,1])
zuihou2 <- rbind(AA0,BB0)
}else{
zuihou1 <- zuihou0[order(Rankres0[,"OverallRank"], decreasing = FALSE),]
zuihou1[,1] <- rank(zuihou1[,1], ties.method = "min")
zuihou2 <- zuihou1
}
return(zuihou2)
###############################################End##########################################################################
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.