从10xGenomics的bam文件分析pre-mRNA数据

当只有 BAM 文件时，那么需要从下载参考基因组和注释文件开始，然后使用这些文件进行比对和注释，最终通过 velocyto 分析生成 .loom 文件进行pre-mRNA分析。以下是详细的步骤：

---

关于配置：

1. 电脑配置

1.1 最低配置

• CPU: 至少8核
• 内存: 至少32GB
• 存储: 至少500GB SSD（用于存储中间文件和结果）
• 操作系统: Linux或macOS（推荐）

1.2 推荐配置

• CPU: 16核或更多
• 内存: 64GB或更多
• 存储: 1TB SSD或更多
• GPU: 可选，但可以加速某些计算任务

2. 软件工具准备

2.1 安装必要的软件和工具

确保你已经安装了以下工具：

• samtools: 用于处理BAM文件。
• STAR: 用于比对和量化。
• velocyto: 用于分析premRNA和mRNA。
• Python: velocyto需要Python环境。

2.2 下载参考基因组和注释文件

从Ensembl或UCSC下载参考基因组和GTF注释文件。

---

1. 下载参考基因组和注释文件

1.1 下载参考基因组

• 小鼠参考基因组：推荐使用 GRCm39（Ensembl）或 mm39（UCSC）。
  • Ensembl:
    • 访问 Ensembl FTP站点
    • 下载文件：Mus_musculus.GRCm39.dna.primary_assembly.fa.gz
  • UCSC:
    • 访问 UCSC Genome Browser
    • 下载文件：mm39.fa.gz

1.2 下载 GTF 注释文件

• 小鼠 GTF 注释文件：确保与参考基因组版本一致。
  • Ensembl:
    • 访问 Ensembl FTP站点
    • 下载文件：Mus_musculus.GRCm39.109.gtf.gz
  • UCSC:
    • 访问 UCSC Genome Browser
    • 下载文件：mm39.refGene.gtf.gz

1.3 解压文件

下载完成后，解压参考基因组和 GTF 注释文件：

gunzip Mus_musculus.GRCm39.dna.primary_assembly.fa.gz
gunzip Mus_musculus.GRCm39.109.gtf.gz

---

2. 使用 velocyto 进行 premRNA 分析

2.1 安装 velocyto

确保已安装 velocyto：

pip install velocyto

2.2 运行 velocyto

假设你已经有了以下文件：

• BAM 文件：your_file.bam
• Barcode 文件：barcodes.tsv（10x Genomics 的 barcode 文件）
• 参考基因组：Mus_musculus.GRCm39.dna.primary_assembly.fa
• GTF 注释文件：Mus_musculus.GRCm39.109.gtf

运行以下命令：

velocyto run -b barcodes.tsv -o output_dir -m repeat_msk.gtf your_file.bam Mus_musculus.GRCm39.109.gtf

• -b barcodes.tsv：10x Genomics 的 barcode 文件。
• -o output_dir：输出目录。
• -m repeat_msk.gtf：重复序列的 GTF 文件（可选，如果没有可以省略）。
• your_file.bam：你的 BAM 文件。
• Mus_musculus.GRCm39.109.gtf：GTF 注释文件。

2.3 输出文件

运行完成后，output_dir 中会生成一个 .loom 文件，例如 your_file.loom。这个文件包含了：

• spliced 和 unspliced 表达矩阵。
• 细胞和基因的元数据。

如果 .loom 文件中包含细胞和基因的元数据，可以对细胞进行分群（聚类），并进一步分析每一群细胞中的 spliced 和 unspliced 表达矩阵。

---

3. 加载 .loom 文件并进行分群

3.1 使用 scanpy 加载 .loom 文件

import scanpy as sc

# 加载 .loom 文件
adata = sc.read_loom("output_dir/your_file.loom")

# 查看数据
print(adata)

3.2 数据预处理

在分群之前，需要对数据进行标准化和筛选高变基因：

# 标准化
sc.pp.normalize_total(adata, target_sum=1e4)  # 将每个细胞的总表达量标准化到 10,000
sc.pp.log1p(adata)  # 对数转换

# 筛选高变基因
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata = adata[:, adata.var.highly_variable]  # 保留高变基因

3.3 降维和聚类

使用 PCA、UMAP 和 Leiden 聚类对细胞进行分群：

# 降维
sc.pp.pca(adata, n_comps=50)  # PCA
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)  # 构建 KNN 图
sc.tl.umap(adata)  # UMAP 降维

# 聚类
sc.tl.leiden(adata, resolution=0.5)  # Leiden 聚类

# 可视化
sc.pl.umap(adata, color='leiden')  # 可视化聚类结果

结果解释

• adata.obs['leiden'] 中存储了每个细胞的聚类标签。
• 聚类结果可以通过 UMAP 或 t-SNE 可视化。

---

4. 提取某一群细胞的 spliced 和 unspliced 表达矩阵

4.1 根据聚类标签提取细胞

假设你想提取聚类标签为 1 的细胞：

# 提取聚类标签为 1 的细胞
cluster_1_cells = adata[adata.obs['leiden'] == '1']

# 提取 spliced 和 unspliced 表达矩阵
spliced_matrix_cluster_1 = cluster_1_cells.layers['spliced']
unspliced_matrix_cluster_1 = cluster_1_cells.layers['unspliced']

4.2 分析特定基因的表达

如果你想分析某一基因（如 GeneA）在聚类 1 中的表达：

# 获取基因的索引
gene_index = adata.var_names.get_loc('GeneA')

# 提取 GeneA 的 spliced 和 unspliced 表达量
spliced_geneA_cluster_1 = spliced_matrix_cluster_1[:, gene_index]
unspliced_geneA_cluster_1 = unspliced_matrix_cluster_1[:, gene_index]

---

5. 分析某一群细胞的基因表达动态

5.1 计算群内基因的平均表达量

import numpy as np

# 计算聚类 1 中所有基因的平均 spliced 和 unspliced 表达量
mean_spliced_cluster_1 = np.array(spliced_matrix_cluster_1.mean(axis=0)).flatten()
mean_unspliced_cluster_1 = np.array(unspliced_matrix_cluster_1.mean(axis=0)).flatten()

5.2 可视化基因表达动态

import matplotlib.pyplot as plt

# 绘制 spliced 和 unspliced 表达量的散点图
plt.scatter(mean_spliced_cluster_1, mean_unspliced_cluster_1, alpha=0.5)
plt.xlabel('Spliced Expression')
plt.ylabel('Unspliced Expression')
plt.title('Spliced vs Unspliced Expression in Cluster 1')
plt.show()

---

6. 进一步分析

6.1 差异表达分析

比较不同聚类之间的基因表达差异：

# 使用 scanpy 进行差异表达分析
sc.tl.rank_genes_groups(adata, groupby='leiden', method='wilcoxon')

# 查看聚类 1 的差异表达基因
sc.pl.rank_genes_groups(adata, groups=['1'], n_genes=20)

6.2 RNA 速度分析

研究某一群细胞的 RNA 速度：

# 计算 RNA 速度
scv.tl.velocity(adata, mode='stochastic')

# 可视化 RNA 速度
scv.pl.velocity_embedding(adata, basis='umap', color='leiden')

7. 总结

• 细胞分群：可以通过降维和聚类方法（如 UMAP + Leiden）对细胞进行分群。
• 提取表达矩阵：根据聚类标签提取某一群细胞的 spliced 和 unspliced 表达矩阵。
• 分析基因表达动态：计算群内基因的平均表达量，并可视化 spliced 和 unspliced 表达关系。
• 进一步分析：可以进行差异表达分析、RNA 速度分析等，深入研究细胞群的特性。

通过这些步骤，你可以对单细胞数据中的细胞进行分群，并分析每一群细胞的转录动态和基因表达特征。

---