Mercurial > repos > iuc > raceid_clustering
comparison scripts/cluster.R @ 10:49776718ae90 draft default tip
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/raceid3 commit 0ffa71ef9f8d020fe7ba94502db8cec26fd8741f
| author | iuc | 
|---|---|
| date | Tue, 05 Nov 2024 16:33:40 +0000 | 
| parents | 0bff0ee0683a | 
| children | 
   comparison
  equal
  deleted
  inserted
  replaced
| 9:0bff0ee0683a | 10:49776718ae90 | 
|---|---|
| 2 VERSION <- "0.5" # nolint | 2 VERSION <- "0.5" # nolint | 
| 3 | 3 | 
| 4 args <- commandArgs(trailingOnly = TRUE) | 4 args <- commandArgs(trailingOnly = TRUE) | 
| 5 | 5 | 
| 6 if (length(args) != 1) { | 6 if (length(args) != 1) { | 
| 7 message(paste("VERSION:", VERSION)) | 7 message(paste("VERSION:", VERSION)) | 
| 8 stop("Please provide the config file") | 8 stop("Please provide the config file") | 
| 9 } | 9 } | 
| 10 | 10 | 
| 11 suppressWarnings(suppressPackageStartupMessages(require(RaceID))) | 11 suppressWarnings(suppressPackageStartupMessages(require(RaceID))) | 
| 12 ## suppressWarnings(suppressPackageStartupMessages(require(scran))) # nolint | 12 ## suppressWarnings(suppressPackageStartupMessages(require(scran))) # nolint | 
| 13 source(args[1]) | 13 source(args[1]) | 
| 15 | 15 | 
| 16 do.filter <- function(sc) { # nolint | 16 do.filter <- function(sc) { # nolint | 
| 17 if (!is.null(filt.lbatch.regexes)) { | 17 if (!is.null(filt.lbatch.regexes)) { | 
| 18 lar <- filt.lbatch.regexes | 18 lar <- filt.lbatch.regexes | 
| 19 nn <- colnames(sc@expdata) | 19 nn <- colnames(sc@expdata) | 
| 20 filt$LBatch <- lapply(1:length(lar), function(m) { # nolint | 20 filt$LBatch <- lapply(1:length(lar), function(m) { # nolint | 
| 21 return(nn[grep(lar[[m]], nn)])}) | 21 return(nn[grep(lar[[m]], nn)]) | 
| 22 }) | |
| 22 } | 23 } | 
| 23 | 24 | 
| 24 sc <- do.call(filterdata, c(sc, filt)) | 25 sc <- do.call(filterdata, c(sc, filt)) | 
| 25 | 26 | 
| 26 ## Get histogram metrics for library size and number of features | 27 ## Get histogram metrics for library size and number of features | 
| 42 print(tmp) | 43 print(tmp) | 
| 43 ## required, for extracting midpoint | 44 ## required, for extracting midpoint | 
| 44 unq <- unique(filt_feat) | 45 unq <- unique(filt_feat) | 
| 45 if (length(unq) == 1) { | 46 if (length(unq) == 1) { | 
| 46 abline(v = unq, col = "red", lw = 2) | 47 abline(v = unq, col = "red", lw = 2) | 
| 47 text(tmp$mids, table(filt_feat)[[1]] - 100, pos = 1, | 48 text(tmp$mids, table(filt_feat)[[1]] - 100, | 
| 48 paste(10^unq, "\nFeatures\nin remaining\nCells", | 49 pos = 1, | 
| 49 sep = ""), cex = 0.8) | 50 paste(10^unq, "\nFeatures\nin remaining\nCells", | 
| 51 sep = "" | |
| 52 ), cex = 0.8 | |
| 53 ) | |
| 50 } | 54 } | 
| 51 | 55 | 
| 52 if (filt.use.ccorrect) { | 56 if (filt.use.ccorrect) { | 
| 53 par(mfrow = c(2, 2)) | 57 par(mfrow = c(2, 2)) | 
| 54 sc <- do.call(CCcorrect, c(sc, filt.ccc)) | 58 sc <- do.call(CCcorrect, c(sc, filt.ccc)) | 
| 78 print(plotsensitivity(sc)) | 82 print(plotsensitivity(sc)) | 
| 79 print(plotoutlierprobs(sc)) | 83 print(plotoutlierprobs(sc)) | 
| 80 ## Heatmaps | 84 ## Heatmaps | 
| 81 test1 <- list() | 85 test1 <- list() | 
| 82 test1$side <- 3 | 86 test1$side <- 3 | 
| 83 test1$line <- 0 #1 #3 | 87 test1$line <- 0 # 1 #3 | 
| 84 | 88 | 
| 85 x <- clustheatmap(sc, final = FALSE) | 89 x <- clustheatmap(sc, final = FALSE) | 
| 86 print(do.call(mtext, c(paste("(Initial)"), test1))) | 90 print(do.call(mtext, c(paste("(Initial)"), test1))) | 
| 87 x <- clustheatmap(sc, final = TRUE) | 91 x <- clustheatmap(sc, final = TRUE) | 
| 88 print(do.call(mtext, c(paste("(Final)"), test1))) | 92 print(do.call(mtext, c(paste("(Final)"), test1))) | 
| 120 buffer <- paste(rep("", 36), collapse = " ") | 124 buffer <- paste(rep("", 36), collapse = " ") | 
| 121 print(do.call(mtext, c(paste(buffer, "Cluster ", n), test))) | 125 print(do.call(mtext, c(paste(buffer, "Cluster ", n), test))) | 
| 122 test$line <- -1 | 126 test$line <- -1 | 
| 123 print(do.call(mtext, c(paste(buffer, "Sig. Genes"), test))) | 127 print(do.call(mtext, c(paste(buffer, "Sig. Genes"), test))) | 
| 124 test$line <- 0 | 128 test$line <- 0 | 
| 125 print(do.call(mtext, c(paste(buffer, "(fc > ", | 129 print(do.call(mtext, c(paste( | 
| 126 genelist.foldchange, ")"), test))) | 130 buffer, "(fc > ", | 
| 131 genelist.foldchange, ")" | |
| 132 ), test))) | |
| 127 }) | 133 }) | 
| 128 write.table(df, file = out.genelist, sep = "\t", quote = FALSE) | 134 write.table(df, file = out.genelist, sep = "\t", quote = FALSE) | 
| 129 } | 135 } | 
| 130 | 136 | 
| 131 | 137 | 
| 132 writecellassignments <- function(sc) { | 138 writecellassignments <- function(sc) { | 
| 133 dat <- sc@cluster$kpart | 139 dat <- sc@cluster$kpart | 
| 134 tab <- data.frame(row.names = NULL, | 140 tab <- data.frame( | 
| 135 cells = names(dat), | 141 row.names = NULL, | 
| 136 cluster.initial = dat, | 142 cells = names(dat), | 
| 137 cluster.final = sc@cpart, | 143 cluster.initial = dat, | 
| 138 is.outlier = names(dat) %in% sc@out$out) | 144 cluster.final = sc@cpart, | 
| 139 | 145 is.outlier = names(dat) %in% sc@out$out | 
| 140 write.table(tab, file = out.assignments, sep = "\t", | 146 ) | 
| 141 quote = FALSE, row.names = FALSE) | 147 | 
| 148 write.table(tab, | |
| 149 file = out.assignments, sep = "\t", | |
| 150 quote = FALSE, row.names = FALSE | |
| 151 ) | |
| 142 } | 152 } | 
| 143 | 153 | 
| 144 | 154 | 
| 145 pdf(out.pdf) | 155 pdf(out.pdf) | 
| 146 | 156 | 
| 147 if (use.filtnormconf) { | 157 if (use.filtnormconf) { | 
| 148 sc <- do.filter(sc) | 158 sc <- do.filter(sc) | 
| 149 message(paste(" - Source:: genes:", nrow(sc@expdata), | 159 message(paste( | 
| 150 ", cells:", ncol(sc@expdata))) | 160 " - Source:: genes:", nrow(sc@expdata), | 
| 151 message(paste(" - Filter:: genes:", nrow(as.matrix(getfdata(sc))), | 161 ", cells:", ncol(sc@expdata) | 
| 152 ", cells:", ncol(as.matrix(getfdata(sc))))) | 162 )) | 
| 153 message(paste(" :: ", | 163 message(paste( | 
| 154 sprintf("%.1f", 100 * nrow(as.matrix( | 164 " - Filter:: genes:", nrow(as.matrix(getfdata(sc))), | 
| 155 getfdata(sc))) / nrow(sc@expdata)), | 165 ", cells:", ncol(as.matrix(getfdata(sc))) | 
| 156 "% of genes remain,", | 166 )) | 
| 157 sprintf("%.1f", 100 * ncol(as.matrix( | 167 message(paste( | 
| 158 getfdata(sc))) / ncol(sc@expdata)), | 168 " :: ", | 
| 159 "% of cells remain")) | 169 sprintf("%.1f", 100 * nrow(as.matrix( | 
| 160 write.table(as.matrix(sc@ndata), file = out.table, col.names = NA, | 170 getfdata(sc) | 
| 161 row.names = TRUE, sep = "\t", quote = FALSE) | 171 )) / nrow(sc@expdata)), | 
| 172 "% of genes remain,", | |
| 173 sprintf("%.1f", 100 * ncol(as.matrix( | |
| 174 getfdata(sc) | |
| 175 )) / ncol(sc@expdata)), | |
| 176 "% of cells remain" | |
| 177 )) | |
| 178 write.table(as.matrix(sc@ndata), | |
| 179 file = out.table, col.names = NA, | |
| 180 row.names = TRUE, sep = "\t", quote = FALSE | |
| 181 ) | |
| 162 } | 182 } | 
| 163 | 183 | 
| 164 if (use.cluster) { | 184 if (use.cluster) { | 
| 165 par(mfrow = c(2, 2)) | 185 par(mfrow = c(2, 2)) | 
| 166 sc <- do.cluster(sc) | 186 sc <- do.cluster(sc) | 
