Mercurial > repos > iuc > raceid_filtnormconf

comparison scripts/cluster.R @ 3:d55e29ac02e3 draft

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 d94b3b8a4c7cf8c604279eb1eea24d32b3868922
author iuc
date Mon, 15 Apr 2019 17:55:17 -0400
parents 56a093c2a3f9
children 5d5b14dbd092
comparison
equal deleted inserted replaced
2:56a093c2a3f9 3:d55e29ac02e3
1 #!/usr/bin/env R 1 #!/usr/bin/env R
2 VERSION = "0.3" 2 VERSION = "0.4"
3 3
4 args = commandArgs(trailingOnly = T) 4 args = commandArgs(trailingOnly = T)
5 5
6 if (length(args) != 1){ 6 if (length(args) != 1){
7 message(paste("VERSION:", VERSION)) 7 message(paste("VERSION:", VERSION))
25 ## Get histogram metrics for library size and number of features 25 ## Get histogram metrics for library size and number of features
26 raw.lib <- log10(colSums(as.matrix(sc@expdata))) 26 raw.lib <- log10(colSums(as.matrix(sc@expdata)))
27 raw.feat <- log10(colSums(as.matrix(sc@expdata)>0)) 27 raw.feat <- log10(colSums(as.matrix(sc@expdata)>0))
28 filt.lib <- log10(colSums(getfdata(sc))) 28 filt.lib <- log10(colSums(getfdata(sc)))
29 filt.feat <- log10(colSums(getfdata(sc)>0)) 29 filt.feat <- log10(colSums(getfdata(sc)>0))
30
31 if (filt.geqone){
32 filt.feat <- log10(colSums(getfdata(sc)>=1))
33 }
30 34
31 br <- 50 35 br <- 50
32 ## Determine limits on plots based on the unfiltered data 36 ## Determine limits on plots based on the unfiltered data
33 ## (doesn't work, R rejects limits and norm data is too different to compare to exp data 37 ## (doesn't work, R rejects limits and norm data is too different to compare to exp data
34 ## so let them keep their own ranges) 38 ## so let them keep their own ranges)
122 dg.goi.table <- head(dg.goi, genelist.tablelim) 126 dg.goi.table <- head(dg.goi, genelist.tablelim)
123 df <<- rbind(df, cbind(n, dg.goi.table)) 127 df <<- rbind(df, cbind(n, dg.goi.table))
124 128
125 goi <- head(rownames(dg.goi.table), genelist.plotlim) 129 goi <- head(rownames(dg.goi.table), genelist.plotlim)
126 print(plotmarkergenes(sc, goi)) 130 print(plotmarkergenes(sc, goi))
127 print(do.call(mtext, c(paste(" Cluster ",n), test))) ## spacing is a hack 131 buffer <- paste(rep("", 36), collapse=" ")
132 print(do.call(mtext, c(paste(buffer, "Cluster ",n), test))) ## spacing is a hack
128 test$line=-1 133 test$line=-1
129 print(do.call(mtext, c(paste(" Sig. Genes"), test))) ## spacing is a hack 134 print(do.call(mtext, c(paste(buffer, "Sig. Genes"), test))) ## spacing is a hack
130 test$line=-2 135 test$line=-2
131 print(do.call(mtext, c(paste(" (fc > ", genelist.foldchange,")"), test))) ## spacing is a hack 136 print(do.call(mtext, c(paste(buffer, "(fc > ", genelist.foldchange,")"), test))) ## spacing is a hack
132 137
133 }) 138 })
134 write.table(df, file=out.genelist, sep="\t", quote=F) 139 write.table(df, file=out.genelist, sep="\t", quote=F)
135 } 140 }
136 141
137 pdf(out.pdf) 142 pdf(out.pdf)
138 143
139 if (use.filtnormconf){ 144 if (use.filtnormconf){
140 sc <- do.filter(sc) 145 sc <- do.filter(sc)
141 message(paste(" - Source:: genes:",nrow(sc@expdata),", cells:",ncol(sc@expdata))) 146 message(paste(" - Source:: genes:",nrow(sc@expdata),", cells:",ncol(sc@expdata)))
142 message(paste(" - Filter:: genes:",nrow(sc@ndata),", cells:",ncol(sc@ndata))) 147 message(paste(" - Filter:: genes:",nrow(getfdata(sc)),", cells:",ncol(getfdata(sc))))
143 message(paste(" :: ", 148 message(paste(" :: ",
144 sprintf("%.1f", 100 * nrow(sc@ndata)/nrow(sc@expdata)), "% of genes remain,", 149 sprintf("%.1f", 100 * nrow(getfdata(sc))/nrow(sc@expdata)), "% of genes remain,",
145 sprintf("%.1f", 100 * ncol(sc@ndata)/ncol(sc@expdata)), "% of cells remain")) 150 sprintf("%.1f", 100 * ncol(getfdata(sc))/ncol(sc@expdata)), "% of cells remain"))
146 } 151 }
147 152
148 if (use.cluster){ 153 if (use.cluster){
149 par(mfrow=c(2,2)) 154 par(mfrow=c(2,2))
150 sc <- do.cluster(sc) 155 sc <- do.cluster(sc)