Closed kfuku52 closed 1 year ago
The problem stems from a wrong integrate
output which looks like this, with no seq num & length info. amalgkit integrate
finished without error, but with a warning zcat: unknown compression format
.
scientific_name | curate_group | run | read1_path | read2_path | is_sampled | is_qualified | exclusion | lib_layout | spot_length | total_spots | total_bases | size | private_file |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Please add in format: Genus species | Please add | DMexp001FlL1 | /Volumes/kfT7/Dropbox/repos/amalgkit/data/getfastq/DMexp0001Fl_R1/DMexp001FlL1_R1.fq | unavailable | yes | yes | no | single | 0 | 0 | 0 | 1.6856E+10 | yes |
@Hego-CCTB Instead of cat or zcat, probably we should always use seqkit seq
, which can automatically handle different compression types of fastq files. I will replace zcat causing the above error with seqkit seq
, but the same improvement may be applied to other parts.
https://bioinf.shenwei.me/seqkit/usage/#seq
Also, seqkit can compress fastq faster than pigz. We should replace it as well at some point.
Seqkit writes gzip files very fast, much faster than the multi-threaded pigz, so there's no need to pipe the result to gzip/pigz.
I now figured out why the fastq processing was so slow that we needed the --accurate_size
in integrate
. It's most likely due to zcat
. Without it, seqkit stat
should finish in seconds. Now, --accurate_size
should be defaulted to 'yes'. I will apply this change as well if this is the case.
Wait...zcat
is not used with --accurate_size yes
, and integrate
finished in 25 sec to process 1 single-end fastq. This is fast enough, but @Hego-CCTB do we have any reasons to keep --accurate_size no
?
--accurate_size no
uses just the first thousand reads to calculate the size, while --accurate_size yes
runs seqkit on everything.
While seqkit
can process a whole fastq.gz
file on its own (when--accurate_size yes
), we need zcat
to extract the first 1000 reads out of a .gz
file --accurate_size no
. Depending on the fastq
files this can make a significant difference in computing time.
Thank you, but please check my above comments. zcat
was not necessary and I will get rid of it. Without zcat
, --accurate_size yes
took 27 secs, and --accurate_size no
took 27 secs as well for a 17-Gbp RNA-seq data. There seems to be no advantage to subsample.
I have a test set of 3 SRA samples, 2 single-end and one paired-end. It takes me 42 seconds to complete integrate --accurate_size no
and 69 seconds to complete integrate --accurate_size yes
. Running integrate --accurate_size yes
on the paired-end sample takes by far the most time. (This is running it on my mac).
...with zcat
, right? The speed is not much different already, but did you see that difference after getting rid of zcat
, which is slower than seqkit seq
?
I'm not sure I understand you correctly, integrate --accurate_size no
(faster for me) always uses zcat
and integrate --accurate_size yes
(slower for me) never uses zcat
. What do you mean by "getting rid of zcat
" ?
Ah, I misunderstood the time. OK, it's faster with no
. The question now is whether we want to keep this functionality for the speedup at <10 sec/sample. I remember we needed --accurate_size
because the first implementation took >1 h per sample, but now this is not the case anymore with only 69 secs for 3 samples.
Yes, now we are on the same page!
Let me run a quick test on some of the RNAseq data I produced. I remember the difference being quite significant for those too.
If the difference is in the realm of seconds, I don't mind scrapping --accurate_size
altogether.
I used all of the Plumbago files for this, i.e. 12 paired-end samples:
--accurate_size no
: 372 sec (ca. 10 min.)
--accurate_size yes
: 3,016 sec (ca. 50 minutes)
I believe the difference between our experience is that I am running integrate
on fq.gz
files. This means they have to be decompressed first (hence why I make use of zcat to just feed the first 1000 sequences to seqkit).
This may very well be the reason for the warning you got:
zcat: unknown compression format.
(zcat doesn't recognise .fq
as a compression format)
since you were running integrate
on already decompressed .fq
files. seqkit
is presumably almost instant with decompressed files, hence you don't experience any difference in speed between --accurate_size no
and --accurate_size yes
, but I do quite significantly.
So I should include a check for file extension. If it's a .fq
file we can just always run --accurate_size yes
and if it's a .fq.gz
file, we can run the faster --accurate_size no
instead (if the user allows it).
I propose the following changes:
--accurate_size to yes
as default--accurate_size no
should only be run for .gz
compressed files and only if the user wants this (i.e. specifying no
)I see. I like your plan and I would like to hand over the above error to you. amalgkit should recognize and distinguish at least these extensions: .fq, .fastq, .fq.gz, .fastq.gz. Thank you in advance.
This should fix the behaviour. Tested with compressed and decompressed files. Also changed the help message to make it clear this option only applies to decompressed files. https://github.com/kfuku52/amalgkit/commit/9d3ebcf942f8efe67a21e0bb2ea2b90663cce100
Thank you. Could you adjust the warning message so there is a species between the file name and "is"? This kind of mistakes may be avoided by using .format().
like this?
warnings.warn("{} is not a fastq file. Skipping.".format(fastq_files[0]))
Yes, thank you.
command
stdout
stderr