Mercurial > repos > iuc > dropletutils

comparison scripts/dropletutils.Rscript @ 1:94ee8160ace3 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 40d4a92c62d75fe931baba8657cde006a26d84cf"
author iuc
date Mon, 26 Aug 2019 05:06:15 -0400
parents 1335bf5bcd0e
children 58105e5efa8d
comparison
equal deleted inserted replaced
0:1335bf5bcd0e 1:94ee8160ace3
4 stop("Please provide the config file") 4 stop("Please provide the config file")
5 } 5 }
6 6
7 suppressWarnings(suppressPackageStartupMessages(require(DropletUtils))) 7 suppressWarnings(suppressPackageStartupMessages(require(DropletUtils)))
8 suppressWarnings(suppressPackageStartupMessages(require(Matrix))) 8 suppressWarnings(suppressPackageStartupMessages(require(Matrix)))
9 suppressWarnings(suppressPackageStartupMessages(require(scater)))
9 10
10 source(args[1]) 11 source(args[1])
11 12
12 ## Helper functions 13 ## Helper functions
13 setSparse <- function(obj){ 14 setSparse <- function(obj){
32 33
33 34
34 read10xFiles <- function(filein, typein){ 35 read10xFiles <- function(filein, typein){
35 sce <- NULL 36 sce <- NULL
36 if (typein == "tsv"){ 37 if (typein == "tsv"){
37 dat <- read.table(filein, header = TRUE, sep='\t', stringsAsFactors=FALSE, quote="", row.names=1) 38 ## Exploding memory problems occured here
38 sce <- SingleCellExperiment(assays = list(counts = as.matrix(dat))) 39 ## - solution is to use the readSparseCounts function from scater
40 sce <- SingleCellExperiment(assays = list(counts = readSparseCounts(filein)))
39 } 41 }
40 else if (typein == "h5ad"){ 42 else if (typein == "h5ad"){
41 sce <- read10xCounts(filein, col.names=T, type="HDF5") # use barcodes.tsv as column names 43 sce <- read10xCounts(filein, col.names=T, type="HDF5") # use barcodes.tsv as column names
42 } 44 }
43 else if (typein == "directory"){ 45 else if (typein == "directory"){
89 91
90 92
91 doDefaultDrops <- function(files, dparams, in.type="directory", out.type="h5ad"){ 93 doDefaultDrops <- function(files, dparams, in.type="directory", out.type="h5ad"){
92 sce <- read10xFiles(files$infile, in.type) 94 sce <- read10xFiles(files$infile, in.type)
93 95
94 dparams$m = as.matrix(counts(sce)) 96 dparams$m = counts(sce)
95 called <- do.call(defaultDrops, c(dparams)) 97 called <- do.call(defaultDrops, c(dparams))
96 print(table(called)) 98 print(table(called))
97 99
98 # Filtered 100 # Filtered
99 sce.filtered <- sce[,called] 101 sce.filtered <- sce[,called]
104 106
105 doBarcodeRankings <- function(files, bparams, in.type="directory"){ 107 doBarcodeRankings <- function(files, bparams, in.type="directory"){
106 sce <- read10xFiles(files$infile, in.type) 108 sce <- read10xFiles(files$infile, in.type)
107 109
108 bparams$... <- NULL ## hack 110 bparams$... <- NULL ## hack
109 bparams$m = as.matrix(counts(sce)) 111 bparams$m = counts(sce)
110 112
111 br.out <- do.call(barcodeRanks, c(bparams)) 113 br.out <- do.call(barcodeRanks, c(bparams))
112 114
113 png(files$plot) 115 png(files$plot)
114 plot(br.out$rank, br.out$total, log="xy", xlab="(log) Rank", ylab="(log) Total Number of Barcodes") 116 plot(br.out$rank, br.out$total, log="xy", xlab="(log) Rank", ylab="(log) Total Number of Barcodes")