Skip to content

Commit 5d152c6

Browse files
Gyuuul2Gyuuul2
authored
Merge pull request #1120 from Gyuuul2/linclust2_dev
Fix linclust OOM, inconsistency & Add switch consensus rep
2 parents cb12a2d + 9418280 commit 5d152c6

18 files changed

Lines changed: 1224 additions & 542 deletions

‎data/workflow/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,7 @@ set(GENERATED_WORKFLOWS
2727
workflow/linsearch.sh
2828
workflow/databases.sh
2929
workflow/pickconsensusrep.sh
30+
workflow/pickconsensusrepfast.sh
3031
workflow/nucleotide_clustering.sh
3132
workflow/iterativepp.sh
3233
workflow/tsv2exprofiledb.sh

‎data/workflow/cascaded_clustering.sh‎

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -66,6 +66,43 @@ if [ "$CLUSTER_MODULE" = "cluster2" ]; then
6666
"$MMSEQS" mergeclusters "$SOURCE" "$2" "${TMP_PATH}/clu_redundancy" "${TMP_PATH}/clu2_rep" $MERGECLU_PAR \
6767
|| fail "Merging of clusters has died"
6868
fi
69+
70+
# Expose alignment results (only produced when --include-align-files is set).
71+
# Two alignment sources exist: the redundancy pre-clustering (internal linclust,
72+
# clu_redundancy_aln) and the cluster2 step (clu2_rep_aln). Union them keyed by the
73+
# final representatives ($2), keeping only each final cluster's member lines, so the
74+
# result has exactly one entry per cluster with that cluster's rep->member alignments.
75+
if [ -f "${TMP_PATH}/clu2_rep_aln.dbtype" ]; then
76+
if [ -f "${TMP_PATH}/clu_redundancy_aln.dbtype" ]; then
77+
# shellcheck disable=SC2086
78+
"$MMSEQS" mergedbs "${2}" "${2}_aln" \
79+
"${TMP_PATH}/clu_redundancy_aln" "${TMP_PATH}/clu2_rep_aln" \
80+
--merge-filter-target 1 ${VERBOSITY} \
81+
|| fail "mergedbs clu2_rep_aln died"
82+
else
83+
# shellcheck disable=SC2086
84+
"$MMSEQS" mvdb "${TMP_PATH}/clu2_rep_aln" "${2}_aln" ${VERBOSITY} \
85+
|| fail "mvdb clu2_rep_aln died"
86+
fi
87+
fi
88+
89+
# Optionally replace representatives by the most profile-consistent observed member,
90+
# reusing the alignments in ${2}_aln (no profile-vs-member realignment).
91+
if [ -n "$SWITCH_CONSENSUS_REP" ]; then
92+
# shellcheck disable=SC2086
93+
"$MMSEQS" pickconsensusrepfast "$1" "$2" "${TMP_PATH}/clu_switched" "${TMP_PATH}/switch_tmp" ${PICKREP_PAR} \
94+
|| fail "pickconsensusrepfast (switch representatives) died"
95+
# shellcheck disable=SC2086
96+
"$MMSEQS" rmdb "$2" ${VERBOSITY}
97+
# shellcheck disable=SC2086
98+
"$MMSEQS" mvdb "${TMP_PATH}/clu_switched" "$2" ${VERBOSITY} \
99+
|| fail "mvdb switched clustering died"
100+
if [ -z "$KEEP_SWITCH_ALN" ]; then
101+
# shellcheck disable=SC2086
102+
"$MMSEQS" rmdb "${2}_aln" ${VERBOSITY}
103+
fi
104+
rm -rf "${TMP_PATH}/switch_tmp"
105+
fi
69106
elif [ "$CLUSTER_MODULE" = "cluster1" ]; then
70107
INPUT="${TMP_PATH}/input_step_redundancy"
71108
STEP=0
@@ -281,12 +318,22 @@ if [ -n "$REMOVE_TMP" ]; then
281318
"$MMSEQS" rmdb "${TMP_PATH}/input_step_redundancy" ${VERBOSITY}
282319
# shellcheck disable=SC2086
283320
"$MMSEQS" rmdb "${TMP_PATH}/input_step_redundancy_h" ${VERBOSITY}
321+
# inner linclust scratch dir and its merged align output
322+
rm -rf "${TMP_PATH}/linclust"
323+
if [ -f "${TMP_PATH}/clu_redundancy_aln.dbtype" ]; then
324+
# shellcheck disable=SC2086
325+
"$MMSEQS" rmdb "${TMP_PATH}/clu_redundancy_aln" ${VERBOSITY}
326+
fi
284327

285328
if [ "$CLUSTER_MODULE" = "cluster2" ]; then
286329
# shellcheck disable=SC2086
287330
"$MMSEQS" rmdb "${TMP_PATH}/pref_rep" ${VERBOSITY}
288331
# shellcheck disable=SC2086
289332
"$MMSEQS" rmdb "${TMP_PATH}/clu2_rep" ${VERBOSITY}
333+
if [ -f "${TMP_PATH}/clu2_rep_aln.dbtype" ]; then
334+
# shellcheck disable=SC2086
335+
"$MMSEQS" rmdb "${TMP_PATH}/clu2_rep_aln" ${VERBOSITY}
336+
fi
290337
elif [ "$CLUSTER_MODULE" = "cluster1" ]; then
291338

292339
STEP=0

‎data/workflow/linclust.sh‎

Lines changed: 57 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,6 +102,43 @@ if [ "$LINCLUST_MODULE" = "linclust2" ]; then
102102
|| fail "mergeclusters died"
103103
fi
104104

105+
# Expose alignment results (only produced when --include-align-files is set).
106+
# The two align2clust passes each emit their own alignments; union them keyed by the
107+
# final representatives ($2) and keep only lines whose target is a member of that
108+
# final cluster (--merge-filter-target), so the result has exactly one entry per
109+
# cluster containing exactly that cluster's rep->member alignments.
110+
if [ -f "${TMP_PATH}/clu_aln.dbtype" ]; then
111+
if [ -f "${TMP_PATH}/clu_rep_aln.dbtype" ]; then
112+
# shellcheck disable=SC2086
113+
"$MMSEQS" mergedbs "${2}" "${2}_aln" \
114+
"${TMP_PATH}/clu_aln" "${TMP_PATH}/clu_rep_aln" \
115+
--merge-filter-target 1 ${VERBOSITY} \
116+
|| fail "mergedbs clu_aln died"
117+
else
118+
# shellcheck disable=SC2086
119+
"$MMSEQS" mvdb "${TMP_PATH}/clu_aln" "${2}_aln" ${VERBOSITY} \
120+
|| fail "mvdb clu_aln died"
121+
fi
122+
fi
123+
124+
# Optionally replace representatives by the most profile-consistent observed member,
125+
# reusing the alignments in ${2}_aln (no profile-vs-member realignment).
126+
if [ -n "$SWITCH_CONSENSUS_REP" ]; then
127+
# shellcheck disable=SC2086
128+
"$MMSEQS" pickconsensusrepfast "$1" "$2" "${TMP_PATH}/clu_switched" "${TMP_PATH}/switch_tmp" ${PICKREP_PAR} \
129+
|| fail "pickconsensusrepfast (switch representatives) died"
130+
# shellcheck disable=SC2086
131+
"$MMSEQS" rmdb "$2" ${VERBOSITY}
132+
# shellcheck disable=SC2086
133+
"$MMSEQS" mvdb "${TMP_PATH}/clu_switched" "$2" ${VERBOSITY} \
134+
|| fail "mvdb switched clustering died"
135+
if [ -z "$KEEP_SWITCH_ALN" ]; then
136+
# shellcheck disable=SC2086
137+
"$MMSEQS" rmdb "${2}_aln" ${VERBOSITY}
138+
fi
139+
rm -rf "${TMP_PATH}/switch_tmp"
140+
fi
141+
105142
elif [ "$LINCLUST_MODULE" = "linclust1" ]; then
106143
# 0. clusthash
107144
if [ -n "$CLUSTHASH" ]; then
@@ -217,6 +254,25 @@ if [ -n "$REMOVE_TMP" ]; then
217254
if [ "$LINCLUST_MODULE" = "linclust2" ]; then
218255
# shellcheck disable=SC2086
219256
"$MMSEQS" rmdb "${TMP_PATH}/pref" ${VERBOSITY}
257+
# shellcheck disable=SC2086
258+
"$MMSEQS" rmdb "${TMP_PATH}/clu" ${VERBOSITY}
259+
# shellcheck disable=SC2086
260+
"$MMSEQS" rmdb "${TMP_PATH}/input_rep" ${VERBOSITY}
261+
# shellcheck disable=SC2086
262+
"$MMSEQS" rmdb "${TMP_PATH}/input_rep_h" ${VERBOSITY}
263+
# shellcheck disable=SC2086
264+
"$MMSEQS" rmdb "${TMP_PATH}/pref_rep" ${VERBOSITY}
265+
# shellcheck disable=SC2086
266+
"$MMSEQS" rmdb "${TMP_PATH}/clu_rep" ${VERBOSITY}
267+
# align intermediates (only present with --include-align-files)
268+
if [ -f "${TMP_PATH}/clu_aln.dbtype" ]; then
269+
# shellcheck disable=SC2086
270+
"$MMSEQS" rmdb "${TMP_PATH}/clu_aln" ${VERBOSITY}
271+
fi
272+
if [ -f "${TMP_PATH}/clu_rep_aln.dbtype" ]; then
273+
# shellcheck disable=SC2086
274+
"$MMSEQS" rmdb "${TMP_PATH}/clu_rep_aln" ${VERBOSITY}
275+
fi
220276
if [ -n "$CLUSTHASH" ]; then
221277
# shellcheck disable=SC2086
222278
"$MMSEQS" rmdb "${TMP_PATH}/input_clusthash" ${VERBOSITY}
@@ -225,7 +281,7 @@ if [ -n "$REMOVE_TMP" ]; then
225281
# shellcheck disable=SC2086
226282
"$MMSEQS" rmdb "${TMP_PATH}/input_clusthash_redundancy" ${VERBOSITY}
227283
# shellcheck disable=SC2086
228-
"$MMSEQS" rmdb "${TMP_PATH}/clu" ${VERBOSITY}
284+
"$MMSEQS" rmdb "${TMP_PATH}/clu_merged" ${VERBOSITY}
229285
rm -f "${TMP_PATH}/order_clusthash_redundancy"
230286
fi
231287
rm -f "${TMP_PATH}/linclust.sh"

‎data/workflow/pickconsensusrepfast.sh‎

Lines changed: 67 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,67 @@
1+
#!/bin/sh -e
2+
# Fast profile-guided representative selection.
3+
# Reuses the representative-to-member alignments produced by clustering (run with
4+
# --include-align-files) instead of realigning a profile against every member, then
5+
# rewrites the cluster DB with the best-scoring observed member as the new rep.
6+
# Output is 1:1 with pickconsensusrep (a new cluster DB), only the core is faster.
7+
8+
fail() {
9+
echo "Error: $1"
10+
exit 1
11+
}
12+
13+
notExists() {
14+
[ ! -f "$1" ]
15+
}
16+
17+
[ -z "$MMSEQS" ] && echo "Please set the environment variable \$MMSEQS to your MMSEQS binary." && exit 1;
18+
[ "$#" -ne 4 ] && echo "Please provide <seqDB> <clusterDB> <outClusterDB> <tmpDir>" && exit 1;
19+
[ ! -f "$1.dbtype" ] && echo "$1.dbtype not found!" && exit 1;
20+
[ ! -f "$2.dbtype" ] && echo "$2.dbtype not found!" && exit 1;
21+
[ ! -d "$4" ] && echo "tmp directory $4 not found!" && mkdir -p "$4";
22+
TMP_PATH="$4"
23+
24+
# Prerequisite: the clustering alignment DB (${clusterDB}_aln), produced by linclust/cluster
25+
# when run with '--include-align-files 1 -a'. Without it there is nothing to reuse.
26+
if notExists "${2}_aln.dbtype"; then
27+
fail "${2}_aln not found. Re-run linclust/cluster with '--include-align-files 1 -a' so the representative-to-member alignments are available for reuse."
28+
fi
29+
30+
# 1. Score every observed member against its cluster profile and pick the best one.
31+
if notExists "${TMP_PATH}/rep_map.dbtype"; then
32+
# shellcheck disable=SC2086
33+
"$MMSEQS" pickrepprofile "$1" "${2}_aln" "${TMP_PATH}/rep_map" ${PICKREP_PAR} \
34+
|| fail "pickrepprofile failed"
35+
fi
36+
37+
# 2. Extract the oldRep -> newRep mapping (key space).
38+
if notExists "${TMP_PATH}/rep_map.tsv"; then
39+
# shellcheck disable=SC2086
40+
"$MMSEQS" prefixid "${TMP_PATH}/rep_map" "${TMP_PATH}/rep_map.tsv" --tsv ${VERBOSITY} \
41+
|| fail "prefixid rep_map failed"
42+
fi
43+
awk '{ print $1 "\t" $2 }' "${TMP_PATH}/rep_map.tsv" > "${TMP_PATH}/rep_mapping.txt"
44+
45+
# 3. Old clustering as a flat oldRep -> member table.
46+
if notExists "${TMP_PATH}/clu.tsv"; then
47+
# shellcheck disable=SC2086
48+
"$MMSEQS" prefixid "$2" "${TMP_PATH}/clu.tsv" --tsv ${VERBOSITY} \
49+
|| fail "prefixid clu failed"
50+
fi
51+
52+
# 4. Rewrite the clustering with the new representatives.
53+
awk 'FNR == NR{ f[$1] = $2; next }
54+
$1 != prev { print f[$1] "\t" f[$1]; prev = $1; }
55+
$1 in f && $2 != f[$1]{ print f[$1] "\t" $2 }' \
56+
"${TMP_PATH}/rep_mapping.txt" "${TMP_PATH}/clu.tsv" > "${TMP_PATH}/updated_clu.tsv"
57+
58+
# shellcheck disable=SC2086
59+
"$MMSEQS" tsv2db "${TMP_PATH}/updated_clu.tsv" "$3" --output-dbtype 6 ${VERBOSITY} \
60+
|| fail "tsv2db failed"
61+
62+
if [ -n "$REMOVE_TMP" ]; then
63+
# shellcheck disable=SC2086
64+
"$MMSEQS" rmdb "${TMP_PATH}/rep_map" ${VERBOSITY}
65+
rm -f "${TMP_PATH}/rep_map.tsv" "${TMP_PATH}/rep_mapping.txt" "${TMP_PATH}/clu.tsv" "${TMP_PATH}/updated_clu.tsv"
66+
rm -f "${TMP_PATH}/pickconsensusrepfast.sh"
67+
fi

‎src/CommandDeclarations.h‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -93,6 +93,8 @@ extern int offsetalignment(int argc, const char **argv, const Command& command);
9393
extern int orftocontig(int argc, const char **argv, const Command& command);
9494
extern int touchdb(int argc, const char **argv, const Command& command);
9595
extern int pickconsensusrep(int argc, const char **argv, const Command& command);
96+
extern int pickconsensusrepfast(int argc, const char **argv, const Command& command);
97+
extern int pickrepprofile(int argc, const char **argv, const Command& command);
9698
extern int prefilter(int argc, const char **argv, const Command& command);
9799
extern int prefixid(int argc, const char **argv, const Command& command);
98100
extern int profile2cs(int argc, const char **argv, const Command& command);

‎src/MMseqsBase.cpp‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -618,6 +618,23 @@ std::vector<Command> baseCommands = {
618618
{"clusterDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::clusterDb },
619619
{"clusterDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::clusterDb },
620620
{"tmpDir", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory }}},
621+
{"pickconsensusrepfast", pickconsensusrepfast, &par.pickconsensusrepfast, COMMAND_CLUSTER,
622+
"Fast profile-guided representative selection reusing clustering alignments (workflow)",
623+
"Reuses the representative-to-member alignments produced by clustering (run with --include-align-files) to score observed members against each cluster profile and rewrite the cluster DB with the best-scoring member as the new representative, without profile-vs-member realignment.",
624+
"Gyuri Kim <gyuribio@snu.ac.kr>",
625+
"<i:seqDB> <i:clusterDB> <o:clusterDB> <tmpDir>",
626+
CITATION_MMSEQS2, {{"seqDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::sequenceDb },
627+
{"clusterDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::clusterDb },
628+
{"clusterDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::clusterDb },
629+
{"tmpDir", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::directory }}},
630+
{"pickrepprofile", pickrepprofile, &par.pickrepprofile, COMMAND_CLUSTER,
631+
"Score observed cluster members against the cluster profile and pick the best (core of pickconsensusrepfast)",
632+
"Reuses representative-to-member alignments to build the center-star MSA/profile and selects the best-scoring observed member, without profile-vs-member realignment. Emits a per-cluster mapping oldRepKey -> newRepKey score coverage clusterSize.",
633+
"Gyuri Kim <gyuribio@snu.ac.kr>",
634+
"<i:seqDB> <i:alnResultDB> <o:repMappingDB>",
635+
CITATION_MMSEQS2, {{"seqDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::sequenceDb },
636+
{"alnResultDB", DbType::ACCESS_MODE_INPUT, DbType::NEED_DATA, &DbValidator::resultDb },
637+
{"repMappingDB", DbType::ACCESS_MODE_OUTPUT, DbType::NEED_DATA, &DbValidator::genericDb }}},
621638
{"prefilter", prefilter, &par.prefilter, COMMAND_PREFILTER,
622639
"Double consecutive diagonal k-mer search",
623640
NULL,

0 commit comments

Comments
 (0)