IDEIS/Generate_final_mtx.R at main · Lab-of-Adaptive-Immunity/IDEIS · 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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
########################################################################
#
# License:
#
# IDEIS (c) by Lab of Adaptive Immunity from Institute of Molecular Genetics of the Czech Academy of Sciences
#
# IDEIS and all its parts are licensed under a
# Creative Commons Attribution 4.0 International License.
#
# You should have received a copy of the license along with this
# work. If not, see <https://creativecommons.org/licenses/by/4.0/>.
#
########################################################################
#
# Acknowledgements:
#
# This project was supported by the National Institute of Virology and
# Bacteriology (Programme EXCELES, LX22NPO5103 to Ondrej Stepanek) -
# funded by the European Union - Next Generation EU.
#
########################################################################
#
# Description:
#
# Processes final matrices generated by IDEIS to format usable directly
# with Seurat as .rds
#
########################################################################

library(Seurat)
library(data.table)
library(Matrix)
library(Matrix.utils)
library(magrittr)

args = commandArgs(trailingOnly=TRUE)
if (length(args) < 6 || length(args) > 6) {
  stop('Incorrect number of arguments detected. A path to input directory, to output directory and sequencing specs are needed.')
}

input.dir <- args[1]
output.path <- args[2]
sequencing.type <- args[3]
whitelist.path <- args[4]
data.set.input.path <- args[5]
data.set.output.path <- args[6]

all.barcodes <- ''
# no whitelist nor .rds data-set with Seurat Object supplied
if(whitelist.path == 'None' && data.set.input.path == 'None'){
  if(sequencing.type == '3-prime'){
    all.barcodes <- read.csv('cr_whitelists/3M-february-2018.txt', h=F)[,1]
  }else if(sequencing.type == '5-prime'){
    all.barcodes <- read.csv('cr_whitelists/737K-august-2016.txt', h=F)[,1]
  }
# whitelist was given
}else if(whitelist.path != 'None'){
  all.barcodes <- read.csv(whitelist.path, h=F)[,1]
# .rds data-set with Seurat Object was given
}else{
  all.barcodes <- read.csv(file.path(input.dir, 'whitelist.csv'), h=F)[,1]
}

# Now create final matrix - either whitelisted that can be directly imported to data or complete 10X-like that must be subset first.
loaded.matrix <- readMM(file.path(input.dir, 'alevin', 'alevin', 'quants_mat.mtx.gz'))
rownames(loaded.matrix) <- read.csv(file.path(input.dir, 'alevin', 'alevin', 'quants_mat_rows.txt'), h=F)[,1]
colnames(loaded.matrix) <- read.csv(file.path(input.dir, 'alevin', 'alevin', 'quants_mat_cols.txt'), h=F)[,1]
loaded.matrix <- t(loaded.matrix)

# fill in missing barcodes for compatibility with 10X
missing.barcodes <- all.barcodes[!(all.barcodes %in% colnames(loaded.matrix))]
filler.matrix <- sparseMatrix(i = integer(0), j = integer(0),  dims = c(nrow(loaded.matrix), length(missing.barcodes)))
colnames(filler.matrix) <- missing.barcodes

# create complete matrix
complete.matrix <- cbind(loaded.matrix, filler.matrix)
complete.matrix <- complete.matrix[,order(colnames(complete.matrix))]

# subset to only barcodes that appear in whitelist
# it is possible that alevin gets data for non-existent cell due to sequencing error; it will get removed here
complete.matrix <- complete.matrix[,colnames(complete.matrix) %in% all.barcodes]

# Add -1 to barcodes for compatibility with other 10X formatted data
colnames(complete.matrix) <- paste0(colnames(complete.matrix), '-1')

qu = Matrix::writeMM(complete.matrix, file.path(output.path, 'raw_feature_matrix', 'matrix.mtx'))
write(rownames(complete.matrix), file.path(output.path, 'raw_feature_matrix', 'genes.tsv'))
write(colnames(complete.matrix), file.path(output.path, 'raw_feature_matrix', 'barcodes.tsv'))

saveRDS(complete.matrix, file.path(output.path, 'raw_matrix.rds'))

# if SeuratObject was given, then add complete matrix to it as CD45isoforms assay
if(data.set.input.path != 'None'){
  cat('The .rds file of SeuratObject provided, creating Assay now.\n')
  data.set <- readRDS(data.set.input.path)
  data.set[['CD45isoforms']] <- CreateAssayObject(complete.matrix)
  saveRDS(data.set, data.set.output.path)
}