blastreduce code should not remove self-alignment lines

[rbdavid]

The self-alignment lines in the all-by-all blast result files are important to maintain, specifically for when the sequence set is compressed/reduced before running the alignments. The below line needs to be removed:

EST/pipelines/est/templates/reduce-template.sql

Line 10 in 2c470f5

WHERE NOT qseqid = sseqid

In the same PR, the logic of the blastreduce's reduce-template.sql code should be verified to be replicating the original perl code's results. Comments in that sql template file says that it does not exactly match the original code.

EST/pipelines/est/templates/reduce-template.sql

Lines 14 to 17 in 2c470f5

	-- The original blastreduce uses a sort + picking the first occurrence of a
	-- (qseqid, sseqid) pair to deduplicate. That is replicated here by paritioning
	-- on (qseqid, sseqid) and assigning a rank to a specific sorted order. It does
	-- not exactly match the previous method but it is extremly close

I haven't done any sort of speed testing but I highly doubt processing blast result file(s) is the bottle neck for the EST workflow. Is this duckdb strategy worth the slight inaccuracies?

[rbdavid]

Okay, here's the new (as of PR #255) breakdown of the est.nf workflow, focusing down on the all_by_all.nf subworkflow:

1. If params.multiplex, reduce the input fasta file down to only unique, representative sequences using CD-HIT. Otherwise, just work with the input fasta file.

2. Create the blast-db files for the sequences contained in the blast_input_file. Files like database.phr, etc are written, containing binary encodings for all the sequences used in the formatdb call. These database files are optimized for the all-to-all alignments about to be done.

3. blastreduce_transcode_fasta (a terribly named process) is run to gather sequence lengths for all sequences in the original, input fasta file.

4. The blast_input_file is split up into shards to run embarrassingly parallel calculations.

5. The all_by_all_blast process block does a lot of stuff.

   • blastall -p blastp is run where the fasta file shard contains the query sequences that are aligned against all sequences in the database.* files (these are the subject sequences). This will inherently include self-alignments because sequence A was present in the fasta file used to create the database files. I will share paths in slack to files where you can see the self alignment results, written to the .tab file output from the blastall call.
   • transcode_blast.py script is run to convert the tab separated blast output file to a parquet file. During this process, columns in the blast out file are discarded leaving just these as columns: ["qseqid", "sseqid", "pident", "alignment_length", "bitscore"]. This script writes a .tab.parquet file for speedy data analysis/manipulation processes.
   • This .tab.parquet file is input to render_prereduce_sql_template.py along with the --sql-template EST/pipelines/est/templates/prereduce-template.sql to render an sql script. That template file includes the first instance of the guilty line WHERE NOT qseqid = sseqid (

     EST/pipelines/est/templates/prereduce-template.sql

     Line 8 in 42faa7f

     WHERE NOT qseqid = sseqid

     ). This template sql code also does the alphabetical sorting of the qseqid and sseqid columns.
   • Duckdb is called using the templated sql code, producing a .tab.sorted.parquet file of the all-to-all blast alignment results (BUT already missing the self-alignments).

6. The blastreduce process block does another round of sql-templating and duckdb parquet processing. This time, it uses the EST/pipelines/est/templates/reduce-template.sql template file, wherein the second instance of the guilty line resides (

   EST/pipelines/est/templates/reduce-template.sql

   Line 10 in 42faa7f

   WHERE NOT qseqid = sseqid

   ) that removes self-alignments. Also, if I understand this sql block (

   EST/pipelines/est/templates/reduce-template.sql

   Lines 18 to 26 in 42faa7f

   	CREATE TABLE reduced AS
   	SELECT qseqid, sseqid, pident, alignment_length, bitscore
   	FROM (
   	SELECT qseqid, sseqid, pident, alignment_length, bitscore, ROW_NUMBER()
   	OVER (PARTITION BY qseqid, sseqid
   	ORDER BY bitscore DESC, pident ASC, alignment_length ASC) ranked_order
   	FROM blast_results
   	) t
   	WHERE t.ranked_order = 1;

   ) correctly, this is where the top ranked row for qseqid and sseqid pairs (no matter the ordering) is kept while all remaining instances are dropped.

7. Finally, if params.multiplex, the clustered/reduced sequence set is expanded to include non-representative sequences. But! since the blastreduce sql templates include the guilty line, self-alignments won't be included and so cluster children will never be given intra-cluster alignment scores.

[nilsoberg]

I haven't done any sort of speed testing but I highly doubt processing blast result file(s) is the bottle neck for the EST workflow. Is this duckdb strategy worth the slight inaccuracies?

Historically, the blastreduce step has taken longer than the all-by-all and dominated the computation pipeline. With the refactoring that dynamic has changed, and the all-by-all is now the dominant factor.