python做生物信息学分析_Python从零开始第五章生物信息学①提取差异基因

目前来说，做生物信息学的人越来越多，但是我觉得目前而言做生信的主要有三类人：

老本行是做实验的，做生信可能是为了辅助研究或者是为了发paper(有非常多的临床生选择趟生信这波水)

主要是做生信的，主要涵盖高通量测序数据分析，组学数据分析等等，专门从事生物学数据分析的这群人，其大部分也是本科生物狗作为强大的生力军，以调包写R，python为主。那么这群人就要熟悉看各种包的tutorial以及如何进行常规的数据的处理和分析等等。那么其实这群人，我个人认为对python的编程要求不是很高，在coursera先上两门课程，然后照着参考书在项目中练练就熟络了。

也是做生信的，但是主要以写包为主，什么意思？就是以开发算法供给第二类人用，做这部分工作的人需要具有良好的数理基础，所以一般大部分都是本科属于物理，数学等理工科的人在做这部分工作。当然对编程也提出较高要求。

我入门生物信息学是通过R语言入门的，但是接触到了python，这个也是目前用户数量数一数二的语言。python去做生信得优点是①过程更加直观，因为常见的R包功能一般已经封装好了，直接应用就可，虽然足够简单友好，但是不利于长期学习②基因组数据一般比较大，python速度一般比R快。

下面我就记录一个通过python做生信分析的流程。

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

导入包

%clear

%reset -f

# In[*]

import seaborn as sns

import matplotlib.pyplot as plt

%matplotlib inline

import os

import numpy as np

import pandas as pd

os.chdir('D:\\train')

from scipy import stats

导入数据

# In[*]

data = pd.read_table("http://ahmedmoustafa.io/notebooks/GSE5583.txt",

header = 0,

index_col = 0)

data.head()

Out[27]:

WT.GSM130365 WT.GSM130366 ... KO.GSM130369 KO.GSM130370

100001_at 11.5 5.6 ... 36.0 42.0

100002_at 20.5 32.4 ... 14.4 22.9

100003_at 72.4 89.0 ... 130.1 86.7

100004_at 261.0 226.2 ... 447.3 288.1

100005_at 1086.2 1555.6 ... 1365.9 1436.2

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

查看数据维度

# In[*]

# Check the dimension of the loaded data (rows & columns)

data.shape

(12488, 6)

# In[*]

# Number of rows

number_of_genes = len(data.index)

print(number_of_genes)\

12488

总共12488个基因，6个样本。

标准化

这一步是标准化，使用的是常见的log2()标准化

# In[*]

data2 = np.log2(data+0.0001)

data2.head()

Out[30]:

WT.GSM130365 WT.GSM130366 ... KO.GSM130369 KO.GSM130370

100001_at 3.523575 2.485453 ... 5.169929 5.392321

100002_at 4.357559 5.017926 ... 3.848007 4.517282

100003_at 6.177920 6.475735 ... 7.023478 6.437962

100004_at 8.027907 7.821456 ... 8.805099 8.170426

100005_at 10.085074 10.603256 ... 10.415636 10.488041

绘制盒须图

# In[*]

# Boxplot of each microarray

plt.show(data2.plot(kind = 'box', title = 'GSE5583 Boxplot', rot = 90))

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

# Density

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

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

# In[*]

# The mean of expression of the wt samples for each gene (row)

wt = data2.loc[:, 'WT.GSM130365' : 'WT.GSM130367'].mean(axis = 1)

wt.head()

# In[*]

# The mean of expression of the ko samples for each gene (row)

ko = data2.loc[:,'KO.GSM130368':'KO.GSM130370'].mean(axis = 1)

ko.head()

查看基因差异的差异倍数fold分布

# In[*]

fold = ko - wt

# Histogram of the fold-change

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.ix[i,0:3], data2.ix[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()