forked from digitalcytometry/ecotyper
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathEcoTyper_recovery_visium.R
180 lines (152 loc) · 7.11 KB
/
EcoTyper_recovery_visium.R
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
suppressPackageStartupMessages({
library(config)
library(argparse)
source("pipeline/lib/config.R")
source("pipeline/lib/misc.R")
source("pipeline/lib/multithreading.R")
})
parser <- ArgumentParser(add_help = F)
arguments = parser$add_argument_group('Arguments')
arguments$add_argument("-c", "--config", type = "character", metavar="<PATH>",
help="Path to the config files [required].")
arguments$add_argument("-h", "--help", action='store_true', help="Print help message.")
args <- parser$parse_args()
if(args$h || is.null(args$config))
{
parser$print_help()
quit()
}
config_file = abspath(args$config)
config <- config::get(file = config_file)
discovery = config$Input$"Discovery dataset name"
recovery = config$Input$"Recovery dataset name"
input_path = config$Input$"Input Visium directory"
fractions = config$Input$"Recovery cell type fractions"
coo = config$Input$"Background cell type"
if(is.null(coo))
{
coo = config$Input$"Malignant cell of origin"
}
coo = config$Input$"Malignant cell of origin"
CSx_username = config$Input$"CIBERSORTx username"
CSx_token = config$Input$"CIBERSORTx token"
n_threads = config$"Pipeline settings"$"Number of threads"
CSx_singularity_path_fractions = config$"Pipeline settings"$"CIBERSORTx fractions Singularity path"
final_output = config$"Output"$"Output folder"
suppressWarnings({
input_path = abspath(input_path)
final_output = abspath(file.path(final_output, recovery))
})
suppressWarnings({
fractions_path = abspath(fractions)
})
if(!is.na(CSx_singularity_path_fractions) && !is.null(CSx_singularity_path_fractions))
{
if(!file.exists(CSx_singularity_path_fractions))
{
stop(paste0("CIBERSORTx fractions Singularity path provided does not exist:", CSx_singularity_path_fractions))
}
}else{
CSx_singularity_path_fractions = "NULL"
}
if(!discovery %in% c("Carcinoma", "Lymphoma"))
{
discovery_config_file = file.path("EcoTyper", discovery, "config_used.yml")
if(!file.exists(discovery_config_file))
{
stop("Error: Cannot read the configuration file used for the discovery of cell states and ecotypes. It should be in the following path: '", config_file, "'. Please make sure that the '--discovery (-d)' argument provided is correct!")
}
config <- config::get(file = discovery_config_file)
discovery_fractions = config$"Pipeline settings"$"Cell type fractions"
if(!is.null(discovery_fractions) && discovery_fractions %in% c("Carcinoma_Fractions", "Lymphoma_Fractions"))
{
fractions = discovery_fractions
}else{
if(!is.null(config$"Pipeline settings"$"Filter non cell type specific genes"))
{
if(config$"Pipeline settings"$"Filter non cell type specific genes")
{
fractions = "Cell_type_specific_genes"
}else{
fractions = "All_genes"
}
}else{
fractions = "Custom"
}
}
}else{
fractions = paste0(discovery, "_Fractions")
}
#Starting EcoTyper
setwd("pipeline")
start = Sys.time()
cat("\nLoading visium data...\n")
PushToJobQueue(paste("Rscript spatial_load_visium_data.R", recovery, input_path))
RunJobQueue()
if(fractions %in% c("Carcinoma_Fractions", "Lymphoma_Fractions") && !file.exists(fractions_path))
{
cat("\nRunning CIBERSORTxFractions on the visium dataset...\n")
print(paste("Rscript csx_fractions.R", "visium", recovery, "LM22", "B_mode", CSx_username, CSx_token, CSx_singularity_path_fractions, TRUE))
PushToJobQueue(paste("Rscript csx_fractions.R", "visium", recovery, "LM22", "B_mode", CSx_username, CSx_token, CSx_singularity_path_fractions, TRUE))
RunJobQueue()
PushToJobQueue(paste("Rscript csx_fractions.R", "visium", recovery, "TR4", "B_mode", CSx_username, CSx_token, CSx_singularity_path_fractions, TRUE))
RunJobQueue()
PushToJobQueue(paste("Rscript csx_fractions_two_tiered.R", "visium", recovery, "TR4", "B_mode", "LM22", "B_mode", fractions))
RunJobQueue()
coo = "Epithelial.cells"
if(fractions %in% "Lymphoma_Fractions")
{
coo = "B.cells"
}
}else{
cat("\nLoading user-provided cell type fractions...\n")
dir.create(file.path("../CIBERSORTx/fractions/visium", recovery, fractions), recursive = T, showWarnings = F)
PushToJobQueue(paste("cp -f", fractions_path, file.path("../CIBERSORTx/fractions/visium", recovery, fractions, "CIBERSORTx_Adjusted.txt")))
RunJobQueue()
}
cat("\nRunning cell state recovery on the visium dataset...\n")
key = read.delim(file.path("../EcoTyper", discovery, fractions, "Analysis", "rank_selection", "rank_data.txt"))
for(cell_type in key[,1])
{
n_states = key[key[,1] == cell_type, 2]
PushToJobQueue((paste("Rscript state_recovery_visium.R", discovery, fractions, cell_type, n_states, recovery, "FALSE")))
}
RunJobQueue()
cat("\nCalculating cell state abundances...\n")
print(paste("Rscript spatial_states.R", discovery, recovery, fractions, coo))
PushToJobQueue((paste("Rscript spatial_states.R", discovery, recovery, fractions, coo)))
RunJobQueue()
cat("\nCalculating ecotype abundances...\n")
PushToJobQueue((paste("Rscript spatial_ecotypes.R", discovery, recovery, fractions, coo)))
RunJobQueue()
cat("\nPlotting cell state heatmaps...\n")
key = read.delim(file.path("../EcoTyper", discovery, fractions, "Analysis", "rank_selection", "rank_data.txt"))
for(cell_type in key[,1])
{
n_states = key[key[,1] == cell_type, 2]
PushToJobQueue((paste("Rscript spatial_plot_states.R", discovery, recovery, fractions, coo, cell_type)))
}
RunJobQueue()
PushToJobQueue((paste("Rscript spatial_plot_ecotypes.R", discovery, recovery, fractions, coo)))
RunJobQueue()
cat("\nCopying EcoTyper results to the output folder!\n")
if(file.exists(final_output) && length(list.files(final_output)) > 0)
{
old_results_folder = paste0(final_output, format(Sys.time(), " %a %b %d %X %Y"))
dir.create(old_results_folder, recursive = T, showWarnings = F)
warning(paste0("The output folder contains files from a previous run. Moving those files to: '", old_results_folder, "'"))
system(paste0("mv -f ", final_output, "/* '", old_results_folder, "'"))
}
dir.create(final_output, recursive = T, showWarnings = F)
key = read.delim(file.path("../EcoTyper", discovery, fractions, "Analysis", "rank_selection", "rank_data.txt"))
for(cell_type in key[,1])
{
system(paste("cp -f", file.path("../EcoTyper", discovery, fractions, "Cell_States", "recovery", recovery, paste0(cell_type, "_spatial_heatmaps.pdf")), final_output))
system(paste("cp -f", file.path("../EcoTyper", discovery, fractions, "Cell_States", "recovery", recovery, paste0(cell_type, "_spatial_heatmaps.png")), final_output))
}
system(paste("cp -f", file.path("../EcoTyper", discovery, fractions, "Cell_States", "recovery", recovery, "state_abundances.txt"), final_output))
system(paste("cp -f", file.path("../EcoTyper", discovery, fractions, "Cell_States", "recovery", recovery, "ecotype_abundances.txt"), final_output))
system(paste("cp -f", file.path("../EcoTyper", discovery, fractions, "Cell_States", "recovery", recovery, "Ecotype_spatial_heatmaps.pdf"), final_output))
system(paste("cp -f", file.path("../EcoTyper", discovery, fractions, "Cell_States", "recovery", recovery, "Ecotype_spatial_heatmaps.png"), final_output))
end = Sys.time()
cat(paste0("\nEcoTyper finished succesfully! Please find the results in: '", final_output, "'.\nRun time: ", format(end - start, digits = 1), "\n"))