# scanpy_pipeline

**Repository Path**: mo_xinwu/scanpy_pipeline

## Basic Information

- **Project Name**: scanpy_pipeline
- **Description**: use scanpy and scvi-tools to analysis single cell RNA data
- **Primary Language**: Unknown
- **License**: GPL-3.0
- **Default Branch**: master
- **Homepage**: None
- **GVP Project**: No

## Statistics

- **Stars**: 1
- **Forks**: 0
- **Created**: 2023-06-28
- **Last Updated**: 2025-03-13

## Categories & Tags

**Categories**: Uncategorized

**Tags**: None

## README

# README
## 环境安装
使用 `conda_env/environment.yml` 即可完成环境配置。
```bash
mamba env create -f conda_env/environment.yml
```

## 路过踩坑
### 1. 创建 scanpy adata 对象
从 `gene x cell` 的表达矩阵创建单细胞分析对象(Seurat 或 scanpy), 都是非常重要且不可缺少的一种方式, 总是会有各种原因得不到 `cellranger` 的 `filtered_feature_bc_matrix` 或 `filtered_feature_bc_matrix.h5` 文件, 如：多样本数据无法通过 `cellranger aggr` 进行整合, 只能合并矩阵; 公共数据仅提供了 `gene x cell` 矩阵等原因。
|Symbol   |CAGCTGGCAAGCGTAG-2|TAGACCATCTATCCCG-2|GTAGTCATCTCCCTGA-1|GACCTGGTCCGCATAA-2|
|---------|------------------|------------------|------------------|------------------|
|LINC01128|                 0|                 0|                 0|                 0|
|HES4     |                 0|                 0|                 0|                 0|
|ISG15    |                 0|                 0|                 0|                 0|
|TNFRSF4  |                 0|                 0|                 0|                 0|

目前为止踩得最深的一个坑, 当时看别人教程时看到从 `gene x cell` 表格创建 `adata` 对象的操作是这样的, 如下：
```python
data = pd.read_table("gene_cell_exprs_table_symbol.xls",index_col = 0)
data.index.name = ''
data = pd.DataFrame(data.values.T,index = data.columns,columns = data.index)
adata = sc.AnnData(data)
```
这样创建的 `adata` 对象, 对于后续的操作(filter、doublet、write), 似乎没什么影响, 但最后使用 `adata.write('project.h5ad')` 保存后, 再使用 `sc.read_h5ad('project.h5ad')`, 就会报错, `anndata._io.utils.AnnDataReadError: Above error raised while reading key '/var' of type <class 'h5py._hl.dataset.Dataset'> from /.`, 无法读取。
这就导致分析结果无法保存的窘境, 那这个流程也就是无意义的。一度怀疑是不是 `scanpy` 版本的原因, 试了别的版本, 也查了下类似的报错, 都没有解决。

目前看起来一切正常的从 `gene x cell` 的表达矩阵创建 `adata` 的方式, 一句话就搞定, 整这么多花里胡哨的操作。
```python
adata = sc.read_csv("gene_cell_exprs_table_symbol.xls", delimiter = '\t').T
```

## scanpy pipeline 设计规范
1. 参数设计
 - 相对固定且复杂的参数, 如：物种的线粒体基因、核糖体基因、红细胞基因 ... 。
 - 常用且需要根据项目进行设置的参数。
 - 灵活调整的参数, 这种参数一般不进行调整, 但是特殊场景下会有需求, 可以赋予 pipeline 更多的可能性。如：默认根据 Cluster 进行 marker 基因分析, 但有时候需要调整到按照 CellType 进行分析。
解决方案:
 - 相对固定且复杂的参数 和 灵活调整的参数, 一律保存在 `config_scanpy_pipeline.json` 中, 通过解析 `json` 文件获取。
 - 需要根据项目灵活调整的参数, 在使用脚本时进行传递。

2. 数据过滤策略
 - 无论是单样本还是多样本在执行数据过滤, 即 根据 min_genes、max_genes、min_cells、MT、rRNA、redcell 是一致的, 不做区分。

3. 对双细胞的预测目前已经实现两种方式
 - scanpy 集成的 external.pp.scrublet: filter doublet 本质上是需要单个样本(批次)进行, 但是 scanpy 集成的 `external.pp.scrublet` 只需要指定 `batch_key` 就可按单个批次进行双细胞预测。
 - 使用 scvi-tools 进行双细胞预测, 概率模型, 运行时间比 external.pp.scrublet 上, 具体代码见 `external_function.py` 中的 `scvi_doublet_predict`, 为了节约 model train 的时间, 使用高变基因进行模型训练。

4. 批次效应去除策略
 - 默认不进行批次效应的去除
 - 去批次效应可以使用的方式有: combat、bbknn、harmony、scanorama
 - 使用 scvi-tools 去批次效应, scvi 是一个快速的多功能方法，整合了来自机器学习、统计模型、概率图模型等领域的思想和模型，应用到了单细胞数据的多种下游分析上来，相比其他针对单一下游分析的方法来说，从生物信息学的角度来讲是一个很大的进步。

5. 固定字段, 在 adata.obs 中使用如下固定字段
 - Sample：区分样本名称
 - CellType: 自定义的细胞类型

6. 在创建并进行细胞过滤后的 `Anndata` 对象, 使用 `adata.layers["counts"] = adata.X.copy()` 对原始表达数据(counts)进行保存。

## Demo data
demo data 使用新冠的研究文章 `A molecular single-cell lung atlas of lethal COVID-19` 中的数据, 该文章一共 27 个样本 `https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE171524` 链接中可以下载到表达矩阵文件 和 meta data 文件, 稍微处理一下, 将 27 个样本的数据处理成一个基础的 adata 文件(仅合并样本和附加 obs ,不进行任何处理), 代码见 `demo_data_create_adata.py`
```python
import scanpy as sc
adata = sc.read_h5ad('project_COVID.h5ad')
adata
# AnnData object with n_obs × n_vars = 116313 × 34546
#     obs: 'Sample', 'CellType', 'group'
adata.obs.head(5)
#                       Sample          CellType   group
# ATTCACTGTAACAGGC-1_1  C51ctr  Epithelial cells  normal
# TAACTTCCAACCACGC-1_1  C51ctr           Myeloid  normal
# TTGGGTACACGACAAG-1_1  C51ctr  Epithelial cells  normal
# AGGCCACAGAGTCACG-1_1  C51ctr  Epithelial cells  normal
# CACTGAAGTCGAAGCA-1_1  C51ctr  Epithelial cells  normal
```


## tips
1. `adata.obs.keys()` 与 `adata.obs.columns` 是一致的, 但`keys()`用途更广, 在使用数据前可以使用 `if` 或 `assert` 判断某些字段是否包含在 `adata` 对象中。
```python
# 使用 assert 断言, 如果 Sample 字段在 adata.obs.columns 中则不报错, 否则报 AssertionError, 来提前结束脚本。
batch_key = "Sample"
try:
    assert batch_key in adata.obs.columns
except AssertionError:
    print(f'{batch_key} not in adata.obs.columns, please check!')

# 使用 if raise 来判断并提前结束脚本。
if batch_key not in adata.obs.keys():
    raise ValueError('`batch_key` must be a column of .obs in the input annData object.')
```
| Command            | description              |
|--------------------|--------------------------|
| adata.obs.keys()   | 获取 obs 的列名           |
| adata.var.keys()   | 获取 var 的列名           |
| adata.uns.keys()   | 获取分析方法及参数         |
| adata.obsm.keys()  | 获取降维使用的方法         |
| adata.varm.keys()  | 获取 varm 中保存的数据名   |
| adata.layers.keys()| 获取 layers 中保存的数据名 |
| adata.obsp.keys()  | 获取 obsp 中保存的数据名   |

2. 矩阵 -> 稀疏矩阵
```python
from scipy.sparse import csr_matrix
adata.X = csr_matrix(adata.X) # 矩阵 -> 稀疏矩阵
adata.X = adata.X.todense()   # 稀疏矩阵 -> 矩阵
```

3. 从 `adata.raw` 取回数据
```python
adata.raw = adata  # freeze the state in `.raw`
bdata = adata.raw.to_adata()
# 基于我们最终能够取回到 过滤细胞后的全部细胞、全部基因 的原始数据 所以脚本中将 adata.raw = adata 操作提前至 细胞过滤后, normalize_total 之前进行。
```

## scanpy 不完美的地方
### Marker gene鉴定及分组差异分析
这一点个人认为还是挺伤的, 输出结果不像 Seurat 的 Marker gene 鉴定的输出结果, 且没有类似 Seurat 中 `FindMarkers` 的功能, 可用性较差。
1. Marker gene鉴定输出结果的差异
- Seurat：`FindAllMarkers`

|p_val|avg_log2FC      |pct.1|pct.2|p_val_adj|cluster|gene   |
|-----|----------------|-----|-----|---------|-------|-------|
|0    |1.17389378811916|0.65 |0.434|0        |1      |Stx1b  |
|0    |1.17349122562029|0.537|0.374|0        |1      |Ctnnb1 |
|0    |1.16742833651545|0.75 |0.493|0        |1      |Serinc1|
|0    |1.15207941822329|0.66 |0.398|0        |1      |Gm1673 |
|0    |1.13982337741822|0.351|0.277|0        |1      |Lin7a  |

pct.1: 在对应 cluster 细胞中检测到该基因表达的细胞比例
pct.2: 在其它 cluster 细胞中检测到该基因表达的细胞比例
- Scanpy: `rank_genes_groups`

Scanpy输出的结果, 很容易出现 `logfoldchanges` 为 `NaN` 的情况, 所以使用 cluster/cellType 的 marker gene 做 GO/KEGG 富集分析、火山图 之类的分析, 该怎么筛选差异显著的基因呢？
scanpy好像是根据 `scores` 来确定的 marker gene, 标准该怎么定呢？ Seurat 中筛选标准为 `Pvalue <= 0.05, 并且 >=1.5 倍表达差异的基因`, 在 scanpy 的输出结果下怎么定？

|group|names  |scores    |logfoldchanges |pvals              |pvals_adj         |
|-----|-------|----------|---------------|-------------------|------------------|
|1    |SFTPB  |167.24551 |               |0.0                |0.0               |
|1    |LRRK2  |156.55324 |               |0.0                |0.0               |
|1    |PDE4D  |149.79997 |               |0.0                |0.0               |
|1    |ZNF385B|149.71375 |               |0.0                |0.0               |
|1    |ANK3   |145.11653 |               |0.0                |0.0               |
|1    |SRSF12 |-0.832998 |3.25954        |0.404845856603561  |0.9999996480546481|
|1    |RAD54L |-0.8325566|2.187088       |0.4050948392834366 |0.9999996480546481|
|1    |FBXO15 |-0.8316275|2.2213638      |0.40561922944151585|0.9999996480546481|
|1    |EPG5   |0.8316002 |               |0.40563467102711626|0.9999996480546481|
|1    |ANKRD31|0.83153677|               |0.40567045368540255|0.9999996480546481|

<img src='http://xinwu-markdown.obs.cn-east-3.myhuaweicloud.com/xinwu-img/scanpy_pipeline/scanpy_cluster1_markergene.png'>

2. 分组差异分析

- Seurat

```R
# [情景1] 按分组比较全部细胞的差异, 比较样本Stim 与 Ctrl 的差异基因
Stim_vs_Ctrl <- FindMarkers(project, ident.1= "Stim", ident.2= "Ctrl",
                            min.pct = 0.1, logfc.threshold = 0.25,
                            group.by="Sample", test.use = "wilcox", only.pos = F)

# [情景2] 比较某一种细胞类型中组间的差异, 比较样本Stim 与 Ctrl 在 cluster 4、9 的差异基因
Bcell_SvsC <- FindMarkers(project, subset.ident = c("4","9"), 
                          ident.1 = "Stim", ident.2 = "Ctrl", 
                          min.pct = 0.1, logfc.threshold = 0.25, 
                          group.by = "Sample", test.use = "wilcox", only.pos = F)

# [情景3] 比较指定的2个细胞类型
cluster4_vs_cluster9 <- FindMarkers(project, ident.1="4", ident.2="9", 
                                    min.pct = 0.1, logfc.threshold = 0.25,
                                    test.use = "wilcox", only.pos = F)

# [情景4] 直接指定2个细胞集(Barcode list)进行比较
cell1 <- sample(colnames(project)[project$Sample == "Ctrl"],200)
cell2 <- sample(colnames(project)[project$Sample == "Stim"],200)
test_CvsS <- FindMarkers(project, ident.1 = cell1, ident.2 = cell2, only.pos = F) 
```
|gene    |p_val               |avg_log2FC        |pct.1|pct.2|p_val_adj            |
|--------|--------------------|------------------|-----|-----|---------------------|
|HLA-B   |0                   |0.666908199601761 |0.998|0.997|0                    |
|RPS29   |1.6140675923033e-123|-0.323522708719915|0.994|0.996|3.92960896022162e-119|
|HLA-DQB1|1.5512025229439e-67 |0.373008627782914 |0.111|0.01 |3.77655766235921e-63 |
|RPL38   |1.95262429735216e-67|-0.270204845113196|0.975|0.989|4.75385911433357e-63 |
|RPL36AL |2.95447468165896e-57|-0.311619976404586|0.91 |0.955|7.1929640599669e-53  |

- Scanpy

Scanpy 压根没有类似的 `FindMarkers` 方法, 所以要进行这样的操作或许只能将数据进行 `subset` 操作, 使用 `rank_genes_groups` 进行操作。
```python
# 寻找 0 vs 1 的差异基因, 局限性大, 没 seurat 好用, 要实现复杂点的还要自己凑, 或另寻他法
sc.tl.rank_genes_groups(adata, groupby='leiden', groups=['0'], reference='1', use_raw = False, method = 'wilcoxon')
```

### scvi-tools marker gene 鉴定及分组差异分析
鉴于 scanpy 的 marker gene 鉴定结果, 个人认为差强人意, 于是探索了一下 scvi-tools 做 marker gene 鉴定及分组差异分析。要使用 scvi-tools 做marker gene 鉴定及分组差异分析, 就需要进行 model train, 代码参考 `external_function.py` 中的 `scvi_doublet_predict` 或者 `scvi_DE` 。
- marker gene
```python
bdata = adata.copy()
scvi.model.SCVI.setup_anndata(adata = bdata, batch_key = 'Sample')
model = scvi.model.SCVI(bdata)
model.train()
# 注意 经过 model.train 后, 只要不改变 adata 的 shape, 可以对 adata 做任何操作 如 neighbors、tl.umap、tl.tsne、tl.leiden ...
# groupby 可以是在 model 构建之后加上去的 adata.obs 的 key, 没有问题的。
result = model.differential_expression(groupby = 'leiden')
result = result.rename_axis('gene').reset_index()
result.to_csv('scvi_marker_gene.xls', sep = "\t", header = True, index = False)
```
- 分组差异分析
```python
# [情景1] 按分组比较全部细胞的差异, 比较样本Stim 与 Ctrl 的差异基因
result = model.differential_expression(groupby='Sample', group1 = "Stim" , group2 = "Ctrl")
# [情景2] 比较某一种细胞类型中组间的差异, 比较样本Stim 与 Ctrl 在 cluster 4、9 的差异基因
result = model.differential_expression(idx1 = [(adata.obs.Sample == 'Stim') & (adata.obs['leiden_res_0.9'].isin([4,9]))],
                                       idx2 = [(adata.obs.Sample == 'Ctrl') & (adata.obs['leiden_res_0.9'].isin([4,9])])
# [情景3] 比较指定的2个细胞类型
result = model.differential_expression(groupby='leiden_res_0.9', group1 = "4" , group2 = "9")
# one v all
result = model.differential_expression(groupby='CellType', group1 = "CD4 T")
```
注： SCVI.differential_expression 的输出结果列数较多, 将近 20 列, 每列的具体含义后面要查一查。

## 更新
1. 2023-6-30
    - 基于我们想最终能够取回 过滤细胞后的全部细胞、全部基因 的原始数据 所以脚本中将 `adata.raw = adata` 操作提前至 `细胞过滤后, normalize_total 之前` 进行。
2. 2023-7-4
    - 在使用 scvi-tools 进行`双细胞预测`时, 依旧使用高变基因进行, 可以节约模型训练的时间。
    - 在使用 scvi-tools 进行`去批次效应`时, 不再使用 `高变基因` 而是使用 `全部基因` 进行, 并将 `model` 进行保存。原因是：去批次效应时训练出来的模型可以直接用于后续的 scvi-tools 进行差异分析(marker gene), 如果使用高变基因, 那使用 scvi-tools 进行差异分析时, 还要训练出一个全部基因的 model, 就显得没必要了。
    - 新增使用 scvi-tools 进行marker gene鉴定 和 分组差异分析

## 未完成
1. marker gene 的可视化, 热图 气泡图 小提琴图
2. 细胞类型(cluster/celltype)占比图
3. 细胞周期评分部分
4. 待更新...