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

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

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 = pd.read_csv("normal_data.csv", index_col=0, header=0)
tdata = pd.read_csv("tumor_data.csv", index_col=0, header=0)

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()

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)
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 = []

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')
plt.hist(result['log2FC'], bins=10, color='blue', alpha=0.6, edgecolor='black')
plt.title("fold change")
plt.show()

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()