hecatomb's same output files from the same dataset are with different sizes about hecatomb HOT 10 CLOSED

Hi Leran,
The crash I think relates to #34 and the logging should be addressed in the next build. Let's revisit it once the pipeline is finishing correctly.

from hecatomb.

Some additional data:
134 samples run on Pathogen:
assembly.fasta:24280
seqtable.fasta: 35,306,586

Same 134 samples running on HTCF (MMSEQS_AA_PRIMARY step):
assembly.fasta:24568
seqtable.fasta: 9,344,202

I don't think the second seqtable will be added to after this step

from hecatomb.

It seems like the difference is mainly arising during clustering:

Pathogen:
more p09_cluster_similar_sequences.M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10.log
Size of the sequence database: 341461
Size of the alignment database: 341461
Number of clusters: 325347

more p08_remove_exact_dups.M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10.log *****
Input: 971106 reads 223037112 bases.
Duplicates: 65529 reads (6.75%) 15441315 bases (6.92%) 0 collisions.
Result: 905577 reads (93.25%) 207595797 bases (93.08%)

HTCF (slurm):
more p09_cluster_similar_sequences.M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10.log *****
Size of the sequence database: 147516
Size of the alignment database: 147516
Number of clusters: 77134

more p08_remove_exact_dups.M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10.log
Input: 971504 reads 223125309 bases.
Duplicates: 65541 reads (6.75%) 15444031 bases (6.92%) 0 collisions.
Result: 905963 reads (93.25%) 207681278 bases (93.08%)

from hecatomb.

This very well may be due to some changes we made to the clustering parameters for linclust. There was some toying with these settings over the past month, so if you ran the Pathogen run a while back and the HTCF run more recently you would likely get different results.

Can you check the cluster settings for each run (should be in your config file, but may be in the rule file depending on which version and when Mike made updates). That large of a difference is more likely explained by a major setting change than a compute error.

from hecatomb.

Here are the differences:
Pathogen:
mmseqs easy-linclust hecatomb_out/PROCESSING/TMP/p08/M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10_R1.deduped.out.fastq hecatomb_out/PROCESSING/TMP/p09/M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10_R1 hecatomb_out/PROCESSING/TMP/p09/M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10_TMP --kmer-per-seq-scale 0.3 -c 0.97 --cov-mode 1 --threads 16 &> hecatomb_out/STDERR/p09_cluster_similar_sequences.M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10.log;

HTCF:
mmseqs easy-linclust hecatomb_out/PROCESSING/TMP/p08/M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10_R1.deduped.out.fastq
hecatomb_out/PROCESSING/TMP/p09/M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10_R1 hecatomb_out/PROCESSING/TMP/p09/M721_I9060_34
044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10_TMP --kmer-per-seq-scale 0.3 -c 0.7 --cov-mode 1 --min-seq-id 0.95 --alignment-mode 3 --threads 24 &> hecatomb_out/STDERR/p09_cluster_similar_sequences.M721_I9060_34044_Parkes_IBD_2034A_25_11_20_NEB_46_TCCCGAAT_S10.log;

from hecatomb.

NOTE:
-c = minimum coverage originally 0.97 to 0.7
--min-seq-id + --alignment-mode 3 (number of identical aligned residues divided by the number of aligned columns including internal gap columns) -> now 0.95

• not clear to me what happens without these params

from hecatomb.

Well dropping from 0.97 to 0.7 is going to have the greatest effect and is likely way lower than what we intend. This should be set back to 0.97 or 0.95.

You can read about alignment mode here: https://github.com/soedinglab/MMseqs2/wiki#how-to-set-the-right-alignment-coverage-to-cluster

@mroach-awri can we go back to my original settings as the default for now?

from hecatomb.

Yes I think I pushed my settings to the last build. The original default (-c .97 --cov-mode 0) I think specifies 97% residue matches in the longer sequence, which will only cluster sequences that are end-to-end almost identical and most clusters are n=1. The current settings (--cov-mode 1 -c 0.7 --min-seq-id 0.95) specify 70% alignment coverage of the member sequence by the rep sequence at 95% identity. These are the setting I'll probably be using as I want to maximize runtime performance.

I'll be pushing a new build soon; did you want the original defaults or did you have a specific coverage and seq identity in mind?

from hecatomb.

from hecatomb.

Running:
grep \> seqtable.fasta | sed 's/.*:\(.*\):.*/\1/' | awk '{n+=$1}END{print n}'
gives consistent number of reads on seqtable.fasta file.

Leran

from hecatomb.