mmseqs as a replacement for cd-hit-est

[hjruscheweyh]

Context

Hi!

We're currently looking for a replacement for CD-HIT-EST for clustering genes in nucleotide space. We have a few conditions when clustering. The CD-HIT-COMMAND is straightforward:

cd-hit-est -i genes.fa -o cluster.fa -c 0.95 -T 64 -M 0 -G 0 -aS 0.9 -g 1 -r 0 -d 0

Basically, we want to perform greedy clustering at 95% identity and at least 90% of the shorter gene should be aligned. The longest sequence of a cluster should become the representative and a member of a cluster should not align against any other representative with a better score/percid.

We have now come up with the following command(s):

mmseqs createdb genes.fa genes.mmseqs.db --dbtype 2 --shuffle 0

mmseqs cluster genes.mmseqs.db genes.mmseqs.cluster tmp --kmer-per-seq-scale 0  --kmer-per-seq 1000  \
-s 4  --max-seq-len 80000 --remove-tmp-files 0 --cluster-mode 2 --min-seq-id 0.95 --threads 64  --cov-mode 1 \
-c 0.9  --spaced-kmer-mode 0  --alignment-mode 3 --cluster-reassign 1

The key parameter seems to be --spaced-kmer-mode 0 as this prevents very similar sequences being placed in different clusters. as reported by others (#831 #489).

Using these parameters the clustering compares quite well against the clusters produced by cd-hit:

mmseqs clusters = 13'471'787
cdhit clusters  = 13'640'333
runtime(mmseqs) / runtime(cdhit) = 0.09
Rand Score: 0.99
Adjusted Rand Score: 0.91
V-measure: 0.99

So it seems that mmseqs is slightly overclustering (might be due to local vs global alignment) and the cluster consistency metrics indicate that even the cluster membership (and not only the representatives) seem to be mostly correct.

I would like to get your feedback on the command(s) before I make this our standard (nucleotide) gene clustering pipeline. Do you think there are other parameters that could be tweaked?

Your Environment

version: 6f45232ac8daca14e354ae320a4359056ec524c2
Downloaded static binary from GitHub

Best and thanks,
Hans

[martin-steinegger]

Your parameters are correct and should result in our desired clustering; however, the options -s 4 and —cluster-reassign 1 will not have any effect. The variations you're seeing may be a result from our transitive clustering steps, which our nucleotide clustering cannot correct (no support for —cluster-reassign). The reason why reassign is not working is because we lose the alignment anchor information needed to realign members.

Also —alignment-mode 3 and —cluster-mode 2 are optional.

The nucleotide cluster is still work in progress but for your use case it should work given the issues mentioned above.

[hjruscheweyh]

Thank you very much. I will adapt the command and will explore the differences between cd-hit and mmseqs further.

Hans

[bayegy]

According to the MMseqs2 help, --cov-mode 1 uses target coverage and --cov-mode 5 uses the coverage of the shorter sequence. Wouldn't --cov-mode 5 then be a better fit for the requirement that "at least 90% of the shorter gene should be aligned"?