編者注：在上篇文章《沒有自己的伺服器如何學習生物資料分析》上篇，我們對 IBM 雲端計算平臺有了基本瞭解，也學習瞭如何對資料進行下載上傳以及基本的預處理。 在《沒有自己的伺服器如何學習生物資料分析》下篇，我們將繼續跟隨作者的腳步學習如何利用IBM雲端計算平臺處理實際的生物學資料分析問題。題目來自生信技能樹論壇，論壇網址：http://biotrainee.com/forum.php/ 如果你沒有看過上篇內容，建議你先去閱讀沒有自己的伺服器如何學習生物資料分析（上篇） 祝閱讀愉快，下面是文章正文！

首先思考一下提問者的幾個問題：

• 每條染色體基因個數的分佈？
• 所有基因平均有多少個轉錄本？
• 所有轉錄本平均有多個exon和intron？

那如何將這幾句話翻譯成 SQL 語句呢

• 每條染色體基因個數的分佈？

思考一下，問的其實是：

每個 Chrom 值，對應 幾種、不重複的Gene？

    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:

    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語句分析了。

更多內容參考文件http://spark.apache.org/docs/2.0.2/sql-programming-guide.html

程式碼塊【6】：

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()

結果：

執行過程時間有點長，請耐心等待。

因為 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】：

pd_genesInEachChr = sqlDF_genesInEachChr.toPandas()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 的同義語句：

rm input.chrSort.txtfor chr in {1..22} X Y MTdo    grep -w ${chr} input.txt >>input.chrSort.txtdone

程式碼塊【8】：

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()

結果：

程式碼塊【9】：

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

結果：

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

大家看是不是實現了 ggplot2 的效果？更多例子請檢視 seaborn文件。

OK，快速解決剩下來的問題：

• 所有基因平均有多少個轉錄本？

程式碼塊【10】：

sqlDF_transInEachGene = spark.sql("""    SELECT Gene, COUNT(DISTINCT(Tran)) AS Cnt FROM GTF    WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")    GROUP BY Gene""")pd_transInEachGene = sqlDF_transInEachGene.toPandas()sns.distplot(pd_transInEachGene['Cnt'])

結果：

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

• 圖片展示效果調整

畫好了，拿給老闆看，這個肯定得捱罵，不好看啊，長尾效應太明顯。Python 作圖微調效果如何？好用的話你畫一個 0~25 的柱狀分佈唄？

既然要微調，我就用原始的python 作圖 matplotlib 庫了，他和 seaborn 的關係如同 R 的 plot 與 ggolot2。

matplotlib 庫有非常精美的 gallery，程式碼拿來就能在jupyter上畫，再次強調，如果不顯示影象請像我這裡最開始import 一樣加 %matplotlib inline 魔法。

畫之前先簡單看下資料分佈，類似 R 的 summary

程式碼塊【11】：

pd_transInEachGene['Cnt'].describe()

結果：

    count    57992.000000    mean         3.413143    std          5.103533    min          1.000000    25%          1.000000    50%          1.000000    75%          4.000000    max        170.000000    Name: Cnt, dtype: float64

程式碼塊【12】：

plt.hist(pd_transInEachGene['Cnt'], max(pd_transInEachGene['Cnt']), histtype='bar')plt.xlim(0, 25)

結果：

    (0, 25)

OK，這個應該拿給老闆不會捱罵了，下一個問題：

• 所有轉錄本平均有多個exon？

程式碼塊【13】：

sqlDF_exonsInEachTran = spark.sql("""        SELECT Tran, COUNT(DISTINCT(Exon))  AS Cnt FROM GTF        WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")        GROUP BY Tran""")pd_exonsInEachTran = sqlDF_exonsInEachTran.toPandas()pd_exonsInEachTran.head()

結果：

程式碼塊【14】：

print("Median Value %d " % (pd_exonsInEachTran.median(0)))plt.hist(pd_exonsInEachTran['Cnt'], max(pd_exonsInEachTran['Cnt']), histtype='bar')plt.xlim(0, 15)

結果：

    Median Value 4     (0, 15)

• 老闆覺得似乎不對，想起了什麼……

“這裡看的是所有基因，我要你看的是編碼基因。那麼多非編碼的看他幹嘛！除了把中位數往下帶跑偏影響結果，有什麼意義？！”

此時，聽話的人就直接 grep protein_coding 去了。而對此，我認為，如果長期以往，只能一直做菜鳥。我們要多長一個心眼，裡面不還有 lincRNA 嘛，也挺重要，萬一老闆哪天讓我比一下lincRNA 和編碼的，我是不是還得再算一次？萬一又要看其他的呢？

防止這種情況，很簡單，把基因型別那一列加進去，分不同基因類別，全算出來放那裡就好了。

如果是用 perl 的 hash表做這件事，就會出來個似乎是（原諒我幾年不寫perl全忘光了）這樣的資料架構：

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

相信有過這種噩夢般經歷的讀者此時會是懵逼的。哪地方該有括號，用 $ @ 還是%，小駱駝根本就沒有，寫錯一個就報錯，想深入學習，要麼去看大神的程式碼，要麼就得去看一本叫做 《Perl高階程式設計》的書，京東購買連結 在這裡，點開發現無貨的別急，這本書我幾年前學這個的時候，就早已斷貨了。

Python 就沒有這麼多規矩，我最早就為的這個轉的 python。

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

當然我們現在有了 pyspark，更不用去折騰 Hash 結構去了，直接在 SQL 裡，說人話。

程式碼塊【15】：

pat_gene = '''gene_ids+"(S+)";'''pat_tran = '''transcript_ids+"(S+)";'''pat_exon = '''exon_numbers+"*(w+)"*'''pat_type = '''gene_biotypes+"*(w+)"*'''pattern_gene = re.compile( pat_gene )pattern_tran = re.compile( pat_tran )pattern_exon = re.compile( pat_exon )pattern_type = re.compile( pat_type )def parseEachLineV2(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] )    match_type = pattern_type.search( f_line[-1] )    gene = "NULL"    tran = "NULL"    exon = "NULL"    gtype = "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)    if match_type:        gtype = match_type.group(1)    return [gene, tran, exon, gtype,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: parseEachLineV2(x))rdd.take(5)

結果：

    [[u'ENSG00000223972',      'NULL',      'NULL',      u'transcribed_unprocessed_pseudogene',      u'1'],     [u'ENSG00000223972',      u'ENST00000456328',      'NULL',      u'transcribed_unprocessed_pseudogene',      u'1'],     [u'ENSG00000223972',      u'ENST00000456328',      u'1',      u'transcribed_unprocessed_pseudogene',      u'1'],     [u'ENSG00000223972',      u'ENST00000456328',      u'2',      u'transcribed_unprocessed_pseudogene',      u'1'],     [u'ENSG00000223972',      u'ENST00000456328',      u'3',      u'transcribed_unprocessed_pseudogene',      u'1']]

程式碼塊【16】：

from pyspark.sql.types import *schema2=StructType(    [StructField("Gene",  StringType())] +     [StructField("Tran",  StringType())] +     [StructField("Exon",  StringType())] +    [StructField("Type",  StringType())] +    [StructField("Chrom", StringType())])df2 = sqlCtx.createDataFrame(rdd, schema2)df2.show()

結果：

    +---------------+---------------+----+--------------------+-----+    |           Gene|           Tran|Exon|        Type|Chrom|    +---------------+---------------+----+--------------------+-----+    |ENSG00000223972|           NULL|NULL|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000456328|NULL|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000456328|   1|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000456328|   2|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000456328|   3|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|NULL|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|   1|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|   2|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|   3|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|   4|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|   5|transcribed_unpro...|    1|    |ENSG00000223972|ENST00000450305|   6|transcribed_unpro...|    1|    |ENSG00000227232|           NULL|NULL|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|NULL|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|   1|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|   2|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|   3|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|   4|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|   5|unprocessed_pseud...|    1|    |ENSG00000227232|ENST00000488147|   6|unprocessed_pseud...|    1|    +---------------+---------------+----+--------------------+-----+    only showing top 20 rows

OK, 新的一列成果加入表格，然後寫SQL 分析資料。既然要看各種基因型別、每個轉錄本有幾種外顯子，那麼 GROUP BY 就加一個 Type 列，SELECT 也加一個 Type 列顯示出來。

程式碼塊【17】：

df2.registerTempTable("GTF2")sqlDF_exonsInEachTran = spark.sql("""        SELECT Tran, Type, COUNT(DISTINCT(Exon))  AS Cnt FROM GTF2        WHERE (Gene != "NULL")  AND (Tran != "NULL") AND (Exon != "NULL")        GROUP BY Tran, Type""")pd_exonsInEachTran = sqlDF_exonsInEachTran.toPandas()pd_exonsInEachTran.head()

結果：

Pandas 也可以進行分組資訊統計，如同 R 的 ddply。

程式碼塊【18】：

pd_sort = pd_exonsInEachTran[['Type', 'Cnt']].groupby(['Type'])                    .agg([len,np.median, np.mean, np.std])['Cnt']                    .sort_values(['median'], ascending=False)pd_sort.head()

結果：

再排序畫圖看看

程式碼塊【19】：

pd_exonsInEachTran_sort = reduce(lambda x,y: x.append(y),     map(lambda x: pd_exonsInEachTran[pd_exonsInEachTran['Type']==x], pd_sort.index[0:10])                        )pd_exonsInEachTran_sort.head()

結果：

最後畫一個複雜點的啊

程式碼塊【20】：

fig = plt.figure(figsize=(8, 8))ax1 = fig.add_subplot(121)ax2 = fig.add_subplot(122)sns.boxplot(data=pd_exonsInEachTran_sort, y="Type", x="Cnt", ax=ax1)sns.violinplot(data=pd_exonsInEachTran_sort, y="Type", x="Cnt", ax=ax2, bw=0.5)ax1.set_xlim(0, 20)ax2.set_xlim(-5, 60)ax2.set_yticklabels([])ax2.set_xticks([0, 20, 40, 60])ax2.set_xticklabels([0, 20, 40, 60])ax1.set_title("Boxplot")ax2.set_title("Violinplot")

結果：

    <matplotlib.text.Text at 0x7fb19e203050>

IBM data science 上程式碼：

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

編者寫在最後：

通過《沒有自己的伺服器如何學習生物資料分析》（點選連結閱讀上篇）上下兩篇文章，我們為大家介紹了IBM大資料計算平臺相關知識，同時也用一個簡單的例項告訴大家如何上手進行分析。

作為初學者的你，看到這樣的長篇難免一頭霧水。

但是萬事最難在開頭，如果通過這樣一篇文章能夠激發你瞭解大資料相關知識的熱情，我們的目的也就達到了。

隨著各種雲端計算平臺的快速發展以及網際網路巨頭公司的大力投入，包括生物資訊從業者在內的各種資料分析工程師現在幾乎可以不受本地工作環境和條件的制約。

只要你有心，只要你願意，就可以學習自己想學的東西。

共勉！