搜索
查看: 2206|回复: 0

使用 IBM datascience 平台统计 hg38每条染色体基因,转录本的分...

[复制链接]

2

主题

3

帖子

26

积分

新手上路

Rank: 1

积分
26
发表于 2017-2-16 17:13:06 | 显示全部楼层 |阅读模式
本帖最后由 hubq 于 2017-2-16 17:26 编辑

由于论坛不支持markdown编写,粘贴比较麻烦,只展示了文章一部分,完整版本见我的个人博客 http://huboqiang.cn/2017/02/16/AnalyseTheAnnotationOfHumanGenome

同时所有代码在我的IBM data science 平台上的仓库公开分享


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

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

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

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

注释文件一般以gtf/gff格式记录着!

基础作业,就是对这个文件 ftp://ftp.ensembl.org/pub/releas ... RCh38.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 平台对注册用户首月免费,默认提供一个 2核 CPU,预装 Rstudio, Jupyter。同时 IBM 公司作为 Spark 项目的重要商业支持者,对自家的Spark 大数据计算完美支持。

如果不使用 IBM data science 平台,也可以自己下载 anaconda 安装科学计算版 Python。但这样 spark 需要自己预装,特别是Spark 通常需要预装 Hadoop,对生信菜鸟比较费事,大型机上管理员也不一定让你装。

所以我这里建议使用 IBM data science ,作为生信菜鸟,也来体验一把高大上的大数据+云计算。

默认环境配置如下:
Instance
DataSciX


Language
Python 2.7 Spark as a Service
Apache Spark 2.0

Preinstalled Libraries

biopython-1.66
bitarray-0.8.1
brunel-1.1
iso8601-0.1.11
jsonschema-2.5.1
lxml-3.5.0
matplotlib-1.5.0
networkx-1.10
nose-1.3.7
numexpr-2.4.6
numpy-1.10.4
pandas-0.17.1
Pillow-3.0.0
pip-8.1.0
pyparsing-2.0.6
pytz-2015.7
requests-2.9.1
scikit-learn-0.17
scipy-0.17.0
simplejson-3.8.1



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

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

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

代码块【1】:
[Python] 纯文本查看 复制代码
!pip install --user seaborn


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

代码块【2】:

[Python] 纯文本查看 复制代码
# -*- coding: utf-8 -*-
import pandas  as pd               
import matplotlib.pyplot as plt    
import seaborn as sns              
import requests, StringIO, json    
import re                         
import numpy as np                


%matplotlib inline                
sns.set_style("white")            

简单解释下:
import pandas  as pd               # 数据分析包
import matplotlib.pyplot as plt    # 作图包
import seaborn as sns              # 类似 ggplot2 的高级作图包
import requests, StringIO, json    # 在IBM datascience 的存储系统里读取数据必须
import re                          # 正则表达式
import numpy as np                 # 矩阵运算,这里用在最后算中位数


%matplotlib inline                 # 允许作图函数执行后,在代码块下面直接显示所画的图
sns.set_style("white")             # 类似 ggplot2 的 theme_bw()

2. 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 分布式存储系统中,命令也很简单:

[Python] 纯文本查看 复制代码
/hadoop/bin/hdfs dfs -put 外星人.GTF hdfs://[HDFS系统IP]:[HDFS系统端口]:[文件路径/外星人.GTF]

继续说正事,IBM有简单的接口实现spark读取输入文件,会自动生成如下代码:

代码块【3】:

[Python] 纯文本查看 复制代码
from pyspark.sql import SparkSession


# @hidden_cell
# This function is used to setup the access of Spark to your Object Storage. The definition contains your credentials.
# You might want to remove those credentials before you share your notebook.
def set_hadoop_config_with_credentials_4fadd3c3d7a24e3b9b78d895f084bb84(name):
    """This function sets the Hadoop configuration so it is possible to
    access data from Bluemix Object Storage using Spark"""


    prefix = 'fs.swift.service.' + name
    hconf = sc._jsc.hadoopConfiguration()
    hconf.set(prefix + '.auth.url', 'https://identity.open.softlayer.com'+'/v3/auth/tokens')
    hconf.set(prefix + '.auth.endpoint.prefix', 'endpoints')
    hconf.set(prefix + '.tenant', '个人保密')
    hconf.set(prefix + '.username', '个人保密')
    hconf.set(prefix + '.password', '个人保密')
    hconf.setInt(prefix + '.http.port', 8080)
    hconf.set(prefix + '.region', 'dallas')
    hconf.setBoolean(prefix + '.public', False)


# you can choose any name
name = 'keystone'
set_hadoop_config_with_credentials_4fadd3c3d7a24e3b9b78d895f084bb84(name)


spark = SparkSession.builder.getOrCreate()


# Please read the documentation of PySpark to learn more about the possibilities to load data files.
# PySpark documentation: [url=https://spark.apache.org/docs/2.]https://spark.apache.org/docs/2.[/url] ... rk.sql.SparkSession
# The SparkSession object is already initalized for you.
# The following variable contains the path to your file on your Object Storage.
path_1 = "swift://mydemo." + name + "/Homo_sapiens.GRCh38.87.chr.gtf.gz"


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


结果:
[AppleScript] 纯文本查看 复制代码
+--------------------+
|               value|
+--------------------+
|#!genome-build GR...|
|#!genome-version ...|
|#!genome-date 201...|
|#!genome-build-ac...|
|#!genebuild-last-...|
|1        havana        gene        118...|
|1        havana        transcri...|
|1        havana        exon        118...|
|1        havana        exon        126...|
|1        havana        exon        132...|
|1        havana        transcri...|
|1        havana        exon        120...|
|1        havana        exon        121...|
|1        havana        exon        126...|
|1        havana        exon        129...|
|1        havana        exon        132...|
|1        havana        exon        134...|
|1        havana        gene        144...|
|1        havana        transcri...|
|1        havana        exon        295...|
+--------------------+
only showing top 20 rows

直接读,发现其实行是识别出来了,是个 DataFrame,但是列不对。

首先是前几行注释需要扔掉,其次是我们需要的基因名称、外显子名称这些内容需要单独被分出一列。

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

Q: 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住海量数据的真实原因,数据不需要全扔进内存。


代码块【4】:
[Python] 纯文本查看 复制代码
pat_gene = '''gene_id\s+\"(\S+)\";'''
pat_tran = '''transcript_id\s+\"(\S+)\";'''
pat_exon = '''exon_number\s+\"*(\w+)\"*'''

pattern_gene = re.compile( pat_gene )
pattern_tran = re.compile( pat_tran )
pattern_exon = re.compile( pat_exon )

def parseEachLine(f_line):
    match_gene = pattern_gene.search( f_line[-1] )
    match_tran = pattern_tran.search( f_line[-1] )
    match_exon = pattern_exon.search( f_line[-1] )
         
    gene = "NULL"
    tran = "NULL"
    exon = "NULL"
    if match_gene:
        gene = match_gene.group(1)
    if match_tran:
        tran = match_tran.group(1)
    if match_exon:
        exon = match_exon.group(1)
    
    return [gene, tran, exon, f_line[0]]

rdd = spark.read.text(path_1).rdd\
          .filter(lambda x: x.value[0]!= "#")\
          .map(lambda x: x.value.split("\t"))\
          .map(lambda x: parseEachLine(x))

rdd.take(5)

简单解释下:

第一部分,设定正则表达式:
[Python] 纯文本查看 复制代码
pat_gene = '''gene_id\s+\"(\S+)\";'''
pat_tran = '''transcript_id\s+\"(\S+)\";'''
pat_exon = '''exon_number\s+\"*(\w+)\"*'''

pattern_gene = re.compile( pat_gene )
pattern_tran = re.compile( pat_tran )
pattern_exon = re.compile( pat_exon )

def parseEachLine(f_line):
    match_gene = pattern_gene.search( f_line[-1] )
    match_tran = pattern_tran.search( f_line[-1] )
    match_exon = pattern_exon.search( f_line[-1] )
         
    gene = "NULL"
    tran = "NULL"
    exon = "NULL"
    if match_gene:
    gene = match_gene.group(1)
    if match_tran:
    tran = match_tran.group(1)
    if match_exon:
    exon = match_exon.group(1)
    
    return [gene, tran, exon, f_line[0]]

这部分是正则表达式。前几行规定我们从 gene_id transcript_id exon_id 这几个字段后面抓数据,并且抓引号里面的内容。
第二部分,正则表达式函数用于rdd 的各行,提取我们需要的信息
这个被进一步写成了一个函数。这个函数是可以直接用于逐行读取结构的,如下:

[Python] 纯文本查看 复制代码
with gzip.open("input.gtf.gz", "rb") as f_gtf:   
    for line in f_gtf:        
        if line[0] == "#":            
            continue        
         f_line = f_gtf.split("\t")      
         [gene, tran, exon, chrom] = parseEachLine(f_line)
也可以被直接用在进过类似 split 处理的 RDD 字段上。


[Python] 纯文本查看 复制代码
rdd = spark.read.text(path_1).rdd\ 
           .filter(lambda x: x.value[0]!= "#")\ 
           .map(lambda x: x.value.split("\t"))\
           .map(lambda x: parseEachLine(x))

处理后的 rdd 长这样:

[Python] 纯文本查看 复制代码
[[u'ENSG00000223972', 'NULL', 'NULL', u'1'],
 [u'ENSG00000223972', u'ENST00000456328', 'NULL', u'1'],
 [u'ENSG00000223972', u'ENST00000456328', u'1', u'1'],
 [u'ENSG00000223972', u'ENST00000456328', u'2', u'1'],
 [u'ENSG00000223972', u'ENST00000456328', u'3', u'1']]



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

[Python] 纯文本查看 复制代码
from pyspark.sql.types import *


schema=StructType(
    [StructField("Gene",  StringType())] + 
    [StructField("Tran",  StringType())] + 
    [StructField("Exon",  StringType())] +
    [StructField("Chrom", StringType())]
)


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

结果:
[AppleScript] 纯文本查看 复制代码
+---------------+---------------+----+-----+
|           Gene|           Tran|Exon|Chrom|
+---------------+---------------+----+-----+
|ENSG00000223972|           NULL|NULL|    1|
|ENSG00000223972|ENST00000456328|NULL|    1|
|ENSG00000223972|ENST00000456328|   1|    1|
|ENSG00000223972|ENST00000456328|   2|    1|
|ENSG00000223972|ENST00000456328|   3|    1|
|ENSG00000223972|ENST00000450305|NULL|    1|
|ENSG00000223972|ENST00000450305|   1|    1|
|ENSG00000223972|ENST00000450305|   2|    1|
|ENSG00000223972|ENST00000450305|   3|    1|
|ENSG00000223972|ENST00000450305|   4|    1|
|ENSG00000223972|ENST00000450305|   5|    1|
|ENSG00000223972|ENST00000450305|   6|    1|
|ENSG00000227232|           NULL|NULL|    1|
|ENSG00000227232|ENST00000488147|NULL|    1|
|ENSG00000227232|ENST00000488147|   1|    1|
|ENSG00000227232|ENST00000488147|   2|    1|
|ENSG00000227232|ENST00000488147|   3|    1|
|ENSG00000227232|ENST00000488147|   4|    1|
|ENSG00000227232|ENST00000488147|   5|    1|
|ENSG00000227232|ENST00000488147|   6|    1|
+---------------+---------------+----+-----+
only showing top 20 rows

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

我们思考一下提问者的几个问题:
  • 每条染色体基因个数的分布?
  • 所有基因平均有多少个转录本?
  • 所有转录本平均有多个exon和intron?


我们思考一下怎么翻译这几句话成 SQL 语句:

每条染色体基因个数的分布?
思考一下,问的其实是:
每个 Chrom 值,对应 几种、不重复的Gene?

[SQL] 纯文本查看 复制代码
    SELECT Chrom, COUNT(DISTINCT(Gene)) FROM GTF    GROUP BY Chrom


“每个” Chrom 意味着 GROUP BY Chrom, 与此同时,前面SELECT 就得加上 Chrom,这样最后的结果才会显示是哪个染色体。
“几种” 翻译成 COUNT,不重复翻译成 “DISTINCT”,于是合并后就是 COUNT(DISTINC(Gene))
然后我们要保证只看 Gene Tran Exon 都非空的行,即 WHERE (Gene != “NULL”) AND (Tran != “NULL”) AND (Exon != “NULL”)。
于是写出SQL:
   
[Python] 纯文本查看 复制代码
SELECT Chrom, COUNT(DISTINCT(Gene)) FROM GTF    WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")    GROUP BY Chrom

了解更多 SQL 可以在这个网站 https://www.w3schools.com/sql/ 在线学习基本的SQL语句。
SQL 语句在调取 UCSC 数据集中同样作用巨大,具体这里不表。
这句话怎么在 DataFrame 执行?需要先 registerTempTable("GTF") 把 df 这个 dataFrame 给 SparkSQL,取一个表名,叫做 “GTF”。这样 df 就可以直接用 SQL语句分析了。
代码块【6】:

[Python] 纯文本查看 复制代码
df.registerTempTable("GTF")
sqlDF_genesInEachChr = spark.sql("""
    SELECT Chrom, COUNT(DISTINCT(Gene)) AS Cnt FROM GTF
    WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")
    GROUP BY Chrom
""")
sqlDF_genesInEachChr.show()

结果:

[AppleScript] 纯文本查看 复制代码
+-----+----+
|Chrom| Cnt|
+-----+----+
|    7|2867|
|   15|2152|
|   11|3235|
|    3|3010|
|    8|2353|
|   22|1318|
|   16|2511|
|    5|2868|
|   18|1170|
|    Y| 523|
|   17|2995|
|   MT|  37|
|    6|2863|
|   19|2926|
|    X|2359|
|    9|2242|
|    1|5194|
|   20|1386|
|   10|2204|
|    4|2505|
+-----+----+
only showing top 20 rows

好啦,SparkSQL 的结果已经只有20+行了,此时可以收进内存里面了。


不过 SparkSQL 的结果是个 DataFrame, R 语言倒是能直接收进去,Python 默认的数据类型,没有这个,怎么办?来,我们先抑制住重复造轮子、准备自己写一个的冲动,由于我们最开始 Import 了 pandas,这个包引入后, Python 也就支持 DataFrame 了。这里直接用SparkSQL 的 toPandas 方法,就可以得到Pandas 的 DataFrame 了:

代码块【7】:
[Python] 纯文本查看 复制代码
pd_genesInEachChr = sqlDF_genesInEachChr.toPandas()
pd_genesInEachChr.head()


结果:
[AppleScript] 纯文本查看 复制代码
Chrom        Cnt
0        7        2867
1        15        2152
2        11        3235
3        3        3010
4        8        2353




得到表了,有人要说,你最开始就讲 Jupyter 能画图,有个包叫做 seaborn 的还跟 ggplot2 一样简单,记忆力强的还念叨着 set_style(‘white’) 相当于 theme_bw(),现场画一个呗?

没问题。首先,Pandas 的DataFrame 没有R语言的 factor 这种让人又爱又恨的东西(掉过这个坑的在下面举手)。所以如果要调整顺序,得自己想办法。我就用了高阶函数做这个事情。具体大家参考 廖雪峰大神的Python 教程之 匿名函数篇高阶函数篇。简单说, 下面的 lambda 属于匿名函数,对我这种懒人而言不用写 def 定义函数了。map 是对一个列表每个值执行一个函数, reduce 把返回结果一个接在另一个尾巴上。有Python基础的注意,由于 map 返回的是 pandas 的 DataFrame 而不是 Python 默认的list,实际上 reduce 的 append 是 Pandas的append 而不是系统 append。

还不清楚的,这里写一个 shell 的同义语句:

[Bash shell] 纯文本查看 复制代码
rm input.chrSort.txt
for chr in {1..22} X Y MT
do
    grep -w ${chr} input.txt >>input.chrSort.txt
done


代码块【8】:

[Python] 纯文本查看 复制代码
l_plotChrOrder = map(lambda x: str(x), range(1, 23)) + ['X', 'Y', 'MT']
pd_genesInEachChrSort = reduce(lambda x,y: x.append(y), 
                               map(lambda x: pd_genesInEachChr[pd_genesInEachChr['Chrom'] == x], l_plotChrOrder)
                        )
pd_genesInEachChrSort.head()


结果

[AppleScript] 纯文本查看 复制代码
Chrom        Cnt
16        1        5194
24        2        3971
3        3        3010
19        4        2505
7        5        2868


代码块【9】:


[Python] 纯文本查看 复制代码
sns.barplot(data=pd_genesInEachChrSort, x="Cnt", y="Chrom")


结果:





大家看是不是实现了 ggplot2 的效果?更多例子请查看 seaborn文档



由于论坛不支持markdown编写,粘贴比较麻烦,只展示了文章一部分,完整版本见我的个人博客 http://huboqiang.cn/2017/02/16/AnalyseTheAnnotationOfHumanGenome












上一篇:直接分析肿瘤转录组数据的网站-Cancer RNA-Seq Nexus
下一篇:直接在线用TCGA的数据来做预后分析OncoLnc
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-11-21 17:44 , Processed in 0.053580 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.