Mercurial > repos > iuc > raceid_clustering

comparison scripts/cluster.R @ 2:528a43b1cbcf 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 1d6e79ba92ce98c7c91f0c4076c9ca5e4e3f3a20
author iuc
date Thu, 28 Feb 2019 17:40:05 -0500
parents 89ee61bcc310
children d4646ce659b4
comparison
equal deleted inserted replaced
1:89ee61bcc310 2:528a43b1cbcf
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 <- log10(colSums(as.matrix(sc@expdata))) 26 raw.lib <- log10(colSums(as.matrix(sc@expdata)))
27 raw.feat <- log10(rowSums(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(rowSums(getfdata(sc)>0)) 29 filt.feat <- log10(colSums(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="RawData Log10(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="RawData Log10(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 Log10(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 tmp <- hist(filt.feat, breaks=br, main="FiltData Log10(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 54 print(tmp)
55 ## required, for extracting midpoint
55 unq <- unique(filt.feat) 56 unq <- unique(filt.feat)
56 if (length(unq) == 1){ 57 if (length(unq) == 1){
57 text(tmp$mids, table(filt.feat)[[1]] - 100, pos=1, paste(format(10^unq, scientific=T, digits=3), 58 abline(v=unq, col="red", lw=2)
58 " Features in all Cells", sep=""), cex=0.8) 59 text(tmp$mids, table(filt.feat)[[1]] - 100, pos=1, paste(10^unq, "\nFeatures\nin remaining\nCells", sep=""), cex=0.8)
59 } 60 }
60
61 61
62 if (filt.use.ccorrect){ 62 if (filt.use.ccorrect){
63 par(mfrow=c(2,2)) 63 par(mfrow=c(2,2))
64 sc <- do.call(CCcorrect, c(sc, filt.ccc)) 64 sc <- do.call(CCcorrect, c(sc, filt.ccc))
65 print(plotdimsat(sc, change=T)) 65 print(plotdimsat(sc, change=T))