纳速健身

 找回密码
 注册

扫一扫,访问微社区

查看: 17|回复: 0

[Python] [python生物信息学]Python从零开始第五章生物信息学一、提取差异基因

[复制链接]
发表于 2021-1-10 11:26:36 | 显示全部楼层 |阅读模式
[python生物信息学]Python从零开始第五章生物信息学一、提取差异基因

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

  • 老本行是做实验的,做生信可能是为了辅助研究或者是为了发paper(有非常多的临床生选择趟生信这波水)
  • 主要是做生信的,主要涵盖高通量测序数据分析,组学数据分析等等,专门从事生物学数据分析的这群人,其大部分也是本科生物狗作为强大的生力军,以调包写R,python为主。那么这群人就要熟悉看各种包的tutorial以及如何进行常规的数据的处理和分析等等。那么其实这群人,我个人认为对python的编程要求不是很高,在coursera先上两门课程,然后照着参考书在项目中练练就熟络了。
  • 也是做生信的,但是主要以写包为主,什么意思?就是以开发算法供给第二类人用,做这部分工作的人需要具有良好的数理基础,所以一般大部分都是本科属于物理,数学等理工科的人在做这部分工作。当然对编程也提出较高要求。
  • 我入门生物信息学是通过R语言入门的,但是接触到了python,这个也是目前用户数量数一数二的语言。python去做生信得优点是①过程更加直观,因为常见的R包功能一般已经封装好了,直接应用就可,虽然足够简单友好,但是不利于长期学习②基因组数据一般比较大,python速度一般比R快。


下面我就记录一个通过python做生信分析的流程。
使用的数据集是GSE5583,来自于2006年的基因芯片结果,该芯片目的是提取野生型和HDAC1小鼠胚胎干细胞用于Affymetrix微阵列上的差异RNA。

导入包

  1. %clear
  2. %reset -f
  3. # In[*]
  4. import seaborn as sns
  5. import matplotlib.pyplot as plt
  6. %matplotlib inline
  7. import os
  8. import numpy as np
  9. import pandas as pd
  10. os.chdir('D:\\train')
  11. from scipy import stats
复制代码

导入数据
  1. # In[*]
  2. data = pd.read_table("http://ahmedmoustafa.io/notebooks/GSE5583.txt",
  3.                      header = 0,
  4.                      index_col = 0)
  5. 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

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

查看数据维度
  1. # In[*]
  2. # Check the dimension of the loaded data (rows & columns)
  3. data.shape
  4. (12488, 6)
  5. # In[*]
  6. # Number of rows
  7. number_of_genes = len(data.index)
  8. print(number_of_genes)\
复制代码
12488

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

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

  1. # In[*]
  2. data2 = np.log2(data+0.0001)
  3. 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

绘制盒须图
  1. # In[*]
  2. # Boxplot of each microarray
  3. plt.show(data2.plot(kind = 'box', title = 'GSE5583 Boxplot', rot = 90))
复制代码


01.webp.jpg


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

  1. # Density
  2. plt.show(data2.plot(kind = 'density', title = 'GSE5583 Density'))
复制代码



02.webp.jpg


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

  1. # In[*]
  2. # The mean of expression of the wt samples for each gene (row)
  3. wt = data2.loc[:, 'WT.GSM130365' : 'WT.GSM130367'].mean(axis = 1)
  4. wt.head()

  5. # In[*]
  6. # The mean of expression of the ko samples for each gene (row)
  7. ko = data2.loc[:,'KO.GSM130368':'KO.GSM130370'].mean(axis = 1)
  8. ko.head()
复制代码


查看基因差异的差异倍数fold分布
  1. # In[*]
  2. fold = ko - wt


  3. # Histogram of the fold-change
  4. plt.hist(fold)
  5. plt.title("Histogram of fold-change")
  6. plt.show()
复制代码



03.webp.jpg


查看基因差异的P值分布
  1. from scipy import stats

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

  6. # Histogram of the p-values
  7. plt.hist(-np.log(pvalue))
  8. plt.title("Histogram of p-value")
  9. plt.show()
复制代码


04.webp.jpg
您需要登录后才可以回帖 登录 | 注册

本版积分规则

快速回复 返回顶部 返回列表