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

编者注：在上篇文章《没有自己的服务器如何学习生物数据分析》上篇，我们对 IBM 云计算平台有了基本了解，也学习了如何对数据进行下载上传以及基本的预处理。

在《没有自己的服务器如何学习生物数据分析》下篇，我们将继续跟随作者的脚步学习如何利用IBM云计算平台处理实际的生物学数据分析问题。题目来自生信技能树论坛，论坛网址：http://biotrainee.com/forum.php/

如果你没有看过上篇内容，建议你先去阅读没有自己的服务器如何学习生物数据分析（上篇）

祝阅读愉快，下面是文章正文！

首先思考一下提问者的几个问题：

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

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

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

那如何将这几句话翻译成 SQL 语句呢

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

思考一下，问的其实是：

每个 Chrom 值，对应 几种、不重复的Gene？

  1.    SELECT Chrom, COUNT(DISTINCT(Gene)) FROM GTF

  2.    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:

  1.    SELECT Chrom, COUNT(DISTINCT(Gene)) FROM GTF

  2.    WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")

  3.    GROUP BY Chrom

了解更多 SQL 可以在这个网站 https://www.w3schools.com/sql/ 在线学习基本的SQL语句。

SQL 语句在调取 UCSC 数据集中同样作用巨大，具体这里不表。

这句话怎么在 DataFrame 执行？需要先 registerTempTable("GTF")把 df 这个 dataFrame 给 SparkSQL，取一个表名，叫做 “GTF”。这样 df 就可以直接用 SQL语句分析了。

更多内容参考文档http://spark.apache.org/docs/2.0.2/sql-programming-guide.html

代码块【6】：

  1. df.registerTempTable("GTF")

  2. sqlDF_genesInEachChr = spark.sql("""

  3.    SELECT Chrom, COUNT(DISTINCT(Gene)) AS Cnt FROM GTF

  4.    WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")

  5.    GROUP BY Chrom

  6. """)

  7. sqlDF_genesInEachChr.show()

结果：

运行过程时间有点长，请耐心等待。

因为 IBM 的免费机器是 2 核心单机模式，体现不出 Spark 大数据分析的威力。如果你在Spark集群模式下，几台 48 线程的机器上对一个大文件执行SparkSQL（前提是没人使用 + 满CPU使用），在等待的过程中去后台 top 一下，会看见计算节点上全部都是恐怖的 4800% 的 CPU 使用率，共同执行同一个任务。

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

不过 SparkSQL 的结果是个 DataFrame, R 语言倒是能直接收进去，Python 默认的数据类型，没有这个，怎么办？

来，我们先抑制住重复造轮子、准备自己写一个的冲动，由于我们最开始 Import 了 pandas，这个包引入后， Python 也就支持 DataFrame 了。这里直接用SparkSQL 的 toPandas 方法，就可以得到Pandas 的 DataFrame 了：

代码块【7】：

  1. pd_genesInEachChr = sqlDF_genesInEachChr.toPandas()

  2. pd_genesInEachChr.head()

结果：

  • 画图

得到表了，有人要说，你最开始就讲* Jupyter 能画图*，有个包叫做 seaborn 的还跟 ggplot2 一样简单，记忆力强的还念叨着 setstyle('white') 相当于 themebw()，现场画一个呗？

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

还不清楚的，这里写一个 shell 的同义语句：

  1. rm input.chrSort.txt

  2. for chr in {1..22} X Y MT

  3. do

  4.    grep -w ${chr} input.txt >>input.chrSort.txt

  5. done

代码块【8】：

  1. l_plotChrOrder = map(lambda x: str(x), range(1, 23)) + ['X', 'Y', 'MT']

  2. pd_genesInEachChrSort = reduce(lambda x,y: x.append(y),

  3.      map(lambda x: pd_genesInEachChr[pd_genesInEachChr['Chrom'] == x], l_plotChrOrder)

  4.          )

  5. pd_genesInEachChrSort.head()

结果：

代码块【9】：

  1. sns.barplot(data=pd_genesInEachChrSort, x="Cnt", y="Chrom")

结果：

  1.    <matplotlib.axes._subplots.AxesSubplot at 0x7fb1a31e9110>

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

OK，快速解决剩下来的问题：

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

代码块【10】：

  1. sqlDF_transInEachGene = spark.sql("""

  2.    SELECT Gene, COUNT(DISTINCT(Tran)) AS Cnt FROM GTF

  3.    WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")

  4.    GROUP BY Gene

  5. """)

  6. pd_transInEachGene = sqlDF_transInEachGene.toPandas()

  7. sns.distplot(pd_transInEachGene['Cnt'])

结果：

  1.    <matplotlib.axes._subplots.AxesSubplot at 0x7fb1a09a6bd0>

  • 图片展示效果调整

画好了，拿给老板看，这个肯定得挨骂，不好看啊，长尾效应太明显。Python 作图微调效果如何？好用的话你画一个 0~25 的柱状分布呗？

既然要微调，我就用原始的python 作图 matplotlib 库了，他和 seaborn 的关系如同 R 的 plot 与 ggolot2。

matplotlib 库有非常精美的 gallery，代码拿来就能在jupyter上画，再次强调，如果不显示图像请像我这里最开始import 一样加 %matplotlib inline 魔法。

画之前先简单看下数据分布，类似 R 的 summary

代码块【11】：

  1. pd_transInEachGene['Cnt'].describe()

结果：

  1.    count    57992.000000

  2.    mean         3.413143

  3.    std          5.103533

  4.    min          1.000000

  5.    25%          1.000000

  6.    50%          1.000000

  7.    75%          4.000000

  8.    max        170.000000

  9.    Name: Cnt, dtype: float64

代码块【12】：

  1. plt.hist(pd_transInEachGene['Cnt'], max(pd_transInEachGene['Cnt']), histtype='bar')

  2. plt.xlim(0, 25)

结果：

  1.    (0, 25)

OK，这个应该拿给老板不会挨骂了，下一个问题：

  • 所有转录本平均有多个exon？

代码块【13】：

  1. sqlDF_exonsInEachTran = spark.sql("""

  2.        SELECT Tran, COUNT(DISTINCT(Exon))  AS Cnt FROM GTF

  3.        WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")

  4.        GROUP BY Tran

  5. """)

  6. pd_exonsInEachTran = sqlDF_exonsInEachTran.toPandas()

  7. pd_exonsInEachTran.head()

结果：

代码块【14】：

  1. print("Median Value %d " % (pd_exonsInEachTran.median(0)))

  2. plt.hist(pd_exonsInEachTran['Cnt'], max(pd_exonsInEachTran['Cnt']), histtype='bar')

  3. plt.xlim(0, 15)

结果：

  1.    Median Value 4

  2.    (0, 15)

  • 老板觉得似乎不对，想起了什么……

“这里看的是所有基因，我要你看的是编码基因。那么多非编码的看他干嘛！除了把中位数往下带跑偏影响结果，有什么意义？！”

此时，听话的人就直接 grep protein_coding 去了。而对此，我认为，如果长期以往，只能一直做菜鸟。我们要多长一个心眼，里面不还有 lincRNA 嘛，也挺重要，万一老板哪天让我比一下lincRNA 和编码的，我是不是还得再算一次？万一又要看其他的呢？

防止这种情况，很简单，把基因类型那一列加进去，分不同基因类别，全算出来放那里就好了。

如果是用 perl 的 hash表做这件事，就会出来个似乎是（原谅我几年不写perl全忘光了）这样的数据架构：

  1. push(@{$TypeTranExons{$gtype}{$tran}}, $exon);

相信有过这种噩梦般经历的读者此时会是懵逼的。哪地方该有括号，用 $ @ 还是%，小骆驼根本就没有，写错一个就报错，想深入学习，要么去看大神的代码，要么就得去看一本叫做 《Perl高级编程》的书，京东购买链接 在这里，点开发现无货的别急，这本书我几年前学这个的时候，就早已断货了。

Python 就没有这么多规矩，我最早就为的这个转的 python。

  1. TypeTranExons[gtype][tran].append(exon)

当然我们现在有了 pyspark，更不用去折腾 Hash 结构去了，直接在 SQL 里，说人话。

代码块【15】：

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

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

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

  4. pat_type = '''gene_biotype\s+\"*(\w+)\"*'''

  5. pattern_gene = re.compile( pat_gene )

  6. pattern_tran = re.compile( pat_tran )

  7. pattern_exon = re.compile( pat_exon )

  8. pattern_type = re.compile( pat_type )

  9. def parseEachLineV2(f_line):

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

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

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

  13.    match_type = pattern_type.search( f_line[-1] )

  14.    gene = "NULL"

  15.    tran = "NULL"

  16.    exon = "NULL"

  17.    gtype = "NULL"

  18.    if match_gene:

  19.        gene = match_gene.group(1)

  20.    if match_tran:

  21.        tran = match_tran.group(1)

  22.    if match_exon:

  23.        exon = match_exon.group(1)

  24.    if match_type:

  25.        gtype = match_type.group(1)

  26.    return [gene, tran, exon, gtype,f_line[0]]

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

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

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

  30.            .map(lambda x: parseEachLineV2(x))

  31. rdd.take(5)

结果：

  1.    [[u'ENSG00000223972',

  2.      'NULL',

  3.      'NULL',

  4.      u'transcribed_unprocessed_pseudogene',

  5.      u'1'],

  6.     [u'ENSG00000223972',

  7.      u'ENST00000456328',

  8.      'NULL',

  9.      u'transcribed_unprocessed_pseudogene',

  10.      u'1'],

  11.     [u'ENSG00000223972',

  12.      u'ENST00000456328',

  13.      u'1',

  14.      u'transcribed_unprocessed_pseudogene',

  15.      u'1'],

  16.     [u'ENSG00000223972',

  17.      u'ENST00000456328',

  18.      u'2',

  19.      u'transcribed_unprocessed_pseudogene',

  20.      u'1'],

  21.     [u'ENSG00000223972',

  22.      u'ENST00000456328',

  23.      u'3',

  24.      u'transcribed_unprocessed_pseudogene',

  25.      u'1']]

代码块【16】：

  1. from pyspark.sql.types import *

  2. schema2=StructType(

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

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

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

  6.    [StructField("Type",  StringType())] +

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

  8. )

  9. df2 = sqlCtx.createDataFrame(rdd, schema2)

  10. df2.show()

结果：

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

  2.    |           Gene|           Tran|Exon|        Type|Chrom|

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

  4.    |ENSG00000223972|           NULL|NULL|transcribed_unpro...|    1|

  5.    |ENSG00000223972|ENST00000456328|NULL|transcribed_unpro...|    1|

  6.    |ENSG00000223972|ENST00000456328|   1|transcribed_unpro...|    1|

  7.    |ENSG00000223972|ENST00000456328|   2|transcribed_unpro...|    1|

  8.    |ENSG00000223972|ENST00000456328|   3|transcribed_unpro...|    1|

  9.    |ENSG00000223972|ENST00000450305|NULL|transcribed_unpro...|    1|

  10.    |ENSG00000223972|ENST00000450305|   1|transcribed_unpro...|    1|

  11.    |ENSG00000223972|ENST00000450305|   2|transcribed_unpro...|    1|

  12.    |ENSG00000223972|ENST00000450305|   3|transcribed_unpro...|    1|

  13.    |ENSG00000223972|ENST00000450305|   4|transcribed_unpro...|    1|

  14.    |ENSG00000223972|ENST00000450305|   5|transcribed_unpro...|    1|

  15.    |ENSG00000223972|ENST00000450305|   6|transcribed_unpro...|    1|

  16.    |ENSG00000227232|           NULL|NULL|unprocessed_pseud...|    1|

  17.    |ENSG00000227232|ENST00000488147|NULL|unprocessed_pseud...|    1|

  18.    |ENSG00000227232|ENST00000488147|   1|unprocessed_pseud...|    1|

  19.    |ENSG00000227232|ENST00000488147|   2|unprocessed_pseud...|    1|

  20.    |ENSG00000227232|ENST00000488147|   3|unprocessed_pseud...|    1|

  21.    |ENSG00000227232|ENST00000488147|   4|unprocessed_pseud...|    1|

  22.    |ENSG00000227232|ENST00000488147|   5|unprocessed_pseud...|    1|

  23.    |ENSG00000227232|ENST00000488147|   6|unprocessed_pseud...|    1|

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

  25.    only showing top 20 rows

OK, 新的一列成果加入表格，然后写SQL 分析数据。既然要看各种基因类型、每个转录本有几种外显子，那么 GROUP BY 就加一个 Type 列，SELECT 也加一个 Type 列显示出来。

代码块【17】：

  1. df2.registerTempTable("GTF2")

  2. sqlDF_exonsInEachTran = spark.sql("""

  3.        SELECT Tran, Type, COUNT(DISTINCT(Exon))  AS Cnt FROM GTF2

  4.        WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")

  5.        GROUP BY Tran, Type

  6. """)

  7. pd_exonsInEachTran = sqlDF_exonsInEachTran.toPandas()

  8. pd_exonsInEachTran.head()

结果：

Pandas 也可以进行分组信息统计，如同 R 的 ddply。

代码块【18】：

  1. pd_sort = pd_exonsInEachTran[['Type', 'Cnt']].groupby(['Type'])\

  2.                    .agg([len,np.median, np.mean, np.std])['Cnt']\

  3.                    .sort_values(['median'], ascending=False)

  4. pd_sort.head()

结果：

再排序画图看看

代码块【19】：

  1. pd_exonsInEachTran_sort = reduce(lambda x,y: x.append(y),

  2.     map(lambda x: pd_exonsInEachTran[pd_exonsInEachTran['Type']==x], pd_sort.index[0:10])

  3.                        )

  4. pd_exonsInEachTran_sort.head()

结果：

最后画一个复杂点的啊

代码块【20】：

  1. fig = plt.figure(figsize=(8, 8))

  2. ax1 = fig.add_subplot(121)

  3. ax2 = fig.add_subplot(122)

  4. sns.boxplot(data=pd_exonsInEachTran_sort, y="Type", x="Cnt", ax=ax1)

  5. sns.violinplot(data=pd_exonsInEachTran_sort, y="Type", x="Cnt", ax=ax2, bw=0.5)

  6. ax1.set_xlim(0, 20)

  7. ax2.set_xlim(-5, 60)

  8. ax2.set_yticklabels([])

  9. ax2.set_xticks([0, 20, 40, 60])

  10. ax2.set_xticklabels([0, 20, 40, 60])

  11. ax1.set_title("Boxplot")

  12. ax2.set_title("Violinplot")

结果：

  1.    <matplotlib.text.Text at 0x7fb19e203050>

IBM data science 上代码：

https://apsportal.ibm.com/analytics/notebooks/d3400624-fd7f-483b-96f3-b9d07876f455/view?access_token=499996f6a4e6f93e448907bf219bae6310975c0d02521c7c67ef02b79b1ccf77

赞 (0)

相关推荐

  • GFF和GTF的异同及相互转换

    GFF(gff)全称为:general feature format GTF(gtf)全称为:gene transfer format 前者用来注释基因组,后者用来注释基因. 异同点: GTF文件和G ...

  • 7 对比对和变异结果用IGV进行可视化

    BRCA1基因为例 1 找到BRCA1在gtf文件中的坐标 zcat /mnt/f/kelly/bioTree/server/wesproject/reference/gencode.v25.anno ...

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

    编者注:完整文章首发于作者博客 http://huboqiang.cn/ 在这篇文章中,作者利用大数据平台 IBM data science 对生信技能树论坛的一道生物信息入门题进行了分析. 由于文章 ...

  • 学习生物分类法

    "教.读.写.看"流水(2021年春季第13周) 周一(2021.5.24) 7:10到校.雨. 查了点元代科学家郭守敬的资料.将郭守敬归入科学家,这是李约瑟的做法,我们也乐得其成 ...

  • 高三学生如何有效率的学习生物

    高中生物课是学习生命活动基本规律相关基础知识的一门课程,那么该如何学习才能提高生物成绩呢?有途网小编为大家讲解一下. 学好高中生物最首要的就是要突破难点. 有些知识比较复杂,或是过于抽象,同学们学高中 ...

  • 解读《对比Excel轻松学习Python数据分析》

    内容:鹤鹤 编辑:紫岸 这本书是一本工具书,对于非专业领域的人来说,可能不是特别适合.我之前是看过一些关于写python的一些基础语法知识的书,但是在实际的工作应用中,暂时没有实操过,很多知识就忘记了 ...

  • 【玩转腾讯云】使用云服务器进行生信数据分析 | 服务器文档网

    很多小伙伴手头有生信数据分析,但苦于没有服务器,没法完成自己需要的数据分析,特别是处于学习阶段的同学.这里,向大家推荐一下使用腾讯云CVM服务器,按量计费进行数据分析.一般认为,除了公司,普通人难以负 ...

  • 跟着阿拉丁学习民生银行数据分析

    近年来,数据分析因其自身显而易见的应用价值,从而得到了迅猛发展,与此同时"数据化管理.数据化运营.数据化决策"等的管理理念也渗透到了各行各业.放眼当下,企业间的经营竞争越发剑拔弩张 ...

  • 高中该如何学习生物？生物逆袭学霸的学习技巧

    生物是理综里面最像文科的一门学科了,但是该如何学习生物也是很多学生关注的重点,那么小编就来分享几个学霸学习生物的技巧吧,希望对你学习生物有所帮助. 作者 | 纸盆 1.教材最重要 生物的很多基础知识都 ...

  • 这样学英语，想不成为牛娃都难——《做孩子最好的英语学习规划师》下篇

    自从开始关注孩子的英语启蒙,加入了一些英语学习社区,常常看到许多5.6岁的孩子就已经自主阅读原版书了,一线城市的孩子二.三年级就通过KET.PET考试的例子比比皆是.作为一个西北地区的二线城市的家长, ...

  • 2019雪球嘉年华学习笔记（下篇）

    在下午的活动开场中,雪球CEO李楠先生用数据揭秘了基金投资者亏损的一大原因--短期投资.短视投机.我国公募偏股型基金的投资者绝大部分为短期投资,只有约19%的投资者能够持有偏股型基金一年以上. 我认为 ...