corrected output folder for the filtered count matrix; information for the...

corrected output folder for the filtered count matrix; information for the selected PC are now additionally exported as a R object which can be visualized in the metrics_explorer app.

sc_workflow/seurat_analysis.R

@@ -107,9 +107,6 @@ count_matrix <- count_matrix[,which(!colnames(count_matrix) %in% outlier)]
 pheno_data <- pheno_data[which(!pheno_data$outlier),]
 if (ncol(count_matrix) != nrow(pheno_data)) {stop("Count matrix and design have different sample numbers")}
 
-# export filtered count matrix
-count_matrix %>% as.data.frame() %>% mutate(gene_id = rownames(.)) %>% push_left("gene_id") %>% write_tsv(paste0(str_replace(opts$count_matrix, ".[a-z]+$", ""), "_filtered.txt"))
-
 
 
 ## export gene information ---------------------------------------------------------------------------------------------
@@ -181,6 +178,10 @@ cells_after = ncol(seo@data)
 data.frame(remaining_cells = cells_after, removed_cells = cells_before - cells_after) %>% kable()
 
 
+# export filtered count matrix
+seo@raw.data %>% as.matrix() %>% as.data.frame() %>% mutate(gene_id = rownames(.)) %>% push_left("gene_id") %>% write_tsv(paste0(str_replace(basename(opts$count_matrix), ".[a-z]+$", ""), "_filtered.txt"))
+
+
 
 ## log-normalize count data and extrat highly variable genes -----------------------------------------------------------
 # The following CPP codes are the implementation for the LogNormalize.
@@ -360,12 +361,6 @@ data.frame(standard_deviation = sdev, variance = sdev^2) %>% mutate(pc = 1:n())
 ## create folter for cluster comparison
 dir.create("cluster_comp_data", showWarnings = FALSE)
 
-seo@scale.data %>%
-    as.data.frame() %>%
-    mutate(ensembl_gene_id = rownames(.)) %>%
-    left_join(gene_info, by = "ensembl_gene_id") %>%
-    write_tsv("cluster_comp_data/scaled_counts.txt")
-
 
 pc_num <- str_split(as.numeric(sop[which(sop$parameter == "calculate_pcs"),]$settings), "[+]") %>% unlist()
 n_cores <- ifelse(length(pc_num) >= 11, 10, length(pc_num))
@@ -418,6 +413,25 @@ all_cluster_infos <- mclapply(pc_num, function(x){
 write_tsv(all_cluster_infos, "all_cluster_infos.txt")
 seo@meta.data %>% left_join(all_cluster_infos, by = "cell_id") %>% write_tsv("cell_infos.txt")
 
+
+scaled_data <- seo@scale.data %>%
+    as.data.frame() %>%
+    mutate(ensembl_gene_id = rownames(.)) %>%
+    left_join(gene_info, by = "ensembl_gene_id")
+scaled_data %>%
+    write_tsv("cluster_comp_data/scaled_counts.txt")
+
+
+# combine information for set select_pc results and save it in the working directory
+selp <- paste0(sop[which(sop$parameter == "select_pc"), ]$settings, "pc")
+
+list(cluster = read_tsv(paste0("cluster_comp_data/cluster_info_", selp, ".txt")),
+    marker = read_tsv(paste0("cluster_comp_data/cluster_marker_genes_", selp, ".txt")),
+    tsne = read_tsv(paste0("cluster_comp_data/tsne_info_", selp, ".txt")),
+    scaled_data = scaled_data
+    ) %>% saveRDS(paste0(selp, "_infos.rds"))
+
+
 list(seo = seo, settings = sop) %>% saveRDS("seurat_object_data.rds")