Differentially expressed gene analysis
An important task of bulk rna-seq analysis is the different expression , which we can perform with omicverse. For different expression analysis, ov change the
gene_name of matrix first. When our dataset existed the batch effect, we can use the SizeFactors of DEseq2 to normalize it, and use
wilcoxon to calculate the p-value of genes. Here we demonstrate this pipeline with a matrix from
featureCounts. The same pipeline would generally be used to analyze any collection of RNA-seq tasks.
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 ov.utils.ov_plot_set()
When we need to convert a gene id, we need to prepare a mapping pair file. Here we have pre-processed 6 genome gtf files and generated mapping pairs including
danRer11. If you need to convert other id_mapping, you can generate your own mapping using gtf Place the files in the
Note that this dataset has not been processed in any way and is only exported by
featureCounts, and Sequence alignment was performed from the genome file of CRCm39
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()
We performed the gene_id mapping by the mapping pair file
GRCm39 downloaded before.
Different expression analysis with ov
We can do differential expression analysis very simply by ov, simply by providing an expression matrix. To run DEG, we simply need to:
- Read the raw count by featureCount or any other qualify methods.
- Create an ov DEseq object.
We notes that the gene_name mapping before exist some duplicates, we will process the duplicate indexes to retain only the highest expressed genes
dds.drop_duplicates_index() print('... drop_duplicates_index success')
We also need to remove the batch effect of the expression matrix,
estimateSizeFactors of DEseq2 to be used to normalize our matrix
dds.normalize() print('... estimateSizeFactors and normalize success')
Now we can calculate the different expression gene from matrix, we need to input the treatment and control groups
treatment_groups=['4-3','4-4'] control_groups=['1--1','1--2'] result=dds.deg_analysis(treatment_groups,control_groups,method='ttest') result.head()
One important thing is that we do not filter out low expression genes when processing DEGs, and in future versions I will consider building in the corresponding processing.
print(result.shape) result=result.loc[result['log2(BaseMean)']>1] print(result.shape)
We also need to set the threshold of Foldchange, we prepare a method named
foldchange_set to finish. This function automatically calculates the appropriate threshold based on the log2FC distribution, but you can also enter it manually.
# -1 means automatically calculates dds.foldchange_set(fc_threshold=-1, pval_threshold=0.05, logp_max=6)
Visualize the DEG result and specific genes
To visualize the DEG result, we use
plot_volcano to do it. This fuction can visualize the gene interested or high different expression genes. There are some parameters you need to input:
- title: The title of volcano
- figsize: The size of figure
- plot_genes: The genes you interested
- plot_genes_num: If you don't have interested genes, you can auto plot it.
dds.plot_volcano(title='DEG Analysis',figsize=(4,4), plot_genes_num=8,plot_genes_fontsize=12,)
To visualize the specific genes, we only need to use the
dds.plot_boxplot function to finish it.
dds.plot_boxplot(genes=['Ckap2','Lef1'],treatment_groups=treatment_groups, control_groups=control_groups,figsize=(2,3),fontsize=12, legend_bbox=(2,0.55))
Pathway enrichment analysis by ov
Here we use the
gseapy package, which included the GSEA analysis and Enrichment. We have optimised the output of the package and given some better looking graph drawing functions
Similarly, we need to download the pathway/genesets first. Five genesets we prepare previously, you can use
ov.utils.download_pathway_database() to download automatically. Besides, you can download the pathway you interested from enrichr: https://maayanlab.cloud/Enrichr/#libraries
Note that the
pvalue_type we set to
auto, this is because when the genesets we enrichment if too small, use the
adjusted pvalue we can't get the correct result. So you can set
raw to get the significant geneset.
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')
To visualize the enrichment, we use
geneset_plot to finish it
ov.bulk.geneset_plot(enr,figsize=(2,5),fig_title='Wiki Pathway enrichment', cmap='Reds')