Source code for dnachisel.biotools.bowtie

import subprocess
import tempfile
import os

def run_process(name, parameters):
    process = subprocess.run(
        parameters,
        stdout=subprocess.PIPE,
        stderr=subprocess.PIPE
    )
    if process.returncode:
        error = process.stderr.decode()
        raise OSError("%s failed:\n\n%s" % (name, error))
    return process.stdout

def create_bowtie_index_from_sequences(sequences, path):
    fasta_path = os.path.join(path, 'sequences.fa')
    bowtie_path = os.path.join(path, 'bowtie')
    with open(fasta_path, 'w') as f:
        f.write("\n".join([
            ">%d\n%s" % (i, sequence)
            for i, sequence in enumerate(sequences)
        ]))
    run_process("build-bowtie", [
        "bowtie-build", "-f", fasta_path, bowtie_path, "--quiet"
    ])
    return bowtie_path

def find_all_bowtie_matches(
    sequence, bowtie_index_path, match_length, max_mismatches=0
):
    """Return (short) matches between a sequence and a Bowtie index.
    
    The result is of the form [(start, end), n_mismatches)] where (start, end)
    indicates the position of the match in the sequence, and n_mismatches is
    the number of mismatches with the closest homology in the index.
    
    
    """

    # CREATE THE PARAMETERS

    parameters = ["bowtie"]
    parameters += ["--best", "-k", "1"]  # only return the best alignments
    parameters += ["-v", str(max_mismatches)]  # only allow that N mismatches
    parameters += [bowtie_index_path]
    parameters += ["--quiet", "--suppress", "2,3,4,5,6,7"]  # small output
    k = match_length
    kmers = [sequence[i : i + k] for i in range(len(sequence) - k + 1)]
    if k * len(kmers) < 10000:
        # Input the sequences directly
        tmp_fasta_path = None
        parameters += ["-c", ",".join(kmers)]
    else:
        # Write sequences to a file if too many.
        tmp_fasta_path = tempfile.mktemp(".fa")
        with open(tmp_fasta_path, "w") as f:
            entries = [">%d\n%s" % (i, s) for i, s in enumerate(kmers)]
            f.write("\n\n".join(entries))

        parameters += ["-f", tmp_fasta_path]

    # RUN THE PROCESS
    try:
        output = run_process("BOWTIE", parameters)
    except Exception as err:
        raise err
    finally:
        if tmp_fasta_path is not None:
            os.remove(tmp_fasta_path)
    output_records = [
        line.split("\t") for line in output.decode().split("\n") if len(line)
    ]
    return [
        ((int(index), int(index) + k), edits.count(":"))
        for index, edits in output_records
    ]