Julia 语言在生信上的一个简单应用

前言

前年(2018 年)8 月份知乎曾经小小地热烈地讨论了 Julia 1.0 的发布。

机器之心:MIT发布Julia 1.0:Python、R、C++三合一

Julia 的开发者之一、就职于 MIT 计算机科学与人工智能实验室(CSAIL)的教授 Alan Edelman 表示:“Julia 1.0 的发布证明,该语言已经做好准备,将 Python 和 R 的高效性和易用性与 C++的闪电速度结合在一起,改变技术世界。”

前一段时间本专栏也有同学写了相关文章:

徐寅生:R、Python与Julia基础简介与入门

而今天我主要对比一下我在工作中遇到一个很简单的例子来直观感受下 Julia 的特点。

到底官方的宣传有几分真实。

任务

任务本身特别简单,就是读取两个 bed 文件,将存在 overlap 的项目记录下来。

要求记录的信息包括染色体号、两个 bed 文件中的 ID,overlap 起止点,overlap 长度。

不同语言的实现

由于我日常主要使用 R 语言进行工作这里就对近乎同样算法的 R 和 Julia 程序进行比较。本身两个 bed 文件都很大,这里为了测试所以只各取 1000 个项目进行比对。

R 语言实现

rm(list = ls())
# load data
NONCODE_data <- read.csv(file = "./NONCODE_t.bed", sep = "\t", header = FALSE)
common_p_s_data <- read.csv("./common_p_s.bed", sep = "\t", header = FALSE)
NONCODE_data <- NONCODE_data[1:10000, ]
common_p_s_data <- common_p_s_data[1:10000, ]
get_peak <- function() {
  chr_ids <- vector("character")
  nct_ids <- vector("character")
  peak_ids <- vector("character")
  starts <- vector("integer")
  ends <- vector("integer")
  widths <- vector("integer")
  for (i in 1:10000) {
    ago_sdf <- common_p_s_data[i, ]
    qry_start <- ago_sdf$V2
    qry_end <- ago_sdf$V3
    qry_chr <- as.character(ago_sdf$V1)
    qry_strand <- as.character(ago_sdf$V6)
    tar_df <- NONCODE_data[(NONCODE_data$V1 == qry_chr) &
      (NONCODE_data$V6 == qry_strand), ]
    for (j in seq_along(row.names(tar_df))) {
      tar_start <- tar_df[j, ]$V2
      tar_end <- tar_df[j, ]$V3
      if ((qry_start < tar_end) & (qry_end > tar_start)) {
        vector_tmp <- c(qry_start, qry_end, tar_start, tar_end)
        vector_tmp <- sort(vector_tmp)
        chr_ids <- c(chr_ids, as.character(tar_df[j, ]$V1))
        nct_ids <- c(nct_ids, as.character(tar_df[j, ]$V4))
        peak_ids <- c(peak_ids, as.integer(ago_sdf$V4))
        starts <- c(starts, as.integer(vector_tmp[2]))
        ends <- c(ends, as.integer(vector_tmp[3]))
        widths <- c(widths, as.integer(vector_tmp[3] - vector_tmp[2]))
      }
    }
  }
  peak_df <- data.frame(chr_ids, nct_ids, peak_ids, starts, ends, widths)
  return(peak_df)
}
system.time({
  peak_df <- get_peak()
})
#   user  system elapsed 
#1949.80    9.80 1960.86 

最后耗时 1960.86 秒。

Julia 实现

using Queryverse
noncode_df = load("../NONCODE_t.csv", header_exists=false) |> DataFrame
common_df = load("../common_p_s.tsv", delim = '\t',header_exists = false) |> DataFrame
noncode_df = noncode_df[1:10000,1:6]
common_df = common_df[1:10000,:];
function get_peak_new(nc_df,cm_df)
    chr_ids = String[]
    nct_ids = String[]
    peak_ids = String[]
    starts = Int64[]
    ends = Int64[]
    widths = Int64[]
    for i in range(1,10000,step=1)
        ago_sdf = cm_df[i,:]
        qry_start = ago_sdf[2,]
        qry_end = ago_sdf[3,]
        qry_chr = ago_sdf[1,]
        qry_strand = ago_sdf[6,]
        tar_df = nc_df |> @filter(_.Column1 == qry_chr && _.Column6 == qry_strand) |> DataFrame
        for j in range(1,nrow(tar_df),step=1)
            tar_start = tar_df[j,2]
            tar_end = tar_df[j,3]
            if tar_start < qry_end && tar_end > qry_start
                vector_tmp = Int64[qry_start,qry_end,tar_start,tar_end]
                vector_tmp = sort(vector_tmp)
                push!(chr_ids,qry_chr)
                push!(nct_ids,tar_df[j,4])
                push!(peak_ids,ago_sdf[4])
                push!(starts,vector_tmp[2])
                push!(ends,vector_tmp[3])
                push!(widths,(vector_tmp[3] - vector_tmp[2]))
            end
        end
    end
    peak_df = DataFrame(Chr_ID = chr_ids, Nct_ID = nct_ids, PEAK_ID = peak_ids,
        START = starts, END = ends, WIDTH = widths)
    return(peak_df)
end
@time peak_df = get_peak_new(noncode_df,common_df)
# 8.938335 seconds (240.79 M allocations: 10.647 GiB, 19.22% gc time)

最终耗时 8.938335 秒。几乎一样的算法,Julia 程序的运行速度为 R 程序的 220 倍左右(1960.86/8.94≈219.34)。

Julia 的特点

我们可以很明显地看到,语法上 Julia 程序和 R 的差别其实并不大,差别最大应该就是下面这句:

tar_df = nc_df |> @filter(_.Column1 == qry_chr && _.Column6 == qry_strand) |> DataFrame

但其实只是因为我用了 Queryverse.jl,Julia 中一个模仿 R 的 Tidyverse 的包合集。如果在 R 中使用 Tidyverse,也可以写出类似的代码。
所以 Julia 是一门语法和 R 或者 Python 差别都不大(会更接近 Python 些),运行速度却相当快的语言。当然 Julia 还有很多高级的特性,如果你想,写得比你之前的 R 或者 Python 程序更复杂完全不是问题。但是哪怕你只是写出“Julia 描述的 R 程序”或者“Julia 描述的 Python 程序”也能获得不错的性能上的收益,且这样写需要的学习成本很低。

生信分析人员学习 Julia 的必要性

从上面的例子我们可以很清晰地看到,Julia 语法简单,性能强大。那么是不是所有生信分析人员都有必要学习 Julia 呢?这要分情况。

在读硕士生和本科生

这些人不用花太多精力去学习,可以尝试了解和学习,但没有较强的必要性。在这个阶段,大部分人所接触的数据量,程序运行的速度并不会成为项目的瓶颈。项目本身的思路(哪怕思路完全是老师或者师兄师姐给的,你也需要理解)可能远比运算这一步骤消耗更多的时间。如果遇到需要自己编写程序,又无法接受 R 的慢速的情况,其实反而建议学习一下 Python 科学计算三件套(NumPy、SciPy 和 Pandas)。因为 Python 是必学的东西,再花点精力学习下它的科学计算库也不是什么问题。而且这样也可以相对于 R 程序获得大几十倍的速度提升。

在读博士

博士的项目数据量相对于硕士,遇到需要追求更高的运行速度的可能性就大为增加,所以学习 Julia 的必要性比本科生和硕士生要高不少。

生信从业人员

学!只要你的工作会涉及到自己写程序进行计算。早点下班,保护头发不好吗?完成项目快点,多分成不香吗?

Julia 语言的学习方法

初学者建议阅读《Julia 编程基础》,目前作者刚刚把初稿交给出版社编辑,面世时间未知。不过作者慷慨的在 GitHub 上开源了前 12 章,作为入门可以说绰绰有余。

image.png

hyper0x/JuliaBasics

已经有一定使用经验的情况下,想要进阶的同学建议阅读《Julia 语言程序设计》,对大量细节和高级的应用有着清晰的讲解。

image.png

《Julia语言程序设计》

对了,在这里还是要再次推荐下《利用 Python 进行数据分析》,作者就是 Pandas 创始人。对 Python 在科学计算上的应用讲得很不错。在使用 Python 更方便的场景下,会这本书所讲解的技能绝对大有脾益。

image.png

《利用Python进行数据分析》

可能是废话的个人感言

之前我在很多推荐 Julia 文章下面看到主要两种不赞成学习 Julia 的观点。

一是“我选择 Numba、Cython”。

Numba 并不支持很多 Python 的操作,限制较大。至于 Cython,其实你如果想完全发挥它的威力,是需要写一些 C 的东西。C 和 R、Python 的差异要远大于 Julia 和 R、Python。比如这篇文章里的 Julia 程序,我就是对着原有的 R 程序,边看 Julia 程序文档边写出来的。加上调试,前后花了 10 个小时。而我之前对于 Julia 的了解仅限于怎么输出“Hello World”。如果同样的起点,用 Cython 显然是不可能的。

二是“Julia 好复杂”。

说实话我个人认为这和 Julia 中文社区的规模有很大的关系。Julia 这门语言的核心开发者是搞科学计算的,这门语言最大的优势也在于科学计算。所以初期用户必然是以高阶用户为主。我曾经看 Julia 中文社区的年会直播,讲的东西专业性特别强,涉及大量早就忘却的高等数学知识,听了是一头雾水。但其实,Julia 也可以写得很简单,用来处理没那么复杂的问题。比如本文的例子,其实初中数学知识就足够理解。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 217,734评论 6 505
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,931评论 3 394
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 164,133评论 0 354
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,532评论 1 293
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,585评论 6 392
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,462评论 1 302
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,262评论 3 418
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 39,153评论 0 276
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,587评论 1 314
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,792评论 3 336
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,919评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,635评论 5 345
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,237评论 3 329
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,855评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,983评论 1 269
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 48,048评论 3 370
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,864评论 2 354

推荐阅读更多精彩内容