<commit_before>"""
SAM alignment format. There are other tools that handles better SAM and BAM.
This script simply parses the lines in SAM into human readable fields.

http://samtools.sourceforge.net/SAM1.pdf
"""

import os.path as op
import sys
import logging

from itertools import groupby
from optparse import OptionParser

from Bio import SeqIO
from jcvi.formats.base import LineFile
from jcvi.formats.fasta import Fasta
from jcvi.utils.cbook import fill
from jcvi.assembly.base import Astat
from jcvi.apps.base import ActionDispatcher, need_update, sh, debug, set_outfile
debug()


class SamLine (object):

    def __init__(self, row):

        args = row.strip().split("\t")
        self.qname = args[0]
        self.flag = args[1]
        self.rname = args[2]
        self.pos = args[3]
        self.mapq = args[4]
        self.cigar = args[5]
        self.mrnm = args[6]
        self.mpos = args[7]
        self.isize = args[8]
        self.seq = args[9]
        self.qual = args[10]

    @property
    def pairline(self):
        qpos = self.cigar.split('H', 1)[0]
        return "%s:%s\t%s:%s" % (self.qname, qpos, self.rname, self.pos)


class Sam (LineFile):

    def __init__(self, filename, callback=None):

        fp = open(filename)
        for row in fp:
            if row[0] == '@':
                continue
            s = SamLine(row)
            if callback:
                callback(s)


def main():

    actions = (
        ('pair', 'parse sam file and get pairs'),
        ('pairs', 'print paired-end reads from BAM file'),
        ('chimera', 'parse sam file from `bwasw` and list multi-hit reads'),
        ('ace', 'convert sam file to ace'),
        ('index', 'convert to bam, sort and then index'),
        ('consensus', 'convert bam alignments to consensus FASTA'),
        ('bcf', 'run mpileup on a set of bam files'),
            )

    p = ActionDispatcher(actions)
    p.dispatch(globals())


def pairs(args):
    """
    See __doc__ for set_options_pairs().
    """
    import jcvi.formats.bed
    from jcvi.formats.blast import set_options_pairs
<commit_msg>add formats.sam.fpkm to convert BAM file to FPKM values<commit_after>"""
SAM alignment format. There are other tools that handles better SAM and BAM.
This script simply parses the lines in SAM into human readable fields.

http://samtools.sourceforge.net/SAM1.pdf
"""

import os.path as op
import sys
import logging

from itertools import groupby
from optparse import OptionParser

from Bio import SeqIO
from jcvi.formats.base import LineFile
from jcvi.formats.fasta import Fasta
from jcvi.utils.cbook import fill
from jcvi.assembly.base import Astat
from jcvi.apps.base import ActionDispatcher, need_update, sh, debug, set_outfile
debug()


class SamLine (object):

    def __init__(self, row):

        args = row.strip().split("\t")
        self.qname = args[0]
        self.flag = args[1]
        self.rname = args[2]
        self.pos = args[3]
        self.mapq = args[4]
        self.cigar = args[5]
        self.mrnm = args[6]
        self.mpos = args[7]
        self.isize = args[8]
        self.seq = args[9]
        self.qual = args[10]

    @property
    def pairline(self):
        qpos = self.cigar.split('H', 1)[0]
        return "%s:%s\t%s:%s" % (self.qname, qpos, self.rname, self.pos)


class Sam (LineFile):

    def __init__(self, filename, callback=None):

        fp = open(filename)
        for row in fp:
            if row[0] == '@':
                continue
            s = SamLine(row)
            if callback:
                callback(s)


def main():

    actions = (
        ('pair', 'parse sam file and get pairs'),
        ('pairs', 'print paired-end reads from BAM file'),
        ('chimera', 'parse sam file from `bwasw` and list multi-hit reads'),
        ('ace', 'convert sam file to ace'),
        ('index', 'convert to bam, sort and then index'),
        ('consensus', 'convert bam alignments to consensus FASTA'),
        ('fpkm', 'calculate FPKM values from BAM file'),
        ('bcf', 'run mpileup on a set of bam files'),
            )

    p = ActionDispatcher(actions)
    p.dispatch(globals())


def fpkm(args):
    """
    %prog fpkm fastafile *.bam

    Calculate FPKM values from BAM file.
    """
    from jcvi.formats.fasta import Fasta

    p = OptionParser(fpkm.__doc__)
    opts, args = p.parse_args(args)

    if len(args) < 2:
        sys.exit(not p.print_help())

    fastafile = args[0]
    bamfiles = args[1:]
    # Create a DUMMY gff file for cufdiff
    gffile = fastafile.rsplit(".", 1)[0] + ".gff"
    if need_update(fastafile, gffile):
        fw = open(gffile, "w")
        f = Fasta(fastafile)
        for key, size in f.itersizes_ordered():
            print >> fw, "\t".join(str(x) for x in (key, "dummy", "transcript",\
                1, size, ".", ".", ".", "ID=" + key))
        fw.close()
        logging.debug("Dummy GFF created: {0}".format(gffile))

    cmd = "cuffdiff {0} {1}".format(gffile, " ".join(bamfiles))
    sh(cmd)


def pairs(args):
    """
    See __doc__ for set_options_pairs().
    """
    import jcvi.formats.bed
    from jcvi.formats.blast import set_options_pairs

{'fraction_non_alphanumeric': 0.12599495243642012, 'fraction_numerical': 0.007959619491360901, 'mean_word_length': 2.4096624751819986, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 5, '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': '2170595', 'n_tokens_mistral': 1752, 'n_tokens_neox': 1680, 'n_words': 507}