Contributions icon indicating copy to clipboard operation
Contributions copied to clipboard

MAPFX

Open HsiaoChiLiao opened this issue 1 year ago • 9 comments

Update the following URL to point to the GitHub repository of the package you wish to submit to Bioconductor

  • Repository: https://github.com/HsiaoChiLiao/MAPFX

Confirm the following by editing each check box to '[x]'

  • [x] I understand that by submitting my package to Bioconductor, the package source and all review commentary are visible to the general public.

  • [x] I have read the Bioconductor Package Submission instructions. My package is consistent with the Bioconductor Package Guidelines.

  • [x] I understand Bioconductor Package Naming Policy and acknowledge Bioconductor may retain use of package name.

  • [x] I understand that a minimum requirement for package acceptance is to pass R CMD check and R CMD BiocCheck with no ERROR or WARNINGS. Passing these checks does not result in automatic acceptance. The package will then undergo a formal review and recommendations for acceptance regarding other Bioconductor standards will be addressed.

  • [x] My package addresses statistical or bioinformatic issues related to the analysis and comprehension of high throughput genomic data.

  • [x] I am committed to the long-term maintenance of my package. This includes monitoring the support site for issues that users may have, subscribing to the bioc-devel mailing list to stay aware of developments in the Bioconductor community, responding promptly to requests for updates from the Core team in response to changes in R or underlying software.

  • [x] I am familiar with the Bioconductor code of conduct and agree to abide by it.

I am familiar with the essential aspects of Bioconductor software management, including:

  • [x] The 'devel' branch for new packages and features.
  • [x] The stable 'release' branch, made available every six months, for bug fixes.
  • [x] Bioconductor version control using Git (optionally via GitHub).

For questions/help about the submission process, including questions about the output of the automatic reports generated by the SPB (Single Package Builder), please use the #package-submission channel of our Community Slack. Follow the link on the home page of the Bioconductor website to sign up.

HsiaoChiLiao avatar Mar 18 '24 12:03 HsiaoChiLiao

Hi @HsiaoChiLiao

Thanks for submitting your package. We are taking a quick look at it and you will hear back from us soon.

The DESCRIPTION file for this package is:

Package: MAPFX
Title: MAssively Parallel Flow cytometry Xplorer (MAPFX): A Toolbox for
        Analysing Data from the Massively-Parallel Cytometry
        Experiments
Version: 0.99.0
Authors@R:
	c(person("Hsiao-Chi", "Liao", email = "[email protected]", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-9586-1246")),
	  person("Agus", "Salim", role = c("ctb")), person("infinityFlow", role = c("ctb")))
Description: MAPFX is an end-to-end toolbox that pre-processes the raw data from MPC experiments (e.g., BioLegend's LEGENDScreen and BD Lyoplates assays), and further imputes the ‘missing’ infinity markers in the wells without those measurements. The pipeline starts by performing background correction on raw intensities to remove the noise from electronic baseline restoration and fluorescence compensation by adapting a normal-exponential convolution model. Unwanted technical variation, from sources such as well effects, is then removed using a log-normal model with plate, column, and row factors, after which infinity markers are imputed using the informative backbone markers as predictors. The completed dataset can then be used for clustering and other statistical analyses. Additionally, MAPFX can be used to normalise data from FFC assays as well.
Depends: R (>= 4.4.0)
License: GPL-2
Encoding: UTF-8
RoxygenNote: 7.3.1
Roxygen: list(markdown = TRUE)
biocViews: Software, FlowCytometry, CellBasedAssays, SingleCell, Proteomics, Clustering
Imports: flowCore, Biobase, stringr, uwot, iCellR, igraph,
        ggplot2, RColorBrewer, Rfast, ggrepel, ComplexHeatmap,
        circlize, glmnetUtils, e1071, xgboost, foreach, doParallel,
        parallel, pbapply, reshape2, gtools, utils, stats, cowplot, methods, grDevices, graphics
Suggests: knitr, rmarkdown
VignetteBuilder: knitr
NeedsCompilation: no
URL: https://github.com/HsiaoChiLiao/MAPFX
BugReports: https://github.com/HsiaoChiLiao/MAPFX/issues

bioc-issue-bot avatar Mar 18 '24 12:03 bioc-issue-bot

Please remove the compiled MAPFX_Vignette.html from the vignettes directory. This file should be build automatically through R CMD build .

Also

I'm getting an ERROR when I try to build locally

PreReview$ R CMD build MAPFX 
* checking for file 'MAPFX/DESCRIPTION' ... OK
* preparing 'MAPFX':
* checking DESCRIPTION meta-information ... OK
* installing the package to build vignettes
* creating vignettes ... ERROR
--- re-building ‘MAPFX_Vignette.Rmd’ using rmarkdown
pandoc: /Volumes/Transcend_4T/PhD_desktop/ProteinExpression/InFlow/Programming/package/2024/Vignette/graphs_vignettes/F1.exp_compu_pipeline.png: openBinaryFile: does not exist (No such file or directory)
Error: processing vignette 'MAPFX_Vignette.Rmd' failed with diagnostics:
pandoc document conversion failed with error 1
--- failed re-building ‘MAPFX_Vignette.Rmd’

SUMMARY: processing the following file failed:
  ‘MAPFX_Vignette.Rmd’

Error: Vignette re-building failed.
Execution halted

lshep avatar Mar 21 '24 14:03 lshep

@lshep Thank you! I've removed the .html file from the vignettes directory and also resolved the above error.

HsiaoChiLiao avatar Mar 22 '24 02:03 HsiaoChiLiao

You might considering hiding the start up messages when you load all the library calls in the vignette because it hides what libraries are actually being loaded.

None of your main functions are actually executed with any data in vignette, man pages, or tests. In its current state it is unacceptable as we can not see results or that any of your functions are working correctly.

lshep avatar Mar 27 '24 14:03 lshep

@lshep thank you for your message. The start-up messages of the library calls in the vignette have been hidden. Also, I added small example datasets for my two exported functions. Apologise that I did not do this in the first place because I was planning to wait until my another Data package published.

It will take around 3 mins to complete R CMD build MAPFX --resave-data.

Many thanks, Hsiao-Chi

HsiaoChiLiao avatar Mar 30 '24 13:03 HsiaoChiLiao

In your vignette, Please change getwd to use a temporary directory with tempdir() and Could you please add a section in your vignette introduction that provides motivation for inclusion in Bioconductor and when appropriate a review and comparison to existing Bioconductor packages with similar functionality or scope.

The note about --resave-data you should apply yourself. You should resave the objects using the compress=xz option.

lshep avatar Apr 04 '24 13:04 lshep

Hi @lshep I've changed getwd() to tempdir(), added a motivation section in my vignette, and resaved the data objects using the compress='xz' option.

Thank you, Hsiao-Chi

HsiaoChiLiao avatar Apr 07 '24 10:04 HsiaoChiLiao

Your package has been added to git.bioconductor.org to continue the pre-review process. A build report will be posted shortly. Please fix any ERROR and WARNING in the build report before a reviewer is assigned or provide a justification on why you feel the ERROR or WARNING should be granted an exception.

IMPORTANT: Please read this documentation for setting up remotes to push to git.bioconductor.org. All changes should be pushed to git.bioconductor.org moving forward. It is required to push a version bump to git.bioconductor.org to trigger a new build report.

Bioconductor utilized your github ssh-keys for git.bioconductor.org access. To manage keys and future access you may want to active your Bioconductor Git Credentials Account

bioc-issue-bot avatar Apr 10 '24 15:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

On one or more platforms, the build results were: "ERROR". This may mean there is a problem with the package that you need to fix. Or it may mean that there is a problem with the build system itself.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.0.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 10 '24 16:04 bioc-issue-bot

Received a valid push on git.bioconductor.org; starting a build for commit id: 3c18ec7bef8bd4168ef9b681998ed5a7909debc5

bioc-issue-bot avatar Apr 11 '24 23:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

Congratulations! The package built without errors or warnings on all platforms.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.1.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 12 '24 00:04 bioc-issue-bot

A reviewer has been assigned to your package for an indepth review. Please respond accordingly to any further comments from the reviewer.

bioc-issue-bot avatar Apr 12 '24 16:04 bioc-issue-bot

Package 'MAPFX' Review

Thank you for submitting your package to Bioconductor. The package passed check and build. However there are several things need to be fixed. Please try to answer the comments line by line when you are ready for a second review. Code: Note: please consider; Important: must be addressed.

The NAMESPACE file

  • [ ] Important: Selective imports using importFrom instead of import all with import.
    • in line 5 import(ComplexHeatmap)
    • in line 6 import(cowplot)
    • in line 7 import(doParallel)
    • in line 8 import(e1071)
    • in line 9 import(foreach)
    • in line 11 import(ggrepel)
    • in line 12 import(glmnetUtils)
    • in line 13 import(xgboost)

General package development

  • [ ] Important: Consider adding unit tests. We strongly encourage them. See http://bioconductor.org/developers/how-to/unitTesting-guidelines/.

R code

  • [ ] Important: warning, message, stop instead of cat and print outside of show methods.
    • In file R/022.bkc_pe.R:
      • at line 140 found ' message("Could not find enough cells (>=10) when used "mle.mean+3*mle.sd", so estimated alpha with the "largest 10" cells:\n", cat(do.call(c, few.cell.well.ls), sep=", "))'
  • [ ] NOTE: :: is not suggested in source code unless you can make sure all the packages are imported. Some people think it is better to keep ::. However, please be aware that you will need to manually double-check the imported items if you make any changes to the DESCRIPTION file during development. My suggestion is to remove one or two repetitions to trigger the dependency check.
    • In file R/051.imputation_bkb.predictors.R:
      • at line 67 found ' res <- apply(gtools::combinations(n,deg,variables,repeats.allowed=TRUE),1,paste,collapse="*")'
      • at line 124 found ' model <- do.call(function(...){e1071::svm(...,x=x[w,chans],y=x[w,yvar])},params)'
      • at line 139 found ' model <- do.call(xgboost::xgboost, args)'
      • at line 160 found ' x <- xgboost::xgb.Booster.complete(x)'
      • at line 161 found ' xgboost::xgb.parameters(x) <- list(nthread = 1)'
      • at line 220 found ' cl <- parallel::makeCluster(min(cores,length(unique(events.code))))'
      • at line 230 found ' mc.reset.stream <- parallel::mc.reset.stream'
      • at line 328 found ' cl <- parallel::makeCluster(cores)'
      • at line 335 found ' mc.reset.stream <- parallel::mc.reset.stream'
  • [ ] NOTE: Vectorize: for loops present, try to replace them by *apply funcitons.
    • In file R/011.fcs_to_rds.R:
      • at line 38 found ' for(i in seq_along(fs)){'
      • at line 95 found ' for(i in seq_along(fs)){'
    • In file R/012.fcs_to_rds_bkb.R:
      • at line 40 found ' for(i in seq_along(fs)){'
      • at line 95 found ' for(i in seq_along(fs)){'
      • at line 132 found ' for(i in seq_along(fs)){'
    • In file R/021.bkc_bkb.R:
      • at line 75 found ' for(m in seq_len(length(bkb.v))){'
      • at line 77 found ' for(w in seq_len(length(unique(metadata.cell$Well.lab)))){'
      • at line 137 found ' for(m in skip.phy.meas){'
      • at line 145 found ' for(w in seq_len(length(unique(metadata.cell$Well.lab)))){'
      • at line 159 found ' for(i in seq_along(marker.nam)){'
    • In file R/022.bkc_pe.R:
      • at line 78 found ' for(i in seq_along(para.df[,1])){ #1:nrow(para.df)'
      • at line 150 found ' for(i in seq_along(para.df$Well.lab)){'
    • In file R/023.bkc_bkb_ffc.R:
      • at line 74 found ' for(m in seq_len(length(bkb.v))){'
      • at line 76 found ' for(w in seq_len(length(unique(metadata.cell$Batch)))){'
      • at line 135 found ' for(m in skip.phy.meas){'
      • at line 143 found ' for(w in seq_len(length(unique(metadata.cell$Batch)))){'
      • at line 157 found ' for(i in seq_along(marker.nam)){'
    • In file R/031.initM.R:
      • at line 80 found ' for(chan in chans){'
      • at line 92 found ' for(i in seq_along(well.v)){'
      • at line 106 found ' for(i in seq_along(batch.v)){'
    • In file R/032.rmWellEffect.R:
      • at line 49 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 81 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 110 found ' for(i in seq_along(log.adj.data[1,])){ #1:ncol(log.adj.data)'
      • at line 161 found ' for(i in seq_along(pre.post)){ #pre & post'
      • at line 162 found ' for(j in seq_along(bio.sys)){ #bio & batch'
    • In file R/033.rmBatchEffect.R:
      • at line 47 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 80 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 108 found ' for(i in seq_along(log.adj.data[1,])){ #1:ncol(log.adj.data)'
      • at line 158 found ' for(i in seq_along(pre.post)){ #pre & post'
      • at line 159 found ' for(j in seq_along(bio.sys)){ #bio & batch'
    • In file R/051.imputation_bkb.predictors.R:
      • at line 198 found ' for(m in seq_len(length(models.use))){'
      • at line 295 found ' for(i in seq_along(regression_functions)){'
      • at line 381 found ' for(i in seq_along(models)){'
      • at line 392 found ' for(i in seq_along(models)){'
      • at line 468 found ' for(i in seq_along(models)){'
      • at line 479 found ' for(i in seq_along(models)){'
      • at line 553 found ' for(i in seq_along(models)){'
      • at line 564 found ' for(i in seq_along(models)){'
      • at line 605 found ' for(i in seq_along(preds)){'
      • at line 613 found ' for(j in seq_along(uniq.well)){'
    • In file R/060.cluster.analysis.R:
      • at line 43 found ' for(i in seq_along(preds)){'
  • [ ] Important: Use file.path to replace paste. Do not include the '/' in the file path.
    • In file R/011.fcs_to_rds.R:
      • at line 29 found ' fs <- read.flowSet(path = file.path(paste0(paths["input"], "/fcs/")))'
    • In file R/012.fcs_to_rds_bkb.R:
      • at line 29 found ' fs <- read.flowSet(path = file.path(paste0(paths["input"], "/fcs/")))'
    • In file R/060.cluster.analysis.R:
      • at line 108 found ' jpeg(filename = file.path(paths["graph"], paste0("/ClusterStructure_UMAP_", length(bkb.v), "bkb_colPhenog_600x750.jpeg")), height = 600, width = 750, res = 80)'
      • at line 127 found ' jpeg(filename = file.path(paths["graph"], paste0("/ClusterStructure_UMAP_", names(preds)[i], length(bkb.v), "bkb",'
    • In file R/070.cluster.analysis_bkbOnly.R:
      • at line 71 found ' jpeg(filename = file.path(paths["graph"], paste0("/ClusterStructure_UMAP_", length(bkb.v), "bkb_colPhenog_600x750.jpeg")), height = 600, width = 750, res = 80)'
  • [ ] Important: Remove unused code.
    • In file R/000.MapfxFFC.R:
      • at line 82 found ' #fcs_to_rds_bkb('
      • at line 83 found ' #paths=paths, file_meta=file_meta, MPC = FALSE)'
    • In file R/000.MapfxMPC.R:
      • at line 124 found ' # fcs_to_rds(paths=paths, file_meta=file_meta, yvar = yvar)'
      • at line 213 found ' # lapply('
      • at line 215 found ' # function(fun){'
      • at line 216 found ' # fun(x = NULL, params = NULL)'
      • at line 217 found ' # }'
      • at line 218 found ' # )'
      • at line 220 found ' # if(length(extra_args_regression_params) != length(regression_functions)){'
      • at line 221 found ' # stop("extra_args_regression_params and regression_functions should be lists of the same lengths")'
      • at line 222 found ' # }'
      • at line 303 found ' # message("\n\n\nCluster analysis with adjusted backbone markers...")'
      • at line 304 found ' # if(cluster.analysis.bkb == TRUE){'
      • at line 305 found ' # cluster.analysis.bkbOnly('
      • at line 306 found ' # paths=paths,'
      • at line 307 found ' # bkb.v=bkb.v,'
      • at line 308 found ' # plots = plots.cluster.analysis.bkb)'
      • at line 309 found ' # }'
    • In file R/011.fcs_to_rds.R:
      • at line 87 found ' # unique(fcs.raw.meta.df$Well.lab)'
      • at line 89 found ' # unique(ord.fcs.raw.meta.df$Well.lab)'
      • at line 117 found ' # unique(fcs.raw.meta.df$Well.lab)'
      • at line 119 found ' # unique(ord.fcs.raw.meta.df$Well.lab)'
      • at line 126 found ' # tail(ord.fcs.raw.meta.df); unique(ord.fcs.raw.meta.df$Well.lab)'
    • In file R/012.fcs_to_rds_bkb.R:
      • at line 87 found ' # unique(fcs.raw.meta.df$Well.lab)'
      • at line 89 found ' # unique(ord.fcs.raw.meta.df$Well.lab)'
      • at line 117 found ' # unique(fcs.raw.meta.df$Well.lab)'
      • at line 119 found ' # unique(ord.fcs.raw.meta.df$Well.lab)'
      • at line 150 found ' # unique(fcs.raw.meta.df$Well.lab)'
      • at line 152 found ' # unique(ord.fcs.raw.meta.df$Well.lab)'
    • In file R/021.bkc_bkb.R:
      • at line 43 found ' # xx = the observed intensity used for estimating mu and sig'
      • at line 54 found ' # xx = the observed intensity used for estimating alpha'
      • at line 56 found ' # exp(p[1]) = 1/exponential mean = 1/alpha'
      • at line 117 found ' # exp(p[1]) = 1/exponential mean = 1/alpha'
      • at line 121 found ' # print(paste0('well=',w,',alpha=',round(alpha[length(alpha)],1)))'
      • at line 162 found ' # set.seed(123)'
    • In file R/022.bkc_pe.R:
      • at line 27 found ' # pe.mean.sd=3, #mean+-3sd for extracting extremely large values'
      • at line 34 found ' # upper.quantile=0.9, #cells used for estimating parameter of signal'
      • at line 43 found ' # xx = the observed intensity used for estimating mu and sig'
      • at line 54 found ' # xx = the observed intensity used for estimating alpha'
      • at line 56 found ' # exp(p[1]) = 1/exponential mean = 1/alpha'
      • at line 78 found ' for(i in seq_along(para.df[,1])){ #1:nrow(para.df)'
      • at line 129 found ' # exp(p[1]) = 1/exponential mean = 1/alpha'
      • at line 138 found ' # print(paste0('Iteration=',i,', alpha=',round(alpha[length(alpha)],1)))'
    • In file R/023.bkc_bkb_ffc.R:
      • at line 43 found ' # xx = the observed intensity used for estimating mu and sig'
      • at line 54 found ' # xx = the observed intensity used for estimating alpha'
      • at line 56 found ' # exp(p[1]) = 1/exponential mean = 1/alpha'
      • at line 105 found ' # suppressWarnings('
      • at line 113 found ' # )'
      • at line 115 found ' # exp(p[1]) = 1/exponential mean = 1/alpha'
      • at line 119 found ' # print(paste0('batch=',w,',alpha=',round(alpha[length(alpha)],1)))'
      • at line 160 found ' # set.seed(123)'
    • In file R/031.initM.R:
      • at line 136 found ' n_neighbors = 15, min_dist = 0.2, metric = "euclidean", n_epochs = 2000) #, n_neighbors = 15, min_dist = 0.2, metric = "euclidean", n_epochs = 2000'
    • In file R/032.rmWellEffect.R:
      • at line 49 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 56 found ' # ln.sig.v[i] <- sd(lm.result$residuals)'
      • at line 81 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 108 found ' # glm.results.ls <- list()'
      • at line 110 found ' for(i in seq_along(log.adj.data[1,])){ #1:ncol(log.adj.data)'
      • at line 117 found ' # ln.sig.v[i] <- sd(lm.result$residuals)'
      • at line 214 found ' #display.brewer.pal(n = 8, name = 'Dark2')'
    • In file R/033.rmBatchEffect.R:
      • at line 47 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 55 found ' # ln.sig.v[i] <- sd(lm.result$residuals)'
      • at line 80 found ' for(i in seq_along(bkc.bkb[1,])){ #1:ncol(bkc.bkb)'
      • at line 106 found ' # glm.results.ls <- list()'
      • at line 108 found ' for(i in seq_along(log.adj.data[1,])){ #1:ncol(log.adj.data)'
      • at line 115 found ' # ln.sig.v[i] <- sd(lm.result$residuals)'
      • at line 207 found ' # col. anno for heatmap (be careful of the position of each variable)'
      • at line 208 found ' #display.brewer.pal(n = 8, name = 'Dark2')'
    • In file R/051.imputation_bkb.predictors.R:
      • at line 245 found ' # set.seed(123)'
      • at line 397 found ' ##parLapplyLB('
      • at line 484 found ' ##parLapplyLB('
      • at line 569 found ' ##parLapplyLB('
      • at line 678 found ' # legend.key.size = unit(1, 'cm'), #change legend key size'
      • at line 679 found ' # legend.key.height = unit(1, 'cm'), #change legend key height'
      • at line 680 found ' # legend.key.width = unit(1, 'cm'), #change legend key width'
      • at line 681 found ' # legend.title = element_text(size=10, face="bold"), #change legend title font size'
      • at line 682 found ' # legend.text = element_text(size=8, face="bold")'
    • In file R/060.cluster.analysis.R:
      • at line 44 found ' # message("Cluster analysis using imputed values from ", names(preds)[i], "...")'
  • [ ] NOTE: Avoid 'suppressWarnings'/'*Messages' if possible
    • In file R/021.bkc_bkb.R:
      • at line 89 found ' suppressWarnings('
      • at line 107 found ' suppressWarnings('
    • In file R/022.bkc_pe.R:
      • at line 92 found ' suppressWarnings('
      • at line 119 found ' suppressWarnings('
    • In file R/023.bkc_bkb_ffc.R:
      • at line 87 found ' suppressWarnings('
  • [ ] Important: Please consider to remove ';' in the code to avoid mutiple scripts in one line.
    • In file R/022.bkc_pe.R:
      • at line 38 found ' pe.raw <- as.matrix(rawInten[,"Legend"]); colnames(pe.raw) <- "Legend"'
  • [ ] Important: Please consider to add drop=FALSE to avoid the reduction of dimension for matrices and arrays. Ignore this if using datatable.
    • In file R/011.fcs_to_rds.R:
      • at line 88 found ' ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),]'
      • at line 90 found ' ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),]'
      • at line 118 found ' ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),]'
      • at line 120 found ' ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),]'
      • at line 132 found ' ord.fcs.raw.meta.df.out <- ord.fcs.raw.meta.df[,c(8,1,3,5,6,4,7)]'
    • In file R/012.fcs_to_rds_bkb.R:
      • at line 88 found ' ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),]'
      • at line 90 found ' ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),]'
      • at line 118 found ' ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),]'
      • at line 120 found ' ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),]'
      • at line 127 found ' ord.fcs.raw.meta.df.out <- ord.fcs.raw.meta.df[,c(8,1,3,5,6,4,7)]'
      • at line 151 found ' ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Batch),]'
      • at line 153 found ' ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Batch),]'
      • at line 159 found ' ord.fcs.raw.meta.df.out <- ord.fcs.raw.meta.df[,c(4,1,3)]'
    • In file R/021.bkc_bkb.R:
      • at line 38 found ' sel.bkb.raw <- rawInten[,match(bkb.v, colnames(rawInten))]'
    • In file R/023.bkc_bkb_ffc.R:
      • at line 38 found ' sel.bkb.raw <- rawInten[,match(bkb.v, colnames(rawInten))]'
    • In file R/031.initM.R:
      • at line 44 found ' sel.bkb.raw <- rawInten[,match(bkb.v, colnames(rawInten))]'
      • at line 93 found ' dat.here <- sel.bkb.lgc[which(metadata.cell$Well.lab == well.v[i]),]'
      • at line 107 found ' dat.here <- sel.bkb.lgc[which(metadata.cell$Batch == batch.v[i]),]'
      • at line 156 found ' qual_col_pals <- brewer.pal.info[brewer.pal.info$category == 'qual',]'
    • In file R/032.rmWellEffect.R:
      • at line 139 found ' in.dat <- list(quant.bio.pcr.mt[-c(1,nrow(quant.bio.pcr.mt)),], adj.quant.bio.pcr.mt[-c(1,nrow(adj.quant.bio.pcr.mt)),])'
      • at line 142 found ' in.dat[[1]][seq_len(2),], Plate3 = -colsums(in.dat[[1]][seq_len(2),]),'
      • at line 143 found ' in.dat[[1]][3:13,], Column12 = -colsums(in.dat[[1]][3:13,]),'
      • at line 144 found ' in.dat[[1]][14:20,], Row8 = -colsums(in.dat[[1]][14:20,]),'
      • at line 145 found ' in.dat[[1]][21:nrow(in.dat[[1]]),], init.Mlast = -colsums(in.dat[[1]][21:nrow(in.dat[[1]]),])),'
      • at line 147 found ' in.dat[[2]][seq_len(2),], Plate3 = -colsums(in.dat[[2]][seq_len(2),]),'
      • at line 148 found ' in.dat[[2]][3:13,], Column12 = -colsums(in.dat[[2]][3:13,]),'
      • at line 149 found ' in.dat[[2]][14:20,], Row8 = -colsums(in.dat[[2]][14:20,]),'
      • at line 150 found ' in.dat[[2]][21:nrow(in.dat[[2]]),], init.Mlast = -colsums(in.dat[[2]][21:nrow(in.dat[[2]]),]))'
      • at line 167 found ' bio.dat <- dat0[-seq_len(23),]'
      • at line 206 found ' well.dat <- dat0[(seq_len(23)),]'
    • In file R/060.cluster.analysis.R:
      • at line 49 found ' bkb.dat <- impu.mt[, match(names(bkb.v), colnames(impu.mt))]'
      • at line 69 found ' complete.dat <- impu.mt[, -match(c(yvar, control.wells), colnames(impu.mt))]'
    • In file R/070.cluster.analysis_bkbOnly.R:
      • at line 41 found ' bkb.dat <- normalised.bkb[,match(bkb.v,colnames(normalised.bkb))]'
  • [ ] NOTE: Try to check the edge condition when using seq.int or seq_len. For example using seq.int(min(5, nrow(data))) to replace seq.int(5)
    • In file R/032.rmWellEffect.R:
      • at line 142 found ' in.dat[[1]][seq_len(2),], Plate3 = -colsums(in.dat[[1]][seq_len(2),]),'
      • at line 147 found ' in.dat[[2]][seq_len(2),], Plate3 = -colsums(in.dat[[2]][seq_len(2),]),'
      • at line 156 found ' bio.lim <- c(min(in.dat[[1]][-seq_len(23),], in.dat[[2]][-seq_len(23),]), max(in.dat[[1]][-seq_len(23),], in.dat[[2]][-seq_len(23),]))'
      • at line 157 found ' well.lim <- c(min(in.dat[[1]][seq_len(23),], in.dat[[2]][seq_len(23),]), max(in.dat[[1]][seq_len(23),], in.dat[[2]][seq_len(23),]))'
      • at line 167 found ' bio.dat <- dat0[-seq_len(23),]'
      • at line 206 found ' well.dat <- dat0[(seq_len(23)),]'
  • [ ] NOTE: Functional programming: code repetition.
    • repetition in bkc_bkb, bkc_bkb_ffc, bkc_pe, cluster.analysis, cluster.analysis.bkbOnly, imputation_bkb.predictors, initM, mkImputeMT, rmBatchEffect, and rmWellEffect
      • in bkc_bkb
        • line 1: bkb.v, bkb.upper.quantile = 0.9, bkb.lower.quantile = 0.1,
        • line 2: bkb.min.quantile = 0.01, plots)
        • line 3:{
        • line 4: upper.quantile <- bkb.upper.quantile
        • line 5: lower.quantile <- bkb.lower.quantile
        • line 6: min.quantile <- bkb.min.quantile
        • line 7: rawInten <- readRDS(file = file.path(paths["intermediary"],
        • line 8: "fcs_rawInten_mt.rds"))
        • line 9: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 10: "fcs_metadata_df.rds"))
        • line 11: sel.bkb.raw <- rawInten[, match(bkb.v, colnames(rawInten))]
        • line 12: {
        • line 13: nlogl.norm.v2 <- function(p, xx, n, rm.min.quantile = min.quantile) {
        • line 14: xx <- sort(xx)[-seq_len(round(n * rm.min.quantile))]
        • line 15: loglik <- sum(dnorm(xx, mean = p[1], sd = exp(p[2]),
        • line 16: log = TRUE)) + (n - length(xx)) * pnorm(max(xx),
        • line 17: mean = p[1], sd = exp(p[2]), log.p = TRUE, lower.tail = FALSE) +
        • line 18: round(n * rm.min.quantile) * pnorm(min(xx), mean = p[1],
        • line 19: sd = exp(p[2]), log.p = TRUE)
        • line 20: -loglik
        • line 21: }
        • line 22: nlogl.exp <- function(p, xx, n) {
        • line 23: loglik <- sum(dexp(xx, rate = exp(p[1]), log = TRUE)) +
        • line 24: (n - length(xx)) * pexp(min(xx), rate = exp(p[1]),
        • line 25: log.p = TRUE)
        • line 26: -loglik
        • line 27: }
        • line 28: calib <- function(xx, alpha, mu, sig.2) {
        • line 29: mu.sx <- xx - mu - sig.2/alpha
        • line 30: log.2ndterm <- dnorm(0, mean = mu.sx, sd = sqrt(sig.2),
        • line 31: log = TRUE) - pnorm(0, mean = mu.sx, sd = sqrt(sig.2),
        • line 32: lower.tail = FALSE, log.p = TRUE)
        • line 33: calib.xx <- mu.sx + sig.2 * exp(log.2ndterm)
        • line 34: calib.xx
        • line 35: }
        • line 36: }
        • line 37: {
        • line 38: message("\tEstimating parameters for calibration...")
        • line 39: para.ls <- list()
        • line 40: for (m in seq_len(length(bkb.v))) {
        • line 41: para.df <- data.frame(Well.lab = unique(metadata.cell$Well.lab),
        • line 42: mu = NA, sig = NA, alpha = NA)
        • line 43: for (w in seq_len(length(unique(metadata.cell$Well.lab)))) {
        • line 44: dat.here <- sel.bkb.raw[which(metadata.cell$Well.lab ==
        • line 45: para.df$Well.lab[w]), m]
        • line 46: ordering.bkb <- dat.here[order(dat.here)]
        • line 47: {
        • line 48: first.quantile <- ordering.bkb[seq_len(length(ordering.bkb) *
        • line 49: lower.quantile)]
        • line 50: ini.mu <- mean(first.quantile)
        • line 51: ini.log.sd <- log(sd(first.quantile))
        • line 52: p.init <- c(ini.mu, ini.log.sd)
        • line 53: suppressWarnings(out.1 <- tryCatch(optim(par = p.init,
        • line 54: fn = nlogl.norm.v2, xx = first.quantile,
        • line 55: n = length(ordering.bkb)), error = function(e) {
        • line 56: NA
        • line 57: }))
        • line 58: mle.mean <- out.1$par[1]
        • line 59: mle.sd <- exp(out.1$par[2])
        • line 60: para.df$mu[w] <- mle.mean
        • line 61: para.df$sig[w] <- mle.sd
        • line 62: }
        • line 63: {
        • line 64: last.quantile <- ordering.bkb[which(ordering.bkb >
        • line 65: quantile(ordering.bkb, probs = upper.quantile))]
        • line 66: p.init <- ifelse(mean(last.quantile) < 0, log(1e-06),
        • line 67: log(mean(last.quantile)))
        • line 68: suppressWarnings(out.2 <- tryCatch(optim(par = p.init,
        • line 69: fn = nlogl.exp, xx = last.quantile, n = length(ordering.bkb),
        • line 70: method = "Brent", lower = -20, upper = 0),
        • line 71: error = function(e) {
        • line 72: NA
        • line 73: }))
        • line 74: alpha <- (1/exp(out.2$p))
        • line 75: para.df$alpha[w] <- alpha
        • line 76: }
        • line 77: }
        • line 78: para.ls[[m]] <- para.df
        • line 79: message("backbone: ", m)
        • line 80: }
        • line 81: names(para.ls) <- colnames(sel.bkb.raw)
        • line 82: saveRDS(para.ls, file = file.path(paths["intermediary"],
        • line 83: paste0("para.ls_", ncol(sel.bkb.raw), "bkb_", length(unique(metadata.cell$Well.lab)),
        • line 84: "well.rds")))
        • line 85: message("\tEstimation of parameters... Completed!")
        • line 86: }
        • line 87: {
        • line 88: message("\tCalibrating backbone markers (except for physical measurements)...")
        • line 89: calib.dat <- sel.bkb.raw
        • line 90: {
        • line 91: skip.phy.meas <- which(!(colnames(sel.bkb.raw) %in%
        • line 92: c("FSC-A", "FSC-H", "FSC-W", "SSC-A", "SSC-H",
        • line 93: "SSC-W")))
        • line 94: for (m in skip.phy.meas) {
        • line 95: med.para <- apply(para.ls[[m]][, -1], 2, median)
        • line 96: med.mle.mean <- med.para[1]
        • line 97: med.mle.sd <- med.para[2]
        • line 98: med.alpha <- med.para[3]
        • line 99: calib.temp <- list()
        • line 100: for (w in seq_len(length(unique(metadata.cell$Well.lab)))) {
        • line 101: dat.here <- sel.bkb.raw[which(metadata.cell$Well.lab ==
        • line 102: para.df$Well.lab[w]), m]
        • line 103: calib.temp[[w]] <- calib(dat.here, med.alpha,
        • line 104: med.mle.mean, sig.2 = med.mle.sd^2)
        • line 105: }
        • line 106: calib.dat[, m] <- do.call(c, calib.temp)
        • line 107: }
        • line 108: saveRDS(calib.dat, file = file.path(paths["intermediary"],
        • line 109: "medpara_bkc.bkb_no.bkcPhy_mt.rds"))
        • line 110: if (plots == TRUE) {
        • line 111: marker.nam <- colnames(calib.dat)
        • line 112: for (i in seq_along(marker.nam)) {
        • line 113: jpeg(filename = file.path(paths["graph"], paste0("medpara_sub0.01.",
        • line 114: marker.nam[i], "_para.", upper.quantile,
        • line 115: "up.", lower.quantile, "lo_500x500.jpeg")),
        • line 116: width = 500, height = 500, res = 80)
        • line 117: extr.cell <- sample(seq_len(nrow(calib.dat)),
        • line 118: nrow(calib.dat) * 0.01)
        • line 119: if (i %in% skip.phy.meas) {
        • line 120: plot(x = sel.bkb.raw[extr.cell, i], y = calib.dat[extr.cell,
        • line 121: i], xlab = "Raw", ylab = "Calib.", main = marker.nam[i],
        • line 122: col = "blue", sub = paste0("Parameters Estimated from ",
        • line 123: upper.quantile, " (hi); ", lower.quantile,
        • line 124: " (lo)"))
        • line 125: }
        • line 126: else {
        • line 127: plot(x = sel.bkb.raw[extr.cell, i], y = calib.dat[extr.cell,
        • line 128: i], xlab = "Raw", ylab = "Raw.", main = marker.nam[i],
        • line 129: col = "blue", sub = paste0("No need to calibrate the physical measurements."))
        • line 130: }
        • line 131: abline(a = 0, b = 1)
        • line 132: dev.off()
        • line 133: }
        • line 134: }
        • line 135: }
        • line 136: message("\tCalibration of backbone markers... Completed!")
        • line 137: }
      • in bkc_bkb_ffc
        • line 1: bkb.v, bkb.upper.quantile = 0.9, bkb.lower.quantile = 0.1,
        • line 2: bkb.min.quantile = 0.01, plots = TRUE)
        • line 3:{
        • line 4: upper.quantile <- bkb.upper.quantile
        • line 5: lower.quantile <- bkb.lower.quantile
        • line 6: min.quantile <- bkb.min.quantile
        • line 7: rawInten <- readRDS(file = file.path(paths["intermediary"],
        • line 8: "fcs_rawInten_mt.rds"))
        • line 9: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 10: "fcs_metadata_df.rds"))
        • line 11: sel.bkb.raw <- rawInten[, match(bkb.v, colnames(rawInten))]
        • line 12: {
        • line 13: nlogl.norm.v2 <- function(p, xx, n, rm.min.quantile = min.quantile) {
        • line 14: xx <- sort(xx)[-seq_len(round(n * rm.min.quantile))]
        • line 15: loglik <- sum(dnorm(xx, mean = p[1], sd = exp(p[2]),
        • line 16: log = TRUE)) + (n - length(xx)) * pnorm(max(xx),
        • line 17: mean = p[1], sd = exp(p[2]), log.p = TRUE, lower.tail = FALSE) +
        • line 18: round(n * rm.min.quantile) * pnorm(min(xx), mean = p[1],
        • line 19: sd = exp(p[2]), log.p = TRUE)
        • line 20: -loglik
        • line 21: }
        • line 22: nlogl.exp <- function(p, xx, n) {
        • line 23: loglik <- sum(dexp(xx, rate = exp(p[1]), log = TRUE)) +
        • line 24: (n - length(xx)) * pexp(min(xx), rate = exp(p[1]),
        • line 25: log.p = TRUE)
        • line 26: -loglik
        • line 27: }
        • line 28: calib <- function(xx, alpha, mu, sig.2) {
        • line 29: mu.sx <- xx - mu - sig.2/alpha
        • line 30: log.2ndterm <- dnorm(0, mean = mu.sx, sd = sqrt(sig.2),
        • line 31: log = TRUE) - pnorm(0, mean = mu.sx, sd = sqrt(sig.2),
        • line 32: lower.tail = FALSE, log.p = TRUE)
        • line 33: calib.xx <- mu.sx + sig.2 * exp(log.2ndterm)
        • line 34: calib.xx
        • line 35: }
        • line 36: }
        • line 37: {
        • line 38: message("\tEstimating parameters for calibration...")
        • line 39: para.ls <- list()
        • line 40: for (m in seq_len(length(bkb.v))) {
        • line 41: para.df <- data.frame(Batch = unique(metadata.cell$Batch),
        • line 42: mu = NA, sig = NA, alpha = NA)
        • line 43: for (w in seq_len(length(unique(metadata.cell$Batch)))) {
        • line 44: dat.here <- sel.bkb.raw[which(metadata.cell$Batch ==
        • line 45: para.df$Batch[w]), m]
        • line 46: ordering.bkb <- dat.here[order(dat.here)]
        • line 47: {
        • line 48: first.quantile <- ordering.bkb[seq_len(length(ordering.bkb) *
        • line 49: lower.quantile)]
        • line 50: ini.mu <- mean(first.quantile)
        • line 51: ini.log.sd <- log(sd(first.quantile))
        • line 52: p.init <- c(ini.mu, ini.log.sd)
        • line 53: suppressWarnings(out.1 <- tryCatch(optim(par = p.init,
        • line 54: fn = nlogl.norm.v2, xx = first.quantile,
        • line 55: n = length(ordering.bkb)), error = function(e) {
        • line 56: NA
        • line 57: }))
        • line 58: mle.mean <- out.1$par[1]
        • line 59: mle.sd <- exp(out.1$par[2])
        • line 60: para.df$mu[w] <- mle.mean
        • line 61: para.df$sig[w] <- mle.sd
        • line 62: }
        • line 63: {
        • line 64: last.quantile <- ordering.bkb[which(ordering.bkb >
        • line 65: quantile(ordering.bkb, probs = upper.quantile))]
        • line 66: p.init <- ifelse(mean(last.quantile) < 0, log(1e-06),
        • line 67: log(mean(last.quantile)))
        • line 68: out.2 <- tryCatch(optim(par = p.init, fn = nlogl.exp,
        • line 69: xx = last.quantile, n = length(ordering.bkb),
        • line 70: method = "Brent", lower = -20, upper = 0),
        • line 71: error = function(e) {
        • line 72: NA
        • line 73: })
        • line 74: alpha <- (1/exp(out.2$p))
        • line 75: para.df$alpha[w] <- alpha
        • line 76: }
        • line 77: }
        • line 78: para.ls[[m]] <- para.df
        • line 79: message("backbone: ", m)
        • line 80: }
        • line 81: names(para.ls) <- colnames(sel.bkb.raw)
        • line 82: saveRDS(para.ls, file = file.path(paths["intermediary"],
        • line 83: paste0("para.ls_", ncol(sel.bkb.raw), "bkb_", length(unique(metadata.cell$Batch)),
        • line 84: "batch.rds")))
        • line 85: message("\tEstimation of parameters... Completed!")
        • line 86: }
        • line 87: {
        • line 88: message("\tCalibrating backbone markers (except for physical measurements)...")
        • line 89: calib.dat <- sel.bkb.raw
        • line 90: {
        • line 91: skip.phy.meas <- which(!(colnames(sel.bkb.raw) %in%
        • line 92: c("FSC-A", "FSC-H", "FSC-W", "SSC-A", "SSC-H",
        • line 93: "SSC-W")))
        • line 94: for (m in skip.phy.meas) {
        • line 95: med.para <- apply(para.ls[[m]][, -1], 2, median)
        • line 96: med.mle.mean <- med.para[1]
        • line 97: med.mle.sd <- med.para[2]
        • line 98: med.alpha <- med.para[3]
        • line 99: calib.temp <- list()
        • line 100: for (w in seq_len(length(unique(metadata.cell$Batch)))) {
        • line 101: dat.here <- sel.bkb.raw[which(metadata.cell$Batch ==
        • line 102: para.df$Batch[w]), m]
        • line 103: calib.temp[[w]] <- calib(dat.here, med.alpha,
        • line 104: med.mle.mean, sig.2 = med.mle.sd^2)
        • line 105: }
        • line 106: calib.dat[, m] <- do.call(c, calib.temp)
        • line 107: }
        • line 108: saveRDS(calib.dat, file = file.path(paths["intermediary"],
        • line 109: "medpara_bkc.bkb_no.bkcPhy_mt.rds"))
        • line 110: if (plots == TRUE) {
        • line 111: marker.nam <- colnames(calib.dat)
        • line 112: for (i in seq_along(marker.nam)) {
        • line 113: jpeg(filename = file.path(paths["graph"], paste0("medpara_sub0.01.",
        • line 114: marker.nam[i], "_para.", upper.quantile,
        • line 115: "up.", lower.quantile, "lo_500x500.jpeg")),
        • line 116: width = 500, height = 500, res = 80)
        • line 117: extr.cell <- sample(seq_len(nrow(calib.dat)),
        • line 118: nrow(calib.dat) * 0.01)
        • line 119: if (i %in% skip.phy.meas) {
        • line 120: plot(x = sel.bkb.raw[extr.cell, i], y = calib.dat[extr.cell,
        • line 121: i], xlab = "Raw", ylab = "Calib.", main = marker.nam[i],
        • line 122: col = "blue", sub = paste0("Parameters Estimated from ",
        • line 123: upper.quantile, " (hi); ", lower.quantile,
        • line 124: " (lo)"))
        • line 125: }
        • line 126: else {
        • line 127: plot(x = sel.bkb.raw[extr.cell, i], y = calib.dat[extr.cell,
        • line 128: i], xlab = "Raw", ylab = "Raw.", main = marker.nam[i],
        • line 129: col = "blue", sub = paste0("No need to calibrate the physical measurements."))
        • line 130: }
        • line 131: abline(a = 0, b = 1)
        • line 132: dev.off()
        • line 133: }
        • line 134: }
        • line 135: }
        • line 136: message("\tCalibration of backbone markers... Completed!")
        • line 137: }
      • in bkc_pe
        • line 3: min.quantile <- pe.min.quantile
        • line 4: lower.quantile <- pe.lower.quantile
        • line 5: rawInten <- readRDS(file = file.path(paths["intermediary"],
        • line 6: "fcs_rawInten_mt.rds"))
        • line 7: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 8: "fcs_metadata_df.rds"))
        • line 11: {
        • line 12: nlogl.norm.v2 <- function(p, xx, n, rm.min.quantile = min.quantile) {
        • line 13: xx <- sort(xx)[-seq_len(round(n * rm.min.quantile))]
        • line 14: loglik <- sum(dnorm(xx, mean = p[1], sd = exp(p[2]),
        • line 15: log = TRUE)) + (n - length(xx)) * pnorm(max(xx),
        • line 16: mean = p[1], sd = exp(p[2]), log.p = TRUE, lower.tail = FALSE) +
        • line 17: round(n * rm.min.quantile) * pnorm(min(xx), mean = p[1],
        • line 18: sd = exp(p[2]), log.p = TRUE)
        • line 19: -loglik
        • line 20: }
        • line 21: nlogl.exp <- function(p, xx, n) {
        • line 22: loglik <- sum(dexp(xx, rate = exp(p[1]), log = TRUE)) +
        • line 23: (n - length(xx)) * pexp(min(xx), rate = exp(p[1]),
        • line 24: log.p = TRUE)
        • line 25: -loglik
        • line 26: }
        • line 27: calib <- function(xx, alpha, mu, sig.2) {
        • line 28: mu.sx <- xx - mu - sig.2/alpha
        • line 29: log.2ndterm <- dnorm(0, mean = mu.sx, sd = sqrt(sig.2),
        • line 30: log = TRUE) - pnorm(0, mean = mu.sx, sd = sqrt(sig.2),
        • line 31: lower.tail = FALSE, log.p = TRUE)
        • line 32: calib.xx <- mu.sx + sig.2 * exp(log.2ndterm)
        • line 33: calib.xx
        • line 34: }
        • line 35: }
        • line 48: first.quantile <- ordering.legend[seq_len(length(ordering.legend) *
        • line 49: lower.quantile)]
        • line 50: ini.mu <- mean(first.quantile)
        • line 51: ini.log.sd <- log(sd(first.quantile))
        • line 52: p.init <- c(ini.mu, ini.log.sd)
        • line 53: suppressWarnings(out.1 <- tryCatch(optim(par = p.init,
        • line 54: fn = nlogl.norm.v2, xx = first.quantile, n = length(ordering.legend)),
        • line 55: error = function(e) {
        • line 56: NA
        • line 57: }))
        • line 58: mle.mean <- out.1$par[1]
        • line 59: mle.sd <- exp(out.1$par[2])
        • line 60: para.df$mu[i] <- mle.mean
        • line 61: para.df$sig[i] <- mle.sd
        • line 75: suppressWarnings(out.2 <- tryCatch(optim(par = p.init,
        • line 76: fn = nlogl.exp, xx = last.mu.3sd, n = length(ordering.legend),
        • line 77: method = "Brent", lower = -20, upper = 0), error = function(e) {
        • line 78: NA
        • line 79: }))
        • line 80: alpha <- (1/exp(out.2$p))
        • line 81: para.df$alpha[i] <- alpha
        • line 107: abline(a = 0, b = 1)
        • line 108: dev.off()
        • line 109: }
        • line 110: }
      • in cluster.analysis
        • line 13: message("Running UMAP...")
        • line 14: a <- Sys.time()
        • line 15: umap.bkb <- umap(bkb.dat, n_neighbors = 15, min_dist = 0.2,
        • line 16: metric = "euclidean", n_epochs = 2000)
        • line 17: b <- Sys.time()
        • line 18: b - a
        • line 19: saveRDS(umap.bkb, file = file.path(paths["downstream"],
        • line 20: paste0("ClusterAnalysis_", "umap_", length(bkb.v),
        • line 21: "bkb.rds")))
        • line 22: message("Running Phenograph...")
        • line 23: a <- Sys.time()
        • line 24: phenog.bkb <- Rphenograph(bkb.dat, k = 50)
        • line 25: b <- Sys.time()
        • line 26: b - a
        • line 27: saveRDS(phenog.bkb, file = file.path(paths["intermediary"],
        • line 28: paste0("ClusterAnalysis_", "phenog_", length(bkb.v),
        • line 29: "bkb.rds")))
        • line 30: }
        • line 35: message("Running UMAP...")
        • line 36: a <- Sys.time()
        • line 37: umap.bkb.impuInf <- umap(complete.dat, n_neighbors = 15,
        • line 38: min_dist = 0.2, metric = "euclidean", n_epochs = 2000)
        • line 39: b <- Sys.time()
        • line 40: b - a
        • line 45: message("Running Phenograph...")
        • line 46: a <- Sys.time()
        • line 47: phenog.bkb.impuInf <- Rphenograph(complete.dat, k = 50)
        • line 48: b <- Sys.time()
        • line 60: if (plots == TRUE) {
        • line 61: message("\tVisualising clusters...")
        • line 62: {
        • line 63: {
        • line 67: colnames(graph.dat) <- c("UMAP1", "UMAP2",
        • line 68: "Cluster")
        • line 69: n <- length(unique(graph.dat[, 3]))
        • line 70: qual_col_pals <- brewer.pal.info[brewer.pal.info$category ==
        • line 71: "qual", ]
        • line 72: col.coeff <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors,
        • line 73: rownames(qual_col_pals)))[seq_len(n)]
        • line 74: jpeg(filename = file.path(paths["graph"],
        • line 75: paste0("/ClusterStructure_UMAP_", length(bkb.v),
        • line 76: "bkb_colPhenog_600x750.jpeg")), height = 600,
        • line 77: width = 750, res = 80)
        • line 78: p <- ggplot(data = graph.dat, mapping = aes(x = UMAP1,
        • line 79: y = UMAP2, colour = Cluster)) + geom_point(size = 0.05,
        • line 80: alpha = 0.25) + labs(titles = paste0("Normalised ",
        • line 81: length(bkb.v), "bkb. \n(", nrow(graph.dat),
        • line 82: " cells)")) + guides(colour = guide_legend(override.aes = list(size = 8))) +
        • line 83: theme_bw() + theme(text = element_text(size = 22)) +
        • line 84: scale_color_manual(values = col.coeff[seq_len(n)])
        • line 85: plot(p)
        • line 86: dev.off()
        • line 87: }
        • line 88: }
        • line 89: {
        • line 92: colnames(graph.dat) <- c("UMAP1", "UMAP2",
        • line 93: "Cluster")
        • line 94: n <- length(unique(graph.dat[, 3]))
        • line 95: qual_col_pals <- brewer.pal.info[brewer.pal.info$category ==
        • line 96: "qual", ]
        • line 97: col.coeff <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors,
        • line 98: rownames(qual_col_pals)))[seq_len(n)]
        • line 99: jpeg(filename = file.path(paths["graph"], paste0("/ClusterStructure_UMAP_",
        • line 100: names(preds)[i], length(bkb.v), "bkb", (ncol(complete.dat) -
        • line 101: length(bkb.v)), "impuInf_colPhenog_600x750.jpeg")),
        • line 102: height = 600, width = 750, res = 80)
        • line 103: p <- ggplot(data = graph.dat, mapping = aes(x = UMAP1,
        • line 104: y = UMAP2, colour = Cluster)) + geom_point(size = 0.05,
        • line 105: alpha = 0.25) + labs(titles = paste0(names(preds)[i],
        • line 106: " ", length(bkb.v), "bkb. ", (ncol(complete.dat) -
        • line 107: length(bkb.v)), "impuInf \n(", nrow(graph.dat),
        • line 108: " cells)")) + guides(colour = guide_legend(override.aes = list(size = 8))) +
        • line 109: theme_bw() + theme(text = element_text(size = 22)) +
        • line 110: scale_color_manual(values = col.coeff[seq_len(n)])
        • line 111: plot(p)
        • line 112: dev.off()
        • line 113: }
        • line 114: }
        • line 115: }
        • line 116: }
        • line 117: head(metadata.cell)
        • line 118: saveRDS(metadata.cell, file = file.path(paths["downstream"],
        • line 119: "fcs_metadata_df.rds"))
        • line 120: message("\tCompleted!")
      • in cluster.analysis.bkbOnly
        • line 12: message("Running UMAP...")
        • line 13: a <- Sys.time()
        • line 14: umap.bkb <- umap(bkb.dat, n_neighbors = 15, min_dist = 0.2,
        • line 15: metric = "euclidean", n_epochs = 2000)
        • line 16: b <- Sys.time()
        • line 17: b - a
        • line 18: saveRDS(umap.bkb, file = file.path(paths["downstream"],
        • line 19: paste0("ClusterAnalysis_", "umap_", length(bkb.v),
        • line 20: "bkb.rds")))
        • line 21: message("Running Phenograph...")
        • line 22: a <- Sys.time()
        • line 23: phenog.bkb <- Rphenograph(bkb.dat, k = 50)
        • line 24: b <- Sys.time()
        • line 25: b - a
        • line 26: saveRDS(phenog.bkb, file = file.path(paths["intermediary"],
        • line 27: paste0("ClusterAnalysis_", "phenog_", length(bkb.v),
        • line 28: "bkb.rds")))
        • line 30: if (plots == TRUE) {
        • line 31: message("\tVisualising clusters...")
        • line 32: {
        • line 33: {
        • line 35: colnames(graph.dat) <- c("UMAP1", "UMAP2",
        • line 36: "Cluster")
        • line 37: n <- length(unique(graph.dat[, 3]))
        • line 38: qual_col_pals <- brewer.pal.info[brewer.pal.info$category ==
        • line 39: "qual", ]
        • line 40: col.coeff <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors,
        • line 41: rownames(qual_col_pals)))[seq_len(n)]
        • line 42: jpeg(filename = file.path(paths["graph"], paste0("/ClusterStructure_UMAP_",
        • line 43: length(bkb.v), "bkb_colPhenog_600x750.jpeg")),
        • line 44: height = 600, width = 750, res = 80)
        • line 45: p <- ggplot(data = graph.dat, mapping = aes(x = UMAP1,
        • line 46: y = UMAP2, colour = Cluster)) + geom_point(size = 0.05,
        • line 47: alpha = 0.25) + labs(titles = paste0("Normalised ",
        • line 48: length(bkb.v), "bkb. \n(", nrow(graph.dat),
        • line 49: " cells)")) + guides(colour = guide_legend(override.aes = list(size = 8))) +
        • line 50: theme_bw() + theme(text = element_text(size = 22)) +
        • line 51: scale_color_manual(values = col.coeff[seq_len(n)])
        • line 52: plot(p)
        • line 53: dev.off()
        • line 54: }
        • line 55: }
        • line 56: }
        • line 57: head(metadata.cell)
        • line 58: saveRDS(metadata.cell, file = file.path(paths["downstream"],
        • line 59: "fcs_metadata_df.rds"))
        • line 60: message("\tCompleted!")
      • in imputation_bkb.predictors
        • line 6: xp <- readRDS(file = file.path(paths["intermediary"],
        • line 7: "impu.input_log.mt.rds"))
        • line 8: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 9: "fcs_metadata_df.rds"))
        • line 503: plot(p)
        • line 504: dev.off()
        • line 505: }
        • line 506: }
        • line 507: }
      • in initM
        • line 4: trans.dat <- "cent.lgc"
        • line 5: rawInten <- readRDS(file = file.path(paths["intermediary"],
        • line 6: "fcs_rawInten_mt.rds"))
        • line 7: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 8: "fcs_metadata_df.rds"))
        • line 9: sel.bkb.raw <- rawInten[, match(bkb.v, colnames(rawInten))]
        • line 10: if (trans.dat == "cent.lgc") {
        • line 100: n <- length(unique(graph.dat$Pheno.gp))
        • line 101: qual_col_pals <- brewer.pal.info[brewer.pal.info$category ==
        • line 102: "qual", ]
        • line 103: col.coeff <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors,
        • line 104: rownames(qual_col_pals)))[seq_len(n)]
      • in mkImputeMT
        • line 4: bkc.pe <- readRDS(file = file.path(paths["intermediary"],
        • line 5: "bkc.pe_mt.rds"))
        • line 6: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 7: "fcs_metadata_df.rds"))
      • in rmBatchEffect
        • line 1: plots = TRUE)
        • line 2:{
        • line 3: bkc.bkb <- readRDS(file = file.path(paths["intermediary"],
        • line 4: "medpara_bkc.bkb_no.bkcPhy_mt.rds"))
        • line 5: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 6: "fcs_metadata_df.rds"))
        • line 7: {
        • line 8: message("\tEstimating coefficients for removing batch effect (Rfast - pre.adj)...")
        • line 9: a <- Sys.time()
        • line 10: {
        • line 11: metadata.cell[, c("Batch", "init.M")] <- lapply(metadata.cell[,
        • line 14: contrasts = TRUE)
        • line 15: contrasts(metadata.cell[, "init.M"]) <- contr.sum(length(unique(metadata.cell$init.M)),
        • line 16: contrasts = TRUE)
        • line 17: xx <- as.matrix(model.matrix(~Batch + init.M, data = metadata.cell))
        • line 18: lm.results.ls <- list()
        • line 19: ln.sig.v <- c()
        • line 20: res.df <- bkc.bkb
        • line 21: for (i in seq_along(bkc.bkb[1, ])) {
        • line 22: yy.val <- as.matrix(bkc.bkb[, i])
        • line 23: lm.result <- tryCatch(lmfit(x = xx, y = log(yy.val +
        • line 24: 1e-06)), error = function(err) NA)
        • line 25: lm.results.ls[[i]] <- lm.result
        • line 26: res <- lm.result$residuals
        • line 27: ln.sig.v[i] <- sqrt(sum((res - mean(res))^2)/(length(res) -
        • line 28: length(lm.results.ls[[1]]$be)))
        • line 29: res.df[, i] <- res
        • line 30: message("Processing backbone: ", i)
        • line 31: }
        • line 32: names(lm.results.ls) <- colnames(bkc.bkb)
        • line 33: quant.bio.pcr.mt <- do.call(cbind, lapply(lm.results.ls,
        • line 34: "[[", 1))
        • line 35: colnames(quant.bio.pcr.mt) <- names(lm.results.ls)
        • line 36: quant.bio.pcr.mt <- rbind(quant.bio.pcr.mt, ln.sig = ln.sig.v)
        • line 37: save(lm.results.ls, quant.bio.pcr.mt, file = file.path(paths["intermediary"],
        • line 38: "rmBatchEffect_log.normal.reg_coef.mdl_pre.RData"))
        • line 39: save(res.df, file = file.path(paths["intermediary"],
        • line 40: "rmBatchEffect_log.normal.reg_residuals.RData"))
        • line 41: }
        • line 42: b <- Sys.time()
        • line 43: b - a
        • line 44: message("\tEstimation completed!")
        • line 45: }
        • line 46: {
        • line 47: message("\tRemoving batch effect for backbone markers...")
        • line 48: adj.data <- log.adj.data <- bkc.bkb
        • line 49: a <- Sys.time()
        • line 50: for (i in seq_along(bkc.bkb[1, ])) {
        • line 53: i])
        • line 54: log.adj.data[, i] <- (log(bkc.bkb[, i] + 1e-06) -
        • line 55: unwanted.eff)
        • line 56: adj.data[, i] <- (exp(log.adj.data[, i]) - 1e-06)
        • line 57: }
        • line 58: b <- Sys.time()
        • line 59: b - a
        • line 60: saveRDS(log.adj.data, file = file.path(paths["downstream"],
        • line 61: "bkc.adj.bkb_logScale_mt.rds"))
        • line 62: saveRDS(adj.data, file = file.path(paths["downstream"],
        • line 63: "bkc.adj.bkb_linearScale_mt.rds"))
        • line 64: message("\tAdjustment completed!")
        • line 65: }
        • line 66: {
        • line 67: message("\tExamining the existence of batch effect in the adjusted data (Rfast - post.adj)...")
        • line 68: a <- Sys.time()
        • line 69: {
        • line 70: metadata.cell[, c("Batch", "init.M")] <- lapply(metadata.cell[,
        • line 73: contrasts = TRUE)
        • line 74: contrasts(metadata.cell[, "init.M"]) <- contr.sum(length(unique(metadata.cell$init.M)),
        • line 75: contrasts = TRUE)
        • line 76: xx <- as.matrix(model.matrix(~Batch + init.M, data = metadata.cell))
        • line 77: adj.lm.results.ls <- list()
        • line 78: ln.sig.v <- c()
        • line 79: for (i in seq_along(log.adj.data[1, ])) {
        • line 80: yy.val <- as.matrix(log.adj.data[, i])
        • line 81: lm.result <- tryCatch(lmfit(x = xx, y = yy.val),
        • line 82: error = function(err) NA)
        • line 83: adj.lm.results.ls[[i]] <- lm.result
        • line 84: res <- lm.result$residuals
        • line 85: ln.sig.v[i] <- sqrt(sum((res - mean(res))^2)/(length(res) -
        • line 86: length(lm.results.ls[[1]]$be)))
        • line 87: message("Processing backbone: ", i)
        • line 88: }
        • line 89: names(adj.lm.results.ls) <- colnames(log.adj.data)
        • line 90: adj.quant.bio.pcr.mt <- do.call(cbind, lapply(adj.lm.results.ls,
        • line 91: "[[", 1))
        • line 92: colnames(adj.quant.bio.pcr.mt) <- names(adj.lm.results.ls)
        • line 93: adj.quant.bio.pcr.mt <- rbind(adj.quant.bio.pcr.mt,
        • line 94: ln.sig = ln.sig.v)
        • line 95: save(lm.results.ls, quant.bio.pcr.mt, adj.quant.bio.pcr.mt,
        • line 96: file = file.path(paths["intermediary"], "rmBatchEffect_log.normal.reg_coef.mdl_post.RData"))
        • line 97: }
        • line 98: b <- Sys.time()
        • line 99: b - a
        • line 100: }
        • line 101: if (plots == TRUE) {
        • line 102: {
        • line 103: in.dat <- list(quant.bio.pcr.mt[-c(1, nrow(quant.bio.pcr.mt)),
        • line 104: ], adj.quant.bio.pcr.mt[-c(1, nrow(adj.quant.bio.pcr.mt)),
        • line 105: ])
        • line 106: in.dat <- list(rbind(in.dat[[1]][seq_len(length(unique(metadata.cell$Batch)) -
        • line 134: ]))
        • line 135: pre.post <- c("pre-", "post-")
        • line 136: bio.sys <- c("bio.", "batch.")
        • line 137: for (i in seq_along(pre.post)) {
        • line 138: for (j in seq_along(bio.sys)) {
        • line 139: dat0 <- in.dat[[i]]
        • line 140: if (j == 1) {
        • line 142: ]
        • line 143: meta_mat <- data.frame(coeff = rownames(bio.dat))
        • line 144: meta_mat$Coef <- factor(paste0("init.M",
        • line 145: seq_len(nrow(meta_mat))), levels = paste0("init.M",
        • line 146: seq_len(nrow(meta_mat))))
        • line 147: rownames(meta_mat) <- meta_mat$Coef
        • line 148: n <- nrow(meta_mat)
        • line 149: qual_col_pals <- brewer.pal.info[brewer.pal.info$category ==
        • line 150: "qual", ]
        • line 151: col.coeff <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors,
        • line 152: rownames(qual_col_pals)))[seq_len(n)]
        • line 153: names(col.coeff) <- meta_mat$Coef
        • line 154: row_ha <- rowAnnotation(Effect = meta_mat[,
        • line 155: "Coef"], col = list(Effect = col.coeff),
        • line 156: show_legend = FALSE)
        • line 157: coeff_col_func <- colorRamp2(c(bio.lim[1],
        • line 158: 0, bio.lim[2]), c("blue", "gray0", "darkorange1"))
        • line 159: jpeg(file.path(paths["graph"], paste0("Heatmap_coeffs_",
        • line 160: pre.post[i], bio.sys[j], "_2000x1500.jpeg")),
        • line 161: height = 2000, width = 1500, res = 250)
        • line 162: draw(Heatmap(bio.dat, name = paste0("MLE"),
        • line 163: column_title = paste0(pre.post[i], "batch.adj."),
        • line 164: row_title = bio.sys[j], right_annotation = row_ha,
        • line 165: col = coeff_col_func, column_dend_height = unit(1,
        • line 166: "cm"), cluster_rows = FALSE, cluster_columns = FALSE,
        • line 167: clustering_method_rows = "average", clustering_method_columns = "average",
        • line 168: show_column_names = TRUE, show_row_names = TRUE))
        • line 169: dev.off()
        • line 170: }
        • line 171: else {
        • line 176: col.coeff <- c(brewer.pal(9, "Set1")[c(1,
        • line 177: 3, 2)], brewer.pal(8, "Dark2")[c(4, 1,
        • line 178: 3)])
        • line 179: row_ha <- rowAnnotation(Effect = meta_mat[,
        • line 180: "Coef"], col = list(Effect = c(Batch = col.coeff[4])))
        • line 181: coeff_col_func <- colorRamp2(c(batch.lim[1],
        • line 182: 0, batch.lim[2]), c("blue", "gray0", "darkorange1"))
        • line 183: jpeg(file.path(paths["graph"], paste0("Heatmap_coeffs_",
        • line 184: pre.post[i], bio.sys[j], "_2000x1500.jpeg")),
        • line 185: height = 2000, width = 1500, res = 250)
        • line 186: draw(Heatmap(batch.dat, name = paste0("MLE"),
        • line 187: column_title = paste0(pre.post[i], "batch.adj."),
        • line 188: row_title = bio.sys[j], right_annotation = row_ha,
        • line 189: col = coeff_col_func, column_dend_height = unit(1,
        • line 190: "cm"), cluster_rows = FALSE, cluster_columns = FALSE,
        • line 191: clustering_method_rows = "average", clustering_method_columns = "average",
        • line 192: show_column_names = TRUE, show_row_names = TRUE))
        • line 193: dev.off()
        • line 194: }
        • line 195: }
        • line 196: }
      • in rmWellEffect
        • line 1: plots = TRUE)
        • line 2:{
        • line 3: bkc.bkb <- readRDS(file = file.path(paths["intermediary"],
        • line 4: "medpara_bkc.bkb_no.bkcPhy_mt.rds"))
        • line 5: metadata.cell <- readRDS(file = file.path(paths["intermediary"],
        • line 6: "fcs_metadata_df.rds"))
        • line 7: {
        • line 8: message("\tEstimating coefficients for removing well effect (Rfast - pre.adj)...")
        • line 9: a <- Sys.time()
        • line 10: {
        • line 11: metadata.cell[, c("Plate", "Column", "Row", "init.M")] <- lapply(metadata.cell[,
        • line 18: contrasts = TRUE)
        • line 19: contrasts(metadata.cell[, "init.M"]) <- contr.sum(length(unique(metadata.cell$init.M)),
        • line 20: contrasts = TRUE)
        • line 21: xx <- as.matrix(model.matrix(~Plate + Column + Row +
        • line 22: init.M, data = metadata.cell))
        • line 23: lm.results.ls <- list()
        • line 24: ln.sig.v <- c()
        • line 25: res.df <- bkc.bkb
        • line 26: for (i in seq_along(bkc.bkb[1, ])) {
        • line 27: yy.val <- as.matrix(bkc.bkb[, i])
        • line 28: lm.result <- tryCatch(lmfit(x = xx, y = log(yy.val +
        • line 29: 1e-06)), error = function(err) NA)
        • line 30: lm.results.ls[[i]] <- lm.result
        • line 31: res <- lm.result$residuals
        • line 32: ln.sig.v[i] <- sqrt(sum((res - mean(res))^2)/(length(res) -
        • line 33: length(lm.results.ls[[1]]$be)))
        • line 34: res.df[, i] <- res
        • line 35: message("Processing backbone: ", i)
        • line 36: }
        • line 37: names(lm.results.ls) <- colnames(bkc.bkb)
        • line 38: quant.bio.pcr.mt <- do.call(cbind, lapply(lm.results.ls,
        • line 39: "[[", 1))
        • line 40: colnames(quant.bio.pcr.mt) <- names(lm.results.ls)
        • line 41: quant.bio.pcr.mt <- rbind(quant.bio.pcr.mt, ln.sig = ln.sig.v)
        • line 42: save(lm.results.ls, quant.bio.pcr.mt, file = file.path(paths["intermediary"],
        • line 43: "rmWellEffect_log.normal.reg_coef.mdl_pre.RData"))
        • line 44: save(res.df, file = file.path(paths["intermediary"],
        • line 45: "rmWellEffect_log.normal.reg_residuals.RData"))
        • line 46: }
        • line 47: b <- Sys.time()
        • line 48: b - a
        • line 49: message("\tEstimation completed!")
        • line 50: }
        • line 51: {
        • line 52: message("\tRemoving well effect for backbone markers...")
        • line 53: adj.data <- log.adj.data <- bkc.bkb
        • line 54: a <- Sys.time()
        • line 55: for (i in seq_along(bkc.bkb[1, ])) {
        • line 57: i])
        • line 58: log.adj.data[, i] <- (log(bkc.bkb[, i] + 1e-06) -
        • line 59: unwanted.eff)
        • line 60: adj.data[, i] <- (exp(log.adj.data[, i]) - 1e-06)
        • line 61: }
        • line 62: b <- Sys.time()
        • line 63: b - a
        • line 64: saveRDS(log.adj.data, file = file.path(paths["downstream"],
        • line 65: "bkc.adj.bkb_logScale_mt.rds"))
        • line 66: saveRDS(adj.data, file = file.path(paths["downstream"],
        • line 67: "bkc.adj.bkb_linearScale_mt.rds"))
        • line 68: message("\tAdjustment completed!")
        • line 69: }
        • line 70: {
        • line 71: message("\tExamining the existence of well effect in the adjusted data (Rfast - post.adj)...")
        • line 72: a <- Sys.time()
        • line 73: {
        • line 74: metadata.cell[, c("Plate", "Column", "Row", "init.M")] <- lapply(metadata.cell[,
        • line 81: contrasts = TRUE)
        • line 82: contrasts(metadata.cell[, "init.M"]) <- contr.sum(length(unique(metadata.cell$init.M)),
        • line 83: contrasts = TRUE)
        • line 84: xx <- as.matrix(model.matrix(~Plate + Column + Row +
        • line 85: init.M, data = metadata.cell))
        • line 86: adj.lm.results.ls <- list()
        • line 87: ln.sig.v <- c()
        • line 88: for (i in seq_along(log.adj.data[1, ])) {
        • line 89: yy.val <- as.matrix(log.adj.data[, i])
        • line 90: lm.result <- tryCatch(lmfit(x = xx, y = yy.val),
        • line 91: error = function(err) NA)
        • line 92: adj.lm.results.ls[[i]] <- lm.result
        • line 93: res <- lm.result$residuals
        • line 94: ln.sig.v[i] <- sqrt(sum((res - mean(res))^2)/(length(res) -
        • line 95: length(lm.results.ls[[1]]$be)))
        • line 96: message("Processing backbone: ", i)
        • line 97: }
        • line 98: names(adj.lm.results.ls) <- colnames(log.adj.data)
        • line 99: adj.quant.bio.pcr.mt <- do.call(cbind, lapply(adj.lm.results.ls,
        • line 100: "[[", 1))
        • line 101: colnames(adj.quant.bio.pcr.mt) <- names(adj.lm.results.ls)
        • line 102: adj.quant.bio.pcr.mt <- rbind(adj.quant.bio.pcr.mt,
        • line 103: ln.sig = ln.sig.v)
        • line 104: save(lm.results.ls, quant.bio.pcr.mt, adj.quant.bio.pcr.mt,
        • line 105: file = file.path(paths["intermediary"], "rmWellEffect_log.normal.reg_coef.mdl_post.RData"))
        • line 106: }
        • line 107: b <- Sys.time()
        • line 108: b - a
        • line 109: }
        • line 110: if (plots == TRUE) {
        • line 111: {
        • line 112: in.dat <- list(quant.bio.pcr.mt[-c(1, nrow(quant.bio.pcr.mt)),
        • line 113: ], adj.quant.bio.pcr.mt[-c(1, nrow(adj.quant.bio.pcr.mt)),
        • line 114: ])
        • line 115: in.dat <- list(rbind(in.dat[[1]][seq_len(2), ], Plate3 = -colsums(in.dat[[1]][seq_len(2),
        • line 133: ]))
        • line 134: pre.post <- c("pre-", "post-")
        • line 135: bio.sys <- c("bio.", "well.")
        • line 136: for (i in seq_along(pre.post)) {
        • line 137: for (j in seq_along(bio.sys)) {
        • line 138: dat0 <- in.dat[[i]]
        • line 139: if (j == 1) {
        • line 141: meta_mat <- data.frame(coeff = rownames(bio.dat))
        • line 142: meta_mat$Coef <- factor(paste0("init.M",
        • line 143: seq_len(nrow(meta_mat))), levels = paste0("init.M",
        • line 144: seq_len(nrow(meta_mat))))
        • line 145: rownames(meta_mat) <- meta_mat$Coef
        • line 146: n <- nrow(meta_mat)
        • line 147: qual_col_pals <- brewer.pal.info[brewer.pal.info$category ==
        • line 148: "qual", ]
        • line 149: col.coeff <- unlist(mapply(brewer.pal, qual_col_pals$maxcolors,
        • line 150: rownames(qual_col_pals)))[seq_len(n)]
        • line 151: names(col.coeff) <- meta_mat$Coef
        • line 152: row_ha <- rowAnnotation(Effect = meta_mat[,
        • line 153: "Coef"], col = list(Effect = col.coeff),
        • line 154: show_legend = FALSE)
        • line 155: coeff_col_func <- colorRamp2(c(bio.lim[1],
        • line 156: 0, bio.lim[2]), c("blue", "gray0", "darkorange1"))
        • line 157: jpeg(file.path(paths["graph"], paste0("Heatmap_coeffs_",
        • line 158: pre.post[i], bio.sys[j], "_2000x1500.jpeg")),
        • line 159: height = 2000, width = 1500, res = 250)
        • line 160: draw(Heatmap(bio.dat, name = paste0("MLE"),
        • line 161: column_title = paste0(pre.post[i], "well.adj."),
        • line 162: row_title = bio.sys[j], right_annotation = row_ha,
        • line 163: col = coeff_col_func, column_dend_height = unit(1,
        • line 164: "cm"), cluster_rows = FALSE, cluster_columns = FALSE,
        • line 165: clustering_method_rows = "average", clustering_method_columns = "average",
        • line 166: show_column_names = TRUE, show_row_names = TRUE))
        • line 167: dev.off()
        • line 168: }
        • line 169: else {
        • line 176: col.coeff <- c(brewer.pal(9, "Set1")[c(1,
        • line 177: 3, 2)], brewer.pal(8, "Dark2")[c(4, 1,
        • line 178: 3)])
        • line 179: row_ha <- rowAnnotation(Effect = meta_mat[,
        • line 180: "Coef"], col = list(Effect = c(logn.mu.Plate = col.coeff[4],
        • line 182: coeff_col_func <- colorRamp2(c(well.lim[1],
        • line 183: 0, well.lim[2]), c("blue", "gray0", "darkorange1"))
        • line 184: jpeg(file.path(paths["graph"], paste0("Heatmap_coeffs_",
        • line 185: pre.post[i], bio.sys[j], "_2000x1500.jpeg")),
        • line 186: height = 2000, width = 1500, res = 250)
        • line 187: draw(Heatmap(well.dat, name = paste0("MLE"),
        • line 188: column_title = paste0(pre.post[i], "well.adj."),
        • line 189: row_title = bio.sys[j], right_annotation = row_ha,
        • line 190: col = coeff_col_func, column_dend_height = unit(1,
        • line 191: "cm"), cluster_rows = FALSE, cluster_columns = FALSE,
        • line 192: clustering_method_rows = "average", clustering_method_columns = "average",
        • line 193: show_column_names = TRUE, show_row_names = TRUE))
        • line 194: dev.off()
        • line 195: }
        • line 196: }
        • line 197: }
    • repetition in fcs_to_rds, fcs_to_rds_bkb, MapfxFFC, and MapfxMPC
      • in fcs_to_rds
        • line 1: file_meta, yvar)
        • line 2:{
        • line 3: fs <- read.flowSet(path = file.path(paste0(paths["input"],
        • line 4: "/fcs/")))
        • line 5: name.desc.fs <- pData(parameters(fs[[1]]))
        • line 6: name.desc.fs$desc[which(is.na(name.desc.fs$desc))] <- name.desc.fs$name[which(is.na(name.desc.fs$desc))]
        • line 7: fcs.raw.mt.list <- fcs.raw.meta.list <- list()
        • line 8: if (file_meta == "auto") {
        • line 9: for (i in seq_along(fs)) {
        • line 10: fcs.raw.mt.list[[i]] <- exprs(fs[[i]])
        • line 11: fcs.raw.meta.list[[i]] <- cbind(NO.in.well = seq_len(dim(fcs.raw.mt.list[[i]])[1]),
        • line 12: Filenam = unlist(keyword(fs[[i]], "GUID")))
        • line 13: }
        • line 14: fcs.raw.mt <- do.call(rbind, fcs.raw.mt.list)
        • line 15: colnames(fcs.raw.mt) <- name.desc.fs$desc[match(colnames(fcs.raw.mt),
        • line 16: name.desc.fs$name)]
        • line 17: fcs.raw.meta.df <- as.data.frame(do.call(rbind, fcs.raw.meta.list))
        • line 18: fcs.raw.meta.df$NO.in.well <- as.integer(fcs.raw.meta.df$NO.in.well)
        • line 19: fcs.raw.meta.df$Plate <- str_extract(fcs.raw.meta.df$Filenam,
        • line 20: "Plate[0-9]")
        • line 21: fcs.raw.meta.df$Well <- str_extract(fcs.raw.meta.df$Filenam,
        • line 22: "[A-H]\d{1,2}")
        • line 23: orig.lab <- c(paste0(rep(c("A", "B", "C", "D", "E", "F",
        • line 24: "G", "H"), each = 9), c(seq_len(9))), paste0(rep(c("A",
        • line 25: "B", "C", "D", "E", "F", "G", "H"), each = 9), c(10:12)))
        • line 29: map <- setNames(new.lab, orig.lab)
        • line 30: fcs.raw.meta.df$Well <- map[fcs.raw.meta.df$Well]
        • line 31: fcs.raw.meta.df$Column <- paste0("Col.", sub("[A-H]",
        • line 32: "", fcs.raw.meta.df$Well))
        • line 33: fcs.raw.meta.df$Row <- sub("\d{1,2}", "", fcs.raw.meta.df$Well)
        • line 34: orig.lab <- c("A", "B", "C", "D", "E", "F", "G", "H")
        • line 35: new.lab <- paste0("Row.", c("01", "02", "03", "04", "05",
        • line 36: "06", "07", "08"))
        • line 37: map <- setNames(new.lab, orig.lab)
        • line 38: fcs.raw.meta.df$Row <- map[fcs.raw.meta.df$Row]
        • line 39: fcs.raw.meta.df$Well.lab <- paste0(sub("late", "", fcs.raw.meta.df$Plate),
        • line 40: "_", fcs.raw.meta.df$Well)
        • line 41: if (sum(is.na(fcs.raw.meta.df)) != 0) {
        • line 42: message("\tfile_meta='auto' doesn't work on your data.")
        • line 43: stop("\tPlease provide a 'filename_meta.csv' file and use file_meta='usr' instead.")
        • line 44: }
        • line 45: ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),
        • line 46: ]
        • line 47: ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),
        • line 48: ]
        • line 49: }
        • line 50: if (file_meta == "usr") {
        • line 51: meta.inf <- read.csv(file = file.path(paths["input"],
        • line 52: "/meta/filename_meta.csv"))
        • line 53: for (i in seq_along(fs)) {
        • line 54: fcs.raw.mt.list[[i]] <- exprs(fs[[i]])
        • line 55: filenam <- unlist(keyword(fs[[i]], "GUID"))
        • line 56: fcs.raw.meta.list[[i]] <- cbind(NO.in.well = (seq_len(dim(fcs.raw.mt.list[[i]])[1])),
        • line 57: Filenam = filenam, Plate = meta.inf[which(filenam ==
        • line 58: meta.inf$Filenam), "Plate"], Well = meta.inf[which(filenam ==
        • line 59: meta.inf$Filenam), "Well"], Column = meta.inf[which(filenam ==
        • line 60: meta.inf$Filenam), "Column"], Row = meta.inf[which(filenam ==
        • line 61: meta.inf$Filenam), "Row"], Well.lab = meta.inf[which(filenam ==
        • line 62: meta.inf$Filenam), "Well.lab"])
        • line 63: }
        • line 64: fcs.raw.mt <- do.call(rbind, fcs.raw.mt.list)
        • line 65: colnames(fcs.raw.mt) <- name.desc.fs$desc[match(colnames(fcs.raw.mt),
        • line 66: name.desc.fs$name)]
        • line 67: fcs.raw.meta.df <- as.data.frame(do.call(rbind, fcs.raw.meta.list))
        • line 68: fcs.raw.meta.df$NO.in.well <- as.integer(fcs.raw.meta.df$NO.in.well)
        • line 69: ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),
        • line 70: ]
        • line 71: ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),
        • line 72: ]
        • line 73: }
        • line 74: ord.fcs.raw.meta.df$NO.in.all <- seq_len(nrow(ord.fcs.raw.mt))
        • line 78: 3, 5, 6, 4, 7)]
        • line 79: saveRDS(ord.fcs.raw.meta.df.out, file = file.path(paths["intermediary"],
        • line 80: "/fcs_metadata_df.rds"))
        • line 81: saveRDS(ord.fcs.raw.mt, file = file.path(paths["intermediary"],
        • line 82: "/fcs_rawInten_mt.rds"))
      • in fcs_to_rds_bkb
        • line 1: file_meta, MPC)
        • line 2:{
        • line 3: fs <- read.flowSet(path = file.path(paste0(paths["input"],
        • line 4: "/fcs/")))
        • line 5: name.desc.fs <- pData(parameters(fs[[1]]))
        • line 6: name.desc.fs$desc[which(is.na(name.desc.fs$desc))] <- name.desc.fs$name[which(is.na(name.desc.fs$desc))]
        • line 7: fcs.raw.mt.list <- fcs.raw.meta.list <- list()
        • line 8: if (MPC == TRUE) {
        • line 9: if (file_meta == "auto") {
        • line 10: for (i in seq_along(fs)) {
        • line 11: fcs.raw.mt.list[[i]] <- exprs(fs[[i]])
        • line 12: fcs.raw.meta.list[[i]] <- cbind(NO.in.well = (seq_len(dim(fcs.raw.mt.list[[i]])[1])),
        • line 13: Filenam = unlist(keyword(fs[[i]], "GUID")))
        • line 14: }
        • line 15: fcs.raw.mt <- do.call(rbind, fcs.raw.mt.list)
        • line 16: colnames(fcs.raw.mt) <- name.desc.fs$desc[match(colnames(fcs.raw.mt),
        • line 17: name.desc.fs$name)]
        • line 18: fcs.raw.meta.df <- as.data.frame(do.call(rbind, fcs.raw.meta.list))
        • line 19: fcs.raw.meta.df$NO.in.well <- as.integer(fcs.raw.meta.df$NO.in.well)
        • line 20: fcs.raw.meta.df$Plate <- str_extract(fcs.raw.meta.df$Filenam,
        • line 21: "Plate[0-9]")
        • line 22: fcs.raw.meta.df$Well <- str_extract(fcs.raw.meta.df$Filenam,
        • line 23: "[A-H]\d{1,2}")
        • line 24: orig.lab <- c(paste0(rep(c("A", "B", "C", "D", "E",
        • line 25: "F", "G", "H"), each = 9), seq_len(9)), paste0(rep(c("A",
        • line 26: "B", "C", "D", "E", "F", "G", "H"), each = 9),
        • line 27: c(10:12)))
        • line 28: new.lab <- c(paste0(rep(c("A", "B", "C", "D", "E",
        • line 29: "F", "G", "H"), each = 9), 0, seq_len(9)), paste0(rep(c("A",
        • line 30: "B", "C", "D", "E", "F", "G", "H"), each = 9),
        • line 31: c(10:12)))
        • line 32: map <- setNames(new.lab, orig.lab)
        • line 33: fcs.raw.meta.df$Well <- map[fcs.raw.meta.df$Well]
        • line 34: fcs.raw.meta.df$Column <- paste0("Col.", sub("[A-H]",
        • line 35: "", fcs.raw.meta.df$Well))
        • line 36: fcs.raw.meta.df$Row <- sub("\d{1,2}", "", fcs.raw.meta.df$Well)
        • line 37: orig.lab <- c("A", "B", "C", "D", "E", "F", "G",
        • line 38: "H")
        • line 39: new.lab <- paste0("Row.", c("01", "02", "03", "04",
        • line 40: "05", "06", "07", "08"))
        • line 41: map <- setNames(new.lab, orig.lab)
        • line 42: fcs.raw.meta.df$Row <- map[fcs.raw.meta.df$Row]
        • line 43: fcs.raw.meta.df$Well.lab <- paste0(sub("late", "",
        • line 44: fcs.raw.meta.df$Plate), "_", fcs.raw.meta.df$Well)
        • line 45: if (sum(is.na(fcs.raw.meta.df)) != 0) {
        • line 46: message("\tfile_meta='auto' doesn't work on your data.")
        • line 47: stop("\tPlease provide a 'filename_meta.csv' file and use file_meta='usr' instead.")
        • line 48: }
        • line 49: ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),
        • line 50: ]
        • line 51: ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),
        • line 52: ]
        • line 53: }
        • line 54: if (file_meta == "usr") {
        • line 55: meta.inf <- read.csv(file = file.path(paths["input"],
        • line 56: "/meta/filename_meta.csv"))
        • line 57: for (i in seq_along(fs)) {
        • line 58: fcs.raw.mt.list[[i]] <- exprs(fs[[i]])
        • line 59: filenam <- unlist(keyword(fs[[i]], "GUID"))
        • line 60: fcs.raw.meta.list[[i]] <- cbind(NO.in.well = (seq_len(dim(fcs.raw.mt.list[[i]])[1])),
        • line 61: Filenam = filenam, Plate = meta.inf[which(filenam ==
        • line 62: meta.inf$Filenam), "Plate"], Well = meta.inf[which(filenam ==
        • line 63: meta.inf$Filenam), "Well"], Column = meta.inf[which(filenam ==
        • line 64: meta.inf$Filenam), "Column"], Row = meta.inf[which(filenam ==
        • line 65: meta.inf$Filenam), "Row"], Well.lab = meta.inf[which(filenam ==
        • line 66: meta.inf$Filenam), "Well.lab"])
        • line 67: }
        • line 68: fcs.raw.mt <- do.call(rbind, fcs.raw.mt.list)
        • line 69: colnames(fcs.raw.mt) <- name.desc.fs$desc[match(colnames(fcs.raw.mt),
        • line 70: name.desc.fs$name)]
        • line 71: fcs.raw.meta.df <- as.data.frame(do.call(rbind, fcs.raw.meta.list))
        • line 72: fcs.raw.meta.df$NO.in.well <- as.integer(fcs.raw.meta.df$NO.in.well)
        • line 73: ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Well.lab),
        • line 74: ]
        • line 75: ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Well.lab),
        • line 76: ]
        • line 77: }
        • line 78: ord.fcs.raw.meta.df$NO.in.all <- seq_len(nrow(ord.fcs.raw.mt))
        • line 81: }
        • line 82: if (MPC == FALSE & file_meta == "usr") {
        • line 83: meta.inf <- read.csv(file = file.path(paths["input"],
        • line 84: "/meta/filename_meta.csv"))
        • line 85: for (i in seq_along(fs)) {
        • line 86: fcs.raw.mt.list[[i]] <- exprs(fs[[i]])
        • line 87: filenam <- unlist(keyword(fs[[i]], "GUID"))
        • line 88: fcs.raw.meta.list[[i]] <- cbind(NO.in.batch = (seq_len(dim(fcs.raw.mt.list[[i]])[1])),
        • line 89: Filenam = filenam, Batch = meta.inf[which(filenam ==
        • line 90: meta.inf$Filenam), "Batch"])
        • line 91: }
        • line 92: fcs.raw.mt <- do.call(rbind, fcs.raw.mt.list)
        • line 93: colnames(fcs.raw.mt) <- name.desc.fs$desc[match(colnames(fcs.raw.mt),
        • line 94: name.desc.fs$name)]
        • line 95: fcs.raw.meta.df <- as.data.frame(do.call(rbind, fcs.raw.meta.list))
        • line 96: fcs.raw.meta.df$NO.in.batch <- as.integer(fcs.raw.meta.df$NO.in.batch)
        • line 97: ord.fcs.raw.meta.df <- fcs.raw.meta.df[order(fcs.raw.meta.df$Batch),
        • line 98: ]
        • line 99: ord.fcs.raw.mt <- fcs.raw.mt[order(fcs.raw.meta.df$Batch),
        • line 100: ]
        • line 101: ord.fcs.raw.meta.df$NO.in.all <- seq_len(nrow(ord.fcs.raw.mt))
        • line 107: }
        • line 108: saveRDS(ord.fcs.raw.meta.df.out, file = file.path(paths["intermediary"],
        • line 109: "/fcs_metadata_df.rds"))
        • line 110: saveRDS(ord.fcs.raw.mt, file = file.path(paths["intermediary"],
        • line 111: "/fcs_rawInten_mt.rds"))
      • in MapfxFFC
        • line 6:{
        • line 7: if (runVignette == TRUE) {
        • line 8: if (is.null(Outpath)) {
        • line 9: message("\nPlease make sure the output path has been specified!")
        • line 10: }
        • line 11: message("\n\n\nCreating directories for output...")
        • line 12: settings <- initialize(path_to_fcs = FCSpath, path_to_output = Outpath,
        • line 13: verbose = TRUE)
        • line 14: paths <- settings$paths
        • line 15: saveRDS(runVignette_meta, file = file.path(paths["intermediary"],
        • line 16: "/fcs_metadata_df.rds"))
        • line 17: saveRDS(runVignette_rawInten, file = file.path(paths["intermediary"],
        • line 18: "/fcs_rawInten_mt.rds"))
        • line 19: message("\n\n\nBackground correcting proteins...")
        • line 20: bkc_bkb_ffc(paths = paths, bkb.v = protein.v, bkb.upper.quantile = protein.upper.quantile,
        • line 21: bkb.lower.quantile = protein.lower.quantile, bkb.min.quantile = protein.min.quantile,
        • line 23: message("\n\n\nForming a matrix of biology (M) for removal of batch effect...")
        • line 24: initM(paths = paths, assay = "FFC", bkb.v = protein.v,
        • line 25: plots = plots.initM)
        • line 26: message("\n\n\nRemoval of batch effect...")
        • line 27: rmBatchEffect(paths = paths, plots = plots.rmBatchEffect)
        • line 30: cluster.analysis.bkbOnly(paths = paths, bkb.v = protein.v,
        • line 31: plots = plots.cluster.analysis.protein)
        • line 32: }
        • line 33: message("\tCompleted!")
        • line 34: }
        • line 35: if (runVignette == FALSE) {
        • line 36: if (runVignette == FALSE | is.null(FCSpath) | is.null(Outpath) |
        • line 37: is.null(protein.v)) {
        • line 38: message("\nPlease make sure you have specified\n(1) path to the FCS files, and metadata (if file_meta="usr") \n(2) output path \n(3) a list of protein markers")
        • line 39: }
        • line 40: if (!is.null(FCSpath) & !is.null(Outpath) & !is.null(protein.v)) {
        • line 41: message("\n\n\nCreating directories for output...")
        • line 42: settings <- initialize(path_to_fcs = FCSpath, path_to_output = Outpath,
        • line 43: verbose = TRUE)
        • line 44: paths <- settings$paths
        • line 45: message("\n\n\nTransforming FCS to RDS files...")
        • line 46: fcs_to_rds_bkb(paths = paths, file_meta = file_meta,
        • line 47: MPC = FALSE)
        • line 48: message("\n\n\nBackground correcting proteins...")
        • line 51: bkb.min.quantile = protein.min.quantile, plots = plots.bkc.protein)
        • line 52: message("\n\n\nForming a matrix of biology for removal of batch effect...")
        • line 53: initM(paths = paths, assay = "FFC", bkb.v = protein.v,
        • line 54: plots = plots.initM)
        • line 55: message("\n\n\nRemoval of batch effect...")
        • line 56: rmBatchEffect(paths = paths, plots = plots.rmBatchEffect)
        • line 59: cluster.analysis.bkbOnly(paths = paths, bkb.v = protein.v,
        • line 60: plots = plots.cluster.analysis.protein)
        • line 61: }
        • line 62: message("\tCompleted!")
        • line 63: }
        • line 64: }
      • in MapfxMPC
        • line 12:{
        • line 13: if (runVignette == TRUE) {
        • line 14: if (is.null(Outpath)) {
        • line 15: message("\nPlease make sure the output path has been specified!")
        • line 16: }
        • line 17: chans <- bkb.v
        • line 18: message("\n\n\nCreating directories for output...")
        • line 19: settings <- initialize(path_to_fcs = FCSpath, path_to_output = Outpath,
        • line 20: verbose = TRUE)
        • line 21: paths <- settings$paths
        • line 22: saveRDS(runVignette_meta, file = file.path(paths["intermediary"],
        • line 23: "/fcs_metadata_df.rds"))
        • line 24: saveRDS(runVignette_rawInten, file = file.path(paths["intermediary"],
        • line 25: "/fcs_rawInten_mt.rds"))
        • line 26: message("\n\n\nBackground correcting backbone markers...")
        • line 27: bkc_bkb(paths = paths, bkb.v = bkb.v, bkb.upper.quantile = bkb.upper.quantile,
        • line 28: bkb.lower.quantile = bkb.lower.quantile, bkb.min.quantile = bkb.min.quantile,
        • line 54: cluster.analysis.bkbOnly(paths = paths, bkb.v = bkb.v,
        • line 55: plots = plots.cluster.analysis.bkb)
        • line 56: }
        • line 57: message("\tCompleted!")
        • line 58: }
        • line 59: if (runVignette == FALSE) {
        • line 60: if (runVignette == FALSE | is.null(FCSpath) | is.null(Outpath) |
        • line 61: is.null(bkb.v)) {
        • line 62: message("\nPlease make sure you have specified\n(1) path to the FCS files, and metadata (if file_meta="usr") \n(2) output path \n(3) a list of backbone markers")
        • line 63: }
        • line 64: if (!is.null(FCSpath) & !is.null(Outpath) & !is.null(bkb.v)) {
        • line 71: }
        • line 72: message("\n\n\nCreating directories for output...")
        • line 73: settings <- initialize(path_to_fcs = FCSpath, path_to_output = Outpath,
        • line 74: verbose = TRUE)
        • line 75: paths <- settings$paths
        • line 76: message("\n\n\nTransforming FCS to RDS files...")
        • line 77: fcs_to_rds(paths = paths, file_meta = file_meta,
        • line 78: yvar = yvar)
        • line 79: message("\n\n\nBackground correcting backbone markers...")
        • line 85: pe.min.quantile = inf.min.quantile, plots = plots.bkc.inf)
        • line 86: message("\n\n\nForming a matrix of biology (M) for removal of well effect...")
        • line 87: initM(paths = paths, assay = "MPC", bkb.v = bkb.v,
        • line 88: plots = plots.initM)
        • line 89: message("\n\n\nRemoval of well effect...")
        • line 90: rmWellEffect(paths = paths, plots = plots.rmWellEffect)
        • line 103: plots = plots.cluster.analysis.all)
        • line 104: }
        • line 105: }
        • line 106: message("\tCompleted!")
        • line 107: }
        • line 108: }

jianhong avatar Apr 18 '24 21:04 jianhong

Received a valid push on git.bioconductor.org; starting a build for commit id: 7d86518a2e963f45392e941b81f207f4bac0b772

bioc-issue-bot avatar Apr 22 '24 11:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

Congratulations! The package built without errors or warnings on all platforms.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.2.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 22 '24 12:04 bioc-issue-bot

Hi @jianhong! We appreciate your comments. Please see below our answers.

Code: Note: please consider; Important: must be addressed.

The NAMESPACE file

  • [v] Important: Selective imports using importFrom instead of import all with import.

Fixed!

General package development

  • [v] Important: Consider adding unit tests. We strongly encourage them. See http://bioconductor.org/developers/how-to/unitTesting-guidelines/.

We acknowledge that it would be more rigorous to have unit tests implemented, however, the reason why we did not add unit tests is because the two exported functions MapfxMPC and MapfxFFC are meant to be end-to-end pipelines which incorporate several non-exported functions to produce several outputs (normalised data and graphs, etc) and the outputs are directly saved in the specified directory. Both MapfxMPC and MapfxFFC do not return R object(s) which makes unit tests difficult to conduct. We will think about this for the updated version in the future.

R code

  • [v] Important: warning, message, stop instead of cat and print outside of show methods.

Fixed!

  • [v] NOTE: :: is not suggested in source code unless you can make sure all the packages are imported. Some people think it is better to keep ::. However, please be aware that you will need to manually double-check the imported items if you make any changes to the DESCRIPTION file during development. My suggestion is to remove one or two repetitions to trigger the dependency check.

Among the nine lines in your suggestion, I retained most of the :: and have made sure all the packages are imported. Also, I removed repetitions (line 323 and line 330 in my updated R/051.imputation_bkb.predictors.R) to check the dependency. The DESCRIPTION file was modified accordingly.

  • [v] NOTE: Vectorize: for loops present, try to replace them by *apply functions.

Thanks for your suggestion, however, we use for loops for those with complicated calculations. We will try to optimise some of them for the updated version in the future.

  • [v] Important: Use file.path to replace paste. Do not include the '/' in the file path.

Fixed.

  • [v] Important: Remove unused code.

Fixed.

  • [v] NOTE: Avoid 'suppressWarnings'/'*Messages' if possible

These functions are used to estimate parameters for background correction. The warning messages might appear when the amount of cells in each well (or batch) is extremely small and it doesn't happen a lot in real data analysis. We will try to resolve this for the updated version in the future.

  • [v] Important: Please consider to remove ';' in the code to avoid mutiple scripts in one line.

Fixed.

  • [v] Important: Please consider to add drop=FALSE to avoid the reduction of dimension for matrices and arrays. Ignore this if using datatable.

Thank you. Fixed.

  • [v] NOTE: Try to check the edge condition when using seq.int or seq_len. For example using seq.int(min(5, nrow(data))) to replace seq.int(5)

The current code is designed for data from the MPC experiments that have cells from 3 plates, 12 columns, and 8 rows, so the edge condition would not change. However, we appreciate your suggestion as it can make our code more flexible. We've made notes for the updated version in the future.

  • [v] NOTE: Functional programming: code repetition.

Some repetitions are necessary, for example, in different functions, we may read and save files using the same file names and use the same default parameters for data adjustment. Other code repetition is mainly because the functions share the core methodology with some modifications for different purposes (e.g., different experimental data or different types of markers for visualisation). Thank you for the pick-up. We will try to optimise the code for the updated version in the future.

HsiaoChiLiao avatar Apr 22 '24 13:04 HsiaoChiLiao

Package 'MAPFX' Review

The package passed check and build. It improved a lot. However there are several things need to be fixed. Please try to answer the comments line by line when you are ready for a second review. Code: Note: please consider; Important: must be addressed.

General package development

  • [ ] Important: Consider adding unit tests. We strongly encourage them. See http://bioconductor.org/developers/how-to/unitTesting-guidelines/. I mean this is required.

R code

  • [ ] Important: Undefined symbols in input parameter.
    • In file R/00_initialize.R:
      • at line 12 found ' path_to_fcs=path_to_fcs, '
      • at line 13 found ' path_to_output=path_to_output, '
  • [ ] Important: Missing input check for FCSpath, Outpath and file_meta for function MapfxMPC and MapfxFFC.
  • [ ] Important: Avoid to include the '/' in the file path.
    • In file R/00_initialize.R:
      • at line 18 found ' path_to_intermediary <- file.path(path_to_output,"/intermediary")'
      • at line 19 found ' path_to_graph <- file.path(path_to_output,"/graph")'
      • at line 20 found ' path_to_downstream <- file.path(path_to_output,"/downstream")'
      • at line 28 found ' input=file.path(path_to_fcs,"/"),'
      • at line 29 found ' intermediary=file.path(path_to_intermediary,"/"),'
      • at line 30 found ' graph=file.path(path_to_graph,"/"),'
      • at line 31 found ' downstream=file.path(path_to_downstream,"/")'
    • In file R/000.MapfxFFC.R:
      • at line 78 found ' saveRDS(runVignette_meta, file = file.path(paths["intermediary"], "/fcs_metadata_df.rds"))'
      • at line 79 found ' saveRDS(runVignette_rawInten, file = file.path(paths["intermediary"], "/fcs_rawInten_mt.rds"))'
    • In file R/000.MapfxMPC.R:
      • at line 120 found ' saveRDS(runVignette_meta, file = file.path(paths["intermediary"], "/fcs_metadata_df.rds"))'
      • at line 121 found ' saveRDS(runVignette_rawInten, file = file.path(paths["intermediary"], "/fcs_rawInten_mt.rds"))'
    • In file R/011.fcs_to_rds.R:
      • at line 92 found ' meta.inf <- read.csv(file = file.path(paths["input"], "/meta/filename_meta.csv"))'
      • at line 127 found ' saveRDS(ord.fcs.raw.meta.df.out, file = file.path(paths["intermediary"], "/fcs_metadata_df.rds"))'
      • at line 128 found ' saveRDS(ord.fcs.raw.mt, file = file.path(paths["intermediary"], "/fcs_rawInten_mt.rds"))'
    • In file R/012.fcs_to_rds_bkb.R:
      • at line 92 found ' meta.inf <- read.csv(file = file.path(paths["input"], "/meta/filename_meta.csv"))'
      • at line 127 found ' meta.inf <- read.csv(file = file.path(paths["input"], "/meta/filename_meta.csv"))'
      • at line 158 found ' saveRDS(ord.fcs.raw.meta.df.out, file = file.path(paths["intermediary"], "/fcs_metadata_df.rds"))'
      • at line 159 found ' saveRDS(ord.fcs.raw.mt, file = file.path(paths["intermediary"], "/fcs_rawInten_mt.rds"))'
    • In file R/031.initM.R:
      • at line 159 found ' jpeg(filename = file.path(paths["graph"], "/", paste0(trans.dat, "_UMAP_colPhenog_600x750.jpeg")), height = 600, width = 750, res = 80)'
  • [ ] NOTE: Functional programming: code repetition. It will be appriciated to make change during this review if possible.
  • [ ] Important: Remove unused code:
    • In file R/032.rmWellEffect.R:
      • at line 70 found ' b-a'
      • at line 86 found ' b-a'
      • at line 125 found ' b-a'
    • In file R/033.rmBatchEffect.R:
      • at line 69 found ' b-a'
      • at line 86 found ' b-a'
      • at line 123 found ' b-a'
    • In file R/060.cluster.analysis.R:
      • at line 55 found ' b-a'
      • at line 62 found ' b-a # 5.802254 mins'
      • at line 75 found ' b-a # 11.9143 mins'
      • at line 84 found ' b-a # 17.67467 mins'
    • In file R/070.cluster.analysis_bkbOnly.R:
      • at line 46 found ' b-a # 12.46998 mins'
      • at line 53 found ' b-a # 5.802254 mins'
    • In file R/031.initM.R:
      • at line 83 found ' head(sel.bkb.lgc)'
    • In file R/060.cluster.analysis.R:
      • at line 142 found ' head(metadata.cell)'
    • In file R/070.cluster.analysis_bkbOnly.R:
      • at line 85 found ' head(metadata.cell)'

Documentation

  • [ ] Important: Update the return section for most of the functions to real returns and change the current description of return to details.
  • [ ] Important: Put a script or a description file into inst folder to describe how you generate the data.

jianhong avatar Apr 22 '24 14:04 jianhong

Received a valid push on git.bioconductor.org; starting a build for commit id: 2445c9c17622723ad32f52590e5e089d5b822a87

bioc-issue-bot avatar Apr 23 '24 09:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

On one or more platforms, the build results were: "ERROR". This may mean there is a problem with the package that you need to fix. Or it may mean that there is a problem with the build system itself.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.3.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 23 '24 10:04 bioc-issue-bot

Received a valid push on git.bioconductor.org; starting a build for commit id: 21f5b11703578f41e88477b9b0a04abbc92d587f

bioc-issue-bot avatar Apr 23 '24 11:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

On one or more platforms, the build results were: "TIMEOUT". This may mean there is a problem with the package that you need to fix. Or it may mean that there is a problem with the build system itself.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.4.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 23 '24 11:04 bioc-issue-bot

Received a valid push on git.bioconductor.org; starting a build for commit id: 70b2a07b825893c6c9f2fb91d2900f85b3111baf

bioc-issue-bot avatar Apr 23 '24 12:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

Congratulations! The package built without errors or warnings on all platforms.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.5.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 23 '24 12:04 bioc-issue-bot

Received a valid push on git.bioconductor.org; starting a build for commit id: 4a1669f71191ca6deaeaf2c44a8446f72d7121ed

bioc-issue-bot avatar Apr 23 '24 12:04 bioc-issue-bot

Dear Package contributor,

This is the automated single package builder at bioconductor.org.

Your package has been built on the Bioconductor Single Package Builder.

Congratulations! The package built without errors or warnings on all platforms.

Please see the build report for more details.

The following are build products from R CMD build on the Single Package Builder: Linux (Ubuntu 22.04.3 LTS): MAPFX_0.99.6.tar.gz

Links above active for 21 days.

Remember: if you submitted your package after July 7th, 2020, when making changes to your repository push to [email protected]:packages/MAPFX to trigger a new build. A quick tutorial for setting up remotes and pushing to upstream can be found here.

bioc-issue-bot avatar Apr 23 '24 12:04 bioc-issue-bot

Hi @jianhong We appreciate your prompt feedback and comments. We've tried to resolve as many things as possible. Please see below our answers, and let us know if there's anything we should further improve during this review.

Code: Note: please consider; Important: must be addressed.

General package development

  • [v] Important: Consider adding unit tests. We strongly encourage them. See http://bioconductor.org/developers/how-to/unitTesting-guidelines/. I mean this is required.

Unit tests added.

R code

  • [v] Important: Undefined symbols in input parameter.

Fixed.

  • [v] Important: Missing input check for FCSpath, Outpath and file_meta for function MapfxMPC and MapfxFFC.

Input check added: lines 93-95 in R/000.MapfxFFC.R lines 145-153 in R/000.MapfxMPC.R Note that we removed the argument file_meta for the function MapfxFFC as the users must provide their filename_meta.csv file in their FCSpath/meta directory.

  • [v] Important: Avoid to include the '/' in the file path.

Fixed.

  • [v] NOTE: Functional programming: code repetition. It will be appriciated to make change during this review if possible.

To reduce code repetition, we've made below changes:

  1. Modified R/000.MapfxFFC.R and R/000.MapfxMPC.R to remove some repetition.
  2. Removed repeated code in R/012.fcs_to_rds_bkb.R.
  3. Modified R/021.bkc_bkb.R and removed R/023.023.bkc_bkb_ffc.R to reduce repetition.
  • [v] Important: Remove unused code:

Fixed.

Documentation

  • [v] Important: Update the return section for most of the functions to real returns and change the current description of return to details.

Fixed. Both the exported functions MapfxMPC and MapfxFFC now have a returned object. Updated "@return" section and created "@details" to describe the files saved in the output directory in more detail.

  • [v] Important: Put a script or a description file into inst folder to describe how you generate the data.

Added a script inst/GenerateExampleData.R that was used to generate the example data.

HsiaoChiLiao avatar Apr 23 '24 12:04 HsiaoChiLiao

Your package has been accepted. It will be added to the Bioconductor nightly builds.

Thank you for contributing to Bioconductor!

Reviewers for Bioconductor packages are volunteers from the Bioconductor community. If you are interested in becoming a Bioconductor package reviewer, please see Reviewers Expectations.

bioc-issue-bot avatar Apr 23 '24 13:04 bioc-issue-bot

The default branch of your GitHub repository has been added to Bioconductor's git repository as branch devel.

To use the git.bioconductor.org repository, we need an 'ssh' key to associate with your github user name. If your GitHub account already has ssh public keys (https://github.com/HsiaoChiLiao.keys is not empty), then no further steps are required. Otherwise, do the following:

  1. Add an SSH key to your github account
  2. Submit your SSH key to Bioconductor

See further instructions at

https://bioconductor.org/developers/how-to/git/

for working with this repository. See especially

https://bioconductor.org/developers/how-to/git/new-package-workflow/ https://bioconductor.org/developers/how-to/git/sync-existing-repositories/

to keep your GitHub and Bioconductor repositories in sync.

Your package will be included in the next nigthly 'devel' build (check-out from git at about 6 pm Eastern; build completion around 2pm Eastern the next day) at

https://bioconductor.org/checkResults/

(Builds sometimes fail, so ensure that the date stamps on the main landing page are consistent with the addition of your package). Once the package builds successfully, you package will be available for download in the 'Devel' version of Bioconductor using BiocManager::install("MAPFX"). The package 'landing page' will be created at

https://bioconductor.org/packages/MAPFX

If you have any questions, please contact the bioc-devel mailing list (https://stat.ethz.ch/mailman/listinfo/bioc-devel); this issue will not be monitored further.

lshep avatar Apr 25 '24 13:04 lshep