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

文章来源: drugai.blog.csdn.net，作者：DrugAI，版权归原作者所有，如需转载，请联系作者。

原文链接：drugai.blog.csdn.net/article/details/88877273