Bulk RNA-seq 分析的一個重要任務(wù)是分析差異表達基因党饮，我們可以用 omicverse包
來完成這個任務(wù)局雄。在omicverse中评甜，除了最簡單的ttest外，在這里隶校，我們介紹一種類似R語言中的Deseq2等包的模型來計算差異表達基因漏益。
環(huán)境的下載
在這里我們只需要安裝omicverse環(huán)境即可，有兩個方法：
- 一個是使用conda：
conda install omicverse -c conda-forge - 另一個是使用pip：
pip install omicverse -i https://pypi.tuna.tsinghua.edu.cn/simple/惠况。-i的意思是指定清華鏡像源遭庶，在國內(nèi)可能會下載地快一些。
詳細的安裝教程見omicverse官網(wǎng)稠屠。
導(dǎo)入包
我們首先導(dǎo)入分析需要用到的所有包峦睡，包括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
#設(shè)定繪圖格式，分辨率300dpi等
ov.utils.ov_plot_set()
下載基因集
當我們需要轉(zhuǎn)換基因 id 時权埠，我們需要準備一個映射對文件榨了。在這里，我們預(yù)處理了6個基因組 gtf 文件和生成的映射對攘蔽，包括 T2T-CHM13龙屉，GRCh38，GRCh37满俗，GRCm39转捕，danRer7和 danRer11。如果需要轉(zhuǎn)換其他 id唆垃，可以使用 gtf 將文件放在 genesets 目錄中生成自己的映射五芝。
ov.utils.download_geneid_annotation_pair()
讀取數(shù)據(jù)
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()
值得注意的是，我們的數(shù)據(jù)集并沒有經(jīng)過任何處理辕万，featurecounts比對時用的gtf為GRCm39枢步，所以我們這里用GRCm39來做基因id映射
基因id轉(zhuǎn)換
data=ov.bulk.Matrix_ID_mapping(data,'genesets/pair_GRCm39.tsv')
data.head()
差異表達分析
我們可以非常簡單地通過omicverse進行差異表達基因分析，只需要提供一個表達式矩陣渐尿。我們首先創(chuàng)建一個 pyDEG 對象醉途，并使用drop_duplicates_index去除重復(fù)的基因。由于部分基因名相同砖茸，我們的去除保留了表達量最大的基因名隘擎。
dds=ov.bulk.pyDEG(data)
dds.drop_duplicates_index()
print('... drop_duplicates_index success')
現(xiàn)在我們可以從表達矩陣中計算差異表達基因，在計算前我們需要輸入實驗組和對照組凉夯。在這里嵌屎，我們指定 4-3和4-4為實驗組推正，1--1, 1--2為對照組恍涂，我們設(shè)定method為DEseq2也是支持的宝惰，不過流程可能會有一些區(qū)別，我們放到下一期講再沧。
treatment_groups=['4-3','4-4']
control_groups=['1--1','1--2']
result=dds.deg_analysis(treatment_groups,control_groups,
method='DEseq2')
result.head()
Fitting size factors...
... done in 0.00 seconds.
Fitting dispersions...
... done in 1.59 seconds.
Fitting dispersion trend curve...
... done in 2.82 seconds.
logres_prior=1.1538905878789707, sigma_prior=0.25
Fitting MAP dispersions...
... done in 1.57 seconds.
Fitting LFCs...
... done in 1.27 seconds.
Refitting 0 outliers.
Running Wald tests...
... done in 1.33 seconds.
Log2 fold change & Wald test p-value: condition Treatment vs Control
在計算完差異表達基因后尼夺，我們會發(fā)現(xiàn)一個重要的事情，就是低表達基因有很多炒瘸，如果我們不對其進行過濾淤堵，會影響后續(xù)火山圖的繪制，我們設(shè)定基因的平均表達量大于1作為閾值顷扩，將平均表達量低于1的基因全部過濾掉拐邪。
print(result.shape)
result=result.loc[result['log2(BaseMean)']>1]
print(result.shape)
我們還需要設(shè)置 Foldchange 的閾值，我們準備了一個名為 foldchange_set 的方法函數(shù)來完成隘截。此函數(shù)根據(jù) log2FC 分布自動計算適當?shù)拈撝翟祝部梢允謩虞斎腴撝怠Ｔ摵瘮?shù)有三個參數(shù)：
- fc_threshold: 差異表達倍數(shù)的閾值婶芭，-1為自動計算
- pval_threshold: 差異表達基因的p-value過濾值东臀，默認為0.05，在有些情況下可以設(shè)定為0.1犀农，意味著統(tǒng)計學差異不顯著惰赋。
- logp_max: p值的最大值，由于部分p值過小呵哨，甚至為0赁濒，取對數(shù)后火山圖繪制較為困難，我們可以設(shè)定一個上限孟害，高于這個上限的p值全部統(tǒng)一拒炎。
# -1 means automatically calculates
dds.foldchange_set(fc_threshold=-1,
pval_threshold=0.05,
logp_max=6)
差異表達的結(jié)果可視化
omicverse除了有較為完善的分析能力外，還有極強的可視化能力纹坐。首先是火山圖枝冀，我們使用 plot_volcano函數(shù)來實現(xiàn)。該函數(shù)可以繪制你感興趣的基因或高表達的基因耘子。您需要輸入一些參數(shù):
- title: 火山圖的標題
- figsize: 圖像大小
- plot_genes: 需要繪制的基因果漾，格式為list。如
['Gm8925','Snorc'] - plot_genes_num: 需要繪制的基因數(shù)谷誓，該參數(shù)與
plot_genes互斥绒障，如果我們沒有指定需要繪制的基因，可以自動繪制前n個高差異表達倍數(shù)的基因捍歪。
此外户辱，我們還可以指定繪制的顏色等鸵钝，具體的參數(shù)可以使用help(dds.plot_volcano)來查看
dds.plot_volcano(title='DEG Analysis',figsize=(4,4),
plot_genes_num=8,plot_genes_fontsize=12,)
如果我們想繪制特定的基因的箱線圖，我們也可以使用 plot_boxplot 函數(shù)來完成該任務(wù)庐镐。
dds.plot_boxplot(genes=['Ckap2','Lef1'],treatment_groups=treatment_groups,
control_groups=control_groups,figsize=(2,3),fontsize=12,
legend_bbox=(2,0.55))
通路富集分析
在差異表達基因計算出來后恩商，我們需要直接進行的下一步分析往往是看差異表達的基因與哪些通路相關(guān)，這里我們常用的方法是富集分析必逆。omicverse可以一鍵完成富集分析并且可視化怠堪。
我們封裝了gseapy 包進入omicverse，其中包括 GSEA 富集分析的相關(guān)功能名眉。我們優(yōu)化了包的輸出粟矿，并給出了一些更好看的圖形繪制功能
類似地，我們首先需要下載通路數(shù)據(jù)庫损拢。我們已經(jīng)準備好了五個基因集陌粹，可以使用 ov.utils.download_pathway_database()進行自動下載。除此之外福压，你還可以在以下網(wǎng)站找到你感興趣的基因集： https://maayanlab.cloud/enrichr/#libraries
ov.utils.download_pathway_database()
#讀取通路基因集掏秩，我們讀取Wiki通路數(shù)據(jù)庫
pathway_dict=ov.utils.geneset_prepare('genesets/WikiPathways_2019_Mouse.txt',organism='Mouse')
與此前選取差異表達的基因進行通路富集分析不同，在這里隧膏，我們使用基因的差異倍數(shù)和p值進行排序哗讥，進行GSEA分析。我們使用以下公式計算基因的順序胞枕。
rnk=dds.ranking2gsea()
們使用 ov.bulk.pyGSEA 構(gòu)造一個 GSEA 對象來執(zhí)行通路富集分析杆煞。
gsea_obj=ov.bulk.pyGSEA(rnk,pathway_dict)
enrich_res=gsea_obj.enrichment()
通路富集分析的結(jié)果被存放在enrich_res變量中，我們使用.head()函數(shù)來查看
gsea_obj.enrich_res.head()
為了可視化通路富集的結(jié)果腐泻，我們使用.plot_enrichment函數(shù)來完成該目的.
- num: 需要展示的通路富集結(jié)果决乎，默認是前10條.
- node_size: 點的大小標注值. 默認是[5,10,15].
- cax_loc: 圖注的橫坐標. 默認值 2.
- cax_fontsize: 圖注的字體大小，默認值 12.
- fig_title: 富集結(jié)果的標題.
- fig_xlabel: 富集結(jié)果的橫坐標標題派桩，默認值是'Fractions of genes'.
- figsize: 圖像大小构诚，默認值是(2,4).
- cmap: 圖像顏色條，默認值是'YlGnBu'.
gsea_obj.plot_enrichment(num=10,node_size=[10,20,30],
cax_loc=2.5,cax_fontsize=12,
fig_title='Wiki Pathway Enrichment',fig_xlabel='Fractions of genes',
figsize=(2,4),cmap='YlGnBu',
text_knock=2,text_maxsize=30)
不僅僅是基本的富集分析铆惑，pyGSEA 還可以幫助我們用繪制GSEA的排序圖范嘱，我們需要選擇想繪制的通路，給定通路的Term位置员魏，例如0 是Complement and Coagulation Cascades WP449 丑蛤，1 是 Matrix Metalloproteinases WP441
gsea_obj.enrich_res.index[:5]
Index(['Complement and Coagulation Cascades WP449',
'Matrix Metalloproteinases WP441', 'TYROBP Causal Network WP3625',
'PPAR signaling pathway WP2316',
'Metapathway biotransformation WP1251'],
dtype='object', name='Term')
我們可以使用.plot_gsea完成繪制。
- term_num: 需要繪制的通路在index的位置
- gene_set_title: 通路標題撕阎，可自定義受裹，如果原index中的太長
fig=gsea_obj.plot_gsea(term_num=1,
gene_set_title='Matrix Metalloproteinases',
figsize=(3,4),
cmap='RdBu_r',
title_fontsize=14,
title_y=0.95)
