【生物信息学】基因差异分析Deg（数据读取、数据处理、差异分析、结果可视化）

目录

一、实验介绍

二、实验环境

1. 配置虚拟环境

2. 库版本介绍

3. IDE

三、实验内容

0. 导入必要的工具包

1. 定义一些阈值和参数

2. 读取数据

normal_data.csv部分展示

tumor_data.csv部分展示

3. 绘制箱型图

4. 删除表达量低于阈值的基因

5. 计算差异显著的基因

6. 初始化结果数据表格

7. 进行秩和检验和差异倍数计算

8. 可视化分析

9. 代码整合

---

一、实验介绍

        本实验完成了基因差异分析，包括数据读取、数据处理（ 绘制箱型图、删除表达量低于阈值的基因、计算差异显著的基因）、差异分析（进行秩和检验和差异倍数计算）等，成功识别出在正常样本与肿瘤样本之间显著表达差异的基因，并对其进行了进一步的可视化分析（箱型图、差异倍数fold分布图、热力图和散点图）。

        基因差异分析是研究不同条件下基因表达差异的重要手段，能够帮助我们理解生物体内基因调控的变化及其与表型特征的关联。本实验旨在探索正常样本与肿瘤样本之间基因表达的差异，并识别差异显著的基因。

二、实验环境

    本系列实验使用了PyTorch深度学习框架，相关操作如下（基于深度学习系列文章的环境）：

1. 配置虚拟环境

深度学习系列文章的环境

conda create -n DL python=3.7 

conda activate DL

pip install torch==1.8.1+cu102 torchvision==0.9.1+cu102 torchaudio==0.8.1 -f https://download.pytorch.org/whl/torch_stable.html

conda install matplotlib

conda install scikit-learn

新增加

conda install pandas

conda install seaborn

conda install networkx

conda install statsmodels

pip install pyHSICLasso

注：本人的实验环境按照上述顺序安装各种库，若想尝试一起安装（天知道会不会出问题）

2. 库版本介绍

新增

3. IDE

        建议使用Pycharm（其中，pyHSICLasso库在VScode出错，尚未找到解决办法……）

win11 安装 Anaconda（2022.10）+pycharm（2022.3/2023.1.4）+配置虚拟环境_QomolangmaH的博客-CSDN博客https://blog.csdn.net/m0_63834988/article/details/128693741

三、实验内容

0. 导入必要的工具包

import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import zscore
from scipy.stats import ranksums
from pyHSICLasso import HSICLasso
from sklearn.preprocessing import LabelEncoder

1. 定义一些阈值和参数

p_cutoff = 0.001
FC_cutoff = 3
num_cutoff = 10

• p_cutoff 是判断差异显著性的 p 值阈值。
• FC_cutoff 是判断差异的倍数阈值。
• num_cutoff 是基因表达量筛选的阈值。

2. 读取数据

ndata = pd.read_csv("normal_data.csv", index_col=0, header=0)
tdata = pd.read_csv("tumor_data.csv", index_col=0, header=0)

        从文件 "normal_data.csv" 和 "tumor_data.csv" 中读取正常样本和肿瘤样本的基因表达数据。

normal_data.csv部分展示

tumor_data.csv部分展示

3. 绘制箱型图

ndata.iloc[1000:1010, ].transpose().plot(kind='box', title='Normal Sample Gene Boxplot', rot=30)
plt.tick_params(labelsize=10)
plt.savefig('Normal_Sample_Gene_box_plot', bbox_inches='tight')
plt.show()

        绘制正常样本中部分基因的箱型图，并保存为图片文件。

4. 删除表达量低于阈值的基因

ndata = ndata.iloc[29:, :]
tdata = tdata.iloc[29:, :]

        删除正常样本和肿瘤样本中表达量低于阈值的基因。

5. 计算差异显著的基因

gene_sig = ndata[ndata.mean(axis=1) > num_cutoff].index.intersection(tdata[tdata.mean(axis=1) > num_cutoff].index)
Ndata = ndata.loc[gene_sig]
Tdata = tdata.loc[gene_sig]

        选择正常样本和肿瘤样本中表达量高于阈值的基因，并保存为新的数据。

6. 初始化结果数据表格

p_value = [1.] * len(Ndata)
log2FoldChange = []
label = ['0'] * len(Ndata)
result = pd.DataFrame([p_value, log2FoldChange, label])
result = result.transpose()
result.columns = ['p_value', 'log2FC', 'label']
result.index = Ndata.index.tolist()
p = []

        创建一个结果数据表格，包含 p 值、log2 fold change 和标签，并初始化为默认值。

7. 进行秩和检验和差异倍数计算

for i in Ndata.index:
    p1 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='greater')[1]
    p2 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='less')[1]
    p3 = ranksums(Tdata.loc[i], Ndata.loc[i])[1]
    result.loc[i, 'log2FC'] = np.log2(np.average(Tdata.loc[i]) / np.average(Ndata.loc[i]))
    p.append(p3)
    if (p1 <= p_cutoff):
       result.loc[i, 'p_value'] = p1
       if result.loc[i, 'log2FC'] > np.log2(FC_cutoff):
           result.loc[i, 'label'] = '1'
 
    if (p2 <= p_cutoff):
       result.loc[i, 'p_value'] = p2
       if result.loc[i, 'log2FC'] < -np.log2(FC_cutoff):
           result.loc[i, 'label'] = '-1'

        对每个基因进行秩和检验，计算 p 值和差异倍数，并根据设定的阈值确定差异显著的基因，并给予标签。

8. 可视化分析

print('finished')
plt.hist(result['log2FC'], bins=10, color='blue', alpha=0.6, edgecolor='black')
plt.title("fold change")
plt.show()

9. 代码整合

import numpy as np
import pandas as pd
import scipy
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import zscore
from scipy.stats import ranksums
from pyHSICLasso import HSICLasso
from sklearn.preprocessing import LabelEncoder
 
 
p_cutoff = 0.001
FC_cutoff = 3
num_cutoff = 10
 
# 读取数据ndata normal data,tdata tumor data
ndata = pd.read_csv("normal_data.csv", index_col=0, header=0)
tdata = pd.read_csv("tumor_data.csv", index_col=0, header=0)
 
# 箱型图
# 中位数
# QL称为下四分位数，表示全部观察值中有四分之一的数据取值比它小；
# QU称为上四分位数，表示全部观察值中有四分之一的数据取值比它大；
# IQR称为四分位数间距，是上四分位数QU与下四分位数QL之差，其间包含了全部观察值的一半。
# 异常值通常被定义为小于QL－1.5IQR或大于QU＋1.5IQR的值。
ndata.iloc[1000:1010, ].transpose().plot(kind='box', title='Normal Sample Gene Boxplot', rot=30)
plt.tick_params(labelsize=10)
plt.savefig('Normal_Sample_Gene_box_plot', bbox_inches='tight')
plt.show()
 
# delete gene expression < num
ndata = ndata.iloc[29:, :]
tdata = tdata.iloc[29:, :]
gene_sig = ndata[ndata.mean(axis=1) > num_cutoff].index.intersection(tdata[tdata.mean(axis=1) > num_cutoff].index)
Ndata = ndata.loc[gene_sig]
Tdata = tdata.loc[gene_sig]
p_value = [1.] * len(Ndata)
# log2 fold change的意思是取log2，这样可以可以让差异特别大的和差异比较小的数值缩小之间的差距
log2FoldChange = []
label = ['0'] * len(Ndata)
 
result = pd.DataFrame([p_value, log2FoldChange, label])
result = result.transpose()
result.columns = ['p_value', 'log2FC', 'label']
result.index = Ndata.index.tolist()
p = []
 
# 秩和检验(也叫Mann-Whitney U-test)，检验两组数据是否来自具有相同中位数的连续分布，检验它们是否有显著差异。
# ’less‘ 基于 x 的分布小于基于 y 的分布
# ‘greater’ x 的分布大于 y 的分布。
 
for i in Ndata.index:
    p1 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='greater')[1]
    p2 = ranksums(Tdata.loc[i], Ndata.loc[i], alternative='less')[1]
    p3 = ranksums(Tdata.loc[i], Ndata.loc[i])[1]
    result.loc[i, 'log2FC'] = np.log2(np.average(Tdata.loc[i]) / np.average(Ndata.loc[i]))
    p.append(p3)
    if (p1 <= p_cutoff):
       result.loc[i, 'p_value'] = p1
       if result.loc[i, 'log2FC'] > np.log2(FC_cutoff):
           result.loc[i, 'label'] = '1'
 
    if (p2 <= p_cutoff):
       result.loc[i, 'p_value'] = p2
       if result.loc[i, 'log2FC'] < -np.log2(FC_cutoff):
           result.loc[i, 'label'] = '-1'
 
print('finished')
# 查看差异倍数fold分布
plt.hist(result['log2FC'], bins=10, color='blue', alpha=0.6, edgecolor='black')
plt.title("fold change")
plt.show()
 
 
# result.to_csv('result.csv')
up_data_n = Ndata.loc[result.query("(label == '1')").index]
up_data_t = Tdata.loc[result.query("(label == '1')").index]
down_data_n = Ndata.loc[result.query("(label == '-1')").index]
down_data_t = Tdata.loc[result.query("(label == '-1')").index]
data = pd.concat([pd.concat([up_data_n, up_data_t], axis=1), pd.concat([down_data_n, down_data_t], axis=1)], axis=0)
deg_gene = pd.DataFrame(data.index.tolist())
deg_gene.to_csv('deg_gene.csv')
z1 = zscore(np.log2(data+1), axis=0)
 
sns.heatmap(data=z1, cmap="coolwarm", xticklabels=False)
plt.savefig('heatmap_plot', bbox_inches='tight')
plt.show()
 
 
p = -np.log10(p)
ax = sns.scatterplot(x="log2FC", y=p,
                     hue='label',
                     hue_order=('-1', '0', '1'),
                     palette=("#377EB8","grey","#E41A1C"),
                     data=result)
ax.set_ylabel('-log(pvalue)', fontweight='bold')
ax.set_xlabel('FoldChange', fontweight='bold')
 
plt.savefig('deg_p_fc_plot', bbox_inches='tight')
plt.show()