flowchart TD
    A[实验设计与样本信息] --> B[Cell Ranger count]
    B --> C[读取 filtered_feature_bc_matrix]
    C --> D[细胞与基因质控]
    D --> E[双细胞 doublet 评估]
    E --> F[环境 RNA / 低质量细胞检查]
    F --> G[SCTransform 或 NormalizeData]
    G --> H[高变基因 / PCA]
    H --> I{是否多样本或多批次?}
    I -- 是 --> J[IntegrateLayers / Harmony / scVI]
    I -- 否 --> K[FindNeighbors]
    J --> K
    K --> L[FindClusters]
    L --> M[UMAP / t-SNE 可视化]
    M --> N[Marker gene 鉴定]
    N --> O[细胞类型注释]
    O --> P[差异表达 / 富集 / 拟时序 / 通讯]

single_cell_project/
├── 00_metadata/
│   └── sample_info.csv
├── 01_fastq/
│   ├── sampleA/
│   └── sampleB/
├── 02_cellranger/
│   ├── sampleA/outs/
│   └── sampleB/outs/
├── 03_seurat/
│   ├── rds/
│   ├── tables/
│   └── plots/
├── 04_annotation/
├── 05_downstream/
└── reports/

cellranger count   --id=sampleA   --transcriptome=/ref/refdata-gex-GRCh38-2024-A   --fastqs=/data/fastq/sampleA   --sample=sampleA   --localcores=24   --localmem=128

REF=/ref/refdata-gex-GRCh38-2024-A
FASTQ_ROOT=/data/fastq
OUT_ROOT=02_cellranger

for sample in sampleA sampleB sampleC
do
  cellranger count     --id=${sample}     --transcriptome=${REF}     --fastqs=${FASTQ_ROOT}/${sample}     --sample=${sample}     --localcores=24     --localmem=128

  mv ${sample} ${OUT_ROOT}/
done

library(Seurat)
library(dplyr)
library(ggplot2)
library(patchwork)

sample_info <- read.csv("00_metadata/sample_info.csv")

seurat_list <- lapply(seq_len(nrow(sample_info)), function(i) {
  sample_id <- sample_info$sample_id[i]
  matrix_dir <- file.path("02_cellranger", sample_id, "outs", "filtered_feature_bc_matrix")

  counts <- Read10X(data.dir = matrix_dir)
  obj <- CreateSeuratObject(
    counts = counts,
    project = sample_id,
    min.cells = 3,
    min.features = 200
  )

  obj$sample_id <- sample_id
  obj$condition <- sample_info$condition[i]
  obj$batch <- sample_info$batch[i]
  obj$tissue <- sample_info$tissue[i]
  obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
  obj[["percent.ribo"]] <- PercentageFeatureSet(obj, pattern = "^RP[SL]")
  obj
})

names(seurat_list) <- sample_info$sample_id

qc_plot <- VlnPlot(
  seurat_list[[1]],
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
  ncol = 3,
  pt.size = 0.1
)
ggsave("03_seurat/plots/qc_violin_sampleA.pdf", qc_plot, width = 12, height = 4)

filter_one <- function(obj) {
  subset(
    obj,
    subset =
      nFeature_RNA > 300 &
      nFeature_RNA < 6000 &
      nCount_RNA > 500 &
      percent.mt < 15
  )
}

seurat_list <- lapply(seurat_list, filter_one)

# DoubletFinder 参数需要根据数据量和预期 doublet rate 调整。
# MVP 阶段建议先标记 doublet，再比较去除前后的 UMAP 和 marker 表。

seurat_list <- lapply(seurat_list, function(obj) {
  SCTransform(
    obj,
    vars.to.regress = "percent.mt",
    verbose = FALSE
  )
})

obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000)
obj <- ScaleData(obj, vars.to.regress = c("percent.mt"))

combined <- merge(
  seurat_list[[1]],
  y = seurat_list[-1],
  add.cell.ids = names(seurat_list)
)

combined[["RNA"]] <- split(combined[["RNA"]], f = combined$sample_id)
combined <- NormalizeData(combined)
combined <- FindVariableFeatures(combined)
combined <- ScaleData(combined)
combined <- RunPCA(combined)

combined <- IntegrateLayers(
  object = combined,
  method = RPCAIntegration,
  orig.reduction = "pca",
  new.reduction = "integrated.rpca",
  verbose = FALSE
)

combined <- FindNeighbors(combined, reduction = "integrated.rpca", dims = 1:30)
combined <- FindClusters(combined, resolution = 0.4)
combined <- FindClusters(combined, resolution = 0.8)
combined <- RunUMAP(combined, reduction = "integrated.rpca", dims = 1:30)

p1 <- DimPlot(combined, group.by = "sample_id")
p2 <- DimPlot(combined, group.by = "seurat_clusters", label = TRUE)
ggsave("03_seurat/plots/umap_sample_cluster.pdf", p1 + p2, width = 12, height = 5)

combined <- PrepSCTFindMarkers(combined)

markers <- FindAllMarkers(
  combined,
  only.pos = TRUE,
  min.pct = 0.25,
  logfc.threshold = 0.25
)

write.csv(markers, "03_seurat/tables/all_cluster_markers.csv", row.names = FALSE)

new_ids <- c(
  "0" = "CD4 T cells",
  "1" = "CD14 Monocytes",
  "2" = "B cells",
  "3" = "CD8 T cells",
  "4" = "NK cells"
)

combined$celltype <- plyr::mapvalues(
  x = combined$seurat_clusters,
  from = names(new_ids),
  to = as.vector(new_ids)
)

DimPlot(combined, group.by = "celltype", label = TRUE)
ggsave("03_seurat/plots/umap_celltype.pdf", width = 8, height = 6)

counts_by_celltype <- AggregateExpression(
  combined,
  group.by = c("celltype", "sample_id"),
  assays = "RNA",
  slot = "counts",
  return.seurat = FALSE
)

import scanpy as sc

adata = sc.read_10x_mtx(
    "02_cellranger/sampleA/outs/filtered_feature_bc_matrix",
    var_names="gene_symbols",
    cache=True,
)

adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

adata = adata[
    (adata.obs["n_genes_by_counts"] > 300)
    & (adata.obs["n_genes_by_counts"] < 6000)
    & (adata.obs["pct_counts_mt"] < 15),
    :
]

sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
adata = adata[:, adata.var.highly_variable]
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.leiden(adata, resolution=0.8)
sc.tl.umap(adata)
sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon")

adata.write("03_scanpy/sampleA_processed.h5ad")

prop_df <- combined@meta.data |>
  dplyr::count(sample_id, condition, celltype) |>
  dplyr::group_by(sample_id) |>
  dplyr::mutate(prop = n / sum(n))

ggplot(prop_df, aes(x = sample_id, y = prop, fill = celltype)) +
  geom_col(width = 0.8) +
  facet_grid(. ~ condition, scales = "free_x", space = "free_x") +
  theme_bw() +
  labs(x = "Sample", y = "Cell proportion")