mmseq2tax.R

mmseqs2 = function (seqs, db, threads=24, preamble="mamba activate mmseqs2;") {
    id = uuid::UUIDgenerate()
    tmpdir = glue::glue("tmp/{id}")
    dir.create(tmpdir, recursive = T)
    fafile = glue::glue("{tmpdir}/input.fa")
    writeLines(paste(sprintf(">%d", seq_along(seqs)), seqs, sep = "\n"), fafile)
    
    cmd = glue::glue(
        "{preamble} ",
        "mmseqs easy-taxonomy",
        "  {fafile} {db}",
        "  {tmpdir}/output {tmpdir}/tmp",
        "  --tax-lineage 1",
        "  --search-type 3",
        "  --vote-mode 2",
        "  --lca-mode 4",
        "  --report-mode 1",
        "  --threads {threads}",
        "  -s 4")
    print(cmd)
    system(sprintf("bash -c '%s'", cmd))
    readr::read_tsv(glue::glue("{tmpdir}/output_lca.tsv"), col_names=c("asv", "taxid", "rank", "taxon", "lineage"), show_col_types = F) |>
        mutate(asv_seq=seqs[asv])
}

mmres = mmseqs2(getSequences(stabnochim), "tmp/SILVA.seqtaxdb")