R 里使用 python 加速 R 代码运行
随缘
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
4结尾
人生艰难呐...
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群
哦,数据代码已上传至QQ群,欢迎加入下载。
群二维码:
老俊俊微信:
知识星球:
所以今天你学习了吗?
今天的分享就到这里了,敬请期待下一篇!
最后欢迎大家分享转发,您的点赞是对我的鼓励和肯定!
如果觉得对您帮助很大,赏杯快乐水喝喝吧!
往期回顾
◀南京的春
◀pysam 读取 bam 文件准备 Ribo-seq 质控数据
◀...