Mercurial > repos > iuc > raceid_inspectclusters

comparison scripts/cluster.R @ 1:64c5c1bbbdbe 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 71e6b205841c83391ea8fc69e10eac03f212f4d6
author iuc
date Thu, 28 Feb 2019 13:00:52 -0500
parents 9fec5dd8fbb9
children 106718959281
comparison
equal deleted inserted replaced
0:9fec5dd8fbb9 1:64c5c1bbbdbe
1 #!/usr/bin/env R 1 #!/usr/bin/env R
2 VERSION = "0.2" 2 VERSION = "0.3"
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))
21 } 21 }
22 22
23 sc <- do.call(filterdata, c(sc, filt)) 23 sc <- do.call(filterdata, c(sc, filt))
24 24
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 <- log(colSums(as.matrix(sc@expdata))) 26 raw.lib <- log10(colSums(as.matrix(sc@expdata)))
27 raw.feat <- log(rowSums(as.matrix(sc@expdata))) 27 raw.feat <- log10(rowSums(as.matrix(sc@expdata)>0))
28 filt.lib <- log(colSums(getfdata(sc))) 28 filt.lib <- log10(colSums(getfdata(sc)))
29 filt.feat <- log(rowSums(getfdata(sc))) 29 filt.feat <- log10(rowSums(getfdata(sc)>0))
30 30
31 br <- 50 31 br <- 50
32 ## Determine limits on plots based on the unfiltered data 32 ## 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 33 ## (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) 34 ## so let them keep their own ranges)
45 45
46 ## feat.y_lim <- c(0,betterrange(max(tmp.feat$counts))) 46 ## feat.y_lim <- c(0,betterrange(max(tmp.feat$counts)))
47 ## feat.x_lim <- c(0,betterrange(max(tmp.feat$breaks))) 47 ## feat.x_lim <- c(0,betterrange(max(tmp.feat$breaks)))
48 48
49 par(mfrow=c(2,2)) 49 par(mfrow=c(2,2))
50 print(hist(raw.lib, breaks=br, main="ExpData Log(LibSize)")) # , xlim=lib.x_lim, ylim=lib.y_lim) 50 print(hist(raw.lib, breaks=br, main="RawData Log10(LibSize)")) # , xlim=lib.x_lim, ylim=lib.y_lim)
51 print(hist(raw.feat, breaks=br, main="ExpData Log(NumFeat)")) #, xlim=feat.x_lim, ylim=feat.y_lim) 51 print(hist(raw.feat, breaks=br, main="RawData Log10(NumFeat)")) #, xlim=feat.x_lim, ylim=feat.y_lim)
52 print(hist(filt.lib, breaks=br, main="FiltData Log(LibSize)")) # , xlim=lib.x_lim, ylim=lib.y_lim) 52 print(hist(filt.lib, breaks=br, main="FiltData Log10(LibSize)")) # , xlim=lib.x_lim, ylim=lib.y_lim)
53 print(hist(filt.feat, breaks=br, main="FiltData Log(NumFeat)")) # , xlim=feat.x_lim, ylim=feat.y_lim) 53 tmp <- hist(filt.feat, breaks=br, main="FiltData Log10(NumFeat)") # , xlim=feat.x_lim, ylim=feat.y_lim)
54 print(tmp) # required, for extracting midpoint
55 unq <- unique(filt.feat)
56 if (length(unq) == 1){
57 text(tmp$mids, table(filt.feat)[[1]] - 100, pos=1, paste(format(10^unq, scientific=T, digits=3),
58 " Features in all Cells", sep=""), cex=0.8)
59 }
60
54 61
55 if (filt.use.ccorrect){ 62 if (filt.use.ccorrect){
56 par(mfrow=c(2,2)) 63 par(mfrow=c(2,2))
57 sc <- do.call(CCcorrect, c(sc, filt.ccc)) 64 sc <- do.call(CCcorrect, c(sc, filt.ccc))
58 print(plotdimsat(sc, change=T)) 65 print(plotdimsat(sc, change=T))