python, perl 和julia的性能对比

在生物信息学中经常用到的脚本语言主要是python和perl，他们被用来处理文本，大量统计，流程控制等等，其自身也是各有优势。比如说perl天生就为了处理文本而生，但是python确是有名的胶水语言，特别在整合C代码时显示出巨大的优势，其语法简洁易懂，易于维护更让其成为仅次于C和JAVA的第三大语言，但其糟糕的性能在处理大量循环时会让人忍不住抓狂。因此，Julia语言应运而生，其控制了python中没必要的动态性，加之使用JIT技术让其能够保有高性能的同时具备简洁的语法。

说了那么多，在生物信息上我们经常需要处理大量的文本文件，例如Fasta格式的序列文件，那么三者又是谁快呢？

版本控制

• python3 = 3.8.3
• perl = 5.26.2
• julia = 1.5.0-beta
• system = centos 8

计算内容

从UCSC上下载人类参考基因组 hg38.fa.gz 并解压，计算基因组GC含量，N碱基不算在总长中。

代码

perl

#!/usr/bin/perl -w

use strict;

if(@ARGV < 1){
    die "Usage : perl $0 <genome.fa>\n";
}

my $input = shift @ARGV;

my ($sum,$G_num,$C_num,$N_num)=(0,0,0,0);
my $id;
open IN, "< $input" or die $!;
while(my $line = <IN>){
    chomp $line;
    if($line =~ />([^\s]+)/){
        $id = $1;
    }else{
        $sum += length($line);
        $G_num += ($line =~ tr/Gg/Gg/);
        $C_num += ($line =~ tr/Cc/Cc/);
        $N_num += ($line =~ tr/Nn/Nn/);
    }
}
close IN;

my $GC_rate = ($G_num+$C_num)/($sum-$N_num);
printf "GC content: %.3f \n",$GC_rate;

julia

function lineGC(seq::String)
    GCnumber=count(x->(x=='G'||x=='C'||x=='g'||x=='c'),seq)
    lineNum=count(x->(x!='N' && x!='n'),seq)
    (GCnumber,lineNum)
end

function calGC(fs)
    GCnumber=zero(Int)
    lineNum=zero(Int)
    open(fs,"r") do IOstream
        for line in eachline(IOstream)
            if startswith(line,">")
                continue
            else
                GC,all=lineGC(line)
                GCnumber+=GC
                lineNum+=all
            end
        end
    end
    round(GCnumber/lineNum;digits=3)
end

println("GC content: ",calGC(ARGS[1]))

python

import sys

def lineGC(seq):
   tmp=[base for base in seq if base =="G" or base =="g" or base == "C" or base == "c"]
   gcNumber=len(tmp)
   tmp2=[base for base in seq if base !="N" and base !="n"]
   allNumber=len(tmp2)
   return (gcNumber,allNumber)


with open(sys.argv[1],'r') as f:
    gcNum=0
    allNum=0
    for line in f:
       if line.startswith(">"):
           continue
       else:
           gc,alln=lineGC(line.strip("\n"))
           gcNum=gcNum+gc
           allNum=allNum+alln

print("GC content: {:.3f}".format(gcNum/allNum))

运行时间测试

python

python

julia

julia

perl

perl

总结

结果令人咋舌，可以从sys时间看出来python和perl都是立马启动，而julia在函数的即时编辑上花了一点时间（一半时间）。 总体用时上，julia仅比perl快了1秒，而python却用了惊人的9分钟，😭

后记

python 也不是这么不堪，想要提速还是可以有很多办法的，比如切换pypy, 或者也用正则表达式，例如：

import sys
import re


def lineGC(seq):
    pattern_1 = re.compile(r"G|C",re.I)
    pattern_2 = re.compile(r"N",re.I)
    gcNumber=len(pattern_1.findall(seq))
    allNumber=len(seq)-len(pattern_2.findall(seq))
    return (gcNumber,allNumber)


with open(sys.argv[1],'r') as f:
    gcNum=0
    allNum=0
    for line in f:
       if line.startswith(">"):
           continue
       else:
           gc,alln=lineGC(line.strip("\n"))
           gcNum=gcNum+gc
           allNum=allNum+alln

print("GC content: {:.3f}".format(gcNum/allNum))

这样计算下来，大概需要6分20秒，提速很多

另外，我们也可以使用NumPy的向量化运算来提速

import sys
from pyfaidx import Fasta
import numpy as np

def lineGC(seq):
    gc_number = np.where((seq==b'G')|(seq==b'C')|(seq==b'g')|(seq==b'c'))[0].shape[0]
    n_number = np.where((seq==b'N')|(seq==b'n'))[0].shape[0]
    allnumber = seq.shape[0] - n_number
    return (gc_number,allnumber)

def calGC(fs):
    GC = 0
    all = 0
    hg38 = Fasta(fs)
    for record in hg38:
        seq = np.asarray(record)
        gc_number,all_number=lineGC(seq)
        GC = GC + gc_number
        all = all + all_number
    return (GC, all)

if __name__ == "__main__":
    gcNum, allNum = calGC(sys.argv[1])
    print("GC content: {:.3f}".format(gcNum/allNum))

这样的话，就只需要2分22秒了，已经是非常快的了，但是和perl还是有差距的。

最后，难道julia真的速度和perl就相差无几吗?

答案是否定的，因为julia设计是为了科学计算的，但是其字符串的性能并算不上优秀，我们可以调用BioSequence来处理生物序列

using BioSequences
using FASTX

function lineGC(seq)
    GCnumber=count(x->(x==DNA_G||x==DNA_C),seq)
    lineNum=length(seq)-count(isambiguous,seq)
    GCnumber,lineNum
end

function calGC(fs)
    GCnumber=zero(Int)
    lineNum=zero(Int)
    reader=open(FASTA.Reader,fs)
    for record in reader
        GC,all=lineGC(FASTA.sequence(record))
        GCnumber+=GC
        lineNum+=all
    end
    close(reader)
    round(GCnumber/lineNum;digits=3)
end

println("GC content: ",calGC(ARGS[1]))

这样就只需要11秒就可以计算出答案了