384_单细胞核RNA测序snRNA-seq¶

# CellRanger 统计内含子reads（snRNA-seq必须开启）
cellranger count \
  --id sample_nuclear \
  --fastqs /path/to/fastq/ \
  --sample sample_name \
  --transcriptome /path/to/refdata-gex-GRCh38-2024-A \
  --include-introns true \    # 【关键】包含内含子reads（默认false，scRNA-seq才是false）
  --localcores 16 \
  --localmem 64

# STARsolo 等价参数
STARsolo \
  --soloType CB_UMI_Simple \
  --soloCBwhitelist 3M-february-2018.txt \
  --soloFeatures Gene GeneFull \   # GeneFull = 包含内含子
  --readFilesIn R2.fastq.gz R1.fastq.gz \
  --runThreadN 16

# 查看每个核的统计数据（过滤标准与scRNA-seq不同）
echo "snRNA-seq 推荐质控阈值："
echo "  每核基因数：200-5000（scRNA-seq是200-6000）"
echo "  线粒体基因比例：<5%（核内线粒体基因少）"
echo "  双核比例：<15%"

library(SoupX)      # 环境RNA去除
library(Seurat)

# ============================================================
# 方法A：从 CellRanger 输出自动运行（推荐）
# ============================================================
# CellRanger 输出目录结构：
# outs/
#   filtered_feature_bc_matrix/   <- 过滤后的矩阵
#   raw_feature_bc_matrix/        <- 原始矩阵（包含空液滴）
#   molecule_info.h5

sc <- load10X("outs/")   # 自动加载 filtered + raw 两个矩阵

# SoupX 估计"汤"的组成
sc <- autoEstCont(sc)    # 自动估计环境RNA污染比例
cat("估计污染比例:", sc$soupProfile[1], "\n")  # 通常5-20%

# 去除环境RNA（矫正count矩阵）
out_matrix <- adjustCounts(
  sc,
  roundToInt = TRUE    # 输出整数（DESeq2等需要整数counts）
)

# 创建去污染后的 Seurat 对象
seurat_clean <- CreateSeuratObject(
  counts = out_matrix,
  project = "snRNA_clean",
  min.cells = 3,
  min.features = 200
)

# ============================================================
# 方法B：手动指定marker基因辅助估计（当自动估计不准时）
# ============================================================
# 例如，MALAT1 是核RNA的标志（高表达=真核，对外部RNA也高=污染严重）
sc <- setSoupProfile(sc, "MALAT1")   # 用MALAT1辅助校正
sc <- autoEstCont(sc, priorRho = 0.1)   # prior=0.1表示预期10%污染

import scanpy as sc
import scrublet as scr    # 双核检测（pip install scrublet）
import numpy as np
import matplotlib.pyplot as plt

# 读取 CellRanger 输出（GeneFull矩阵，包含内含子）
adata = sc.read_10x_mtx(
    "outs/filtered_feature_bc_matrix/",
    var_names="gene_symbols",
    cache=True
)
adata.var_names_make_unique()

# ============================================================
# 质控（snRNA-seq有特殊阈值）
# ============================================================
# 计算线粒体比例（核内线粒体DNA转录的mRNA很少）
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)

# snRNA-seq 特殊质控：标记核糖体基因（核内rRNA很少）
adata.var["ribo"] = adata.var_names.str.startswith(("RPS", "RPL"))
sc.pp.calculate_qc_metrics(adata, qc_vars=["ribo"], inplace=True)

# 核的质量可视化
sc.pl.violin(adata,
    ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
    jitter=0.4,
    multi_panel=True,
    save="_snRNA_qc.pdf"
)

# 过滤（snRNA-seq阈值宽松一些）
sc.pp.filter_cells(adata, min_genes=200)      # 最少200个基因
sc.pp.filter_genes(adata, min_cells=3)        # 最少3个核表达

adata = adata[
    (adata.obs.pct_counts_mt < 5) &          # 线粒体比例<5%（核内更低）
    (adata.obs.n_genes_by_counts < 5000)     # 上限宽松些（核RNA稍少）
].copy()

print(f"过滤后：{adata.n_obs} 个核，{adata.n_vars} 个基因")

# ============================================================
# 双核检测（Scrublet）
# ============================================================
scrub = scr.Scrublet(adata.X, expected_doublet_rate=0.1)  # snRNA-seq双核率更高
doublet_scores, predicted_doublets = scrub.scrub_doublets(
    min_counts=2,
    min_cells=3,
    min_gene_variability_pctl=85
)
adata.obs["doublet_score"] = doublet_scores         # 双核得分
adata.obs["predicted_doublet"] = predicted_doublets  # 是否双核

# 过滤双核
n_before = adata.n_obs
adata = adata[~adata.obs.predicted_doublet].copy()
print(f"去除双核后：{adata.n_obs} 个核（去掉 {n_before - adata.n_obs} 个双核）")

# 标准化和降维
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, n_top_genes=3000, flavor="seurat_v3")
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata, n_pcs=30)
sc.tl.umap(adata)
sc.tl.leiden(adata, resolution=0.5)
sc.pl.umap(adata, color="leiden", save="_snRNA_clusters.pdf")