OmicVerse是用Python进行多组学（包括Bulk和单细胞分析）的基础框架。前面我们在公众号宣传过一波; Python的转录组学分析框架与生态，因为是需要去github点star后发邮件才能进群交流，所以操作门槛有点高， 我们后续再次开放拉群小助手给大家哈。

欲了解更多信息，请阅读我们的论文：OmicVerse: A single pipeline for exploring the entire transcriptome universe

omicverse最初的名字是Pyomic，但我们希望涵盖整个转录组学的领域，因此将其改名为OmicVerse，旨在解决RNA-seq中的所有任务。目前omicverse已加入scverse生态。您可以在scverse的官网上找到我们。

安装前准备

OmicVerse可以通过conda或pypi进行安装，不过您需要先安装pytorch

为避免潜在的依赖冲突，建议在conda环境中安装。并使用pip install -U omicverse进行更新。

平台

在不同的平台上，最合适的安装方法有所不同。

• Windows：我们建议安装wsl子系统，并在wsl子系统中安装conda以配置omicverse环境。
• Linux：我们可以选择安装anaconda或miniconda，然后使用conda来配置omicverse环境。
• Mac Os：我们建议使用miniforge或mambaforge进行配置。

pip前提条件

• 如果使用conda/mamba，则只需运行conda install -c anaconda pip并跳过此部分。
• 安装Python，我们推荐使用pyenv版本管理系统，以及pyenv-virtualenv。

Apple Silicon前提条件

在装有Apple Silicon的Mac上安装omicverse只能使用本机版本的Python。可以通过使用Apple Silicon版本的mambaforge（可以通过本机版本的homebrew通过brew install --cask mambaforge进行安装）来安装本机版本的Python。

Conda

1. 安装conda。我们通常使用mambaforge发行版。使用python>=3.8，conda考虑使用mamba代替conda。

2. 创建一个新的conda环境：

   conda create -n omicverse python=3.8
   

3. 激活您的环境：

   conda activate omicverse
   

4. 首先安装PyTorch：

   conda install pytorch torchvision torchaudio pytorch-cuda=11.7 -c pytorch -c nvidia
   

5. 安装omicverse：

   conda install omicverse -c conda-forge
   

Pip

可以使用以下命令之一通过pip安装omicverse包：

1. 首先安装PyTorch：有关安装的更多信息可以在PyTorch找到。

   # ROCM 5.2（仅限Linux）
   pip3 install torch torchvision torchaudio --extra-index-url
   pip install torch==1.13.1+rocm5.2 torchvision==0.14.1+rocm5.2 torchaudio==0.13.1 --extra-index-url https://download.pytorch.org/whl/rocm5.2
   # CUDA 11.6
   pip install torch==1.13.1+cu116 torchvision==0.14.1+cu116 torchaudio==0.13.1 --extra-index-url https://download.pytorch.org/whl/cu116
   # CUDA 11.7
   pip install torch==1.13.1+cu117 torchvision==0.14.1+cu117 torchaudio==0.13.1 --extra-index-url https://download.pytorch.org/whl/cu117
   # 仅限CPU
   pip install torch==1.13.1+cpu torchvision==0.14.1+cpu torchaudio==0.13.1 --extra-index-url https://download.pytorch.org/whl/cpu
   

2. 在安装pytorch之后，我们可以通过pip开始安装omicverse

   pip install -U omicverse
   pip install -U numba
   

3. 如果您想使用测试版，请有两种方法供您安装

• 测试版 - 克隆此repo并运行：pip install .
• 使用pip install git+https://github.com/Starlitnightly/omicverse.git

开发者版本

用于开发者测试的版本 - 克隆此repo并运行：

pip install -e ".[dev,docs]"

Python版RNA-seq分析教程：差异表达基因分析

Bulk RNA-seq 分析的一个重要任务是分析差异表达基因，我们可以用 omicverse包 来完成这个任务。对于差异表达分析而言，首先，我们可以先将 gene_id 改为 gene_name。其次，当我们的数据集存在批量效应时，我们可以使用 DEseq2的 SizeFactor 对其进行归一化，从统计学上，使用 wilcoxon的秩和检验或者 t-test来计算基因的 p 值。也可以使用类似edgeR，Deseq2等包的模型来计算p值。在这里，我们用一个从RNA-seq上游的定量包FeatureCounts生成的表达矩阵来演示差异表达分析的流程。我们的流程适用于任何Bulk RNA-seq的差异表达分析。

环境的下载

在这里我们只需要安装omicverse环境即可，有两个方法：

• 一个是使用conda：conda install omicverse -c conda-forge
• 另一个是使用pip：pip install omicverse -i https://pypi.tuna.tsinghua.edu.cn/simple/。-i的意思是指定清华镜像源，在国内可能会下载地快一些。

导入包

我们首先导入分析需要用到的所有包，包括omicverse, pandas, numpy, scanpy matplotlib   和 seaborn.

import omicverse as ov
import pandas as pd
import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt
import seaborn as sns

#设定绘图格式，分辨率300dpi等
ov.utils.ov_plot_set()

下载基因集

当我们需要转换基因 id 时，我们需要准备一个映射对文件。在这里，我们预处理了6个基因组 gtf 文件和生成的映射对，包括 T2T-CHM13，GRCh38，GRCh37，GRCm39，danRer7和 danRer11。如果需要转换其他 id，可以使用 gtf 将文件放在 genesets 目录中生成自己的映射。

ov.utils.download_geneid_annotation_pair()

读取数据

data=pd.read_csv('https://raw.githubusercontent.com/Starlitnightly/ov/master/sample/counts.txt',index_col=0,sep='\t',header=1)
#replace the columns `.bam` to `` 
data.columns=[i.split('/')[-1].replace('.bam','') for i in data.columns]
data.head()

值得注意的是，我们的数据集并没有经过任何处理，featurecounts比对时用的gtf为GRCm39，所以我们这里用GRCm39来做基因id映射

基因id转换

data=ov.bulk.Matrix_ID_mapping(data,'genesets/pair_GRCm39.tsv')
data.head()

差异表达分析

我们可以非常简单地通过omicverse进行差异表达基因分析，只需要提供一个表达式矩阵。我们首先创建一个 pyDEG 对象，并使用drop_duplicates_index去除重复的基因。由于部分基因名相同，我们的去除保留了表达量最大的基因名。

    
dds=ov.bulk.pyDEG(data)
dds.drop_duplicates_index()
print('... drop_duplicates_index success')

我们还需要去除表达矩阵的批次效应 (batch effect)，我们使用DEseq2的的 SizeFactor 来对我们的矩阵计算归一化因子来去除批次效应。

dds.normalize()
print('... estimateSizeFactors and normalize success')

现在我们可以从表达矩阵中计算差异表达基因，在计算前我们需要输入实验组和对照组。在这里，我们指定 4-3和4-4为实验组，1--1, 1--2为对照组，使用ttest进行差异表达分析计算。当然你也可以使用wilcox来计算。此外deseq2也是支持的，不过流程可能会有一些区别，我们放到下一期讲。

treatment_groups=['4-3','4-4']
control_groups=['1--1','1--2']
result=dds.deg_analysis(treatment_groups,control_groups,method='ttest')
result.head()

在计算完差异表达基因后，我们会发现一个重要的事情，就是低表达基因有很多，如果我们不对其进行过滤，会影响后续火山图的绘制，我们设定基因的平均表达量大于1作为阈值，将平均表达量低于1的基因全部过滤掉。

print(result.shape)
result=result.loc[result['log2(BaseMean)']>1]
print(result.shape)

我们还需要设置 Foldchange 的阈值，我们准备了一个名为 foldchange_set 的方法函数来完成。此函数根据 log2FC 分布自动计算适当的阈值，但您也可以手动输入阈值。该函数有三个参数：

• fc_threshold: 差异表达倍数的阈值，-1为自动计算
• pval_threshold: 差异表达基因的p-value过滤值，默认为0.05，在有些情况下可以设定为0.1，意味着统计学差异不显著。
• logp_max: p值的最大值，由于部分p值过小，甚至为0，取对数后火山图绘制较为困难，我们可以设定一个上限，高于这个上限的p值全部统一。

# -1 means automatically calculates
dds.foldchange_set(fc_threshold=-1,
                   pval_threshold=0.05,
                   logp_max=6)

差异表达的结果可视化

omicverse除了有较为完善的分析能力外，还有极强的可视化能力。首先是火山图，我们使用  plot_volcano函数来实现。该函数可以绘制你感兴趣的基因或高表达的基因。您需要输入一些参数:

• title: 火山图的标题
• figsize: 图像大小
• plot_genes: 需要绘制的基因，格式为list。如['Gm8925','Snorc']
• plot_genes_num: 需要绘制的基因数，该参数与plot_genes互斥，如果我们没有指定需要绘制的基因，可以自动绘制前n个高差异表达倍数的基因。

此外，我们还可以指定绘制的颜色等，具体的参数可以使用help(dds.plot_volcano)来查看

dds.plot_volcano(title='DEG Analysis',figsize=(4,4),
                 plot_genes_num=8,plot_genes_fontsize=12,)

如果我们想绘制特定的基因的箱线图，我们也可以使用 plot_boxplot 函数来完成该任务。

dds.plot_boxplot(genes=['Ckap2','Lef1'],treatment_groups=treatment_groups,
                control_groups=control_groups,figsize=(2,3),fontsize=12,
                 legend_bbox=(2,0.55))

通路富集分析

在差异表达基因计算出来后，我们需要直接进行的下一步分析往往是看差异表达的基因与哪些通路相关，这里我们常用的方法是富集分析。omicverse可以一键完成富集分析并且可视化。

我们封装了gseapy 包进入omicverse，其中包括 GSEA 富集分析的相关功能。我们优化了包的输出，并给出了一些更好看的图形绘制功能

类似地，我们首先需要下载通路数据库。我们已经准备好了五个基因集，可以使用 ov.utils.download_pathway_database()进行自动下载。除此之外，你还可以在以下网站找到你感兴趣的基因集：https://maayanlab.cloud/enrichr/#libraries

ov.utils.download_pathway_database()
#读取通路基因集，我们读取Wiki通路数据库
pathway_dict=ov.utils.geneset_prepare('genesets/WikiPathways_2019_Mouse.txt',organism='Mouse')

我们提取前面的差异表达基因来进行通路富集

#差异表达基因提取
deg_genes=dds.result.loc[dds.result['sig']!='normal'].index.tolist()
#通路富集分析
enr=ov.bulk.geneset_enrichment(gene_list=deg_genes,
                                pathways_dict=pathway_dict,
                                pvalue_type='auto',
                                organism='mouse')

我们可以使用geneset_plot来可视化通路富集的结果

ov.bulk.geneset_plot(enr,figsize=(2,5),fig_title='Wiki Pathway enrichment',
                        cmap='Reds')
#如果需要保存的话,使用`plt.savefig`来保存图像
plt.savefig("enr_pathway.png",dpi=300,bbox_inches = 'tight')