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

作者: 契卡 | 来源:发表于2020-02-21 00:39 被阅读0次

    前言

    前年(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 也可以写得很简单,用来处理没那么复杂的问题。比如本文的例子,其实初中数学知识就足够理解。

    相关文章

      网友评论

        本文标题:Julia 语言在生信上的一个简单应用

        本文链接:https://www.haomeiwen.com/subject/fglqqhtx.html