用rpy2修改r对象

Question

我试图使用rpy2在python中使用DESeq2 R/生物导体包。

实际上，我在编写问题时解决了我的问题(使用do_slots允许访问r对象属性)，但我认为这个示例可能对其他人有用，下面是我在R中的操作以及在python中的转换方式：

在R中

我可以根据以下两个数据帧创建一个"DESeqDataSet“：

counts_data <- read.table("long/path/to/file",
                           header=TRUE, row.names="gene")
head(counts_data)
##       WT_RT_1 WT_RT_2 prg1_RT_1 prg1_RT_2
## aap-1     406     311        41        95
## aat-1       5       8         2         0
## aat-2       1       1         0         0
## aat-3      13      12         0         1
## aat-4       6       6         2         3
## aat-5       3       1         1         0

col_data <- DataFrame(lib = c("WT", "WT", "prg1", "prg1"),
                      treat = c("RT", "RT", "RT", "RT"),
                      rep = c("1", "2", "1", "2"), 
                      row.names = colnames(counts_data))
head(col_data)
## DataFrame with 4 rows and 3 columns
##                   lib       treat         rep
##           <character> <character> <character>
## WT_RT_1            WT          RT           1
## WT_RT_2            WT          RT           2
## prg1_RT_1        prg1          RT           1
## prg1_RT_2        prg1          RT           2
dds <- DESeqDataSetFromMatrix(countData = counts_data,
                              colData = col_data,
                              design = ~ lib)
## Warning message:
## In DESeqDataSet(se, design = design, ignoreRank) :
## some variables in design formula are characters, converting to factors

dds
## class: DESeqDataSet 
## dim: 18541 4 
## metadata(1): version
## assays(1): counts
## rownames(18541): aap-1 aat-1 ... WBGene00255550 WBGene00255553
## rowData names(0):
## colnames(4): WT_RT_1 WT_RT_2 prg1_RT_1 prg1_RT_2
## colData names(3): lib treat rep

为了确保分析使用正确的控件，我需要使用“双括号”语法访问一个可以访问的因素：

dds[["lib"]]
## [1] WT   WT   prg1 prg1
## Levels: prg1 WT

dds[["lib"]] <- relevel(dds[["lib"]], ref="WT")
dds[["lib"]]
## [1] WT   WT   prg1 prg1
## Levels: WT prg1

然后我就可以进行分析：

dds <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing

res <- results(dds)

我查看给定基因的结果：

res["his-10",]
## log2 fold change (MAP): lib prg1 vs WT 
## Wald test p-value: lib prg1 vs WT 
## DataFrame with 1 row and 6 columns
##         baseMean log2FoldChange     lfcSE      stat       pvalue        padj
##        <numeric>      <numeric> <numeric> <numeric>    <numeric>   <numeric>
## his-10  586.5464       3.136174 0.2956132  10.60904 2.705026e-26 8.78785e-25

在Python中

现在，我想在python中使用rpy2来做同样的事情。

我似乎成功地从熊猫的数据中创建了这个对象：

import pandas as pd
from rpy2.robjects import r, pandas2ri, Formula
as_df = r("as.data.frame")
from rpy2.robjects.packages import importr
deseq2 = importr("DESeq2")

counts_data = pd.read_table("long/path/to/file", index_col=0)
col_data = pd.DataFrame({
    "cond_names" : counts_data.columns,
    "lib" : ["WT", "WT", "prg1", "prg1"],
    "rep" : ["1", "1", "2", "2"],
    "treat" : ["RT", "RT", "RT", "RT"]})
col_data.set_index("cond_names", inplace=True)

pandas2ri.activate()  # makes some conversions automatic
dds = deseq2.DESeqDataSetFromMatrix(
    countData=counts_data,
    colData=col_data,
    design=Formula("~lib"))

在IPython (我实际上运行了前面的命令)中，我可以使用do_slots查看对象内部，以尝试识别需要发布的因素：

In [229]: tuple(dds.do_slot("colData").slotnames())
Out[229]: ('rownames', 'nrows', 'listData', 'elementType', 'elementMetadata', 'metadata')

In [230]: dds.do_slot("colData").do_slot("listData")
Out[230]: 
R object with classes: ('list',) mapped to:
<ListVector - Python:0x7f2ae2590a08 / R:0x108fcdd0>
[FactorVector, FactorVector, FactorVector]
  lib: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
<FactorVector - Python:0x7f2ae20f1c08 / R:0x136a3920>
[       2,        2,        1,        1]
  rep: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
<FactorVector - Python:0x7f2a9600c948 / R:0x136a30f0>
[       1,        1,        2,        2]
  treat: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
<FactorVector - Python:0x7f2a9600ccc8 / R:0x136a3588>
[       1,        1,        1,        1]

我想，relevel的因素是第一个，因为"lib“是传递给deseq2.DESeqDataSetFromMatrix函数的col_data中的第一列(编辑:我意识到"lib”实际上写在r对象的描述中)。

关于通过relevel访问的属性的do_slots似乎会产生影响：

In [231]: dds.do_slot("colData").do_slot("listData")[0] = r.relevel(dds.do_slot("colData").do_slot("listData")[0], ref="WT")

In [232]: dds.do_slot("colData").do_slot("listData")
Out[232]: 
R object with classes: ('list',) mapped to:
<ListVector - Python:0x7f2a95078508 / R:0x108fcdd0>
[FactorVector, FactorVector, FactorVector]
  lib: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
<FactorVector - Python:0x7f2a9600bb88 / R:0x12a7ff60>
[       1,        1,        2,        2]
  rep: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
<FactorVector - Python:0x7f2ae2568888 / R:0x136a30f0>
[       1,        1,        2,        2]
  treat: <class 'rpy2.robjects.vectors.FactorVector'>
  R object with classes: ('factor',) mapped to:
<FactorVector - Python:0x7f2ae2568848 / R:0x136a3588>
[       1,        1,        1,        1]

然后运行分析部分：

In [233]: dds = deseq2.DESeq(dds)
/home/bli/.local/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: estimating size factors

  warnings.warn(x, RRuntimeWarning)
/home/bli/.local/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: estimating dispersions

  warnings.warn(x, RRuntimeWarning)
/home/bli/.local/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: gene-wise dispersion estimates

  warnings.warn(x, RRuntimeWarning)
/home/bli/.local/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: mean-dispersion relationship

  warnings.warn(x, RRuntimeWarning)
/home/bli/.local/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: final dispersion estimates

  warnings.warn(x, RRuntimeWarning)
/home/bli/.local/lib/python3.6/site-packages/rpy2/rinterface/__init__.py:186: RRuntimeWarning: fitting model and testing

  warnings.warn(x, RRuntimeWarning)

In [234]: res = pandas2ri.ri2py(as_df(deseq2.results(dds)))

In [235]: res.index.names = ["gene"]

dds = deseq2.DESeq(dds)
res = pandas2ri.ri2py(as_df(deseq2.results(dds)))
res.index.names = ["gene"]

现在，检查测试基因的结果：

In [236]: res.loc["his-10"]
Out[236]: 
baseMean          5.865464e+02
log2FoldChange    3.136174e+00
lfcSE             2.956132e-01
stat              1.060904e+01
pvalue            2.705026e-26
padj              8.787850e-25
Name: his-10, dtype: float64

python返回的结果与R.

Answer 1

我在rpy2文档中找到了帮助我解决问题的代码示例：2.8.x/rinterface.html#pass-by-value-paradigm。

可以通过do_slots方法访问r对象的属性，该方法以属性名称为参数。请参阅问题中的完整解决方案。

编辑：

还有一个do_slot_assign方法，例如可以用来更改设计公式：

>>> dds.do_slot("design").r_repr()
'~lib'
>>> dds.do_slot_assign("design", Formula("~ treat"))
>>> dds.do_slot("design").r_repr()
'~treat'