在生物信息学中经常用到的脚本语言主要是
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秒就可以计算出答案了
网友评论