Mercurial > repos > iuc > dropletutils

comparison scripts/dropletutils.Rscript @ 0:4cd9f0008d9c 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/dropletutils/ commit e66ab3d4fc0c1a72523e8f93447cc07cdd6816b7
author iuc
date Tue, 04 Jun 2019 17:19:52 -0400
parents
children cfe1e6c28d95
comparison
equal deleted inserted replaced
-1:000000000000 0:4cd9f0008d9c
1 ## Load in data
2 args = commandArgs(trailingOnly = T)
3 if (length(args) != 1){
4 stop("Please provide the config file")
5 }
6
7 suppressWarnings(suppressPackageStartupMessages(require(DropletUtils)))
8 suppressWarnings(suppressPackageStartupMessages(require(Matrix)))
9
10 source(args[1])
11
12 ## Helper functions
13 setSparse <- function(obj){
14 return(as(obj, "dgCMatrix"))
15 }
16
17 writeTSV <- function(fileout, obj){
18 write.table(as.matrix(obj), file=fileout, col.names=NA, sep='\t', quote=FALSE)
19 }
20
21 writeOut <- function(counts, fileout, typeout){
22 if (typeout == "tsv"){
23 writeTSV(fileout, counts)
24 }
25 else if (typeout == "h5ad"){
26 write10xCounts(fileout, counts, type="HDF5", overwrite=TRUE)
27 }
28 else if (typeout == "directory"){
29 write10xCounts(fileout, counts, type="sparse", overwrite=TRUE)
30 }
31 }
32
33
34 read10xFiles <- function(filein, typein){
35 sce <- NULL
36 if (typein == "tsv"){
37 dat <- read.table(filein, header = TRUE, sep='\t', stringsAsFactors=FALSE, quote="", row.names=1)
38 sce <- SingleCellExperiment(assays = list(counts = as.matrix(dat)))
39 }
40 else if (typein == "h5ad"){
41 sce <- read10xCounts(filein, col.names=T, type="HDF5") # use barcodes.tsv as column names
42 }
43 else if (typein == "directory"){
44 sce <- read10xCounts(filein, col.names=T, type="sparse")
45 }
46 counts(sce) <- setSparse(counts(sce))
47 return(sce)
48 }
49
50
51 ## Methods
52
53
54 doEmptyDrops <- function(files, eparams, in.type="directory", out.type="h5ad", fdr_threshold = 0.01){
55 sce <- read10xFiles(files$infile, in.type)
56
57 eparams$... <- NULL ## hack
58 eparams$m = Matrix(counts(sce), sparse=TRUE)
59
60 e.out <- do.call(emptyDrops, c(eparams))
61
62 bar.names <- colnames(sce)
63 if (length(bar.names) != nrow(e.out)){
64 stop("Length of barcodes and output metrics don't match.")
65 }
66 e.out <- cbind(bar.names, e.out)
67 e.out$is.Cell <- e.out$FDR <= fdr_threshold
68 e.out$is.CellAndLimited <- e.out$is.Cell & e.out$Limited
69
70 # Write to table
71 writeTSV(files$table, e.out)
72
73 # Print to log
74 print(table(Limited=e.out$Limited, Significant=e.out$is.Cell))
75
76 # Write to Plot
77 png(files$plot)
78 plot(e.out$Total, -e.out$LogProb, col=ifelse(e.out$is.Cell, "red", "black"),
79 xlab="Total UMI count", ylab="-Log Probability")
80 dev.off()
81
82 # Filtered
83 called <- e.out$is.CellAndLimited
84 called[is.na(called)] <- FALSE # replace NA's with FALSE
85 sce.filtered <- sce[,called]
86
87 writeOut(counts(sce.filtered), files$out, out.type)
88 }
89
90
91 doDefaultDrops <- function(files, dparams, in.type="directory", out.type="h5ad"){
92 sce <- read10xFiles(files$infile, in.type)
93
94 dparams$m = as.matrix(counts(sce))
95 called <- do.call(defaultDrops, c(dparams))
96 print(table(called))
97
98 # Filtered
99 sce.filtered <- sce[,called]
100
101 writeOut(Matrix(counts(sce.filtered),sparse=TRUE), files$out, out.type)
102 }
103
104
105 doBarcodeRankings <- function(files, bparams, in.type="directory"){
106 sce <- read10xFiles(files$infile, in.type)
107
108 bparams$... <- NULL ## hack
109 bparams$m = as.matrix(counts(sce))
110
111 br.out <- do.call(barcodeRanks, c(bparams))
112
113 png(files$plot)
114 plot(br.out$rank, br.out$total, log="xy", xlab="(log) Rank", ylab="(log) Total Number of Barcodes")
115 o <- order(br.out$rank)
116 lines(br.out$rank[o], br.out$fitted[o], col="red")
117
118 abline(h=br.out$knee, col="dodgerblue", lty=2)
119 abline(h=br.out$inflection, col="forestgreen", lty=2)
120 legend("bottomleft", lty=2, col=c("dodgerblue", "forestgreen"), legend=c("knee", "inflection"))
121 dev.off()
122
123 print(paste("knee =", br.out$knee, ", inflection = ", br.out$inflection))
124 }
125
126 ## Main
127 set.seed(seed.val)
128
129 if (do.method == "barcodeRankings") {
130 doBarcodeRankings(files, bparams, in.type)
131
132 } else if (do.method == "defaultDrops") {
133 doDefaultDrops(files, dparams, in.type, out.type)
134
135 } else if (do.method == "emptyDrops") {
136 doEmptyDrops(files, eparams, in.type, out.type, empty.fdr_threshold)
137 }