Python生物信息学③提取差异基因

python做生信分析的流程

使用的数据集是GSE5583，来自于2006年的基因芯片结果，该芯片目的是提取野生型和HDAC1小鼠胚胎干细胞用于Affymetrix微阵列上的差异RNA。

#导入包
import matplotlib.pyplot as plt
import os
import numpy as np
import pandas as pd
from scipy import stats
import seaborn as sns
%matplotlib inline

#载入数据
data = pd.read_table("GSE5583.txt",header = 0,index_col = 0)
data.head()  #查看前5行

每一行是一个基因，每一列是一个样本，这也是比较经典的芯片数据集

#查看数据维度
data.shape

标准化

常见的log2（）标准化

data2 = np.log2(data+0.0001)
data2.head()

# 每个阵列的箱线图
plt.show(data2.plot(kind = 'box', title = 'GSE5583 Boxplot', rot = 90))

目的是查看不同样本之间是否有总体差异。

# Density
plt.show(data2.plot(kind = 'density', title = 'GSE5583 Density'))

可以看出样本之间没有总体差异，可以做差异分析。

#每个基因（行）wt样本的表达平均值
wt = data2.loc[:, 'WT.GSM130365' : 'WT.GSM130367'].mean(axis = 1)
wt.head()

#每个基因（行）的ko样本的表达平均值
ko = data2.loc[:,'KO.GSM130368':'KO.GSM130370'].mean(axis = 1)
ko.head()

fold = ko - wt

#折叠变化的直方图
plt.hist(fold)
plt.title("Histogram of fold-change")
plt.show()

查看基因差异的P值分布

from scipy import stats

pvalue = []
for i in range(0, number_of_genes):
    ttest = stats.ttest_ind(data2.iloc[i,0:3], data2.iloc[i,3:6])
    pvalue.append(ttest[1])

# Histogram of the p-values
plt.hist(-np.log(pvalue))
plt.title("Histogram of p-value")
plt.show()

---

参考：

https://www.jianshu.com/p/91c98585b79b