Python从零开始第五章生物信息学①提取差异基因

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

    1. 老本行是做实验的,做生信可能是为了辅助研究或者是为了发paper(有非常多的临床生选择趟生信这波水)
    1. 主要是做生信的,主要涵盖高通量测序数据分析,组学数据分析等等,专门从事生物学数据分析的这群人,其大部分也是本科生物狗作为强大的生力军,以调包写R,python为主。那么这群人就要熟悉看各种包的tutorial以及如何进行常规的数据的处理和分析等等。那么其实这群人,我个人认为对python的编程要求不是很高,在coursera先上两门课程,然后照着参考书在项目中练练就熟络了。
    1. 也是做生信的,但是主要以写包为主,什么意思?就是以开发算法供给第二类人用,做这部分工作的人需要具有良好的数理基础,所以一般大部分都是本科属于物理,数学等理工科的人在做这部分工作。当然对编程也提出较高要求。

我入门生物信息学是通过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()

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏大数据文摘

用递归神经网络,撰写一份特朗普式发言稿!

特朗普充满个人特色的语言风格让作者产生了兴趣,如果把他的推文和演讲稿都用于训练数据,再运用递归神经网络能否生成一份有特式风格的发言稿呢?结论是,如果数据和算力足...

11220

视觉信息理论

我喜欢有一个新的思维方式来思考这个世界。我特别喜欢把一些模糊的想法正式化为一个具体的概念。信息理论就是一个很好的例子。

23460
来自专栏量子位

早发arXiv可多获得65%的引用,但……

对研究人员来说,这么做主要有两个好处。一方面尽早占坑,另一方面可以绕过漫长的同行评议时间,加速圈内人交流工作进展。

11840
来自专栏互联网大杂烩

最优化模型 数据挖掘之优化模型

最短路径问题、网络最大流问题、最小费用最大流问题、最小生成树问题(MST)、旅行商问题(TSP)、图的着色问题。

17220
来自专栏斑斓

统计学中的相关性分析

掌握一点儿统计学介绍了统计学中常用到的函数,特别重点介绍了Standard Deviation(标准差)。接下来结合一个案例来谈谈相关性(Correlation...

44270
来自专栏大数据文摘

概率的意义:随机世界与大数法则

30840
来自专栏AI研习社

Mercari Price 比赛分享 —— 语言不仅是算法和公式而已

最近半年一直在忙于各种NLP比赛,除夕因为kaggle的price写到凌晨3点,最后靠rp爬回季军,也算圆了一个solo gold的梦想。这应该是我2017下半...

482120
来自专栏AI2ML人工智能to机器学习

人工智能深度学习人物关系[全]

为庆祝祖国母亲生日和自己终于脱单,给大家献上这篇绝对的干货!也提前祝大家国庆快乐!

15910
来自专栏机器学习AI算法工程

用R语言分析《我是歌手》出场顺序与名次的关系

《我是歌手》吵吵闹闹地落幕了,总决赛这一季是我最关注的一季,很认真的从头看到尾。当然,这篇文章的主旨不在此,我们要看的如题《我是歌手》节目中,出场...

31880
来自专栏PPV课数据科学社区

概率的意义(深度好文)

1987年,是印度传奇数学家拉曼努扬(SrinivasaRamanujan,1887-1920)的百年诞辰。为了纪念他,有一系列的活动。当代著名统计学者, 出生...

35670

扫码关注云+社区

领取腾讯云代金券