當只有 BAM 文件時，那么需要從下載參考基因組和注釋文件開始稍刀，然后使用這些文件進行比對和注釋账月，最終通過 velocyto 分析生成 .loom 文件進行pre-mRNA分析。以下是詳細的步驟：
關(guān)于配置：
1. 電腦配置
1.1 最低配置
- CPU: 至少8核
- 內(nèi)存: 至少32GB
- 存儲: 至少500GB SSD（用于存儲中間文件和結(jié)果）
- 操作系統(tǒng): Linux或macOS（推薦）
1.2 推薦配置
- CPU: 16核或更多
- 內(nèi)存: 64GB或更多
- 存儲: 1TB SSD或更多
- GPU: 可選剧劝，但可以加速某些計算任務(wù)
2. 軟件工具準備
2.1 安裝必要的軟件和工具
確保你已經(jīng)安裝了以下工具：
- samtools: 用于處理BAM文件讥此。
- STAR: 用于比對和量化。
- velocyto: 用于分析premRNA和mRNA卒稳。
- Python: velocyto需要Python環(huán)境展哭。
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
-
Ensembl:
1.2 下載 GTF 注釋文件
-
小鼠 GTF 注釋文件：確保與參考基因組版本一致役衡。
-
Ensembl:
- 訪問 Ensembl FTP站點
- 下載文件：
Mus_musculus.GRCm39.109.gtf.gz
-
UCSC:
- 訪問 UCSC Genome Browser
- 下載文件：
mm39.refGene.gtf.gz
-
Ensembl:
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
假設(shè)你已經(jīng)有了以下文件：
-
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：重復(fù)序列的 GTF 文件（可選吨述，如果沒有可以省略）钞脂。 -
your_file.bam：你的 BAM 文件冰啃。 -
Mus_musculus.GRCm39.109.gtf：GTF 注釋文件。
2.3 輸出文件
運行完成后焚刚，output_dir 中會生成一個 .loom 文件汪榔，例如 your_file.loom肃拜。這個文件包含了：
- spliced 和 unspliced 表達矩陣。
- 細胞和基因的元數(shù)據(jù)士聪。
如果 .loom 文件中包含細胞和基因的元數(shù)據(jù)剥悟，可以對細胞進行分群（聚類），并進一步分析每一群細胞中的 spliced 和 unspliced 表達矩陣区岗。
3. 加載 .loom 文件并進行分群
3.1 使用 scanpy 加載 .loom 文件
import scanpy as sc
# 加載 .loom 文件
adata = sc.read_loom("output_dir/your_file.loom")
# 查看數(shù)據(jù)
print(adata)
3.2 數(shù)據(jù)預(yù)處理
在分群之前慈缔，需要對數(shù)據(jù)進行標準化和篩選高變基因：
# 標準化
sc.pp.normalize_total(adata, target_sum=1e4) # 將每個細胞的總表達量標準化到 10,000
sc.pp.log1p(adata) # 對數(shù)轉(zhuǎn)換
# 篩選高變基因
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) # 構(gòu)建 KNN 圖
sc.tl.umap(adata) # UMAP 降維
# 聚類
sc.tl.leiden(adata, resolution=0.5) # Leiden 聚類
# 可視化
sc.pl.umap(adata, color='leiden') # 可視化聚類結(jié)果
結(jié)果解釋
-
adata.obs['leiden']中存儲了每個細胞的聚類標簽赂韵。 - 聚類結(jié)果可以通過 UMAP 或 t-SNE 可視化祭示。
4. 提取某一群細胞的 spliced 和 unspliced 表達矩陣
4.1 根據(jù)聚類標簽提取細胞
假設(shè)你想提取聚類標簽為 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. 分析某一群細胞的基因表達動態(tài)
5.1 計算群內(nèi)基因的平均表達量
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 可視化基因表達動態(tài)
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. 總結(jié)
- 細胞分群：可以通過降維和聚類方法（如 UMAP + Leiden）對細胞進行分群悄窃。
-
提取表達矩陣：根據(jù)聚類標簽提取某一群細胞的
spliced和unspliced表達矩陣轧抗。 -
分析基因表達動態(tài)：計算群內(nèi)基因的平均表達量瞬测，并可視化
spliced和unspliced表達關(guān)系月趟。 - 進一步分析：可以進行差異表達分析孝宗、RNA 速度分析等，深入研究細胞群的特性问潭。
通過這些步驟，你可以對單細胞數(shù)據(jù)中的細胞進行分群梳虽，并分析每一群細胞的轉(zhuǎn)錄動態(tài)和基因表達特征灾茁。
