pacman::p_load(Seurat, tidyverse, qs2, scDblFinder, BiocParallel)
knitr::opts_knit$set(root.dir = "/home/han/Desktop/guoqiang/stroke_bm/stroke_bm_Rproj")
setwd("/home/han/Desktop/guoqiang/stroke_bm/stroke_bm_Rproj")
source("codes/functions.R")Quality Control
读入 11 例颅骨骨髓 10x 矩阵并做质控：去掉空液滴、双细胞和线粒体过高的细胞。对照按性别、年龄与 ACI 匹配。取样批次和 5′/3′ 文库见首页。
1. Load Data
读入 11 例 10x 矩阵，写上分组、性别、年龄。
all_data <- list()
all_data[["ctrl_1"]] <-
Read10X(data.dir = "data/cr_output/220401M62/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "ctrl_1")
all_data[["ctrl_1"]]@meta.data <- all_data[["ctrl_1"]]@meta.data %>%
mutate(sample = "CTRL-1", group = "CTRL", gender = "Male", age = 62)
all_data[["ctrl_2"]] <-
Read10X(data.dir = "data/cr_output/CYM_SkullBM_240322/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "ctrl_2")
all_data[["ctrl_2"]]@meta.data <- all_data[["ctrl_2"]]@meta.data %>%
mutate(sample = "CTRL-2", group = "CTRL", gender = "Female", age = 59)
all_data[["ctrl_3"]] <-
Read10X(data.dir = "data/cr_output/LWN_SkullBM_240407/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "ctrl_3")
all_data[["ctrl_3"]]@meta.data <- all_data[["ctrl_3"]]@meta.data %>%
mutate(sample = "CTRL-3", group = "CTRL", gender = "Male", age = 77)
all_data[["ctrl_4"]] <-
Read10X(data.dir = "data/cr_output/LXY_SkullBM_240412/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "ctrl_4")
all_data[["ctrl_4"]]@meta.data <- all_data[["ctrl_4"]]@meta.data %>%
mutate(sample = "CTRL-4", group = "CTRL", gender = "Female", age = 60)
all_data[["ctrl_5"]] <-
Read10X(data.dir = "data/cr_output/STX_SkullBM_240229/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "ctrl_5")
all_data[["ctrl_5"]]@meta.data <- all_data[["ctrl_5"]]@meta.data %>%
mutate(sample = "CTRL-5", group = "CTRL", gender = "Male", age = 30)
all_data[["ctrl_6"]] <-
Read10X(data.dir = "data/cr_output/ZM_SkullBM_240620/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "ctrl_6")
all_data[["ctrl_6"]]@meta.data <- all_data[["ctrl_6"]]@meta.data %>%
mutate(sample = "CTRL-6", group = "CTRL", gender = "Female", age = 30)
all_data[["aci_1"]] <-
Read10X(data.dir = "data/cr_output/CBM_3/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "aci_1")
all_data[["aci_1"]]@meta.data <- all_data[["aci_1"]]@meta.data %>%
mutate(sample = "ACI-1", group = "ACI", gender = "Male", age = 54)
all_data[["aci_2"]] <-
Read10X(data.dir = "data/cr_output/CBM_4/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "aci_2")
all_data[["aci_2"]]@meta.data <- all_data[["aci_2"]]@meta.data %>%
mutate(sample = "ACI-2", group = "ACI", gender = "Male", age = 75)
all_data[["aci_3"]] <-
Read10X(data.dir = "data/cr_output/CBM_6/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "aci_3")
all_data[["aci_3"]]@meta.data <- all_data[["aci_3"]]@meta.data %>%
mutate(sample = "ACI-3", group = "ACI", gender = "Male", age = 50)
all_data[["aci_4"]] <-
Read10X(data.dir = "data/cr_output/CBM_5/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "aci_4")
all_data[["aci_4"]]@meta.data <- all_data[["aci_4"]]@meta.data %>%
mutate(sample = "ACI-4", group = "ACI", gender = "Female", age = 65)
all_data[["aci_5"]] <-
Read10X(data.dir = "data/cr_output/CBM_2/filtered_feature_bc_matrix") %>%
CreateSeuratObject(project = "aci_5")
all_data[["aci_5"]]@meta.data <- all_data[["aci_5"]]@meta.data %>%
mutate(sample = "ACI-5", group = "ACI", gender = "Male", age = 35)
qs_save(all_data, file = "data/rdata/preprocess/all_data.qs2", nthreads = 10)2. Merge Layers
合成一个对象，方便后面统一质控。
all_data <- qs_read("data/rdata/preprocess/all_data.qs2", nthreads = 10)all_data_qc <-
merge(all_data[[1]], all_data[-1], add.cell.ids = names(all_data))
Layers(all_data_qc) [1] "counts.ctrl_1" "counts.ctrl_2" "counts.ctrl_3" "counts.ctrl_4"
[5] "counts.ctrl_5" "counts.ctrl_6" "counts.aci_1" "counts.aci_2"
[9] "counts.aci_3" "counts.aci_4" "counts.aci_5" all_data_qc <- JoinLayers(all_data_qc)
rm(all_data)
Layers(all_data_qc)[1] "counts"dim(all_data_qc)[1] 38606 79488dim(all_data_qc@assays$RNA$counts)[1] 38606 79488as.data.frame(all_data_qc@assays$RNA$counts[1:10, 1:2]) ctrl_1_AAACCTGAGCCAGTAG-1 ctrl_1_AAACCTGAGTCCGGTC-1
DDX11L2 0 0
MIR1302-2HG 0 0
FAM138A 0 0
ENSG00000290826 0 0
OR4F5 0 0
ENSG00000238009 0 0
ENSG00000239945 0 0
ENSG00000239906 0 0
ENSG00000241860 0 0
ENSG00000241599 0 0head(all_data_qc@meta.data, 10) orig.ident nCount_RNA nFeature_RNA sample group
ctrl_1_AAACCTGAGCCAGTAG-1 ctrl_1 60591 7154 CTRL-1 CTRL
ctrl_1_AAACCTGAGTCCGGTC-1 ctrl_1 519 436 CTRL-1 CTRL
ctrl_1_AAACCTGCAAGCGCTC-1 ctrl_1 18100 4470 CTRL-1 CTRL
ctrl_1_AAACCTGCACACATGT-1 ctrl_1 2025 1244 CTRL-1 CTRL
ctrl_1_AAACCTGCAGACGCAA-1 ctrl_1 8300 3061 CTRL-1 CTRL
ctrl_1_AAACCTGCAGGTGCCT-1 ctrl_1 556 451 CTRL-1 CTRL
ctrl_1_AAACCTGTCAAGGTAA-1 ctrl_1 769 615 CTRL-1 CTRL
ctrl_1_AAACCTGTCACAATGC-1 ctrl_1 669 387 CTRL-1 CTRL
ctrl_1_AAACCTGTCGCGATCG-1 ctrl_1 580 308 CTRL-1 CTRL
ctrl_1_AAACCTGTCGGAGGTA-1 ctrl_1 31591 6887 CTRL-1 CTRL
gender age
ctrl_1_AAACCTGAGCCAGTAG-1 Male 62
ctrl_1_AAACCTGAGTCCGGTC-1 Male 62
ctrl_1_AAACCTGCAAGCGCTC-1 Male 62
ctrl_1_AAACCTGCACACATGT-1 Male 62
ctrl_1_AAACCTGCAGACGCAA-1 Male 62
ctrl_1_AAACCTGCAGGTGCCT-1 Male 62
ctrl_1_AAACCTGTCAAGGTAA-1 Male 62
ctrl_1_AAACCTGTCACAATGC-1 Male 62
ctrl_1_AAACCTGTCGCGATCG-1 Male 62
ctrl_1_AAACCTGTCGGAGGTA-1 Male 62all_data_qc@meta.data |>
dplyr::select(orig.ident, sample) |>
distinct() orig.ident sample
ctrl_1_AAACCTGAGCCAGTAG-1 ctrl_1 CTRL-1
ctrl_2_AAACCTGAGGGAAACA-1 ctrl_2 CTRL-2
ctrl_3_AAACCTGAGAAACGAG-1 ctrl_3 CTRL-3
ctrl_4_AAACCTGAGACACTAA-1 ctrl_4 CTRL-4
ctrl_5_AAACCTGCATCGGACC-1 ctrl_5 CTRL-5
ctrl_6_AAACCTGAGAAACGAG-1 ctrl_6 CTRL-6
aci_1_AAACCAAAGCAATCTG-1 aci_1 ACI-1
aci_2_AAACCAAAGGTCATCA-1 aci_2 ACI-2
aci_3_AAACCAAAGCAAGGCA-1 aci_3 ACI-3
aci_4_AAACCAAAGATACTGT-1 aci_4 ACI-4
aci_5_AAACCAAAGCAATCGC-1 aci_5 ACI-53. Calcuate QC
算每个细胞的基因数、UMI、线粒体和核糖体比例。线粒体很高常见于濒死或破膜细胞。
grep("^MT-", rownames(all_data_qc), value = TRUE) [1] "MT-ND1" "MT-ND2" "MT-CO1" "MT-CO2" "MT-ATP8" "MT-ATP6" "MT-CO3"
[8] "MT-ND3" "MT-ND4L" "MT-ND4" "MT-ND5" "MT-ND6" "MT-CYB" grep("^RP[SL]", rownames(all_data_qc), value = TRUE) [1] "RPL22" "RPL11" "RPS6KA1" "RPS8" "RPL5"
[6] "RPS27" "RPS6KC1" "RPS7" "RPS27A" "RPL31"
[11] "RPL37A-DT" "RPL37A" "RPL32" "RPL15" "RPSA"
[16] "RPL14" "RPL29" "RPL24" "RPL22L1" "RPL39L"
[21] "RPL35A" "RPL9" "RPL34-DT" "RPL34" "RPS3A"
[26] "RPL37" "RPS23" "RPS14" "RPL26L1-AS1" "RPL26L1"
[31] "RPS18" "RPS10" "RPL10A" "RPL7L1" "RPS12"
[36] "RPS6KA2" "RPS6KA2-IT1" "RPS6KA2-AS1" "RPS20" "RPL7"
[41] "RPL30" "RPL30-AS1" "RPL8" "RPS6" "RPL35"
[46] "RPL12" "RPL7A" "RPS24" "RPLP2" "RPL27A"
[51] "RPS13" "RPS6KA4" "RPS6KB2" "RPS6KB2-AS1" "RPS3"
[56] "RPS25" "RPS26" "RPL41" "RPL6" "RPLP0"
[61] "RPL21" "RPL10L" "RPS29" "RPL36AL" "RPS6KL1"
[66] "RPS6KA5" "RPS27L" "RPL4" "RPLP1" "RPS17"
[71] "RPL3L" "RPS2" "RPS15A" "RPL13" "RPL26"
[76] "RPL23A" "RPL23" "RPL19" "RPL27" "RPS6KB1"
[81] "RPL38" "RPL17" "RPS15" "RPL36" "RPS28"
[86] "RPL18A" "RPSA2" "RPS16" "RPS19" "RPL18"
[91] "RPL13A" "RPS11" "RPS9" "RPL28" "RPS5"
[96] "RPS21" "RPL3" "RPS19BP1" "RPL23AP7" "RPS6KA3"
[101] "RPS4X" "RPS6KA6" "RPL36A" "RPL39" "RPL10"
[106] "RPS4Y1" "RPS4Y2" grep("^HB[ABDEGMQZ]\\d*$", rownames(all_data_qc), value = TRUE) [1] "HBB" "HBD" "HBG1" "HBG2" "HBE1" "HBZ" "HBM" "HBA2" "HBA1" "HBQ1"grep("PF4|PPBP", rownames(all_data_qc), value = TRUE) [1] "NBPF4" "YIPF4" "PRPF40A" "PF4V1" "PF4" "PPBP"
[7] "PPBPP2" "PRPF4B" "PRPF4" "PRPF40B" "APPBP2" "APPBP2-DT"all_data_qc <- PercentageFeatureSet(
all_data_qc,
"^MT-",
assay = "RNA",
col.name = "percent_mt"
)
all_data_qc <- PercentageFeatureSet(
all_data_qc,
"^RP[SL]",
assay = "RNA",
col.name = "percent_rb"
)
all_data_qc <- PercentageFeatureSet(
all_data_qc,
"^HB[ABDEGMQZ]\\d*$",
col.name = "percent_hb"
)
all_data_qc <- PercentageFeatureSet(
all_data_qc,
"PF4|PPBP",
assay = "RNA",
col.name = "percent_pt"
)head(all_data_qc[[]]) orig.ident nCount_RNA nFeature_RNA sample group
ctrl_1_AAACCTGAGCCAGTAG-1 ctrl_1 60591 7154 CTRL-1 CTRL
ctrl_1_AAACCTGAGTCCGGTC-1 ctrl_1 519 436 CTRL-1 CTRL
ctrl_1_AAACCTGCAAGCGCTC-1 ctrl_1 18100 4470 CTRL-1 CTRL
ctrl_1_AAACCTGCACACATGT-1 ctrl_1 2025 1244 CTRL-1 CTRL
ctrl_1_AAACCTGCAGACGCAA-1 ctrl_1 8300 3061 CTRL-1 CTRL
ctrl_1_AAACCTGCAGGTGCCT-1 ctrl_1 556 451 CTRL-1 CTRL
gender age percent_mt percent_rb percent_hb
ctrl_1_AAACCTGAGCCAGTAG-1 Male 62 2.3980459 27.223515 7.28821112
ctrl_1_AAACCTGAGTCCGGTC-1 Male 62 1.5414258 13.872832 0.19267823
ctrl_1_AAACCTGCAAGCGCTC-1 Male 62 0.9337017 42.370166 0.01104972
ctrl_1_AAACCTGCACACATGT-1 Male 62 5.0370370 12.000000 0.00000000
ctrl_1_AAACCTGCAGACGCAA-1 Male 62 2.4939759 22.168675 0.00000000
ctrl_1_AAACCTGCAGGTGCCT-1 Male 62 9.1726619 3.776978 0.71942446
percent_pt
ctrl_1_AAACCTGAGCCAGTAG-1 0.03630902
ctrl_1_AAACCTGAGTCCGGTC-1 0.00000000
ctrl_1_AAACCTGCAAGCGCTC-1 0.01104972
ctrl_1_AAACCTGCACACATGT-1 0.04938272
ctrl_1_AAACCTGCAGACGCAA-1 0.02409639
ctrl_1_AAACCTGCAGGTGCCT-1 0.00000000Violin plot before filtering
vln_qc <- list()
vln_qc$before_sample <-
VlnPlot(all_data_qc, group.by = "orig.ident", layer = "counts",
features = c("nFeature_RNA", "nCount_RNA", "percent_mt",
"percent_rb", "percent_hb", "percent_pt"),
pt.size = 0,
ncol = 3) + NoLegend()https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/tutorials/neutrophils
2^6 = 64
2^7.5 = 181
https://www.nature.com/articles/s41586-022-05400-x, methods section, scRNA-seq data processing
range(all_data_qc$nFeature_RNA)[1] 22 11084vln_qc$before_sample[[1]] <-
vln_qc$before_sample[[1]] +
geom_hline(yintercept = 5000, col = "red") +
geom_hline(yintercept = 100, col = "red")
vln_qc$before_sample[[1]]
range(all_data_qc$nCount_RNA)[1] 500 154080vln_qc$before_sample[[2]] <-
vln_qc$before_sample[[2]] +
geom_hline(yintercept = 20000, col = "red")
vln_qc$before_sample[[2]]
vln_qc$before_sample[[3]] <-
vln_qc$before_sample[[3]] +
geom_hline(yintercept = 5, col = "red")
vln_qc$before_sample[[3]]
# vln_qc$before_sample[[5]] <-
# vln_qc$before_sample[[5]] +
# geom_hline(yintercept = 25, col = "red")
# vln_qc$before_sample[[5]]vln_qc$before_sample
vln_qc$before_group <-
VlnPlot(all_data_qc, group.by = "group", layer = "counts",
features = c("nFeature_RNA", "nCount_RNA", "percent_mt",
"percent_rb", "percent_hb", "percent_pt"),
pt.size = 0,
ncol = 3) + NoLegend()
vln_qc$before_group
FeatureScatter(
all_data_qc,
"nCount_RNA",
"nFeature_RNA",
group.by = "orig.ident",
pt.size = 0.2
)
qs_save(
all_data_qc,
file = "data/rdata/preprocess/all_data_qc.qs2",
nthreads = 10
)
rm(all_data_qc)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 11441727 611.1 20569462 1098.6 20569462 1098.6
Vcells 23551927 179.7 991233340 7562.6 1232058104 9399.94. Detection-based Filtering
按检出基因数和 UMI 去掉空液滴、信息量太低的细胞。看各样本剩下多少，有没有某一例被削得很狠。
all_data_qc <- qs_read("data/rdata/preprocess/all_data_qc.qs2", nthreads = 10)
dim(all_data_qc)[1] 38606 79488selected_c <- WhichCells(
all_data_qc,
slot = "counts",
expression = nFeature_RNA > 100 &
nFeature_RNA < 5000 &
nCount_RNA < 20000
)
length(selected_c)[1] 74403scales::percent(
(dim(all_data_qc)[2] - length(selected_c)) / dim(all_data_qc)[2],
accuracy = 0.01
)[1] "6.40%"selected_f0 <- rownames(all_data_qc)[
Matrix::rowSums(all_data_qc@assays$RNA$counts) > 0
]
length(selected_f0)[1] 33077n_feature_sample <- rowSums(x = GetAssayData(all_data_qc, layer = "counts") > 0)
sum(n_feature_sample > 0)[1] 33077length(n_feature_sample)[1] 38606identical(selected_f0, names(n_feature_sample)[n_feature_sample > 0])[1] TRUEselected_f <- names(n_feature_sample)[n_feature_sample > 100]
length(selected_f)[1] 19681scales::percent(
(length(n_feature_sample) - length(selected_f)) /
length(n_feature_sample),
accuracy = 0.01
)[1] "49.02%"data_filt_qc <- subset(all_data_qc, cells = selected_c)
dim(data_filt_qc)[1] 38606 74403dim(all_data_qc)[1] 38606 79488rm(all_data_qc)counts_matrix <- data_filt_qc@assays$RNA$counts
counts_matrix@x <- counts_matrix@x /
rep.int(colSums(counts_matrix), diff(counts_matrix@p)) * 100
most_expressed <- order(Matrix::rowSums(counts_matrix), decreasing = TRUE)[20:1]
boxplot(as.matrix(t(counts_matrix[most_expressed, ])),
cex = 0.1, las = 1, xlab = "Percent counts per cell",
col = (scales::hue_pal())(20)[20:1], horizontal = TRUE)
rm(counts_matrix, most_expressed)5. Mito/Ribo Filtering
再按线粒体、核糖体比例过滤。阈值是经验性的，目的是去掉明显差的细胞，不是疾病分组。
selected_mito <- WhichCells(
data_filt_qc,
slot = "counts",
expression = percent_mt < 5
)
length(selected_mito)[1] 72346pct_removed <- scales::percent(
(dim(data_filt_qc)[2] - length(selected_mito)) / dim(data_filt_qc)[2],
accuracy = 0.01
)
pct_removed[1] "2.76%"data_filt_qc <- subset(data_filt_qc, cells = selected_mito)
dim(data_filt_qc)[1] 38606 72346table(data_filt_qc$orig.ident) |> sort()
ctrl_1 ctrl_5 ctrl_3 ctrl_2 ctrl_6 ctrl_4 aci_4 aci_2 aci_1 aci_5 aci_3
1953 2096 3671 6085 6828 6986 7871 8037 8064 9717 11038 Do not filter MT and RB genes
https://www.biostars.org/p/355658/
# data_filt_qc <- data_filt_qc[!grepl("^MT-", rownames(data_filt_qc)), ]
# data_filt_qc <- data_filt_qc[!grepl("^RP[SL]", rownames(data_filt_qc)), ]Violin plot after filtering
vln_qc$after_sample <-
VlnPlot(data_filt_qc, group.by = "orig.ident", layer = "counts",
features = c("nFeature_RNA", "nCount_RNA", "percent_mt",
"percent_rb", "percent_hb", "percent_pt"),
pt.size = 0,
ncol = 3) + NoLegend()
vln_qc$after_sample
vln_qc$after_group <-
VlnPlot(data_filt_qc, group.by = "group", layer = "counts",
features = c("nFeature_RNA", "nCount_RNA", "percent_mt",
"percent_rb", "percent_hb", "percent_pt"),
pt.size = 0,
ncol = 3) + NoLegend()
vln_qc$after_group
qs_save(vln_qc, file = "data/rdata/preprocess/violin_qc.qs2", nthreads = 10)
rm(vln_qc)
gc() used (Mb) gc trigger (Mb) max used (Mb)
Ncells 11490351 613.7 20569462 1098.6 20569462 1098.6
Vcells 160602711 1225.4 792986672 6050.1 1232058104 9399.9(Skipped) Check Sample Sex
性别来自临床记录，不从表达推断。
# suppressMessages(require(biomaRt))
# mart <- useMart("ENSEMBL_MART_ENSEMBL", dataset = "hsapiens_gene_ensembl")
#
# # fetch chromosome info plus some other annotations
# genes_table <- try(biomaRt::getBM(attributes = c("ensembl_gene_id",
# "external_gene_name",
# "description",
# "gene_biotype",
# "chromosome_name",
# "start_position"),
# mart = mart, useCache = F))
# write_tsv(genes_table, "data/genes_table.tsv")
#
# genes_table <- read_tsv("data/genes_table.tsv")
#
# chrY_gene <-
# genes_table |>
# filter(!is.na(external_gene_name)) |>
# filter(external_gene_name %in% rownames(data_filt_qc)) |>
# filter(chromosome_name == "Y") |>
# dplyr::select(external_gene_name) |>
# pull
# length(chrY_gene)
#
# data_filt_qc$pct_chrY = colSums(data_filt_qc@assays$RNA$counts[chrY_gene, ]) / colSums(data_filt_qc@assays$RNA$counts)
#
# VlnPlot(data_filt_qc, features = c("XIST", "pct_chrY"), layer = "counts", group.by = "patient")
# sample_info6. Cell Cycle Score
给每个细胞打一个周期分数，后面作图用。没有强制从表达里回归掉周期。
# data_filt_qc <- qs_read("data/rdata/preprocess/data_filt_qc.qs2", nthreads = 10) # nolintdata_filt_qc <- NormalizeData(data_filt_qc, verbose = TRUE)data_filt_qc <- CellCycleScoring(object = data_filt_qc,
g2m.features = cc.genes.updated.2019$g2m.genes,
s.features = cc.genes.updated.2019$s.genes)
colnames(data_filt_qc@meta.data)[12:14] <- c("s_score", "g2m_score", "phase")FeatureScatter(data_filt_qc, "s_score", "g2m_score", group.by = "phase")
data_filt_qc$cc_diff <- data_filt_qc$s_score - data_filt_qc$g2m_score
VlnPlot(data_filt_qc, features = c("s_score", "g2m_score", "cc_diff"),
group.by = "orig.ident", ncol = 3, pt.size = 0)
data_filt_qc[["RNA"]]$data <- NULL
qs_save(
data_filt_qc,
file = "data/rdata/preprocess/data_filt_qc.qs2",
nthreads = 10
)
rm(
data_filt_qc,
n_feature_sample,
selected_c,
selected_f,
selected_f0,
selected_mito
)7. Predict Doublets
预测两个细胞粘在一起的液滴。双细胞会把两种 marker 叠在同一个点上，干扰命名。
data_prep <- qs_read(
"data/rdata/preprocess/data_filt_qc.qs2",
nthreads = 10
)
bp_param <- SnowParam(workers = 15, RNGseed = 123)
sce <- as.SingleCellExperiment(data_prep)
sce <- scDblFinder(sce, samples = "sample", clusters = TRUE, BPPARAM = bp_param)
identical(colnames(data_prep), colnames(sce))[1] TRUEdata_prep$is_doublet <- sce$scDblFinder.class
data_prep$doublet_score <- sce$scDblFinder.scoreqs_save(
data_prep,
file = "data/rdata/preprocess/data_prep.qs2",
nthreads = 10
)
rm(sce, bp_param, data_prep)8. Prepare Final Data
写出质控后的对象，后面聚类都读这份。
data_prep <- qs_read(
"data/rdata/preprocess/data_prep.qs2",
nthreads = 10
)
Assays(data_prep)An object of class "SimpleAssays"
Slot "data":
List of length 1dim(data_prep)[1] 38606 72346data_prep$doublet_score |> range()[1] 3.379262e-06 9.999934e-01doublet_fraction <-
(dim(data_prep)[2] - sum(data_prep$is_doublet == "singlet")) /
dim(data_prep)[2]
scales::percent(doublet_fraction, accuracy = 0.01)[1] "7.33%"head(data_prep@meta.data) orig.ident nCount_RNA nFeature_RNA sample group
ctrl_1_AAACCTGAGTCCGGTC-1 ctrl_1 519 436 CTRL-1 CTRL
ctrl_1_AAACCTGCAAGCGCTC-1 ctrl_1 18100 4470 CTRL-1 CTRL
ctrl_1_AAACCTGCAGACGCAA-1 ctrl_1 8300 3061 CTRL-1 CTRL
ctrl_1_AAACCTGTCAAGGTAA-1 ctrl_1 769 615 CTRL-1 CTRL
ctrl_1_AAACCTGTCACAATGC-1 ctrl_1 669 387 CTRL-1 CTRL
ctrl_1_AAACGGGCAGATGGCA-1 ctrl_1 5152 1881 CTRL-1 CTRL
gender age percent_mt percent_rb percent_hb
ctrl_1_AAACCTGAGTCCGGTC-1 Male 62 1.5414258 13.87283 0.19267823
ctrl_1_AAACCTGCAAGCGCTC-1 Male 62 0.9337017 42.37017 0.01104972
ctrl_1_AAACCTGCAGACGCAA-1 Male 62 2.4939759 22.16867 0.00000000
ctrl_1_AAACCTGTCAAGGTAA-1 Male 62 2.0806242 11.31339 1.69050715
ctrl_1_AAACCTGTCACAATGC-1 Male 62 1.3452915 18.98356 3.73692078
ctrl_1_AAACGGGCAGATGGCA-1 Male 62 1.5333851 39.46040 0.03881988
percent_pt s_score g2m_score phase cc_diff
ctrl_1_AAACCTGAGTCCGGTC-1 0.00000000 0.09752453 -0.03345996 S 0.13098448
ctrl_1_AAACCTGCAAGCGCTC-1 0.01104972 0.59433181 0.05024735 S 0.54408446
ctrl_1_AAACCTGCAGACGCAA-1 0.02409639 -0.05256067 -0.09477101 G1 0.04221034
ctrl_1_AAACCTGTCAAGGTAA-1 0.13003901 0.07119995 -0.07176488 S 0.14296483
ctrl_1_AAACCTGTCACAATGC-1 0.00000000 0.14764930 0.03119341 S 0.11645589
ctrl_1_AAACGGGCAGATGGCA-1 0.09704969 -0.01050198 -0.04324619 G1 0.03274421
is_doublet doublet_score
ctrl_1_AAACCTGAGTCCGGTC-1 singlet 0.003584485
ctrl_1_AAACCTGCAAGCGCTC-1 singlet 0.139647678
ctrl_1_AAACCTGCAGACGCAA-1 singlet 0.154037639
ctrl_1_AAACCTGTCAAGGTAA-1 singlet 0.048985749
ctrl_1_AAACCTGTCACAATGC-1 singlet 0.002627888
ctrl_1_AAACGGGCAGATGGCA-1 singlet 0.023566755data_prep <-
qs_read("data/rdata/preprocess/data_prep.qs2", nthreads = 10)data_prep <-
NormalizeData(data_prep, verbose = TRUE)
data_prep <-
FindVariableFeatures(data_prep, nfeatures = 3000)
data_prep <-
ScaleData(data_prep) # vars.to.regress = c("nFeature_RNA", "percent_mt"))
data_prep <-
RunPCA(data_prep)
data_prep <-
RunUMAP(data_prep, dims = 1:30, reduction = "pca")VlnPlot(
data_prep,
features = "nFeature_RNA",
group.by = "is_doublet",
pt.size = 0
)
DimPlot(
data_prep,
group.by = "is_doublet",
raster = FALSE,
pt.size = 0.1,
alpha = 0.1
) +
NoAxes() +
theme(aspect.ratio = 1)
sessionInfo()R version 4.4.3 (2025-02-28)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 24.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.12.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Asia/Shanghai
tzcode source: system (glibc)
attached base packages:
[1] stats4 stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] BiocParallel_1.40.2 scDblFinder_1.20.2
[3] SingleCellExperiment_1.28.1 SummarizedExperiment_1.36.0
[5] Biobase_2.66.0 GenomicRanges_1.58.0
[7] GenomeInfoDb_1.42.0 IRanges_2.40.0
[9] S4Vectors_0.44.0 BiocGenerics_0.52.0
[11] MatrixGenerics_1.18.0 matrixStats_1.4.1
[13] qs2_0.1.7 lubridate_1.9.3
[15] forcats_1.0.0 stringr_1.5.1
[17] dplyr_1.1.4 purrr_1.0.2
[19] readr_2.1.5 tidyr_1.3.1
[21] tibble_3.2.1 ggplot2_3.5.1
[23] tidyverse_2.0.0 Seurat_5.2.1
[25] SeuratObject_5.0.2 sp_2.1-4
loaded via a namespace (and not attached):
[1] RcppAnnoy_0.0.22 splines_4.4.3 later_1.4.1
[4] BiocIO_1.16.0 bitops_1.0-9 R.oo_1.27.0
[7] polyclip_1.10-7 XML_3.99-0.17 fastDummies_1.7.4
[10] lifecycle_1.0.4 edgeR_4.4.0 globals_0.16.3
[13] lattice_0.22-6 MASS_7.3-65 magrittr_2.0.3
[16] limma_3.62.1 plotly_4.10.4 rmarkdown_2.29
[19] yaml_2.3.10 metapod_1.14.0 httpuv_1.6.15
[22] sctransform_0.4.1 spam_2.11-0 spatstat.sparse_3.1-0
[25] reticulate_1.42.0 cowplot_1.1.3 pbapply_1.7-2
[28] RColorBrewer_1.1-3 abind_1.4-8 zlibbioc_1.52.0
[31] Rtsne_0.17 R.utils_2.12.3 RCurl_1.98-1.16
[34] GenomeInfoDbData_1.2.13 ggrepel_0.9.6 irlba_2.3.5.1
[37] listenv_0.9.1 spatstat.utils_3.1-1 goftest_1.2-3
[40] RSpectra_0.16-2 dqrng_0.4.1 spatstat.random_3.3-2
[43] fitdistrplus_1.2-1 parallelly_1.39.0 codetools_0.2-20
[46] DelayedArray_0.32.0 scuttle_1.16.0 tidyselect_1.2.1
[49] UCSC.utils_1.2.0 farver_2.1.2 viridis_0.6.5
[52] ScaledMatrix_1.14.0 spatstat.explore_3.3-3 GenomicAlignments_1.42.0
[55] jsonlite_1.8.9 BiocNeighbors_2.0.0 progressr_0.15.1
[58] scater_1.34.0 ggridges_0.5.6 survival_3.8-3
[61] tools_4.4.3 snow_0.4-4 ica_1.0-3
[64] Rcpp_1.0.13-1 glue_1.8.0 gridExtra_2.3
[67] SparseArray_1.6.0 xfun_0.49 withr_3.0.2
[70] fastmap_1.2.0 bluster_1.16.0 fansi_1.0.6
[73] digest_0.6.37 rsvd_1.0.5 timechange_0.3.0
[76] R6_2.5.1 mime_0.12 colorspace_2.1-1
[79] scattermore_1.2 tensor_1.5 spatstat.data_3.1-4
[82] R.methodsS3_1.8.2 utf8_1.2.4 generics_0.1.3
[85] data.table_1.16.2 rtracklayer_1.66.0 httr_1.4.7
[88] htmlwidgets_1.6.4 S4Arrays_1.6.0 uwot_0.2.2
[91] pkgconfig_2.0.3 gtable_0.3.6 lmtest_0.9-40
[94] XVector_0.46.0 htmltools_0.5.8.1 dotCall64_1.2
[97] scales_1.3.0 png_0.1-8 spatstat.univar_3.1-1
[100] scran_1.34.0 knitr_1.49 rjson_0.2.23
[103] tzdb_0.4.0 reshape2_1.4.4 nlme_3.1-166
[106] curl_6.0.1 zoo_1.8-12 KernSmooth_2.23-26
[109] vipor_0.4.7 parallel_4.4.3 miniUI_0.1.1.1
[112] ggrastr_1.0.2 restfulr_0.0.15 pillar_1.9.0
[115] grid_4.4.3 vctrs_0.6.5 RANN_2.6.2
[118] promises_1.3.2 stringfish_0.18.0 BiocSingular_1.22.0
[121] beachmat_2.22.0 xtable_1.8-4 cluster_2.1.8.1
[124] beeswarm_0.4.0 evaluate_1.0.1 locfit_1.5-9.10
[127] cli_3.6.3 compiler_4.4.3 Rsamtools_2.22.0
[130] rlang_1.1.4 crayon_1.5.3 future.apply_1.11.3
[133] labeling_0.4.3 ggbeeswarm_0.7.2 plyr_1.8.9
[136] stringi_1.8.4 viridisLite_0.4.2 deldir_2.0-4
[139] munsell_0.5.1 Biostrings_2.74.0 lazyeval_0.2.2
[142] spatstat.geom_3.3-4 pacman_0.5.1 Matrix_1.7-3
[145] RcppHNSW_0.6.0 hms_1.1.3 patchwork_1.3.0
[148] future_1.34.0 statmod_1.5.0 shiny_1.9.1
[151] ROCR_1.0-11 igraph_2.1.2 RcppParallel_5.1.9
[154] xgboost_1.7.9.1