|
Commit message:
planemo upload for repository https://github.com/workflow4metabolomics/tools-metabolomics commit 8d2ca678d973501b60479a8dc3f212eecd56eab8 |
|
added:
diagmfl.R mixmodel.xml mixmodel_script.R mixmodel_wrapper.R test-data/OneFactor_dataMatrix.tsv test-data/OneFactor_sampleMetadata.tsv test-data/OneFactor_variableMetadata.tsv test-data/TwoFactor_dataMatrix.txt test-data/TwoFactor_sampleMetadata.txt test-data/TwoFactor_variableMetadata.txt test-data/mixmodel_TwoFactor_variableMetadata.txt |
| b |
| diff -r 000000000000 -r a4d89d47646f diagmfl.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diagmfl.R Mon May 16 09:25:01 2022 +0000 |
| [ |
| b'@@ -0,0 +1,918 @@\n+#\' Calcul des grandeurs "diagnostiques"\r\n+#\'\r\n+#\' Script adapte de http://www.ime.unicamp.br/~cnaber/residdiag_nlme_v22.R pour fonctionner\r\n+#\' avec un modele lmer (et non lme), des sujets avec des identifiants non numeriques,\r\n+#\' et des observations non ordonnees sujet par sujet (dernier point a verifier.)\r\n+#\'\r\n+#\' @detail Les graphiques, les calculs associ\xc3\xa9s et les notations utilisees dans le script suivent\r\n+#\' l\'article de Singer et al (2016) Graphical Tools for detedcting departures from linear\r\n+#\' mixed model assumptions and some remedial measures, International Statistical Review\r\n+#\' (doi:10.1111/insr.12178)\r\n+#\'\r\n+#\' @param mfl A linear mixed model fitted via lmer or a data frame containing data\r\n+#\' @return A list\r\n+#\' @author Natacha Lenuzza\r\n+#\' @examples\r\n+#\' print("hello !")\r\n+#\'\r\n+#\' @export lmer.computeDiag\r\n+\r\n+\r\n+\r\n+lmer.computeDiag <- function(mfl) {\r\n+\r\n+ ## Check arguments ---------------------------------------------------------\r\n+\r\n+ if (length(mfl@flist) > 1)\r\n+ stop("Several \'grouping level\' for random effect not implemented yet.")\r\n+\r\n+\r\n+ ## extracting information from mfl models -------------------------------------------------------------\r\n+ # data\r\n+ df <- mfl@frame\r\n+ responseC <- names(df)[1]\r\n+ unitC <- names(mfl@flist)[1]\r\n+ # observations\r\n+ yVn <- df[, responseC]\r\n+ nobsN <- length(yVn)\r\n+ # units\r\n+ idunitVc <- levels(mfl@flist[[1]])\r\n+ nunitN <- length(unique(idunitVc))\r\n+ #X\r\n+ xMN <- mfl@pp$X\r\n+ pN <- ncol(xMN)\r\n+ #Z\r\n+ zMN <- t(as.matrix(mfl@pp$Zt))\r\n+ # Estimated covariance matrix of random effects (Gam)\r\n+ aux <- VarCorr(mfl)[[1]] ## assuming only one level of grouping\r\n+ aux2 <- attr(aux, "stddev")\r\n+ gMN <- attr(aux, "correlation") * (aux2 %*% t(aux2))\r\n+ gammaMN <- as.matrix(kronecker(diag(nunitN), gMN))\r\n+ q <- dim(gMN)[1]\r\n+ # Estimated covariance matrix of conditonal error (homoskedastic conditional independance model)\r\n+ sigsqN <- attr(VarCorr(mfl), "sc")^2\r\n+ rMN <- sigsqN * diag(nobsN)\r\n+ # Estimated covariance matrix of Y\r\n+ vMN <- (zMN %*% gammaMN %*% t(zMN)) + rMN\r\n+ invvMN <- MASS::ginv(vMN)\r\n+ # H and Q matrix\r\n+ hMN <- MASS::ginv(t(xMN) %*% invvMN %*% xMN)\r\n+ qMN <- invvMN - invvMN %*% xMN %*% (hMN) %*% t(xMN) %*% invvMN\r\n+ # eblue and eblup\r\n+ eblueVn <- mfl@beta\r\n+ eblupVn <- gammaMN %*% t(zMN) %*% invvMN %*% (yVn - xMN %*% eblueVn) ## equivalent de ranef(mfl)\r\n+ rownames(eblupVn) <- colnames(zMN)\r\n+ ## Calculs of matrices and vectors used in graph diagnosics ---------------------------------------------\r\n+ ## Marginal and individual predictions, residuals and variances\r\n+ marpredVn <- xMN %*% eblueVn\r\n+ marresVn <- yVn - marpredVn\r\n+ marvarMN <- vMN - xMN %*% hMN %*% t(xMN)\r\n+ condpredVn <- marpredVn + zMN %*% eblupVn\r\n+ condresVn <- yVn - condpredVn\r\n+ condvarMN <- rMN %*% qMN %*% rMN\r\n+ ## Analysis of marginal and conditional residuals\r\n+ stmarresVn <- stcondresVn <- rep(0, nobsN)\r\n+ lesverVn <- rep(0, nunitN)\r\n+ names(lesverVn) <- idunitVc\r\n+ for (i in 1:nunitN) {\r\n+ idxiVn <- which(df[, unitC] == idunitVc[i]) ## position des observations du sujet i\r\n+ miN <- length(idxiVn)\r\n+ ## standardization of marginal residual\r\n+ stmarresVn[idxiVn] <- as.vector(solve(sqrtmF(marvarMN[idxiVn, idxiVn])) %*% marresVn[idxiVn])\r\n+ ##Standardized Lessafre and Verbeke\'s measure\r\n+ auxMN <- diag(1, ncol = miN, nrow = miN) - stmarresVn[idxiVn] %*% t(stmarresVn[idxiVn])\r\n+ lesverVn[i] <- sum(diag(auxMN %*% t(auxMN)))\r\n+ ## standardization of conditional residual\r\n+ stcondresVn[idxiVn] <- as.vector(solve(sqrtmF(condvarMN[idxiVn, idxiVn])) %*% condresVn[idxiVn])\r\n+ }\r\n+ lesverVn <- lesverVn / sum(lesverVn)\r\n+ ## Least confounded conditional residuals\r\n+ ## EBLUP analysis (Mahalanobis\' distance)\r\n+ varbMN <- gammaMN %*% t(zMN) %*% qMN %*% zMN %*% gammaMN\r\n+ mdistVn <- rep(0, nunitN)\r\n+ qm <- q - 1\r\n+ for (j in 1:nunitN) {\r\n+ gbi <- varbMN[(q * j - qm):'..b'diff(range(x)) / (2 * IQR(x) / length(x) ^ (1 / 3))\r\n+ )\r\n+ }\r\n+\r\n+\r\n+ daf <- data.frame(x = x)\r\n+ ## graph\r\n+ return(ggplot(data = daf, aes(x)) +\r\n+ geom_histogram(aes(y = ..density..),\r\n+ col = "black", fill = "grey", binwidth = bw) +\r\n+ geom_density(size = 1.2,\r\n+ col = "blue",\r\n+ linetype = "blank",\r\n+ fill = rgb(0, 0, 1, 0.1)) +\r\n+ stat_function(fun = dnorm,\r\n+ args = list(mean = 0, sd = sd_x),\r\n+ col = "blue", size = 1.2) +\r\n+ theme(legend.position = "none") +\r\n+ xlab(""))\r\n+}\r\n+\r\n+\r\n+plot.res.Lmixed <- function(mfl, df, title = "", pvalCutof = 0.05) {\r\n+\r\n+ ## define subscript of the different columns depending if we have only time (ncol(df)=3) or not\r\n+ if (ncol(df) > 3) {\r\n+ varidx <- 4\r\n+ ffidx <- 1\r\n+ timidx <- 2\r\n+ individx <- 3\r\n+ } else {\r\n+ varidx <- 3\r\n+ ffidx <- 1\r\n+ timidx <- 1\r\n+ individx <- 2\r\n+ }\r\n+ nameVar <- colnames(df)[varidx]\r\n+ fflab <- colnames(df)[ffidx]\r\n+ ## Individual time-course\r\n+ rawPlot <-\r\n+ ggplot(data = df, aes(x = df[[timidx]], y = df[[varidx]], colour = df[[ffidx]], group = df[[individx]])) +\r\n+ geom_point() +\r\n+ geom_line() + ggtitle("Individual time-courses (raw data)") +\r\n+ ylab(nameVar) +\r\n+ xlab(label = colnames(df)[2]) +\r\n+ theme(legend.title = element_blank(), legend.position = "none", plot.title = element_text(size = rel(1.2), face = "bold"))\r\n+ ## Boxplot of fixed factor\r\n+ bPlot <-\r\n+ ggplot(data = df, aes(y = df[[varidx]], x = df[[ffidx]], color = df[[ffidx]])) +\r\n+ geom_boxplot(outlier.colour = "red", outlier.shape = 8, outlier.size = 4) +\r\n+ ggtitle(paste("Boxplot by ", fflab, sep = "")) + xlab("") + ylab("") +\r\n+ theme(legend.title = element_blank(), plot.title = element_text(size = rel(1.2), face = "bold"))\r\n+ ## Post-hoc estimates\r\n+ ddlsm1 <- mfl\r\n+ ddlsm1$name <- rownames(ddlsm1)\r\n+ ddlsm1$Significance <- rep("NS", nrow(ddlsm1))\r\n+ ## modif JF pour tenir compte du seuil de pvalues defini par le user\r\n+ options("scipen" = 100, "digits" = 5)\r\n+ pvalCutof <- as.numeric(pvalCutof)\r\n+ bs <- 0.05; bm <- 0.01; bi <- 0.005\r\n+ if (pvalCutof > bm) {\r\n+ bs <- pvalCutof\r\n+ } else\r\n+ if (pvalCutof < bm & pvalCutof > bi) {\r\n+ bm <- pvalCutof; bs <- pvalCutof\r\n+ } else\r\n+ if (pvalCutof < bi) {\r\n+ bi <- pvalCutof; bm <- pvalCutof; bs <- pvalCutof\r\n+ }\r\n+ lbs <- paste("p-value < ", bs, sep = "")\r\n+ lbm <- paste("p-value < ", bm, sep = "")\r\n+ lbi <- paste("p-value < ", bi, sep = "")\r\n+ cols <- paste("p-value < ", bs, sep = "")\r\n+ colm <- paste("p-value < ", bm, sep = "")\r\n+ coli <- paste("p-value < ", bi, sep = "")\r\n+ valcol <- c("grey", "yellow", "orange", "red")\r\n+ names(valcol) <- c("NS", lbs, lbm, lbi)\r\n+ ddlsm1$Significance[which(ddlsm1$p.value <= bs)] <- lbs\r\n+ ddlsm1$Significance[which(ddlsm1$p.value < bs & ddlsm1$p.value >= bm)] <- lbm\r\n+ ddlsm1$Significance[which(ddlsm1$p.value < bi)] <- lbi\r\n+ ddlsm1$levels <- rownames(ddlsm1)\r\n+ ddlsm1$term <- sapply(rownames(ddlsm1), function(namC) {\r\n+ strsplit(namC, split = " ", fixed = TRUE)[[1]][1]\r\n+ })\r\n+\r\n+ phPlot <-\r\n+ ggplot(ddlsm1, aes(x = levels, y = Estimate)) +\r\n+ facet_grid(facets = ~term, ddlsm1, scales = "free", space = "free") +\r\n+ geom_bar(aes(fill = Significance), stat = "identity") +\r\n+ theme(axis.text.x = element_text(angle = 90, hjust = 1)) +\r\n+ scale_fill_manual(\r\n+ values = valcol) +\r\n+ geom_errorbar(aes(ymin = Lower.CI, ymax = Upper.CI), width = 0.25) +\r\n+ ggtitle("Post-hoc estimates ") + xlab("") +\r\n+ theme(plot.title = element_text(size = rel(1.2), face = "bold"))\r\n+\r\n+ ## Final plotting\r\n+ grid.arrange(arrangeGrob(rawPlot, bPlot, ncol = 2),\r\n+ phPlot, nrow = 2,\r\n+ top = textGrob(title, gp = gpar(fontsize = 32, font = 4))\r\n+ )\r\n+\r\n+}\r\n' |
| b |
| diff -r 000000000000 -r a4d89d47646f mixmodel.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel.xml Mon May 16 09:25:01 2022 +0000 |
| [ |
| b'@@ -0,0 +1,264 @@\n+<tool id="mixmodel4repeated_measures" name="mixmodel" version="3.1.0" profile="19.01">\r\n+ <description>ANOVA for repeated measures statistics</description>\r\n+\r\n+ <requirements>\r\n+ <requirement type="package" version="1.1_27">r-lme4</requirement>\r\n+ <requirement type="package" version="3.1_3">r-lmertest</requirement>\r\n+ <requirement type="package" version="3.2_0">r-ggplot2</requirement>\r\n+ <requirement type="package" version="2.3">r-gridExtra</requirement>\r\n+ <requirement type="package" version="2.40.0">bioconductor-multtest</requirement>\r\n+ </requirements>\r\n+\r\n+ <command detect_errors="exit_code"><![CDATA[\r\n+ Rscript \'$__tool_directory__/mixmodel_wrapper.R\'\r\n+\r\n+ dataMatrix_in \'$dataMatrix_in\'\r\n+ sampleMetadata_in \'$sampleMetadata_in\'\r\n+ variableMetadata_in \'$variableMetadata_in\'\r\n+\r\n+ fixfact \'$fixfact\'\r\n+ time \'$time\'\r\n+ subject \'$subject\'\r\n+ adjC \'$adjC\'\r\n+ trf \'$trf\'\r\n+ thrN \'$thrN\'\r\n+ diaR \'$diaR\'\r\n+ dff \'$dff\'\r\n+ rounding \'${roundchoice.rounding}\'\r\n+ #if $roundchoice.rounding == \'yes\' :\r\n+ decplaces \'${roundchoice.decplaces}\'\r\n+ #end if\r\n+\r\n+ variableMetadata_out \'$variableMetadata_out\'\r\n+ out_graph_pdf \'$out_graph_pdf\'\r\n+ out_estim_pdf \'$out_estim_pdf\'\r\n+ information \'$information\'\r\n+\r\n+ ]]></command>\r\n+\r\n+ <inputs>\r\n+ <param name="dataMatrix_in" label="Data matrix file" type="data" format="tabular" help="variable x sample, decimal: \'.\', missing: NA, mode: numerical, sep: tabular" />\r\n+ <param name="sampleMetadata_in" label="Sample metadata file" type="data" format="tabular" help="sample x metadata, decimal: \'.\', missing: NA, mode: character and numerical, sep: tabular" />\r\n+ <param name="variableMetadata_in" label="Variable metadata file" type="data" format="tabular" help="variable x metadata, decimal: \'.\', missing: NA, mode: character and numerical, sep: tabular" />\r\n+ <param name="fixfact" label="Fixed Factor of interest" type="text" help="Name of sample metadata column corresponding to the fixed factor (use none if only time factor"/>\r\n+ <param name="time" label="Repeated factor (time)" type="text" help="Name of the column of the sample metadata table corresponding to the repeated factor"/>\r\n+ <param name="subject" label="Subject factor" type="text" help="Name of the column of the sample metadata table corresponding to the subject factor"/>\r\n+ <param name="dff" label="Degrees of freedom method for the F-tests" type="select" help="">\r\n+ <option value="Satt">Satterthwaite</option>\r\n+ <option value="KenR">Kenward-Roger</option>\r\n+ </param>\r\n+ <param name="adjC" label="Method for multiple testing correction" type="select" help="">\r\n+ <option value="fdr">fdr</option>\r\n+ <option value="BH">BH</option>\r\n+ <option value="bonferroni">bonferroni</option>\r\n+ <option value="BY">BY</option>\r\n+ <option value="hochberg">hochberg</option>\r\n+ <option value="holm">holm</option>\r\n+ <option value="hommel">hommel</option>\r\n+ <option value="none">none</option>\r\n+ </param>\r\n+ <param name="trf" label="Log transform of raw data" type="select" help="Transformation of raw data">\r\n+ <option value="none">none</option>\r\n+ <option value="log10">log10</option>\r\n+ <option value="log2">log2</option>\r\n+ </param>\r\n+ <param name="thrN" type="float" value="0.05" label="(Corrected) p-value significance threshold" help="Must be between 0 and 1"/>\r\n+ <param name="diaR" label="Perform diagnostic of the residuals" type="select" help=" Used to assess the quality of models considering distribution of residuals ">\r\n+ <option value="yes">yes</option>\r\n+ <option value="no">no</option>\r\n'..b'88) ("hommel"), Benjamini and Hochberg (1995) ("BH" or its alias "fdr"), and Benjamini and Yekutieli (2001) ("BY"), respectively. A pass-through option ("none") is also included. The set of methods are contained in the p.adjust.methods vector for the benefit of methods that need to have the method as an option and pass it on to p.adjust. The first four methods are designed to give strong control of the family-wise error rate. There seems no reason to use the unmodified Bonferroni correction because it is dominated by Holm\'s method, which is also valid under arbitrary assumptions. Hochberg\'s and Hommel\'s methods are valid when the hypothesis tests are independent or when they are non-negatively associated (Sarkar, 1998; Sarkar and Chang, 1997). Hommel\'s method is more powerful than Hochberg\'s, but the difference is usually small and the Hochberg p-values are faster to compute. The "BH" (aka "fdr") and "BY" method of Benjamini, Hochberg, and Yekutieli control the false discovery rate, the expected proportion of false discoveries amongst the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so these methods are more powerful than the others."\r\n+\r\n+\r\n+(Corrected) p-value significance threshold\r\n+|\r\n+|\r\n+\r\n+Rounding the result\'s numerical columns\r\n+| This parameter enables the choice of a given number of decimal places to be used to display results\' values in the output columns of the variableMetadata that are generated by the tool.\r\n+| If set to "Yes", you can indicate the number of decimal places wanted, and all the results\' values will be rounded according to this number.\r\n+|\r\n+\r\n+------------\r\n+Output files\r\n+------------\r\n+\r\n+variableMetadata_out.tabular\r\n+| **variableMetadata** file identical to the file given as argument plus\r\n+| pvalue of Shapiro normality test of the residuals\r\n+| pvalues of the main effects and interaction\r\n+| PostHoc test with difference between levels and pvalues of these difference\r\n+|\r\n+\r\n+mixedmodel_Estimates\r\n+| A pdf file is created with a graphical representation of differences among levels of factors with a color code for significance and an error bar.\r\n+| Plots are only provided for features with at least one significant factor. \r\n+\r\n+\r\n+mixedmodel_diagResiduals\r\n+| if Perform diagnostic of the residuals" is set to yes(default) a pdf file is created with a a serie of\r\n+| graphics outputed in order to assess the distribution of residuals to check the adjustment.\r\n+| Plots are only provided for features with at least one factor being significant based on the non-corrected p-values. \r\n+\r\n+information.txt\r\n+| File with all messages and warnings generated during the computation\r\n+| The list of variables with name and a tag if it is significant for at least fixed or repeated factor.\r\n+\r\n+\r\n+\r\n+---------------------------------------------------\r\n+\r\n+Changelog/News\r\n+--------------\r\n+\r\n+**Version 3.1.0 - April 2022**\r\n+\r\n+- NEW FEATURE: the tool now can be used without fixed factor by specify "none" in the concerned parameter field.\r\n+\r\n+- NEW FEATURE: addition of the possibility to choose the degrees of freedom computation method (previously using only the Satterthwaite method).\r\n+\r\n+- NEW FEATURE: addition of the possibility to choose the number of decimal places in the variable-metadata numerical columns\' output.\r\n+\r\n+\r\n+ ]]></help>\r\n+\r\n+ <citations>\r\n+ <citation type="doi">10.18637/jss.v082.i13.</citation>\r\n+ <citation type="bibtex">@ARTICLE{fisher,\r\n+ author = {Benjamini Y. and Hochberg Y.,\r\n+ title = {Controlling the false discovery rate: a practical and powerful approach for multiple testing. Journal of the Royal Statistical Society},\r\n+ journal = {Series B (Methodological)},\r\n+ year = {1995},\r\n+ volume = {57},\r\n+ pages = {289-300}\r\n+ }</citation>\r\n+ <citation type="doi">10.1093/bioinformatics/btu813</citation>\r\n+ </citations>\r\n+\r\n+</tool>\r\n' |
| b |
| diff -r 000000000000 -r a4d89d47646f mixmodel_script.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel_script.R Mon May 16 09:25:01 2022 +0000 |
| [ |
| b'@@ -0,0 +1,481 @@\n+####### R functions to perform linear mixed model for repeated measures\r\n+####### on a multi var dataset using 3 files as used in W4M\r\n+##############################################################################################################\r\n+lmRepeated2FF <- function(ids, ifixfact, itime, isubject, ivd, ndim, nameVar=colnames(ids)[[ivd]],\r\n+ pvalCutof = 0.05, dffOption, visu = visu, tit = "", least.confounded = FALSE, outlier.limit = 3) {\r\n+ ### function to perform linear mixed model with 1 Fixed factor + Time + random factor subject\r\n+ ### based on lmerTest package providing functions giving the same results as SAS proc mixed\r\n+ options(scipen = 50, digits = 5)\r\n+\r\n+ if (!is.numeric(ids[[ivd]])) stop("Dependant variable is not numeric")\r\n+ if (!is.factor(ids[[ifixfact]])) stop("fixed factor is not a factor")\r\n+ if (!is.factor(ids[[itime]])) stop("Repeated factor is not a factor")\r\n+ if (!is.factor(ids[[isubject]])) stop("Random factor is not a factor")\r\n+\r\n+ ## factors\r\n+ time <- ids[[itime]]\r\n+ fixfact <- ids[[ifixfact]]\r\n+ subject <- ids[[isubject]]\r\n+ # dependant variables\r\n+ vd <- ids[[ivd]]\r\n+\r\n+ # argument of the function instead of re re-running ndim <- defColRes(ids, ifixfact, itime)\r\n+ # nfp : number of main factors + model infos (REML, varSubject) + normality test\r\n+ nfp <- ndim[1];\r\n+ # ncff number of comparison of the fixed factor\r\n+ nlff <- ndim[2]; ncff <- ndim[3]\r\n+ # nct number of comparison of the time factor\r\n+ nlt <- ndim[4]; nct <- ndim[5]\r\n+ # nci number of comparison of the interaction\r\n+ nli <- ndim[6]; nci <- ndim[7]\r\n+ # number of all lmer results\r\n+ nresT <- ncff + nct + nci\r\n+ ## initialization of the result vector (1 line)\r\n+ ## 4 * because nresf for : pvalues + Etimates + lower CI + Upper CI\r\n+ res <- data.frame(array(rep(NA, (nfp + 4 * nresT))))\r\n+ colnames(res)[1] <- "resultLM"\r\n+\r\n+ ### if at least one subject have data for only 1 time, mixed model is not possible and variable must be skip\r\n+ ### after excluding NA, table function is used to seek subjects with only 1 data\r\n+ ids <- ids[!is.na(ids[[ivd]]), ]\r\n+ skip <- length(which(table(ids[[isubject]]) == 1))\r\n+\r\n+ if (skip == 0) {\r\n+\r\n+ mfl <- lmer(vd ~ time + fixfact + time:fixfact + (1 | subject), ids) # lmer remix\r\n+\r\n+ rsum <- summary(mfl, ddf = dffOption)\r\n+ ## test Shapiro Wilks on the residus of the model\r\n+ rShapiro <- shapiro.test(rsum$residuals)\r\n+ raov <- anova(mfl, ddf = dffOption)\r\n+ dlsm1 <- data.frame(difflsmeans(mfl, test.effs = NULL))\r\n+ ddlsm1 <- dlsm1\r\n+ ## save rownames and factor names\r\n+ rn <- rownames(ddlsm1)\r\n+ fn <- ddlsm1[, c(1, 2)]\r\n+ ## writing the results on a single line\r\n+ namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesFactPval <- paste("pvalue ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesInter <- rownames(ddlsm1)[-c(1:(nct + ncff))]\r\n+ namesEstimate <- paste("estimate ", namesInter)\r\n+ namespvalues <- paste("pvalue ", namesInter)\r\n+ namesFactprinc <- c("pval_time", "pval_trt", "pval_inter")\r\n+ namesFactEstim <- paste("estimate ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+\r\n+ namesFactLowerCI <- paste("lowerCI ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesLowerCI <- paste("lowerCI ", namesInter, sep = "")\r\n+\r\n+ namesFactUpperCI <- paste("UpperCI ", rownames(ddlsm1)[c(1:(nct + ncff))], sep = "")\r\n+ namesUpperCI <- paste("UpperCI ", namesInter, sep = "")\r\n+\r\n+\r\n+ ### lmer results on 1 vector row\r\n+ # pvalue of shapiro Wilks test of the residuals\r\n+ res[1, ] <- rShapiro$p.value; rownames(res)[1] <- "Shapiro.pvalue.residuals"\r\n+ res[2, ] <- rsum$varcor$subject[1] ;rownames(res)[2] <- "Subject.Variance"\r\n+ res[3, ] <- rsum$devcomp$cmp[7] ; rownames(res)[3] <- "REML"\r\n+ '..b'utof), c((ndim[1] + ndim[5] + 1):(ndim[1] + ndim[5] + ndim[3]))] <- NA\r\n+ ## interaction effect\r\n+ resLM[which(resLM[, firstpval + 2] > pvalCutof), c((lastpval + ndim[5] + ndim[3] + 1):(lastpval + ndim[5] + ndim[3] + ndim[7]))] <- 0\r\n+ resLM[which(resLM[, firstpval + 2] > pvalCutof), c((ndim[1] + ndim[5] + ndim[3] + 1):(ndim[1] + ndim[5] + ndim[3] + ndim[7]))] <- NA\r\n+ } else {\r\n+ ## time effect only\r\n+ resLM[which(resLM[, firstpval] > pvalCutof), c((lastpval + 1):(lastpval + ndim[5]))] <- 0\r\n+ resLM[which(resLM[, firstpval] > pvalCutof), c((firstpval + 1):(firstpval + ndim[5]))] <- NA\r\n+ }\r\n+\r\n+ ## for each variable, estimates plots are performed if at least one factor is significant after p-value correction\r\n+ pdf(pdfE, onefile = TRUE, height = 15, width = 30)\r\n+\r\n+ ## for each variable (in row)\r\n+ for (i in seq_len(nrow(resLM))) {\r\n+\r\n+ ## if any fixed factor + time factor\r\n+ if (ifixfact > 0)\r\n+\r\n+ ## if any main factor after p-value correction is significant -> plot estimates and time course\r\n+ if (length(which(resLM[i, c(4:6)] < pvalCutof)) > 0) {\r\n+\r\n+ ## Plot of time course by fixfact : data prep with factors and quantitative var to be plot\r\n+ subv <- dslm[, colnames(dslm) == rownames(resLM)[i]]\r\n+ subds <- data.frame(dslm[[ifixfact]], dslm[[itime]], dslm[[isubject]], subv)\r\n+ libvar <- c(fixfact, time, subject)\r\n+ colnames(subds) <- c(libvar, rownames(resLM)[i])\r\n+\r\n+ ## Plot of estimates with error bars for all fixed factors and interaction\r\n+ rddlsm1 <- t(resLM[i, ])\r\n+ pval <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "pvalue"]\r\n+ esti <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "estima"]\r\n+ loci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "lowerC"]\r\n+ upci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "UpperC"]\r\n+ rddlsm1 <- data.frame(pval, esti, loci, upci, factorRow)\r\n+ colnames(rddlsm1) <- c("p.value", "Estimate", "Lower.CI", "Upper.CI", colnames(factorRow))\r\n+ rownames(rddlsm1) <- labelRow\r\n+\r\n+ ## function for plotting these 2 graphs\r\n+ plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)\r\n+\r\n+ }\r\n+\r\n+ ## if only a time factor\r\n+ if (ifixfact == 0)\r\n+\r\n+ ## if time factor after p-value correction is significant -> plot time course\r\n+ if (length(which(resLM[i, 4] < pvalCutof)) > 0) {\r\n+\r\n+ ## Plot of time course : data prep with factors and quantitative var to be plot\r\n+ subv <- dslm[, colnames(dslm) == rownames(resLM)[i]]\r\n+ subds <- data.frame(dslm[[itime]], dslm[[isubject]], subv)\r\n+ libvar <- c(time, subject)\r\n+ colnames(subds) <- c(libvar, rownames(resLM)[i])\r\n+\r\n+ ## Plot of estimates with error bars for all fixed factors and interaction\r\n+ rddlsm1 <- t(resLM[i, ])\r\n+ pval <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "pvalue"]\r\n+ esti <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "estima"]\r\n+ loci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "lowerC"]\r\n+ upci <- rddlsm1[substr(rownames(rddlsm1), 1, 6) == "upperC"]\r\n+ rddlsm1 <- data.frame(pval, esti, loci, upci, factorRow)\r\n+ colnames(rddlsm1) <- c("p.value", "Estimate", "Lower.CI", "Upper.CI", colnames(factorRow))\r\n+ rownames(rddlsm1) <- labelRow\r\n+\r\n+ ## function for plotting these 2 graphs\r\n+ plot.res.Lmixed(rddlsm1, subds, title = rownames(resLM)[i], pvalCutof = pvalCutof)\r\n+\r\n+ }\r\n+\r\n+ }\r\n+ dev.off()\r\n+\r\n+ ## return result file with pvalues and estimates (exclude confidence interval used for plotting)\r\n+ iCI <- which(substr(colnames(resLM), 4, 7) == "erCI")\r\n+ resLM <- resLM[, - iCI]\r\n+ resLM <- cbind(varids, resLM)\r\n+ return(resLM)\r\n+}\r\n' |
| b |
| diff -r 000000000000 -r a4d89d47646f mixmodel_wrapper.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mixmodel_wrapper.R Mon May 16 09:25:01 2022 +0000 |
| [ |
| @@ -0,0 +1,200 @@ +#!/usr/bin/env Rscript + +library(lme4) ## mixed model computing +library(Matrix) +library(MASS) +library(lmerTest) ## computing pvalue and lsmeans from results of lme4 package +library(multtest) ## multiple testing +library(ggplot2) +library(gridExtra) +library(grid) + +source_local <- function(fname) { + argv <- commandArgs(trailingOnly = FALSE) + base_dir <- dirname(substring(argv[grep("--file=", argv)], 8)) + source(paste(base_dir, fname, sep = "/")) +} + +source_local("mixmodel_script.R") +source_local("diagmfl.R") + +parse_args <- function() { + args <- commandArgs() + start <- which(args == "--args")[1] + 1 + if (is.na(start)) { + return(list()) + } + seq_by2 <- seq(start, length(args), by = 2) + result <- as.list(args[seq_by2 + 1]) + names(result) <- args[seq_by2] + return(result) +} + +argVc <- unlist(parse_args()) + +##------------------------------ +## Initializing +##------------------------------ + +## options +##-------- + +strAsFacL <- options()$stringsAsFactors +options(stringsAsFactors = FALSE) + +## constants +##---------- + +modNamC <- "mixmodel" ## module name + +topEnvC <- environment() +flagC <- "\n" + +## functions +##---------- + +flgF <- function(tesC, + envC = topEnvC, + txtC = NA) { ## management of warning and error messages + + tesL <- eval(parse(text = tesC), envir = envC) + + if (!tesL) { + + sink(NULL) + stpTxtC <- ifelse(is.na(txtC), + paste0(tesC, " is FALSE"), + txtC) + + stop(stpTxtC, + call. = FALSE) + + } + +} ## flgF + +## log file +##--------- + +sink(argVc["information"]) + +cat("\nStart of the '", modNamC, "' Galaxy module call: ", format(Sys.time(), "%a %d %b %Y %X"), "\n", sep = "") +cat("\nParameters used:\n\n") +print(argVc) +cat("\n\n") + +## loading +##-------- + +datMN <- t(as.matrix(read.table(argVc["dataMatrix_in"], + check.names = FALSE, + header = TRUE, + row.names = 1, + sep = "\t"))) + +samDF <- read.table(argVc["sampleMetadata_in"], + check.names = FALSE, + header = TRUE, + row.names = 1, + sep = "\t") + +varDF <- read.table(argVc["variableMetadata_in"], + check.names = FALSE, + header = TRUE, + row.names = 1, + sep = "\t") + + +## checking +##--------- + +flgF("identical(rownames(datMN), rownames(samDF))", txtC = "Column names of the dataMatrix are not identical to the row names of the sampleMetadata; check your data with the 'Check Format' module in the 'Quality Control' section") +flgF("identical(colnames(datMN), rownames(varDF))", txtC = "Row names of the dataMatrix are not identical to the row names of the variableMetadata; check your data with the 'Check Format' module in the 'Quality Control' section") + +flgF("argVc['time'] %in% colnames(samDF)", txtC = paste0("Required time factor '", argVc["time"], "' could not be found in the column names of the sampleMetadata")) +flgF("argVc['subject'] %in% colnames(samDF)", txtC = paste0("Required subject factor '", argVc["subject"], "' could not be found in the column names of the sampleMetadata")) + +flgF("mode(samDF[, argVc['time']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc["time"], "' column of the sampleMetadata should contain either number only, or character only")) +flgF("mode(samDF[, argVc['subject']]) %in% c('character', 'numeric')", txtC = paste0("The '", argVc["subject"], "' column of the sampleMetadata should contain either number only, or character only")) + +flgF("argVc['adjC'] %in% c('holm', 'hochberg', 'hommel', 'bonferroni', 'BH', 'BY', 'fdr', 'none')") +flgF("argVc['trf'] %in% c('none', 'log10', 'log2')") + +flgF("0 <= as.numeric(argVc['thrN']) && as.numeric(argVc['thrN']) <= 1", txtC = "(corrected) p-value threshold must be between 0 and 1") +flgF("argVc['diaR'] %in% c('no', 'yes')") + + +##------------------------------ +## Formating +##------------------------------ + +if (argVc["dff"] == "Satt") { + dffmeth <- "Satterthwaite" +} else { + dffmeth <- "Kenward-Roger" +} + + +##------------------------------ +## Computation +##------------------------------ + + +varDFout <- lmixedm(datMN = datMN, + samDF = samDF, + varDF = varDF, + fixfact = argVc["fixfact"], + time = argVc["time"], + subject = argVc["subject"], + logtr = argVc["trf"], + pvalCutof = argVc["thrN"], + pvalcorMeth = argVc["adjC"], + dffOption = dffmeth, + visu = argVc["diaR"], + least.confounded = FALSE, + outlier.limit = 3, + pdfC = argVc["out_graph_pdf"], + pdfE = argVc["out_estim_pdf"] + ) + + + + +##------------------------------ +## Rounding +##------------------------------ + +if (argVc["rounding"] == "yes") { + varDFout[, which(!(colnames(varDFout) %in% colnames(varDF)))] <- apply(varDFout[, which(!(colnames(varDFout) %in% colnames(varDF)))], 2, round, digits = as.numeric(argVc["decplaces"])) +} + +##------------------------------ +## Ending +##------------------------------ + + +## saving +##-------- + +varDFout <- cbind.data.frame(variableMetadata = rownames(varDFout), + varDFout) + +write.table(varDFout, + file = argVc["variableMetadata_out"], + quote = FALSE, + row.names = FALSE, + sep = "\t") + +## closing +##-------- + +cat("\n\nEnd of '", modNamC, "' Galaxy module call: ", + as.character(Sys.time()), "\n", sep = "") +cat("\nInformation about R (version, Operating System, attached or loaded packages):\n\n") +sessionInfo() + +sink() + +options(stringsAsFactors = strAsFacL) + +rm(list = ls()) |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/OneFactor_dataMatrix.tsv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/OneFactor_dataMatrix.tsv Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,23 @@ +variableMetadata s1D2 s2D2 s3D2 s1D3 s2D3 s3D3 s1D4 s2D4 s3D4 s1D5 s2D5 s3D5 s1D6 s2D6 s3D6 +V1 1.06E+06 1.07E+06 8.93E+05 9.11E+05 8.07E+05 6.90E+05 7.07E+05 7.66E+05 6.82E+05 6.25E+05 7.25E+05 6.78E+05 7.59E+05 6.70E+05 7.41E+05 +V2 4.67E+06 4.78E+06 4.17E+06 4.84E+06 4.50E+06 4.41E+06 6.96E+06 7.05E+06 6.65E+06 4.12E+06 4.42E+06 4.39E+06 1.26E+07 1.28E+07 1.32E+07 +V3 9.03E+05 8.44E+05 7.31E+05 9.97E+05 7.97E+05 7.47E+05 7.18E+05 8.28E+05 7.81E+05 9.08E+05 8.35E+05 8.37E+05 1.12E+06 1.11E+06 1.16E+06 +V4 1.62E+05 1.27E+05 1.21E+05 1.41E+05 8.48E+04 1.06E+05 9.72E+04 1.52E+05 1.27E+05 1.29E+05 1.42E+05 1.17E+05 1.72E+05 1.67E+05 1.87E+05 +V5 5.22E+04 6.20E+04 5.55E+04 4.78E+04 3.41E+04 5.84E+04 6.30E+04 1.11E+05 1.16E+05 4.59E+04 3.34E+04 4.08E+04 9.14E+04 6.78E+04 1.08E+05 +V6 1.68E+05 1.83E+05 2.76E+05 2.30E+05 1.62E+05 3.13E+05 4.47E+05 5.34E+05 6.59E+05 2.38E+05 2.30E+05 1.81E+05 4.93E+05 3.83E+05 6.30E+05 +V7 8.85E+06 8.39E+06 8.56E+06 5.51E+06 5.07E+06 4.90E+06 2.90E+06 3.10E+06 3.01E+06 2.36E+06 2.40E+06 2.30E+06 5.06E+06 4.94E+06 4.93E+06 +V8 2.53E+05 2.26E+05 1.90E+05 2.75E+05 2.67E+05 2.22E+05 2.12E+05 2.43E+05 1.72E+05 1.51E+05 1.56E+05 1.50E+05 2.49E+05 2.22E+05 2.48E+05 +V9 7.91E+05 7.16E+05 6.57E+05 6.08E+05 4.62E+05 4.88E+05 5.49E+05 6.13E+05 5.54E+05 7.56E+05 7.81E+05 7.59E+05 1.47E+06 1.46E+06 1.60E+06 +V10 1.76E+04 1.79E+04 1.86E+04 4.00E+04 4.27E+04 3.30E+04 5.34E+04 6.26E+04 6.09E+04 1.04E+05 8.24E+04 7.16E+04 2.09E+05 2.07E+05 2.19E+05 +V11 1.03E+05 9.28E+04 9.34E+04 8.70E+04 5.73E+04 8.05E+04 1.82E+06 1.98E+06 1.76E+06 8.41E+04 7.80E+04 7.31E+04 5.62E+04 5.01E+04 7.48E+04 +V12 4.48E+04 6.71E+04 4.27E+04 3.54E+04 3.04E+04 2.73E+04 2.47E+04 2.85E+04 3.26E+04 2.11E+04 1.55E+04 2.71E+04 3.44E+04 4.14E+04 5.87E+04 +V13 4.48E+04 6.71E+04 4.27E+04 3.54E+04 3.04E+04 2.73E+04 2.47E+04 2.85E+04 3.26E+04 2.11E+04 1.55E+04 2.71E+04 3.44E+04 4.14E+04 5.87E+04 +V14 1.16E+04 6.64E+03 1.29E+04 2.59E+04 1.55E+04 2.21E+04 3.15E+04 4.05E+04 3.10E+04 4.37E+04 2.90E+04 2.67E+04 6.62E+04 1.00E+05 8.75E+04 +V15 3.04E+05 3.05E+05 2.69E+05 2.57E+05 1.99E+05 2.16E+05 1.97E+05 1.97E+05 2.03E+05 2.65E+05 2.71E+05 2.60E+05 1.67E+05 1.89E+05 1.99E+05 +V16 2.69E+05 2.45E+05 2.18E+05 1.23E+05 1.27E+05 1.04E+05 2.91E+05 2.82E+05 2.79E+05 2.37E+05 2.11E+05 2.18E+05 5.60E+05 5.89E+05 6.08E+05 +V17 9.30E+04 8.55E+04 7.01E+04 4.84E+04 4.06E+04 3.14E+04 8.77E+04 1.03E+05 9.48E+04 9.88E+04 6.91E+04 7.81E+04 3.19E+05 3.25E+05 3.76E+05 +V18 2.30E+04 3.87E+04 2.01E+04 9.56E+03 1.27E+04 1.69E+04 2.38E+04 2.66E+04 2.99E+04 5.98E+04 4.93E+04 5.13E+04 1.08E+05 1.30E+05 1.36E+05 +V19 6.25E+04 7.43E+04 4.71E+04 1.62E+04 1.05E+04 1.57E+04 3.73E+04 3.08E+04 3.67E+04 3.50E+04 3.90E+04 2.78E+04 4.64E+05 1.60E+05 1.75E+05 +V20 1.72E+04 1.58E+04 1.35E+04 1.20E+04 1.02E+04 8.85E+03 1.35E+04 1.71E+04 1.62E+04 1.64E+04 8.41E+03 6.02E+03 5.07E+04 7.33E+04 6.57E+04 +V21 1.86E+06 1.90E+06 1.90E+06 2.05E+06 2.06E+06 2.07E+06 1.97E+06 1.94E+06 1.90E+06 1.64E+06 1.67E+06 1.78E+06 1.65E+06 1.73E+06 1.77E+06 +V22 1.52E+04 8.17E+03 1.37E+04 1.39E+04 1.34E+04 1.36E+04 1.40E+04 1.67E+04 1.60E+04 1.33E+04 1.69E+04 1.20E+04 1.05E+04 1.15E+04 1.11E+04 |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/OneFactor_sampleMetadata.tsv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/OneFactor_sampleMetadata.tsv Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,16 @@ +sampleMetadata time subject +s1D2 D2 s1 +s2D2 D2 s2 +s3D2 D2 s3 +s1D3 D3 s1 +s2D3 D3 s2 +s3D3 D3 s3 +s1D4 D4 s1 +s2D4 D4 s2 +s3D4 D4 s3 +s1D5 D5 s1 +s2D5 D5 s2 +s3D5 D5 s3 +s1D6 D6 s1 +s2D6 D6 s2 +s3D6 D6 s3 |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/OneFactor_variableMetadata.tsv --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/OneFactor_variableMetadata.tsv Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,23 @@ +variableMetadata name +V1 name1 +V2 name2 +V3 name3 +V4 name4 +V5 name5 +V6 name6 +V7 name7 +V8 name8 +V9 name9 +V10 name10 +V11 name11 +V12 name12 +V13 name13 +V14 name14 +V15 name15 +V16 name16 +V17 name17 +V18 name18 +V19 name19 +V20 name20 +V21 name21 +V22 name22 |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/TwoFactor_dataMatrix.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/TwoFactor_dataMatrix.txt Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,6 @@ +variable S01T0 S01T1 S01T2 S01T3 S02T0 S02T1 S02T2 S02T3 S03T0 S03T1 S03T2 S03T3 S04T0 S04T1 S04T2 S04T3 S05T0 S05T1 S05T2 S05T3 S06T0 S06T1 S06T2 S06T3 S07T0 S07T1 S07T2 S07T3 S08T0 S08T1 S08T2 S08T3 S09T0 S09T1 S09T2 S09T3 S10T0 S10T1 S10T2 S10T3 S11T0 S11T1 S11T2 S11T3 S12T0 S12T1 S12T2 S12T3 S13T0 S13T1 S13T2 S13T3 S14T0 S14T1 S14T2 S14T3 S15T0 S15T1 S15T2 S15T3 S16T0 S16T1 S16T2 S16T3 S17T0 S17T1 S17T2 S17T3 S18T0 S18T1 S18T2 S18T3 +v0.4425_385 126.39 234.11 162.68 115.29 172.44 171.53 180.31 114.94 224.14 195.09 189.45 122.79 212.02 243.48 188.42 138.91 241.59 192.95 214.74 98.69 270.03 176.79 208.73 78.04 247.82 176.51 175.78 98.94 260.63 181.93 178.06 111.58 218.95 176.24 200.29 112.77 140.23 157.38 124.62 114.34 147.23 165.29 125.85 128.41 126.35 139.75 127.95 74.86 124.92 155.12 109.4 86.88 145.57 142 135.29 83.08 116.95 142.28 113.43 118.26 80.81 146.03 130.52 103.4 79.58 168 109.35 88.19 138.12 166.61 115.05 83.05 +v0.4498_625 108.9 133.53 128.66 205.38 121.48 112.23 142.97 199.35 127.44 102.45 173.31 252.14 123.96 96.93 118.7 234.54 146.33 109.1 178.69 233.61 154.27 138.81 157.97 213.18 176.06 108.79 142.11 116.14 138.78 131.42 172.65 214.14 159.19 102.17 137.67 211.04 120.1 69.85 65.75 92.68 121.22 70.04 76.06 91.89 81.18 73.14 80.93 92.69 110.68 66.84 79.72 92.89 93.5 62.98 87.39 98.32 84.59 65.18 99.23 36.94 85.1 80.49 95.17 115.22 113.47 65.32 53.29 113.25 106.68 74.04 66.5 88.96 +v0.4711_463 124.1 166.84 156.18 82.57 158.4 190.09 152.89 100.76 155.75 174.85 153.14 112.76 171.54 150.37 163.49 116.06 198.78 161.53 247.61 122.94 167.38 186.12 194.59 93.05 208.24 150.15 140.97 55.62 200.32 138.59 145.99 89.44 192.49 131.44 145.89 93.38 102.94 67.11 83.99 86.86 110.09 98.61 77.91 100.75 136.32 92.56 86.79 75.62 144.39 65.27 110.24 94.26 136.48 77.62 95.43 60.3 135.98 87.54 65.37 92.04 112.09 58.9 92.35 72.92 132.84 97.12 80.71 79.01 144.75 85.78 73.61 83.68 +v0.4715_659 688.13 739.41 529.62 354.46 828.48 660.76 545.32 318.51 798.21 714.6 617.95 400.94 815.2 696.85 611.31 354.28 909.08 792.05 650.85 386.48 897.32 729.1 638.99 420.19 967.92 765.14 579.41 369.82 929.8 759.62 617.5 339.81 989.48 781.37 606.77 386.24 289.67 359.43 240.13 294.74 314.71 338.43 203.29 264.08 349.01 365.76 245.06 283.41 327.8 354.82 228.03 263.83 300.57 362.51 242.41 282.21 280.17 344.99 253.88 273.4 304.37 325.89 292.93 292.98 317.8 318.16 205.93 245.92 310.23 330.51 230.27 270.3 +v0.4766_757 117.8 110.31 114.56 99.99 144.39 155.33 114.47 122 147.16 132.82 129.78 94.86 128.26 161.38 119.12 96.73 177.26 134.94 149.39 134.12 163.22 146.12 148.26 108.68 159.39 140.75 119.11 65.83 170.45 131.79 163.58 104.63 167.33 112.88 142.58 85.38 80.53 84.3 42.33 62.04 73.74 97.78 58 79.33 91.7 82.39 69.71 83.58 77.86 65.56 74.11 82.9 90.53 75.27 80.46 81.88 92.77 52.94 67.99 84.21 91.27 69.54 77.13 85.27 93.18 83.68 57.01 65.09 92.12 82.63 70.9 75.53 |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/TwoFactor_sampleMetadata.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/TwoFactor_sampleMetadata.txt Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,73 @@ +sample sujet phenotype time +S01T0 S01 WT T0 +S01T1 S01 WT T1 +S01T2 S01 WT T2 +S01T3 S01 WT T3 +S02T0 S02 WT T0 +S02T1 S02 WT T1 +S02T2 S02 WT T2 +S02T3 S02 WT T3 +S03T0 S03 WT T0 +S03T1 S03 WT T1 +S03T2 S03 WT T2 +S03T3 S03 WT T3 +S04T0 S04 WT T0 +S04T1 S04 WT T1 +S04T2 S04 WT T2 +S04T3 S04 WT T3 +S05T0 S05 WT T0 +S05T1 S05 WT T1 +S05T2 S05 WT T2 +S05T3 S05 WT T3 +S06T0 S06 WT T0 +S06T1 S06 WT T1 +S06T2 S06 WT T2 +S06T3 S06 WT T3 +S07T0 S07 WT T0 +S07T1 S07 WT T1 +S07T2 S07 WT T2 +S07T3 S07 WT T3 +S08T0 S08 WT T0 +S08T1 S08 WT T1 +S08T2 S08 WT T2 +S08T3 S08 WT T3 +S09T0 S09 WT T0 +S09T1 S09 WT T1 +S09T2 S09 WT T2 +S09T3 S09 WT T3 +S10T0 S10 MT T0 +S10T1 S10 MT T1 +S10T2 S10 MT T2 +S10T3 S10 MT T3 +S11T0 S11 MT T0 +S11T1 S11 MT T1 +S11T2 S11 MT T2 +S11T3 S11 MT T3 +S12T0 S12 MT T0 +S12T1 S12 MT T1 +S12T2 S12 MT T2 +S12T3 S12 MT T3 +S13T0 S13 MT T0 +S13T1 S13 MT T1 +S13T2 S13 MT T2 +S13T3 S13 MT T3 +S14T0 S14 MT T0 +S14T1 S14 MT T1 +S14T2 S14 MT T2 +S14T3 S14 MT T3 +S15T0 S15 MT T0 +S15T1 S15 MT T1 +S15T2 S15 MT T2 +S15T3 S15 MT T3 +S16T0 S16 MT T0 +S16T1 S16 MT T1 +S16T2 S16 MT T2 +S16T3 S16 MT T3 +S17T0 S17 MT T0 +S17T1 S17 MT T1 +S17T2 S17 MT T2 +S17T3 S17 MT T3 +S18T0 S18 MT T0 +S18T1 S18 MT T1 +S18T2 S18 MT T2 +S18T3 S18 MT T3 |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/TwoFactor_variableMetadata.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/TwoFactor_variableMetadata.txt Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,6 @@ +variable num +v0.4425_385 1 +v0.4498_625 2 +v0.4711_463 3 +v0.4715_659 4 +v0.4766_757 5 |
| b |
| diff -r 000000000000 -r a4d89d47646f test-data/mixmodel_TwoFactor_variableMetadata.txt --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test-data/mixmodel_TwoFactor_variableMetadata.txt Mon May 16 09:25:01 2022 +0000 |
| b |
| @@ -0,0 +1,6 @@ +variableMetadata num Shapiro.pvalue.residuals Subject.Variance REML pval_time pval_trt pval_inter pvalue.timeT0...timeT1 pvalue.timeT0...timeT2 pvalue.timeT0...timeT3 pvalue.timeT1...timeT2 pvalue.timeT1...timeT3 pvalue.timeT2...timeT3 pvalue.fixfactMT...fixfactWT pvalue..timeT0.fixfactMT...timeT1.fixfactMT pvalue..timeT0.fixfactMT...timeT2.fixfactMT pvalue..timeT0.fixfactMT...timeT3.fixfactMT pvalue..timeT0.fixfactMT...timeT0.fixfactWT pvalue..timeT0.fixfactMT...timeT1.fixfactWT pvalue..timeT0.fixfactMT...timeT2.fixfactWT pvalue..timeT0.fixfactMT...timeT3.fixfactWT pvalue..timeT1.fixfactMT...timeT2.fixfactMT pvalue..timeT1.fixfactMT...timeT3.fixfactMT pvalue..timeT1.fixfactMT...timeT0.fixfactWT pvalue..timeT1.fixfactMT...timeT1.fixfactWT pvalue..timeT1.fixfactMT...timeT2.fixfactWT pvalue..timeT1.fixfactMT...timeT3.fixfactWT pvalue..timeT2.fixfactMT...timeT3.fixfactMT pvalue..timeT2.fixfactMT...timeT0.fixfactWT pvalue..timeT2.fixfactMT...timeT1.fixfactWT pvalue..timeT2.fixfactMT...timeT2.fixfactWT pvalue..timeT2.fixfactMT...timeT3.fixfactWT pvalue..timeT3.fixfactMT...timeT0.fixfactWT pvalue..timeT3.fixfactMT...timeT1.fixfactWT pvalue..timeT3.fixfactMT...timeT2.fixfactWT pvalue..timeT3.fixfactMT...timeT3.fixfactWT pvalue..timeT0.fixfactWT...timeT1.fixfactWT pvalue..timeT0.fixfactWT...timeT2.fixfactWT pvalue..timeT0.fixfactWT...timeT3.fixfactWT pvalue..timeT1.fixfactWT...timeT2.fixfactWT pvalue..timeT1.fixfactWT...timeT3.fixfactWT pvalue..timeT2.fixfactWT...timeT3.fixfactWT estimate.timeT0...timeT1 estimate.timeT0...timeT2 estimate.timeT0...timeT3 estimate.timeT1...timeT2 estimate.timeT1...timeT3 estimate.timeT2...timeT3 estimate.fixfactMT...fixfactWT estimate..timeT0.fixfactMT...timeT1.fixfactMT estimate..timeT0.fixfactMT...timeT2.fixfactMT estimate..timeT0.fixfactMT...timeT3.fixfactMT estimate..timeT0.fixfactMT...timeT0.fixfactWT estimate..timeT0.fixfactMT...timeT1.fixfactWT estimate..timeT0.fixfactMT...timeT2.fixfactWT estimate..timeT0.fixfactMT...timeT3.fixfactWT estimate..timeT1.fixfactMT...timeT2.fixfactMT estimate..timeT1.fixfactMT...timeT3.fixfactMT estimate..timeT1.fixfactMT...timeT0.fixfactWT estimate..timeT1.fixfactMT...timeT1.fixfactWT estimate..timeT1.fixfactMT...timeT2.fixfactWT estimate..timeT1.fixfactMT...timeT3.fixfactWT estimate..timeT2.fixfactMT...timeT3.fixfactMT estimate..timeT2.fixfactMT...timeT0.fixfactWT estimate..timeT2.fixfactMT...timeT1.fixfactWT estimate..timeT2.fixfactMT...timeT2.fixfactWT estimate..timeT2.fixfactMT...timeT3.fixfactWT estimate..timeT3.fixfactMT...timeT0.fixfactWT estimate..timeT3.fixfactMT...timeT1.fixfactWT estimate..timeT3.fixfactMT...timeT2.fixfactWT estimate..timeT3.fixfactMT...timeT3.fixfactWT estimate..timeT0.fixfactWT...timeT1.fixfactWT estimate..timeT0.fixfactWT...timeT2.fixfactWT estimate..timeT0.fixfactWT...timeT3.fixfactWT estimate..timeT1.fixfactWT...timeT2.fixfactWT estimate..timeT1.fixfactWT...timeT3.fixfactWT estimate..timeT2.fixfactWT...timeT3.fixfactWT +v0.4425_385 1 0.004001 0 -140.931331 0 0 0.000264 0.191909 0.277993 0 0.018695 0 0 0 0.001835 0.869646 0.006728 0 0 0 0.238471 0.002992 0 0.000044 0.003836 0.009486 0.000036 0.004238 0 0 0 0.180279 0 0 0 0.112058 0.170546 0.091717 0 0.745427 0 0 -0.030835 0.02558 0.192247 0.056415 0.223081 0.166666 -0.148875 -0.107496 -0.005448 0.092614 -0.252536 -0.206709 -0.195927 0.039343 0.102048 0.200111 -0.14504 -0.099213 -0.088431 0.146839 0.098062 -0.247088 -0.201261 -0.19048 0.044791 -0.34515 -0.299323 -0.288542 -0.053272 0.045827 0.056608 0.291879 0.010781 0.246052 0.23527 +v0.4498_625 2 0.000002 0 -120.084702 0 0 0.000801 0.000045 0.135374 0.047823 0.005607 0 0.000776 0 0.000128 0.004174 0.135488 0.000855 0.170794 0.000053 0 0.272071 0.012599 0 0.000001 0 0 0.149351 0 0.000049 0 0 0.000005 0.00515 0 0 0.038427 0.408218 0.000047 0.00448 0 0.000769 0.120489 0.041611 -0.055515 -0.078879 -0.176005 -0.097126 -0.249448 0.158724 0.115616 0.058829 -0.136156 -0.053901 -0.16855 -0.306015 -0.043108 -0.099895 -0.29488 -0.212625 -0.327274 -0.46474 -0.056787 -0.251772 -0.169517 -0.284166 -0.421632 -0.194985 -0.11273 -0.227379 -0.364845 0.082255 -0.032394 -0.16986 -0.114649 -0.252114 -0.137466 +v0.4711_463 3 0.474695 0.000093 -135.57425 0 0 0.000006 0.000013 0.000128 0 0.488485 0.000043 0.000004 0 0 0.000004 0.000001 0.000295 0.006152 0.002433 0.000335 0.51173 0.755875 0 0 0 0.046325 0.729076 0 0 0 0.173445 0 0 0 0.089918 0.319884 0.50045 0 0.745698 0 0 0.11757 0.100696 0.22644 -0.016874 0.108869 0.125743 -0.194248 0.200783 0.178185 0.190094 -0.13207 -0.097712 -0.108863 0.130715 -0.022598 -0.010689 -0.332853 -0.298495 -0.309645 -0.070068 0.011909 -0.310255 -0.275897 -0.287048 -0.04747 -0.322164 -0.287806 -0.298957 -0.059379 0.034358 0.023208 0.262785 -0.01115 0.228427 0.239578 +v0.4715_659 4 0.623926 0.000349 -238.084202 0 0 0 0.219999 0 0 0 0 0 0 0.001773 0 0.000326 0 0 0 0.000025 0 0 0 0 0 0.070602 0.000027 0 0 0 0 0 0 0 0 0.000006 0 0 0 0 0 0.012039 0.138015 0.211119 0.125976 0.19908 0.073104 -0.326171 -0.045354 0.116577 0.053052 -0.44462 -0.375188 -0.285168 -0.075434 0.161931 0.098406 -0.399266 -0.329834 -0.239814 -0.03008 -0.063525 -0.561197 -0.491765 -0.401744 -0.192011 -0.497672 -0.42824 -0.33822 -0.128486 0.069432 0.159453 0.369186 0.09002 0.299754 0.209734 +v0.4766_757 5 0.005552 0.000365 -149.179654 0.00001 0 0.00028 0.015826 0.000092 0.000002 0.083272 0.004665 0.23638 0 0.064552 0.000169 0.102895 0 0 0 0.062636 0.033569 0.819707 0 0 0 0.000447 0.019489 0 0 0 0 0 0 0 0.000898 0.106307 0.055879 0 0.755242 0.000055 0.000152 0.052746 0.090042 0.115324 0.037296 0.062579 0.025283 -0.224878 0.056412 0.121655 0.049578 -0.240112 -0.191032 -0.181683 -0.059041 0.065243 -0.006834 -0.296524 -0.247444 -0.238095 -0.115453 -0.072077 -0.361767 -0.312687 -0.303338 -0.180696 -0.28969 -0.24061 -0.231262 -0.10862 0.04908 0.058428 0.181071 0.009349 0.131991 0.122642 |