Closed ixxmu closed 1 year ago
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
进行配置。conda install -c anaconda pip
并跳过此部分。在装有Apple Silicon的Mac上安装omicverse只能使用本机版本的Python。可以通过使用Apple Silicon版本的mambaforge(可以通过本机版本的homebrew通过brew install --cask mambaforge
进行安装)来安装本机版本的Python。
安装conda。我们通常使用mambaforge
发行版。使用python>=3.8,conda考虑使用mamba代替conda。
创建一个新的conda环境:
conda create -n omicverse python=3.8
激活您的环境:
conda activate omicverse
首先安装PyTorch:
conda install pytorch torchvision torchaudio pytorch-cuda=11.7 -c pytorch -c nvidia
安装omicverse
:
conda install omicverse -c conda-forge
可以使用以下命令之一通过pip安装omicverse
包:
首先安装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
在安装pytorch之后,我们可以通过pip
开始安装omicverse
pip install -U omicverse
pip install -U numba
如果您想使用测试版,请有两种方法供您安装
pip install .
pip install git+https://github.com/Starlitnightly/omicverse.git
用于开发者测试的版本 - 克隆此repo并运行:
pip install -e ".[dev,docs]"
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 install omicverse -c conda-forge
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映射
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 分布自动计算适当的阈值,但您也可以手动输入阈值。该函数有三个参数:
# -1 means automatically calculates
dds.foldchange_set(fc_threshold=-1,
pval_threshold=0.05,
logp_max=6)
omicverse除了有较为完善的分析能力外,还有极强的可视化能力。首先是火山图,我们使用 plot_volcano
函数来实现。该函数可以绘制你感兴趣的基因或高表达的基因。您需要输入一些参数:
['Gm8925','Snorc']
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')
因为前面我们在<生信技能树>公众号宣传过一波; Python的转录组学分析框架与生态,因为是需要去github点star后发邮件才能进群交流,所以操作门槛有点高, 我们后续再次开放拉群小助手给大家哈。
因为交流群早就满200人,所以没办法通过二维码进群了哈。这个时候需要我们生信技能树的官方拉群小助手帮忙拉群哦!!!
(免费,但是名额有限,500人,先到先得!!!另外,因为每次人数太多, 所以是每天上午十点准时拉群,其他时间不予回复,望见谅)
长按识别二维码
烦请备注姓名学校单位信息以及Python
再造一遍轮子有何意义,缝合怪还不如R来的干脆
https://mp.weixin.qq.com/s/fIRVMNHFcffCuCzh_Qozhw