190_单细胞线粒体突变追踪¶

# 安装mgatk
pip install mgatk

# 从10x CellRanger BAM文件提取线粒体变异
mgatk tenx \
    -i possorted_genome_bam.bam \
    -o mgatk_output \
    -c filtered_barcodes.tsv \
    -bt CB \
    -b chrM \
    --nsamples 1000 \
    --ncores 16

# 输出文件
# mgatk_output/
#   final/
#     chrM.variant_stats.tsv.gz     # 变异统计
#     chrM.coverage.txt.gz           # 覆盖度
#     mgatk.rds                       # R对象

library(Signac)
library(Seurat)
library(ggplot2)

# 加载mgatk结果
mgatk <- ReadMGATK("mgatk_output/final/")

# 创建Assay
seurat_obj[["mito"]] <- CreateAssayObject(counts = mgatk$counts)
seurat_obj <- AddMetaData(seurat_obj, metadata = mgatk$depth, col.name = "mtDNA_depth")

# 质控：过滤低覆盖度细胞
VlnPlot(seurat_obj, features = "mtDNA_depth", pt.size = 0)
seurat_obj <- subset(seurat_obj, mtDNA_depth > 10)  # 最低10x线粒体覆盖度

# 识别信息性变异位点（高变异、真实的体细胞突变）
# 使用变异等位基因频率（VAF）
variable_sites <- IdentifyVariants(seurat_obj, assay = "mito", refallele = mgatk$refallele)

# 过滤信息性位点
# strand_correlation > 0.65: 正负链一致性
# vmr > 0.01: 变异足够大
# mean_coverage > 2: 足够的覆盖度
high_conf <- subset(variable_sites,
                     subset = n_cells_conf_detected >= 5 &
                       strand_correlation >= 0.65 &
                       vmr > 0.01)

informative_variants <- high_conf$variant
cat("信息性变异位点数:", length(informative_variants), "\n")

# 提取变异等位基因频率矩阵
DefaultAssay(seurat_obj) <- "mito"
vaf_matrix <- GetAssayData(seurat_obj, slot = "data")[informative_variants, ]

# 基于VAF进行聚类
# 方法1：简单层次聚类
library(pheatmap)
pheatmap(as.matrix(vaf_matrix),
         show_colnames = FALSE,
         clustering_distance_rows = "correlation",
         clustering_distance_cols = "correlation",
         main = "Mitochondrial Mutation VAF")

# 方法2：使用Seurat标准流程
seurat_obj <- ScaleData(seurat_obj, features = informative_variants)
seurat_obj <- RunPCA(seurat_obj, features = informative_variants, reduction.name = "mito_pca")
seurat_obj <- RunUMAP(seurat_obj, reduction = "mito_pca", dims = 1:10,
                       reduction.name = "mito_umap")
seurat_obj <- FindNeighbors(seurat_obj, reduction = "mito_pca", dims = 1:10)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)

# 可视化
DimPlot(seurat_obj, reduction = "mito_umap", label = TRUE) +
  ggtitle("Clonal Structure (mtDNA mutations)")

# 方法3：构建系统发育树
library(ape)
dist_matrix <- dist(t(as.matrix(vaf_matrix)))
tree <- nj(dist_matrix)
plot(tree, type = "unrooted", show.tip.label = FALSE,
     main = "Cell Lineage Tree (mtDNA)")

# 将线粒体克隆信息与细胞类型注释整合
# 检查不同克隆中的细胞类型分布
table(seurat_obj$mito_clusters, seurat_obj$cell_type)

# 可视化特定突变在UMAP上的分布
FeaturePlot(seurat_obj, features = c("3243A>G", "8344A>G"),
            reduction = "umap", min.cutoff = 0.01)

# 克隆-细胞类型桑基图
library(ggalluvial)
clone_celltype <- as.data.frame(table(
  Clone = seurat_obj$mito_clusters,
  CellType = seurat_obj$cell_type
))

ggplot(clone_celltype, aes(axis1 = Clone, axis2 = CellType, y = Freq)) +
  geom_alluvium(aes(fill = Clone), width = 0.3) +
  geom_stratum(width = 0.3) +
  geom_text(stat = "stratum", aes(label = after_stat(stratum)), size = 3) +
  theme_minimal() +
  ggtitle("Clone-CellType Relationship")

# mgatk从10x BAM提取
mgatk tenx -i possorted.bam -o output/ -c barcodes.tsv -bt CB -b chrM --ncores 16

# 检查BAM中线粒体reads
samtools view -c -F 4 possorted.bam chrM

# 每个细胞的线粒体覆盖度
samtools view -F 4 possorted.bam chrM | grep -oP 'CB:Z:\K[^ ]+' | sort | uniq -c | sort -rn