Mercurial > repos > iuc > raceid_clustering

comparison scripts/cluster.R @ 10:49776718ae90 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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)