如果在Python庫中有一些需要的功能找不到贴唇,那么第一反應(yīng)應(yīng)該是檢查下R中有沒有已經(jīng)實現(xiàn)的包搀绣。對于統(tǒng)計方法來說,R仍然是最完整的框架滤蝠。此外豌熄,一些生信工具也只存在于R中授嘀,大概率作為屬于Bioconductor項目的包提供物咳。
rpy2提供了從Python到R的申明式接口。如你所見蹄皱,你可以編寫非常優(yōu)雅的Python代碼實現(xiàn)與R的交互過程览闰。為了展示這種交互(并且嘗試最常見的R數(shù)據(jù)結(jié)構(gòu)DataFrame和最流行的R包ggplot2),我們將下載一些Human 1000Genomes項目的元信息巷折。這不是本關(guān)于R的書压鉴,但我們想提供有趣和有用的示例。
Getting ready
你將需要1000Genomes項目的元信息锻拘。請檢查git倉庫的Datasets.ipynb 和下載sequence.index文件油吭。如果使用Jupyter Notebook,打開Chapter01/Interfacing_R.ipynb執(zhí)行wget命令署拟。
這個文件含有關(guān)于這個項目的fastq文件信息(我們將在接下來的章節(jié)使用來自1000Genomes項目的數(shù)據(jù))婉宰。其中包含fastq文件,樣本ID推穷,群體來源和每個lane的重要統(tǒng)計信息(比如測序reads數(shù)據(jù)和DNA堿基數(shù))心包。
How to do it
- 加載需要的庫
import os
from IPython.display import Image
import rpy2.robjects as robjects
import rpy2.robjects.lib.ggplot2 as ggplot2
from rpy2.robjects.functions import SignatureTranslatedFunction
import pandas as pd
from rpy2.robjects import pandas2ri
from rpy2.robjects import default_converter
from rpy2.robjects.conversion import localconverter
需要注意的是直接conda安裝的rpy2似乎比較舊,不支持高版本的ggplot2馒铃,可以卸載重裝
- 我們將使用R的 read.delim功能從文件中讀取數(shù)據(jù)
read_delim = robjects.r('read.delim')
seq_data = read_delim('sequence.index', header=True, stringsAsFactors=False)
#In R:
# seq.data <- read.delim('sequence.index', header=TRUE, stringsAsFactors=FALSE)
我們導(dǎo)入了所需包之后就可以使用R的read.delim功能蟹腾,該功能允許使用者讀取文件。R語言允許對象的名字中加點区宇。因此娃殖,我們需要把函數(shù)名字改為read_delim。然后议谷,我們說函數(shù)名是正確的珊随;注意一下高度申明的特性。首先,大部分原子對象叶洞,如字符串鲫凶,可以不轉(zhuǎn)換進行傳遞。其次衩辟,參數(shù)名是可以無縫轉(zhuǎn)換的(除了點的問題)螟炫。最后,對象在 Python 命名空間中可用(但對象實際上在 R 命名空間中不可用艺晴;稍后會詳細介紹)昼钻。
作為參考,作者給出了相關(guān)的R代碼封寞。作者希望這種簡單的轉(zhuǎn)換很清晰然评。這個seq_data對象是DataFrame。如果你知道基礎(chǔ)的R或者pandas包狈究,你可能知道這種類型的數(shù)據(jù)結(jié)構(gòu)碗淌; 如果不了解的這些內(nèi)容的話,那么這本質(zhì)上是一個表:一系列行抖锥,其中每一列都具有相同的類型亿眠。
讓我們對這個DataFrame做一個簡單的檢查,如下所示:
print('This data frame has %d columns and %d rows' % (seq_data.ncol, seq_data.nrow))
print(seq_data.colnames)
#In R:
# print(colnames(seq.data))
# print(nrow(seq.data))
# print(ncol(seq.data))
print(seq_data.nrow)
my_cols = robjects.r.ncol(seq_data)
print(my_cols)
你可以直接調(diào)用R的函數(shù)磅废。在這個例子里纳像,我們在它們名字中沒有點情況下調(diào)用了ncol。然而拯勉,需要注意的是竟趾,它會顯示的輸出為含有26的向量[26]而不是數(shù)字26。這是由于默認情況下大部分的R操作返回的是向量宫峦。另外岔帽,R的數(shù)組索引是從1開始的,而Python是從0開始斗遏。
現(xiàn)在山卦,我們需要進行一些數(shù)據(jù)清理。例如诵次,某些列應(yīng)該被解釋為數(shù)字账蓉,但是按照字符串讀了進來:
as_integer = robjects.r('as.integer')
match = robjects.r.match
my_col = match('READ_COUNT', seq_data.colnames)[0] # Vector returned
print('Type of read count before as.integer: %s' % seq_data[my_col - 1].rclass[0])
seq_data[my_col - 1] = as_integer(seq_data[my_col - 1])
print('Type of read count after as.integer: %s' % seq_data[my_col - 1].rclass[0])
match 函數(shù)有點類似于 Python 列表中的 index 方法。 正如預(yù)期的那樣逾一,它返回一個向量铸本,以便我們可以提取。它的索引同樣是從1開始的遵堵,所以我們在Python中使用的時候進行了減一箱玷。這個as_integer函數(shù)可以將列轉(zhuǎn)換為整數(shù)怨规。
關(guān)于我們基于人類1000 Genomes 項目的具體例子,我們將首先通過直方圖展示測序中心的名字锡足,我們使用ggplot來作圖波丰。
ggplot2.theme = SignatureTranslatedFunction(ggplot2.theme,
init_prm_translate = {'axis_text_x': 'axis.text.x'})
bar = ggplot2.ggplot(seq_data) + ggplot2.geom_bar() + ggplot2.aes_string(x='CENTER_NAME') + ggplot2.theme(axis_text_x=ggplot2.element_text(angle=90, hjust=1))
robjects.r.png('out.png')
bar.plot()
dev_off = robjects.r('dev.off')
dev_off()
Image(filename='out.png')