read-normalization/07_compute_qumis_census_sim2.R at main · berenslab/read-normalization · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
suppressPackageStartupMessages(library(SingleCellExperiment))
suppressPackageStartupMessages(library(EDASeq))
suppressPackageStartupMessages(library(Matrix))
library(scRNAseq)
library(quminorm)
library(monocle)
library(scran)

print('reading data')
mtx_tpm = readMM('data//tasic//simulations//sim_only_marker_genes_DE_tpm_transposed.mtx')

print('make sce object')
sce = SingleCellExperiment(list(tpm=mtx_tpm))
cellnames = read.csv('data//tasic//simulations//sim_only_marker_genes_DE_tpm_cells.csv')$cells
genenames = read.csv('data//tasic//simulations//sim_only_marker_genes_DE_tpm_genes.csv')$gene_name
rownames(sce) = genenames
colnames(sce) = cellnames

print('compute qUMIs')
system.time(sce<-quminorm(sce,assayName="tpm",shape=2,mc.cores=16))
qumi=assay(sce,"qumi_poilog_2")

print('write out qUMIs')
writeMM(qumi,'data//tasic//simulations//sim_only_marker_genes_DE_qumi_transposed.mtx')


#monocle/census in batches
print('preparing batches')
batchsize_max = 1000
batch_id <- ceiling(seq(nrow(colData(sce)))/batchsize_max)
batches <- lapply(unique(batch_id), function(x) sce[, batch_id == x])

print('compute census counts')
print('number of batches:')
print(length(unique(batch_id)))

system.time(
for (i in 1:length(batches)) {
    print('working on batch')
    print(i)
    sce_batch = batches[[i]]
    cds_batch = convertTo(sce_batch,"monocle",assay.type='tpm')
    census_batch = relative2abs(cds_batch,cores=16,verbose=TRUE,return_all=TRUE,method="num_genes")
    census_counts_batch = census_batch$norm_cds
    if (i==1){
        results = census_counts_batch
    }
    else {
        results = cbind(results,census_counts_batch)
    }
}
)

print('write out census counts')
writeMM(results,'data//tasic//simulations//sim_only_marker_genes_DE_census_transposed.mtx')