没有自己的服务器如何学习生物数据分析（上篇）

编者注：完整文章首发于作者博客 http://huboqiang.cn/

在这篇文章中，作者利用大数据平台 IBM data science 对生信技能树论坛的一道生物信息入门题进行了分析。

由于文章篇幅较长，我们将分为上篇和下篇分别进行推送。

其中上篇部分主要为大家介绍IBM data science 平台相关知识；下篇则为大家具体展示如何通过该平台运用pySpark来解决我们具体的问题。

希望对那些苦于没有自己的服务器而无法进行生物数据分析学习的朋友有所启发。同时，这篇文章也是非常好的大数据处理平台入门级介绍。

祝阅读愉快！以下为文章正文。

---

使用 IBM data science 平台统计hg38每条染色体转录本分布（上）

• 前言

这是一篇以生物信息学入门习题为例的大数据教程。具体而言，就是在 IBM 云计算平台，使用 pySpark 完成一个很简单的任务。任务描述如下：

每条染色体基因个数的分布？

所有基因平均有多少个转录本？

所有转录本平均有多个exon和intron？

注释文件一般以gtf/gff格式记录着！ 基础作业，就是对这个文件 ftp://ftp.ensembl.org/pub/release-87/gtf/homosapiens/Homosapiens.GRCh38.87.chr.gtf.gz 进行统计，普通人只需要关心第1，3列就好了。

源地址来自生信技能树 http://www.biotrainee.com/thread-626-1-1.html

这些代码可以使用 IBM data science 平台( http://datascience.ibm.com/ )的 Notebook 功能直接执行。

• IBM data science 平台介绍

IBM data science 平台对注册用户首月免费，默认提供一个 2核 CPU，预装 Rstudio, Jupyter。同时 IBM 公司作为 Spark 项目的重要商业支持者，对自家的Spark 大数据计算完美支持。

更重要的是，这个平台提供了很好的数据科学家相互交流的机会。编写的代码可以轻松在技术人员之间直接传阅，写完代码，最后的结果可以直接发给老板。

如果需要使用，首先需要在网站完成注册：

注册完成后，选择 DataHub

然后建立 Notebook，建立后的 Notebook 会在下面列出。

如果希望写个 HelloWorld，推荐你去简单画个图，源码位于 matplotlib 官方 gallery，我们开个选择 Python，生成Jupyter Notebook 画一下：

然后把 matplotlib 上的这个例子的代码复制粘贴下来：

1. import numpy as np

2. import matplotlib.pyplot as plt

3. n_groups = 5

4. means_men = (20, 35, 30, 35, 27)

5. std_men = (2, 3, 4, 1, 2)

6. means_women = (25, 32, 34, 20, 25)

7. std_women = (3, 5, 2, 3, 3)

8. fig, ax = plt.subplots()

9. index = np.arange(n_groups)

10. bar_width = 0.35

11. opacity = 0.4

12. error_config = {'ecolor': '0.3'}

13. rects1 = plt.bar(index, means_men, bar_width,

14.                 alpha=opacity,

15.                 color='b',

16.                 yerr=std_men,

17.                 error_kw=error_config,

18.                 label='Men')

19. rects2 = plt.bar(index + bar_width, means_women, bar_width,

20.                 alpha=opacity,

21.                 color='r',

22.                 yerr=std_women,

23.                 error_kw=error_config,

24.                 label='Women')

25. plt.xlabel('Group')

26. plt.ylabel('Scores')

27. plt.title('Scores by group and gender')

28. plt.xticks(index + bar_width / 2, ('A', 'B', 'C', 'D', 'E'))

29. plt.legend()

30. plt.tight_layout()

31. plt.show()

如图操作，就可以得到 matplotlib 官网上的图。

神马？没有出图像？来，这里有个特殊的地方，需要在 import 完所有库之后，加一行 %matplotlibinline 魔法，允许直接在代码块下面显示，就像我图中写的那样。

Jupyter Notebook 和 Rstudio 一样有很多方便好用的快捷键。具体可以点击 Help => ShortCut 看一下。Help 位于箭头挡住的 run 又上方。

如果不使用 IBM data science 平台，也可以自己下载 anaconda 安装科学计算版 Python。但这样 spark 需要自己预装，特别是Spark 通常需要预装 Hadoop，对生信菜鸟比较费事，大型机上管理员也不一定让你装。不过 anaconda 本身不使用 spark 加成，开 Jupyter Notebook 就已经十分强大了，建议大家试一试。

我在我们的大型机的一个计算节点装好 anaconda 后，根据 Jupyter Notebook 官方文档，设定集群访问http://jupyter-notebook.readthedocs.io/en/latest/public_server.html，需要分析项目，会首先 cd 到项目所在的分析文件夹(鄙视放进 /home 目录里的人)， 接着 cmd 输入 jupyter notebook，这样jupyter 会在后端挂起，然后访问 https://IP:PORT，IP 是该集群的内网 IP，端口在上一步指定，默认 8888，注意是这里是 https 不是 http，然后允许打开网页，输入集群访问密码，就会进入管理页面，如图。

新建Python Notebook 后，会直接进入该文件，管理页面里面会出现一个绿色的 ipynb 文件，linux shell 里面也能看见。

这个文件就是Jupyter Notebook所在的文件，用法与 IBM datascience 的完全相同，大家也可以照着上图 HelloWorld 一下。

我这里建议，如果想体验一把 PySpark，使用 IBM data science ，即使是菜鸟，也可以来体验一把高大上的大数据+云计算。

默认环境配置如下：

1. Instance

2. DataSciX

3. Language

4. Python 2.7

5. Spark as a Service

6. Apache Spark 2.0

7. Preinstalled Libraries

8. biopython-1.66

9. bitarray-0.8.1

10. brunel-1.1

11. iso8601-0.1.11

12. jsonschema-2.5.1

13. lxml-3.5.0

14. matplotlib-1.5.0

15. networkx-1.10

16. nose-1.3.7

17. numexpr-2.4.6

18. numpy-1.10.4

19. pandas-0.17.1

20. Pillow-3.0.0

21. pip-8.1.0

22. pyparsing-2.0.6

23. pytz-2015.7

24. requests-2.9.1

25. scikit-learn-0.17

26. scipy-0.17.0

27. simplejson-3.8.1

当然，在运行程序前，需要预装 python 的类似 ggplot2 画图包，方便作图。

最前头的感叹号的意思是这一行执行 shell 脚本。也就是说这个命令本应在 linux shell 里面执行，但由于 jupyter 把 shell 也给完美的集成了进来，所以在 notebook 中写就 OK。

代码块【1】：

1. !pip install --user seaborn

结果是：

1. Requirement already satisfied (use --upgrade to upgrade): seaborn in /gpfs/global_fs01/sym_shared/YPProdSpark/user/sa9e-127e054d347dc8-cc04bd2554cc/.local/lib/python2.7/site-packages

这里我是因为已经装过，否则会显示下载进度以及提示安装完成。

接下来的程序不以 ！ 开头，写的就是 python 代码了。这几行是载入需要用到的包。

代码块【2】：

1. # -*- coding: utf-8 -*-

2. import pandas  as pd              

3. import matplotlib.pyplot as plt    

4. import seaborn as sns              

5. import requests, StringIO, json    

6. import re                        

7. import numpy as np                

8. %matplotlib inline                

9. sns.set_style("white")            

简单解释下：

1. import pandas  as pd               # 数据分析包

2. import matplotlib.pyplot as plt    # 作图包

3. import seaborn as sns              # 类似 ggplot2 的高级作图包

4. import requests, StringIO, json    # 在IBM datascience 的存储系统里读取数据必须

5. import re                          # 正则表达式

6. import numpy as np                 # 矩阵运算，这里用在最后算中位数

7. %matplotlib inline                 # 允许作图函数执行后，在代码块下面直接显示所画的图

8. sns.set_style("white")             # 类似 ggplot2 的 theme_bw()

• Jupyter 借助 Spark 轻松实现 Python 多核心编程

看起来 Jupyter 既可以像 Rstudio 一样轻松的写 python 代码，又可以在代码块的上下写各种注释，一边写程序，一边出图，然后直接写图注，文章的逻辑写在旁边，是非常完美的轻量级科研利器。

Jupyter + pyspark 虽然轻量，但其实力气一点都不小。写出来的性能，在某种意义上甚至高于 C++ Java 这样的低级语言。我说某种意义，指的是单核运算方面的瓶颈。因为自从本世纪初我们耳熟能详的奔腾4处理器的主频，已经高达1.4GHz了，大家可以去看看自己电脑的主频，会发现这么些年来虽然 CPU 各种进步，主频的提升实际上是非常小的。CPU 的摩尔定律，主要还是在 核心数以及线程数 的提升。家用笔记本现在很多都是2核4线程，而服务器的单 CPU 线程数一般也都在 10 个以上。如果还是通过 for 循环读取各行，就会导致“一核干活，十核围观”，对计算资源造成闲置浪费。

所以，为了进一步跟上时代潮流，重要的软件程序，我们都使用多核心编程技术。我们生物信息领域很多耳熟能详的软件，如比对用的 bwa bowtie 的参数，都有使用几个核心的选项。

那么我们能不能也轻松写一个多核心程序出来呢？C++ 的套路是，用多核心就要调 pthread 库，通常是先在 IO 读取数据，缓存1,000,000行，然后通过 pthread 把缓存里面的各行扔给需要调用的函数。具体把哪一行扔给函数，也需要自己指定，比如当前的行数取余数，余几就扔给几号CPU。然后还需要预留一块内存接各个CPU 执行函数的输出结果，不能直接输出。。。可能菜鸟已经听晕了，不知道在说什么，而听懂的人想必是清楚其中的麻烦是我这几行远远没有说明白的。

这一问题在 Python 和 R 中也或多或少的存在。这两种语言系统支持的多线程技术调起来也只是稍微简单一些，而性能却没法和C++比。于是乎，在这个大数据的时代背景下，他们抱上了 Hadoop Spark 这些最新的大数据工具的大腿。特别是 Spark。

Spark 源码是通过一种叫做 Scala 的语言编写的。Scala 是脱胎于 java 的一种更高效的编程语言，一般人还真不会用，于是 Spark 项目就打通了 Python R 的使用接口。然而为了保证版本升级的进度，Spark 的新功能一般是首先 Java Scala 能用，然后轮到 Python，最后才到 R。比如 Spark 的机器学习库，目前 Python 已经能很好支持了，而 R语言得等到 2.2.0（16年11月 IBM 的 Spark机器学习库编写人员亲口所说）。

虽然 PySpark 用的是一种不完整的 Spark，但用它对列式数据（R 中的 dataframe 类型）搞分组求和、文件清洗，已经足够了。更重要的是，这里由于是和数据科学界接轨，强烈推荐把数据简单处理后（抓取信息，规定每一列的名称，扔掉某些行），放进 SparkSQL中，用 SQL 语句，用 人话 而不是代码，去人机交互，分析数据。

最后，多说一句，就是实际上 Spark 能做的已经不是单机多核心了，哪怕是上百台电脑，处理和本文一模一样的一个 TB 级的基因注释GTF文件（就当是外星人的），代码怎么写？一模一样，只要 Spark 指挥的 Hadoop 集群被合理的配置好，PySpark 代码方面一模一样，上百台电脑，上千个 CPU 核心，共同处理同一文件。当然这个文件需要被放入 HDFS 分布式存储系统中，命令也很简单：

1. /hadoop/bin/hdfs dfs -put 外星人.GTF hdfs://[HDFS系统IP]:[HDFS系统端口]:[文件路径/外星人.GTF]

• 下载和上传数据

继续说正事，分析数据之前，必须做的一件事，就是上传数据。而上传数据的第一步，是得把数据先给下载下来。

我们的数据，就是从 ftp://ftp.ensembl.org/pub/releas ... RCh38.87.chr.gtf.gz 下载的压缩状态的gtf 文件，不解压缩，直接上传到 IBM data 平台。

方法如下：

选择 Insert SparkSession Step 后，系统会自动生成一系列代码，如下，除 spark.read.text(path_1).show()：

代码块【3】：

1. from pyspark.sql import SparkSession

2. # @hidden_cell

3. # This function is used to setup the access of Spark to your Object Storage. The definition contains your credentials.

4. # You might want to remove those credentials before you share your notebook.

5. def set_hadoop_config_with_credentials_4fadd3c3d7a24e3b9b78d895f084bb84(name):

6.    """This function sets the Hadoop configuration so it is possible to

7.    access data from Bluemix Object Storage using Spark"""

8.    prefix = 'fs.swift.service.' + name

9.    hconf = sc._jsc.hadoopConfiguration()

10.    hconf.set(prefix + '.auth.url', 'https://identity.open.softlayer.com'+'/v3/auth/tokens')

11.    hconf.set(prefix + '.auth.endpoint.prefix', 'endpoints')

12.    hconf.set(prefix + '.tenant', '个人保密')

13.    hconf.set(prefix + '.username', '个人保密')

14.    hconf.set(prefix + '.password', '个人保密')

15.    hconf.setInt(prefix + '.http.port', 8080)

16.    hconf.set(prefix + '.region', 'dallas')

17.    hconf.setBoolean(prefix + '.public', False)

18. # you can choose any name

19. name = 'keystone'

20. set_hadoop_config_with_credentials_4fadd3c3d7a24e3b9b78d895f084bb84(name)

21. spark = SparkSession.builder.getOrCreate()

22. # Please read the documentation of PySpark to learn more about the possibilities to load data files.

23. # PySpark documentation: https://spark.apache.org/docs/2.0.1/api/python/pyspark.sql.html#pyspark.sql.SparkSession

24. # The SparkSession object is already initalized for you.

25. # The following variable contains the path to your file on your Object Storage.

26. path_1 = "swift://mydemo." + name + "/Homo_sapiens.GRCh38.87.chr.gtf.gz"

27. spark.read.text(path_1).show()

结果：

1. +--------------------+

2. |               value|

3. +--------------------+

4. |#!genome-build GR...|

5. |#!genome-version ...|

6. |#!genome-date 201...|

7. |#!genome-build-ac...|

8. |#!genebuild-last-...|

9. |1    havana  gene    118...|

10. |1    havana  transcri...|

11. |1    havana  exon    118...|

12. |1    havana  exon    126...|

13. |1    havana  exon    132...|

14. |1    havana  transcri...|

15. |1    havana  exon    120...|

16. |1    havana  exon    121...|

17. |1    havana  exon    126...|

18. |1    havana  exon    129...|

19. |1    havana  exon    132...|

20. |1    havana  exon    134...|

21. |1    havana  gene    144...|

22. |1    havana  transcri...|

23. |1    havana  exon    295...|

24. +--------------------+

25. only showing top 20 rows

直接读，发现其实行是识别出来了，是个 DataFrame，但是列不对。首先是前几行注释需要扔掉，其次是我们需要的基因名称、外显子名称这些内容需要单独被分出一列。

于是我们通过 Python 的正则表达式 re 包，配合 PySpark 的 RDD 相关操作，做数据清洗以及特征提取。

• Spark 的 RDD和 DataFrame 是什么

简单说，RDD 可以理解成我们以前开文件后逐行读取每行信息，不直接处理，好多行给缓存成了个列表，这个列表就类似是 RDD。而 DataFrame 则类似是R 中的 DataFrame，RDD + 表头。

但是 这里的 RDD 虽然类似列表，DataFrame 虽然也跟 R 很像，却都不支持行列操作。只可以显示最上面的几行， 如 rdd.take(5) 或者 DataFrame.show(5)显示最上面的5行，却不支持显示例如第250行这样的命令。原因是， RDD DataFrame 这里并不位于内存，而是位于前面提到的 HDFS，在硬盘里！所以一个行下标是不能直接调出数据的。内存只是存了指针指向了硬盘，多个CPU来要数据时，内存的指针快速给他们在分布式的存储系统给他们分配任务。这也是为什么 Spark 可以Hold住海量数据的真实原因，数据不需要全扔进内存。

更多他们之间的区别联系，推荐阅读：https://databricks.com/blog/2016/07/14/a-tale-of-three-apache-spark-apis-rdds-dataframes-and-datasets.html

代码块【4】：

1. pat_gene = '''gene_id\s+\"(\S+)\";'''

2. pat_tran = '''transcript_id\s+\"(\S+)\";'''

3. pat_exon = '''exon_number\s+\"*(\w+)\"*'''

4. pattern_gene = re.compile( pat_gene )

5. pattern_tran = re.compile( pat_tran )

6. pattern_exon = re.compile( pat_exon )

7. def parseEachLine(f_line):

8.    match_gene = pattern_gene.search( f_line[-1] )

9.    match_tran = pattern_tran.search( f_line[-1] )

10.    match_exon = pattern_exon.search( f_line[-1] )

11.    gene = "NULL"

12.    tran = "NULL"

13.    exon = "NULL"

14.    if match_gene:

15.        gene = match_gene.group(1)

16.    if match_tran:

17.        tran = match_tran.group(1)

18.    if match_exon:

19.        exon = match_exon.group(1)

20.    return [gene, tran, exon, f_line[0]]

21. rdd = spark.read.text(path_1).rdd\

22.            .filter(lambda x: x.value[0]!= "#")\

23.            .map(lambda x: x.value.split("\t"))\

24.            .map(lambda x: parseEachLine(x))

25. rdd.take(5)

结果：

1. [[u'ENSG00000223972', 'NULL', 'NULL', u'1']

2. [u'ENSG00000223972', u'ENST00000456328', 'NULL', u'1']

3. [u'ENSG00000223972', u'ENST00000456328', u'1', u'1'],

4. [u'ENSG00000223972', u'ENST00000456328', u'2', u'1'],

5. [u'ENSG00000223972', u'ENST00000456328', u'3', u'1']]

简单解释下：

•     第一部分，设定正则表达式：

1. pat_gene = '''gene_id\s+\"(\S+)\";'''

2. pat_tran = '''transcript_id\s+\"(\S+)\";'''

3. pat_exon = '''exon_number\s+\"*(\w+)\"*'''

4. pattern_gene = re.compile( pat_gene )

5. pattern_tran = re.compile( pat_tran )

6. pattern_exon = re.compile( pat_exon )

7. def parseEachLine(f_line):

8.    match_gene = pattern_gene.search( f_line[-1] )

9.    match_tran = pattern_tran.search( f_line[-1] )

10.    match_exon = pattern_exon.search( f_line[-1] )

11.    gene = "NULL"

12.    tran = "NULL"

13.    exon = "NULL"

14.    if match_gene:

15.        gene = match_gene.group(1)

16.    if match_tran:

17.        tran = match_tran.group(1)

18.    if match_exon:

19.        exon = match_exon.group(1)

20.    return [gene, tran, exon, f_line[0]]

这部分是正则表达式。前几行规定我们从 geneid transcriptid exon_id 这几个字段后面抓数据，并且抓引号里面的内容。

• 第二部分，正则表达式函数用于rdd 的各行，提取我们需要的信息

这个被进一步写成了一个函数。这个函数是可以直接用于逐行读取结构的，如下：

1. with gzip.open("input.gtf.gz", "rb") as f_gtf:

2.    for line in f_gtf:

3.        if line[0] == "#":

4.            continue

5.        f_line = f_gtf.split("\t")

6.        [gene, tran, exon, chrom] = parseEachLine(f_line)

也可以被直接用在进过类似 split 处理的 RDD 字段上。

1. rdd = spark.read.text(path_1).rdd\

2.            .filter(lambda x: x.value[0]!= "#")\

3.            .map(lambda x: x.value.split("\t"))\

4.            .map(lambda x: parseEachLine(x))

处理后的 rdd 长这样：

1. [[u'ENSG00000223972', 'NULL', 'NULL', u'1'],

2. [u'ENSG00000223972', u'ENST00000456328', 'NULL', u'1'],

3. [u'ENSG00000223972', u'ENST00000456328', u'1', u'1'],

4. [u'ENSG00000223972', u'ENST00000456328', u'2', u'1'],

5. [u'ENSG00000223972', u'ENST00000456328', u'3', u'1']]

得到这个以后，后续处理想必大家都跃跃欲试了。传统的 Hadoop 使用的 MapReduce 结构，有这个就够了。但写出的代码终归不太好看。而我们需要的，是 说人话，怎么说人话呢，就是给RDD加一个表头，变成 DataFrame，然后通过 SparkSQL ，用 SQL 语句进行人机交互。代码如下：

代码块【5】：

1. from pyspark.sql.types import *

2. schema=StructType(

3.    [StructField("Gene",  StringType())] +

4.    [StructField("Tran",  StringType())] +

5.    [StructField("Exon",  StringType())] +

6.    [StructField("Chrom", StringType())]

7. )

8. df = sqlCtx.createDataFrame(rdd, schema)

9. df.show()

结果：

1. +---------------+---------------+----+-----+

2. |           Gene|           Tran|Exon|Chrom|

3. +---------------+---------------+----+-----+

4. |ENSG00000223972|           NULL|NULL|    1|

5. |ENSG00000223972|ENST00000456328|NULL|    1|

6. |ENSG00000223972|ENST00000456328|   1|    1|

7. |ENSG00000223972|ENST00000456328|   2|    1|

8. |ENSG00000223972|ENST00000456328|   3|    1|

9. |ENSG00000223972|ENST00000450305|NULL|    1|

10. |ENSG00000223972|ENST00000450305|   1|    1|

11. |ENSG00000223972|ENST00000450305|   2|    1|

12. |ENSG00000223972|ENST00000450305|   3|    1|

13. |ENSG00000223972|ENST00000450305|   4|    1|

14. |ENSG00000223972|ENST00000450305|   5|    1|

15. |ENSG00000223972|ENST00000450305|   6|    1|

16. |ENSG00000227232|           NULL|NULL|    1|

17. |ENSG00000227232|ENST00000488147|NULL|    1|

18. |ENSG00000227232|ENST00000488147|   1|    1|

19. |ENSG00000227232|ENST00000488147|   2|    1|

20. |ENSG00000227232|ENST00000488147|   3|    1|

21. |ENSG00000227232|ENST00000488147|   4|    1|

22. |ENSG00000227232|ENST00000488147|   5|    1|

23. |ENSG00000227232|ENST00000488147|   6|    1|

24. +---------------+---------------+----+-----+

25. only showing top 20 rows

这里表头已经加上去了。懂得 R 语言 melt ddply dcast 套路的人一定知道该怎么做了。但其实更多的数据从业者不懂 R，而是 SQL 专家，所以 SparkSQL 作为通用处理框架，提供的是 SQL 语句作为解决思路。

（未完待续……）

---

上篇的内容就介绍到这里，如果你耐心读完一定觉得知识量很大。虽然文章较长，但如果你是一个想入门大数据的有心人，仔细阅读体会并且实践 一定会有很大的收获。

再下篇中，我们将介绍如何利用该平台和PySpark具体解决我们的生物信息数据分析问题。

敬请期待！