Clarify when chunking files

[olgabot]

Hello! I ran sourmash compute on ~50k single-cell (SmartSeq2 library prep) RNA-seq samples (they are here if you would like to see them: aws s3 ls s3://olgabot-maca/facs/sourmash/) and wanted to index/compare them all vs our cell-cell distances/clusters/annotations using gene count tables

At first, I thought sourmash compare was broken because it said it was only loading 3444 signatures out of the 50k:

 Wed 13 Jun - 18:34  /mnt/data/maca/facs/sourmash 
 ubuntu@olgabot-reflow  ls | xargs sourmash compare --ksize 21 --csv ../ksize21.csv
loaded 3444 signatures total.

But there's 51,446 files here!!

 Wed 13 Jun - 18:54  /mnt/data/maca/facs/sourmash 
 ubuntu@olgabot-reflow  ls | wc -l
51446

But then sourmash index was more explicit in showing that it was chunking the data:

 Wed 13 Jun - 18:05  /mnt/data/maca/facs/sourmash 
 ubuntu@olgabot-reflow  ls | xargs sourmash index --ksize 21 ../ksize21db
loading 3444 files into SBT
loaded 3444 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3441 files into SBT
loaded 3441 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3441 files into SBT
loaded 3441 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3439 files into SBT
loaded 3439 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3699 files into SBT
loaded 3699 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3437 files into SBT
loaded 3437 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3434 files into SBT
loaded 3434 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3449 files into SBT
loaded 3449 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3447 files into SBT
loaded 3447 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3436 files into SBT
loaded 3436 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3435 files into SBT
loaded 3435 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3435 files into SBT
loaded 3435 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3439 files into SBT
loaded 3439 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3442 files into SBT
loaded 3442 sigs; saving SBT under "../ksize21db"

Finished saving nodes, now saving SBT json file.
loading 3028 files into SBT

So it seems that sourmash compute was not broken after all, but just taking its time through all the samples.

Here are my questions:

1. Can there be more explicit descriptions of how the data gets chunked by sourmash, either as a flag or a description?
2. Is the chunking on a per-cpu basis? If so, that would be helpful to know when starting up an EC instance to know how much to allocate.
3. Could this output, e.g. on chunk 1/20 be output to the stdout?

Thank you!

EDIT: This was run on an AWS EC2 m4.large

[luizirber]

I was surprised to see chunking happening in sourmash... Turns out it's coming from xargs!
(you can control how big the chunk is with -n).

I'm guessing you're using xargs to be able to pass more args in the command line (50k would break the shell limits for args). In this case, may I suggest using --traverse-directory instead? If you're in /mnt/data/maca/facs/sourmash

$ sourmash index --ksize 21 --traverse-directory ../ksize21db .

will create the index in /mnt/data/maca/facs/ksize21db.sbt.json

(I need to document this better, but the genbank and refseq databases use --traverse-directory to create the index, see https://github.com/dib-lab/sourmash_databases/blob/75343ca39543c2f17b29ce44b665480c60890abe/Snakefile#L65 for the command that generates them)

[luizirber]

And you can also use --traverse-directory with sourmash compare:

$ sourmash compare --ksize 21 --csv ../ksize21.csv --traverse-directory .

[luizirber]

(pinging @taylorreiter for more large scale compare tips and tricks =] )

[taylorreiter]

Other than --traverse-directory, the only other thing I have found helpful is --singleton and --name-from-first, so that if you ever have a multifasta you don't need to split it, you can just calculate one signature. compare knows you're comparing all of those signatures still, so the syntax would be:

sourmash compare -o output sig

I have a few other tricks for dealing with the very large compare matrix that is output by this, but those are pretty use-case specific.

[olgabot]

Oh interesting! I'll use --traverse-directory for index.

For compare, I didn't see a --traverse-directory option so I'll have to try something else. In this case, each signature is from one or more pairs of fastq files so I'm not able to use --singleton, but I could definitely cat all the signatures into one huge json. What do you think?

[olgabot]

For the large compare matrix, do you end up saving it as a sparse matrix?

[brooksph]

Hi @olgabot! I added '--traverse-directory' as an option for 'compare' last week if you'd like to give it a try.

[luizirber]

I think --traverse-directory as an option for compare didn't get in 2.0.0a7, but I'm about to release 2.0.0a8 (which will include it)

[taylorreiter]

re sparse matrix, no, I have not experimented with that. I was comparing tetramernucleotide frequency of each contig across a euk genome that was highly fragmented. I ended up doing compare on a subset instead, so 1/4 at a time. Saving it as a csv is also a bad idea when dealing many comparisons. Lastly, @luizirber has recommended dask before.

[olgabot]

Thanks, everyone!