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
| library(Seurat) library(dplyr) library(patchwork)
root_path = "~/zlliu/R_data/hBLA"
output_path = "~/zlliu/R_output/21.07.15.10x.hM1.se/"
setwd(output_path) print(getwd())
tp_samples <- list.files(file.path(root_path, "10x")) tp_dir <- file.path(root_path, "10x", tp_samples) names(tp_dir) <- tp_samples counts <- Read10X(data.dir = tp_dir) scRNA <- CreateSeuratObject(counts, project = 'hBLA', min.cells = 3, min.features = 200) rm(counts) rm(tp_dir) gc()
scRNA[["percent.mt"]] <- PercentageFeatureSet(scRNA, pattern = "^MT-")
x10s <- SplitObject(scRNA, split.by = 'ident') rm(scRNA) gc()
x10s$BA213 <- subset(x10s$BA213, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 10) x10s$BA04 <- subset(x10s$BA04, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mt < 10) x10s$BA09 <- subset(x10s$BA09, subset = nFeature_RNA > 200 & nFeature_RNA < 3500 & percent.mt < 10) x10s$BA21 <- subset(x10s$BA21, subset = nFeature_RNA > 200 & nFeature_RNA < 4500 & percent.mt < 10) x10s$BA22 <- subset(x10s$BA22, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mt < 10) x10s$BA39 <- subset(x10s$BA39, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 10) x10s$BA40 <- subset(x10s$BA40, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 10) x10s$BA44 <- subset(x10s$BA44, subset = nFeature_RNA > 200 & nFeature_RNA < 4000 & percent.mt < 10) x10s$BA45 <- subset(x10s$BA45, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 10)
library(SingleR) library(SummarizedExperiment) library(scater) library(BiocParallel)
if (!("hM1.se" %in% ls())){ hM1.se <- readRDS("/home/rqzhang/zlliu/R_data/human_M1_10x/hM1.se.rds") }
f_get_pred <- function(scRNA, lc_i){ if(file.exists(sprintf( "%02d_hM1_subclass.txt", lc_i))){ result_main_hpca <- read.table(sprintf( "%02d_hM1_subclass.txt", lc_i), stringsAsFactors = F) return (result_main_hpca) } test.count=as.data.frame(scRNA[["RNA"]]@counts) common_hpca <- intersect(rownames(test.count), rownames(hM1.se)) lc_hM1.se <- hM1.se[common_hpca,] test.count_forhpca <- test.count[common_hpca,] gc() test.count_forhpca.se <- SummarizedExperiment(assays=list(counts=test.count_forhpca)) test.count_forhpca.se <- logNormCounts(test.count_forhpca.se) pred.main.hpca <- SingleR(test = test.count_forhpca.se, ref = lc_hM1.se, labels = hM1.se$meta$subclass_label, BPPARAM=MulticoreParam(32)) gc() result_main_hpca <- as.data.frame(pred.main.hpca$labels) result_main_hpca$CB <- rownames(pred.main.hpca) colnames(result_main_hpca) <- c('hM1_subclass', 'CB') write.table(result_main_hpca, sprintf( "%02d_hM1_subclass.txt", lc_i)) result_main_hpca }
tp_sc <- x10s
gc() tp_sc_st <- 1 tp_sc_len <- length(tp_sc) for (lc_i in tp_sc_st:tp_sc_len) { print(lc_i) f_get_pred(tp_sc[[lc_i]], lc_i) }
|