<commit_before>#!/usr/bin/env python


import os
import sys
import subprocess

from bamsurgeon.common import *
from shutil import move
from re import sub
from uuid import uuid4

import logging
FORMAT = '%(levelname)s %(asctime)s %(message)s'
logging.basicConfig(format=FORMAT)
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

#
# Remapping functions bam --> bam
#


supported_aligners_bam   = ['backtrack', 'mem', 'novoalign', 'gsnap', 'STAR', 'bowtie2', 'tmap', 'bwakit', 'minimap2']
supported_aligners_fastq = ['backtrack', 'mem', 'novoalign']

def checkoptions(name, options, picardjar, sv=False):
    ''' checks if necessary options have been specified '''

    if sv:
        if name not in supported_aligners_fastq:
            raise ValueError("ERROR unsupported aligner: " + name)
    else:
        if name not in supported_aligners_bam:
            raise ValueError("ERROR unsupported aligner: " + name)

    if name != 'backtrack' and not sv:
        if picardjar is None:
            raise ValueError("ERROR '--aligner " + name + "' requires '--picardjar' to be specified\n")

    if name == 'novoalign':
        if 'novoref' not in options:
            raise ValueError("ERROR '--aligner novoalign' requires '--alignopts novoref:path_to_novoalign_reference' to be set\n")

    if name == 'gsnap':
        if 'gsnaprefdir' not in options or 'gsnaprefname' not in options:
            raise ValueError("ERROR '--aligner gsnap' requires '--alignopts gsnaprefdir:GSNAP_reference_dir,gsnaprefname:GSNAP_ref_name\n")

    if name == 'STAR':
        if 'STARrefdir' not in options:
            raise ValueError("ERROR '--aligner STAR' requires '--alignopts STARrefdir:/path/to/STAR_reference_dir\n")

    if name == 'bowtie2':
        if 'bowtie2ref' not in options:
            raise ValueError("ERROR '--aligner bowtie2' requires '--alignopts bowtie2ref:bowtie2_ref_basepath\n")

    if name == 'minimap2':
        if 'x' not in options:
            raise ValueError("ERROR '--aligner minimap2' requires '--alignopts x:map-pb or x:map-ont\n")


def remap_bam(name, bamfn, fastaref, options, mutid='null', threads=1, paired=True, picardjar=None, insane=False):
    ''' remap bam file with supported alignment method. "options" param is a dict of aligner-specific required options '''

    checkoptions(name, options, picardjar)

    assert os.path.exists(bamfn), 'cannot locate bam: %s' % bamfn

    if picardjar is not None:
        assert os.path.exists(picardjar), 'cannot locate picard.jar: %s' % picardjar

    if name == 'backtrack':
        remap_backtrack_bam(bamfn, threads, fastaref, mutid=mutid, paired=paired)

    if name == 'mem':
        remap_bwamem_bam(bamfn, threads, fastaref, picardjar, mutid=mutid, paired=paired, insane=insane)

    if name == 'minimap2':
        remap_mm2_bam(bamfn, threads, fastaref, picardjar, options['x'], mutid=mutid, paired=paired, insane=insane)
        if paired:
            raise ValueError("ERROR --aligner minimap2 requires --single option to be set\n")

    if name == 'bwakit':
        remap_bwakit_bam(bamfn, threads, fastaref, picardjar, mutid=mutid, paired=paired, insane=insane)

    if name == 'novoalign':
        remap_novoalign_bam(bamfn, threads, fastaref, picardjar, options['novoref'], mutid=mutid, paired=paired, insane=insane)

    if name == 'gsnap':
        remap_gsnap_bam(bamfn, threads, fastaref, picardjar, options['gsnaprefdir'], options['gsnaprefname'], mutid=mutid, paired=paired, insane=insane)

    if name == 'STAR':
        remap_STAR_bam(bamfn, threads, fastaref, picardjar, options['STARrefdir'], mutid=mutid, paired=paired, insane=insane)

    if name == 'bowtie2':
        remap_bowtie2_bam(bamfn, threads, fastaref, picardjar, options['bowtie2ref'], mutid=mutid, paired=paired, insane=insane)

    if name == 'tmap':
        remap_tmap_bam(bamfn, threads, fastaref, picardjar, mutid=mutid, paired=paired, insane=insane)


def remap_backtrack_bam(bamfn, threads, fastaref, mutid='null', paired=True):
    """ call bwa/samtools to remap .bam
    """
    if paired:
        sai1fn = bamfn + ".1.sai"
        sai2fn = bamfn + ".2.sai"
        samfn  = bamfn + ".sam"
        refidx = fastaref + ".fai"

        sai1args = ['bwa', 'aln', '-q', '5', '-l', '32', '-k', '3', '-t', str(threads), '-o', '1', '-f', sai1fn, '-b1', fastaref, bamfn]
        sai2args = ['bwa', 'aln', '-q', '5', '-l', '32', '-k', '3', '-t', str(threads), '-o', '1', '-f', sai2fn, '-b2', fastaref, bamfn]
        samargs  = ['bwa', 'sampe', '-P', '-f', samfn, fastaref, sai1fn, sai2fn, bamfn, bamfn]
        bamargs  = ['samtools', 'view', '-bt', refidx, '-o', bamfn, samfn] 

        logger.info("%s mapping 1st end, cmd: %s" % (mutid, " ".join(sai1args)))
        subprocess.check_call(sai1args)
        logger.info("%s mapping 2nd end, cmd: %s" % (mutid, " ".join(sai2args)))
        subprocess.check_call(sai2args)
        logger.info("%s pairing ends, building .sam, cmd: %s" % (mutid, " ".join(samargs)))
        subprocess.check_call(samargs)
        logger.info("%s sam --> bam, cmd: %s" % (mutid, " ".join(bamargs)))
        subprocess.check_call(bamargs)

        sortfn = bamfn + ".sort.bam"
        sortargs = ['samtools','sort', '-T', sortfn ,'-o',sortfn, bamfn]
        logger.info("%s sorting, cmd: %s" % (mutid, " ".join(sortargs)))
        subprocess.check_call(sortargs)
        os.rename(sortfn, bamfn)

        indexargs = ['samtools','index',bamfn]
        logger.info("%s indexing, cmd: %s" % (mutid, " ".join(indexargs)))
        subprocess.check_call(indexargs)

        # cleanup
        os.remove(sai1fn)
        os.remove(sai2fn)
        os.remove(samfn)

    else:
        saifn = bamfn + ".sai"
        samfn  = bamfn + ".sam"
        refidx = fastaref + ".fai"

        saiargs = ['bwa', 'aln', '-q', '5', '-l', '32', '-k', '3', '-t', str(threads), '-o', '1', '-f', saifn, '-b1', fastaref, bamfn]
        samargs  = ['bwa', 'samse', '-f', samfn, fastaref, saifn, bamfn]
        bamargs  = ['samtools', 'view', '-bt', refidx, '-o', bamfn, samfn] 

        logger.info("%s mapping, cmd: %s" % (mutid, " ".join(saiargs)))
        subprocess.check_call(saiargs)
        logger.info("%s pairing ends, building .sam, cmd: %s" % (mutid, " ".join(samargs)))
        subprocess.check_call(samargs)
        logger.info("%s sam --> bam, cmd: %s" % (mutid, " ".join(bamargs)))
        subprocess.check_call(bamargs)

        sortfn = bamfn + ".sort.bam"
        sortargs = ['samtools','sort', '-T', sortfn, '-o',sortfn, bamfn]
        logger.info("%s sorting, cmd: %s" % (mutid, " ".join(sortargs)))
        subprocess.check_call(sortargs)
        os.rename(sortfn,bamfn)

        indexargs = ['samtools','index',bamfn]
        logger.info("%s indexing, cmd: %s" % (mutid, " ".join(indexargs)))
        subprocess.check_call(indexargs)

        # cleanup
        os.remove(saifn)
        os.remove(samfn)


def remap_bwamem_bam(bamfn, threads, fastaref, picardjar, mutid='null', paired=True, deltmp=True, insane=False):
    """ call bwa mem and samtools to remap .bam
    """
    assert os.path.exists(picardjar)
    assert bamreadcount(bamfn) > 0
    if paired:
        assert bamreadcount(bamfn) > 1 

    sam_out  = bamfn + '.realign.sam'
    sort_out = bamfn + '.realign.sorted.bam'

    logger.info("%s converting %s to fastq" % (mutid, bamfn))
    fastq = bamtofastq(bamfn, picardjar, threads=threads, paired=paired)[0]

    sam_cmd = []

    if paired:
        sam_cmd  = ['bwa', 'mem', '-t', str(threads), '-M', '-p', fastaref, fastq] # interleaved
    else:
        sam_cmd  = ['bwa', 'mem', '-t', str(threads), '-M', fastaref, fastq] # single-end

    assert len(sam_cmd) > 0

    bam_cmd  = ['samtools', 'view', '-bt', fastaref + '.fai', '-o', bamfn, sam_out]<commit_msg>fix: Add softclipping (-Y) to bwa mem command<commit_after>#!/usr/bin/env python


import os
import sys
import subprocess

from bamsurgeon.common import *
from shutil import move
from re import sub
from uuid import uuid4

import logging
FORMAT = '%(levelname)s %(asctime)s %(message)s'
logging.basicConfig(format=FORMAT)
logger = logging.getLogger(__name__)
logger.setLevel(logging.INFO)

#
# Remapping functions bam --> bam
#


supported_aligners_bam   = ['backtrack', 'mem', 'novoalign', 'gsnap', 'STAR', 'bowtie2', 'tmap', 'bwakit', 'minimap2']
supported_aligners_fastq = ['backtrack', 'mem', 'novoalign']

def checkoptions(name, options, picardjar, sv=False):
    ''' checks if necessary options have been specified '''

    if sv:
        if name not in supported_aligners_fastq:
            raise ValueError("ERROR unsupported aligner: " + name)
    else:
        if name not in supported_aligners_bam:
            raise ValueError("ERROR unsupported aligner: " + name)

    if name != 'backtrack' and not sv:
        if picardjar is None:
            raise ValueError("ERROR '--aligner " + name + "' requires '--picardjar' to be specified\n")

    if name == 'novoalign':
        if 'novoref' not in options:
            raise ValueError("ERROR '--aligner novoalign' requires '--alignopts novoref:path_to_novoalign_reference' to be set\n")

    if name == 'gsnap':
        if 'gsnaprefdir' not in options or 'gsnaprefname' not in options:
            raise ValueError("ERROR '--aligner gsnap' requires '--alignopts gsnaprefdir:GSNAP_reference_dir,gsnaprefname:GSNAP_ref_name\n")

    if name == 'STAR':
        if 'STARrefdir' not in options:
            raise ValueError("ERROR '--aligner STAR' requires '--alignopts STARrefdir:/path/to/STAR_reference_dir\n")

    if name == 'bowtie2':
        if 'bowtie2ref' not in options:
            raise ValueError("ERROR '--aligner bowtie2' requires '--alignopts bowtie2ref:bowtie2_ref_basepath\n")

    if name == 'minimap2':
        if 'x' not in options:
            raise ValueError("ERROR '--aligner minimap2' requires '--alignopts x:map-pb or x:map-ont\n")


def remap_bam(name, bamfn, fastaref, options, mutid='null', threads=1, paired=True, picardjar=None, insane=False):
    ''' remap bam file with supported alignment method. "options" param is a dict of aligner-specific required options '''

    checkoptions(name, options, picardjar)

    assert os.path.exists(bamfn), 'cannot locate bam: %s' % bamfn

    if picardjar is not None:
        assert os.path.exists(picardjar), 'cannot locate picard.jar: %s' % picardjar

    if name == 'backtrack':
        remap_backtrack_bam(bamfn, threads, fastaref, mutid=mutid, paired=paired)

    if name == 'mem':
        remap_bwamem_bam(bamfn, threads, fastaref, picardjar, mutid=mutid, paired=paired, insane=insane)

    if name == 'minimap2':
        remap_mm2_bam(bamfn, threads, fastaref, picardjar, options['x'], mutid=mutid, paired=paired, insane=insane)
        if paired:
            raise ValueError("ERROR --aligner minimap2 requires --single option to be set\n")

    if name == 'bwakit':
        remap_bwakit_bam(bamfn, threads, fastaref, picardjar, mutid=mutid, paired=paired, insane=insane)

    if name == 'novoalign':
        remap_novoalign_bam(bamfn, threads, fastaref, picardjar, options['novoref'], mutid=mutid, paired=paired, insane=insane)

    if name == 'gsnap':
        remap_gsnap_bam(bamfn, threads, fastaref, picardjar, options['gsnaprefdir'], options['gsnaprefname'], mutid=mutid, paired=paired, insane=insane)

    if name == 'STAR':
        remap_STAR_bam(bamfn, threads, fastaref, picardjar, options['STARrefdir'], mutid=mutid, paired=paired, insane=insane)

    if name == 'bowtie2':
        remap_bowtie2_bam(bamfn, threads, fastaref, picardjar, options['bowtie2ref'], mutid=mutid, paired=paired, insane=insane)

    if name == 'tmap':
        remap_tmap_bam(bamfn, threads, fastaref, picardjar, mutid=mutid, paired=paired, insane=insane)


def remap_backtrack_bam(bamfn, threads, fastaref, mutid='null', paired=True):
    """ call bwa/samtools to remap .bam
    """
    if paired:
        sai1fn = bamfn + ".1.sai"
        sai2fn = bamfn + ".2.sai"
        samfn  = bamfn + ".sam"
        refidx = fastaref + ".fai"

        sai1args = ['bwa', 'aln', '-q', '5', '-l', '32', '-k', '3', '-t', str(threads), '-o', '1', '-f', sai1fn, '-b1', fastaref, bamfn]
        sai2args = ['bwa', 'aln', '-q', '5', '-l', '32', '-k', '3', '-t', str(threads), '-o', '1', '-f', sai2fn, '-b2', fastaref, bamfn]
        samargs  = ['bwa', 'sampe', '-P', '-f', samfn, fastaref, sai1fn, sai2fn, bamfn, bamfn]
        bamargs  = ['samtools', 'view', '-bt', refidx, '-o', bamfn, samfn] 

        logger.info("%s mapping 1st end, cmd: %s" % (mutid, " ".join(sai1args)))
        subprocess.check_call(sai1args)
        logger.info("%s mapping 2nd end, cmd: %s" % (mutid, " ".join(sai2args)))
        subprocess.check_call(sai2args)
        logger.info("%s pairing ends, building .sam, cmd: %s" % (mutid, " ".join(samargs)))
        subprocess.check_call(samargs)
        logger.info("%s sam --> bam, cmd: %s" % (mutid, " ".join(bamargs)))
        subprocess.check_call(bamargs)

        sortfn = bamfn + ".sort.bam"
        sortargs = ['samtools','sort', '-T', sortfn ,'-o',sortfn, bamfn]
        logger.info("%s sorting, cmd: %s" % (mutid, " ".join(sortargs)))
        subprocess.check_call(sortargs)
        os.rename(sortfn, bamfn)

        indexargs = ['samtools','index',bamfn]
        logger.info("%s indexing, cmd: %s" % (mutid, " ".join(indexargs)))
        subprocess.check_call(indexargs)

        # cleanup
        os.remove(sai1fn)
        os.remove(sai2fn)
        os.remove(samfn)

    else:
        saifn = bamfn + ".sai"
        samfn  = bamfn + ".sam"
        refidx = fastaref + ".fai"

        saiargs = ['bwa', 'aln', '-q', '5', '-l', '32', '-k', '3', '-t', str(threads), '-o', '1', '-f', saifn, '-b1', fastaref, bamfn]
        samargs  = ['bwa', 'samse', '-f', samfn, fastaref, saifn, bamfn]
        bamargs  = ['samtools', 'view', '-bt', refidx, '-o', bamfn, samfn] 

        logger.info("%s mapping, cmd: %s" % (mutid, " ".join(saiargs)))
        subprocess.check_call(saiargs)
        logger.info("%s pairing ends, building .sam, cmd: %s" % (mutid, " ".join(samargs)))
        subprocess.check_call(samargs)
        logger.info("%s sam --> bam, cmd: %s" % (mutid, " ".join(bamargs)))
        subprocess.check_call(bamargs)

        sortfn = bamfn + ".sort.bam"
        sortargs = ['samtools','sort', '-T', sortfn, '-o',sortfn, bamfn]
        logger.info("%s sorting, cmd: %s" % (mutid, " ".join(sortargs)))
        subprocess.check_call(sortargs)
        os.rename(sortfn,bamfn)

        indexargs = ['samtools','index',bamfn]
        logger.info("%s indexing, cmd: %s" % (mutid, " ".join(indexargs)))
        subprocess.check_call(indexargs)

        # cleanup
        os.remove(saifn)
        os.remove(samfn)


def remap_bwamem_bam(bamfn, threads, fastaref, picardjar, mutid='null', paired=True, deltmp=True, insane=False):
    """ call bwa mem and samtools to remap .bam
    """
    assert os.path.exists(picardjar)
    assert bamreadcount(bamfn) > 0
    if paired:
        assert bamreadcount(bamfn) > 1 

    sam_out  = bamfn + '.realign.sam'
    sort_out = bamfn + '.realign.sorted.bam'

    logger.info("%s converting %s to fastq" % (mutid, bamfn))
    fastq = bamtofastq(bamfn, picardjar, threads=threads, paired=paired)[0]

    sam_cmd = []

    if paired:
        sam_cmd  = ['bwa', 'mem', '-Y', '-t', str(threads), '-M', '-p', fastaref, fastq] # interleaved
    else:
        sam_cmd  = ['bwa', 'mem', '-Y', '-t', str(threads), '-M', fastaref, fastq] # single-end

    assert len(sam_cmd) > 0

    bam_cmd  = ['samtools', 'view', '-bt', fastaref + '.fai', '-o', bamfn, sam_out]
{'fraction_non_alphanumeric': 0.16241507499038585, 'fraction_numerical': 0.007306755544161005, 'mean_word_length': 3.164131305044035, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '7020501', 'n_tokens_mistral': 5809, 'n_tokens_neox': 5590, 'n_words': 1350}