Reading very large text files

[rsidd]

Apologies if this is covered previously, but I am unable to find an answer after extensive reading. I did find one reference, discussed below.

Let's say I am reading the human genome chromosome 1, approx 0.25GB unzipped. The zipped file is available here
http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr1.fa.gz

In this file the first line is a "header" which I don't care about, and the rest is the sequence. I want the output to be one long string representing the sequence, concatenating all lines after throwing away the first line and removing trailing newline characters. A simple-minded Python function that does this is as follows. On my machine it takes less than 700ms to run.

def readfasta(filename):
    f = open(filename)
    l0 = f.readline()
    l0 = f.readline()[0:-1]
    for l in f:
        l0 += l[0:-1]
    return l0

A direct translation to Julia takes -- well, it doesn't seem to finish in any reasonable time, but it takes 11 seconds on a file 1/200 the size and the time seems to increase more than linearly with file size.

I found this post apparently addressing the matter: https://groups.google.com/forum/#!topic/julia-dev/UDllYRfm64w
The OP said reading the file into an array of strings helped. I tried that, and the following runs in 20 seconds -- still about 30 times slower than python.

function readfasta(filename::String)
    f = open(filename)
    l0 = readline(f)
    strarr = map(x::String->x[1:end-1], readlines(f))
    return join(strarr,"")
end

A comment on that post suggested using IOBuffer, but I am not clear on how to do this (joining strings after dropping the last newline characters) more efficiently than what is being done above. Any help would be very welcome.

[bicycle1885]

I think issue is not the best place to ask questions. You can ask at the julia-users group or stackoverflow. BTW, if you'd like to load sequences in a FASTA file, Bio.jl has a FASTA parser.

[JeffBezanson]

I get ~0.8 seconds in python and ~5 seconds in julia, so it seems like there is a performance issue to fix here. I'm not sure why you're seeing 30x slower; what version of julia are you using?

[rsidd]

Thanks for the reopen. This is julia 0.5.0 via Ubuntu PPA (staticfloat).

Did you get 5 secs with exactly the julia code I posted? I get 16-24 secs but generally about 20 secs each time I try. Consistently about 600-700ms with the python code I posted.

[rsidd]

Answering @bicycle1885 above. Using Bio.seq,

@time for record in open(FASTAReader,"chr1.fa")
           print(length(record.seq))
end

takes 19 seconds to run, similar to my code.

[bicycle1885]

Your Julia seems to be much slower than mine (note that the first call includes JIT compiling time):

julia> @time open(collect, FASTAReader, "chr1.fa");
  2.541710 seconds (1.34 M allocations: 548.084 MB, 9.27% gc time)

julia> @time open(collect, FASTAReader, "chr1.fa");
  1.478964 seconds (81 allocations: 494.553 MB, 15.05% gc time)

Also note that what your Python code does is just reading and concatenating strings while Bio.jl does checking validity and encoding data into BioSequence objects.

As a reference, I ran your Python code on my computer (Python 3.4.1):

In [5]: %time s = readfasta("chr1.fa")
CPU times: user 2.18 s, sys: 199 ms, total: 2.38 s
Wall time: 2.41 s

and your Julia code (Julia 0.5.0):

julia> @time readfasta("chr1.fa");
  4.506265 seconds (19.94 M allocations: 1.942 GB, 23.07% gc time)

[rsidd]

I get

julia> @time open(collect, FASTAReader, "chr1.fa");
 20.262316 seconds (25.87 M allocations: 690.482 MB, 0.84% gc time)

julia> @time open(collect, FASTAReader, "chr1.fa");
 21.800200 seconds (24.87 M allocations: 874.410 MB, 1.03% gc time)

julia> @time readfasta("chr1.fa");
 22.814977 seconds (19.94 M allocations: 1.978 GB, 4.46% gc time)

Note the vastly bigger number of allocations for the Bio.jl code compared to what you see. But for my readfasta function there is no major difference in allocations. This is with a fresh update of julia to the latest in the PPA (still 0.5.0), and after a Pkg.update(). Should I be trying a Julia install from github directly?

[rsidd]

I installed julia 0.5.1-pre+4 (2016-10-28 19:39 UTC) from github. Pretty much the same timings as above. Baffled why @JeffBezanson sees faster-but-still-slow performance while @bicycle1885 sees performance comparable to python. Can it be something OS-related? This is Ubuntu Xenial, Intel i7-4770 CPU @ 3.40GHz, 32GB RAM. ZFS filesystem, but I tried reading from an ext4 filesystem and the results are the same.

On my laptop (i5-5200U CPU @2.20 GHz, 4GB RAM, Ubuntu Xenial, ZFS) both of these are much faster: the readfasta function takes 6-16 seconds (very variable, and later calls aren't always faster) and the FASTAReader call takes about 5-6 seconds. I don't understand the difference, since it's a slower machine with less RAM running the same OS. Even so, python 2.7.12 takes 700ms on the function I posted.

If I translate my julia function to python (read a list of lines, truncate them, join them) it takes about 10 seconds, not too different from julia. But if I translate my python function to julia it takes a very very long time apparently because concatenating strings is slow in julia.

[rsidd]

This is very mysterious: after a reboot (previous uptime had been about 2 months) the exact same readfasta function as above takes 1.8 seconds, 10x faster than before. The Bio.jl function takes 1.8 seconds too, again 10x faster. The python function takes 0.6 seconds, almost unchanged from earlier.

So I am wondering if it is an OS issue, but if so, why was it affecting julia and not python (or, as far as I can tell, any other software that read/wrote large files, such as bedtools)?

[jiahao]

@rsidd what happens when you time the result after running gc_enable(false)? I wonder if the performance problem you are seeing has been visited before in #10428

[rsidd]

With gc_enable(false) the readfasta time goes down from avg 1.6-1.8 seconds to 1.1-1.2 seconds. If performance slows in the future to the extent I saw previously, I will give it a try and report back.