R 里使用 python 加速 R 代码运行

Original JunJunLab 老俊俊的生信笔记 2022-08-15

收录于合集 #R操作小技巧（Tips） 22个

随缘

1引言

R 里面运行大数据实在是非常鸡肋, 当你使用 lappy, future.apply, llply, purrr 等函数依然无法提高你的 R 代码时,是不是非常绝望?

你可以把数据导出去,然后使用 python 代码进行分析,但是我又想用 R 来可视化呢?要再导入进来吗?

这样来回切换确实有点麻烦,不过 reticulate 这个包给 Rstudio 兼容了 python 代码,方便我们在 Rstudio 里进行切换,今天尝试了一下,可以加速我们的数据分析过程。

大概就是把 R 里面的变量从 python 里读入。
然后写 python 代码处理得到结果。
在 R 里导入 python 分析的结果。

因为变量都在 R 里面,所以可以比较方便的调用!

2起因

下面是我的代码:

# normalize by meanNorm for every gene
system.time({llply(unique(mean_density_st$type),function(x){
  tmp = mean_density_st[type %in% x]
  llply(unique(tmp$rname)[1:5],function(z){
    tmp1 = tmp[rname %in% z]
    meanNorm <- end5_gene_meanNorm[(type %in% x & rname %in% z),.(meanNorm)] %>% as.numeric()
    tmp1$norm_explen <- tmp1$rpm/meanNorm
    return(tmp1)
  }) %>% do.call('rbind',.) -> res
  return(res)
}) %>% do.call('rbind',.) %>% data.table() %>%
    .[,.(mean_read_density = sum(norm_explen),ngenes = .N),by = .(type,rel_st)] -> mean_read_density_st
})

需要运行 70003.44s(116 分钟), 我肯定是不想等这么久的,本身数据可能也比较大,有 3372079 行。

输入的两个数据为两个表格:

end5_gene_meanNorm

mean_read_density_st

3在 R 里写个 python 代码

加载这个包:

library(reticulate)

新建一个 python 脚本,把这两个表格变量加载到 python 里面,使用 r.变量:

# load into python
end5_gene_meanNorm = r.end5_gene_meanNorm
mean_density_st = r.mean_density_st

然后写一些和 R 一样完成任务的代码:

# 1. end5_gene_meanNorm save in dict
gene_meanNorm = {}
for type,rname,meanNorm in zip(list(end5_gene_meanNorm['type']),list(end5_gene_meanNorm['rname']),list(end5_gene_meanNorm['meanNorm'])):
  # print(type,rname,meanNorm)
  key = '|'.join([type,rname])
  gene_meanNorm[key] = meanNorm

# 2. mean_density_st save in dict
density_st = {}
for type,rel_st,rname,rpm in zip(list(mean_density_st['type']),list(mean_density_st['rel_st']),list(mean_density_st['rname']),list(mean_density_st['rpm'])):
  key = '|'.join([type,str(rel_st),rname])
  density_st[key] = rpm

# 3. save normaled value and gene counts
norm_expleninfo = {}
count = {}
for key,val in density_st.items():
  filed = key.split(sep='|')
  new_id = '|'.join([filed[0],filed[2]])
  if new_id in gene_meanNorm:
    norm_explen = val/gene_meanNorm[new_id]
    fin_id = '|'.join([filed[0],filed[1]])
    norm_expleninfo.setdefault(fin_id,0)
    norm_expleninfo[fin_id] += norm_explen
    count.setdefault(fin_id,0)
    count[fin_id] += 1

# 4. final results
final_res = {}
for key,val in norm_expleninfo.items():
  mean_density = val/count[key]
  final_res[key] = [val,count[key],mean_density]

# 5.make it to dataframe format
type = [];rel_st = [];density = [];ngenes =[];normed_density = []

for key,val in final_res.items():
  fileds = key.split(sep='|')
  type.append(fileds[0])
  rel_st.append(fileds[1])
  density.append(val[0])
  ngenes.append(val[1])
  normed_density.append(val[2])

大概四五分钟就结束了, 最后将字典转为数据框:

# 6. to dataframe
dict_info = {'type':type,'rel_st':rel_st,'density':density,'ngenes':ngenes,'normed_density':normed_density}

import pandas as pd
final_res = pd.DataFrame.from_dict(dict_info)

点击可以直接打开:

注意此时我们是在 python 环境里,R 里面环境变量并没有这个变量

所以我们需要加其加载到 R 里面,使用 py$final_res即可:

# 7. load into R
final_res_st <- py$final_res

这样就可以继续做后面的分析了,画个图看看:

# plot meata-gene
pst <- ggplot(final_res_st,
              aes(x = as.numeric(rel_st),y = normed_density)) +
  geom_line(aes(color = type)) +
  theme_classic(base_size = 16) +
  ylab('Mean read density [AU]') +
  xlab('Distance from start codon (nt)') +
  scale_color_d3(name = '') +
  theme(legend.position = c(0.85,0.85),
        legend.background = element_blank())

pst