<commit_before>import sys
import os
import os.path as op
import argparse
from toolshed import nopen

def faseq(fa, chrom, start, end, cache=[None]):
    """
    this is called by pileup which is ordered by chrom
    so we can speed things up by reading in a chrom at
    a time into memory
    """
    if cache[0] is None or cache[0][0] != chrom:
        seq = "".join(x.strip() for i, x in
            enumerate(nopen("|samtools faidx %s %s" % (fa, chrom))) if i >
            0).upper()
        cache[0] = (chrom, seq)
    chrom, seq = cache[0]
    return seq[start - 1: end]

def get_context(seq5, forward):
    """
    >>> get_context('GACGG', True)
    'CG+'                  
    """
    if forward:
        assert seq5[2] == "C", seq5
        if seq5[3] == "G": return "CG+"
        if seq5[4] == "G": return "CHG+"
        return "CHH+"
    else: # reverse complement
        assert seq5[2] == "G", seq5
        if seq5[1] == "C": return "CG-"
        if seq5[0] == "C": return "CHG-"
        return "CHH-"

def tabulate_main(args):
    __doc__ = """
    tabulate methylation from bwa-meth.py call
    """
    p = argparse.ArgumentParser(__doc__)
    p.add_argument("--reference", help="reference fasta", required=True)
    p.add_argument("--region", help="specific region in which to run mpileup",
            default="")
    p.add_argument("--g-only", default=False, action="store_true",
        help="only use information from the reverse strand (G->A conversions)")
    p.add_argument("--map-q", "-q", type=int, default=10, help="only tabulate "
                   "methylation for reads with at least this mapping quality")
    p.add_argument("--base-q", "-Q", type=int, default=10, help="only tabulate "
                   "methylation for bases with at least this quality")
    p.add_argument("--skip-left", type=int, default=3, help="don't use the "
            "first N bases from the start of the read to avoid bias")
    p.add_argument("--skip-right", type=int, default=3, help="don't use the "
            "last N bases from the end of the read to avoid bias")
    p.add_argument("--read-length", type=int, help="length of reads"
            " before trimming for used with skipping", required=True)
    p.add_argument("bams", nargs="+")

    a = p.parse_args(args)
    assert os.path.exists(a.reference)
    if a.region:
        a.region = ("-l " if op.exists(a.region) else "-r ") + a.region
    cmd = "|samtools mpileup -Of {reference} {region} -d100000 -BQ {base_q} -q {map_q} {bams}"
    cmd = cmd.format(reference=a.reference, map_q=a.map_q, base_q=a.base_q,
                     bams=" ".join(a.bams), region=a.region)
    print >>sys.stderr, "generating pileup with command:", cmd
    samples = [op.basename(b)[:-4] for b in a.bams]
    tabulate_methylation(cmd, a.reference, samples, a.g_only, a.skip_left,
            a.read_length - a.skip_right)

def tabulate_methylation(fpileup, reference, samples, g_only=False,
        skip_left=1, skip_right=99):

    conversion = {"C": "T", "G": "A"}
    fhs = []
    for sample in samples:
        fhs.append(open('%s.methylation.txt' % sample, 'w'))
        fhs[-1].write("#chrom\tpos0\tpos1\tpct\tn_same\tn_converted\n")

    import re
    # regexp to remove the ^, $ stuff
    end_re = re.compile("\^.|\$")
    deletion_re = re.compile("(-\d+[ACTGNatcgn]+)")

    insertion_re = re.compile("\+(\d+)([ACTGNatcgn]+)")
    deletion_re = re.compile("-(\d+)([ACTGNatcgn]+)")

    def del_sub(m):
        """ with a base string like "-1ACTG" just return CTG
        for -3CCCTG return TG
        """
        n = int(m.groups()[0])
        bases = m.groups()[1]
        return bases[n:]
<commit_msg>add filtering/required flags<commit_after>import sys
import os
import os.path as op
import argparse
from toolshed import nopen

def faseq(fa, chrom, start, end, cache=[None]):
    """
    this is called by pileup which is ordered by chrom
    so we can speed things up by reading in a chrom at
    a time into memory
    """
    if cache[0] is None or cache[0][0] != chrom:
        seq = "".join(x.strip() for i, x in
            enumerate(nopen("|samtools faidx %s %s" % (fa, chrom))) if i >
            0).upper()
        cache[0] = (chrom, seq)
    chrom, seq = cache[0]
    return seq[start - 1: end]

def get_context(seq5, forward):
    """
    >>> get_context('GACGG', True)
    'CG+'                  
    """
    if forward:
        assert seq5[2] == "C", seq5
        if seq5[3] == "G": return "CG+"
        if seq5[4] == "G": return "CHG+"
        return "CHH+"
    else: # reverse complement
        assert seq5[2] == "G", seq5
        if seq5[1] == "C": return "CG-"
        if seq5[0] == "C": return "CHG-"
        return "CHH-"

def tabulate_main(args):
    __doc__ = """
    tabulate methylation from bwa-meth.py call
    """
    p = argparse.ArgumentParser(__doc__)
    p.add_argument("--reference", help="reference fasta", required=True)
    p.add_argument("--region", help="specific region in which to run mpileup",
            default="")
    p.add_argument("--g-only", default=False, action="store_true",
        help="only use information from the reverse strand (G->A conversions)")
    p.add_argument("--map-q", "-q", type=int, default=10, help="only tabulate "
                   "methylation for reads with at least this mapping quality")
    p.add_argument("--base-q", "-Q", type=int, default=10, help="only tabulate "
                   "methylation for bases with at least this quality")
    p.add_argument("--skip-left", type=int, default=3, help="don't use the "
            "first N bases from the start of the read to avoid bias")
    p.add_argument("--skip-right", type=int, default=3, help="don't use the "
            "last N bases from the end of the read to avoid bias")
    p.add_argument("--read-length", type=int, help="length of reads"
            " before trimming for used with skipping", required=True)
    p.add_argument("--rf", action='append', default=[0x2], type=int)
    p.add_argument("--ff", action='append', default=[0x200, 0x400, 0x800],
            type=int)
    p.add_argument("bams", nargs="+")

    a = p.parse_args(args)
    assert os.path.exists(a.reference)
    if a.region:
        a.region = ("-l " if op.exists(a.region) else "-r ") + a.region
    ff = 0
    for flag in a.ff: ff |= flag
    rf = 0
    for flag in a.rf: rf |= flag
    cmd = ("|samtools mpileup --rf {rf} --ff {ff} -Of {reference} {region}"
           " -d100000 -BQ {base_q} -q {map_q} {bams}")

    cmd = cmd.format(reference=a.reference, map_q=a.map_q, base_q=a.base_q,
                     bams=" ".join(a.bams), region=a.region, rf=rf, ff=ff)
    print >>sys.stderr, "generating pileup with command:", cmd
    samples = [op.basename(b)[:-4] for b in a.bams]
    tabulate_methylation(cmd, a.reference, samples, a.g_only, a.skip_left,
            a.read_length - a.skip_right)

def tabulate_methylation(fpileup, reference, samples, g_only=False,
        skip_left=1, skip_right=99):

    conversion = {"C": "T", "G": "A"}
    fhs = []
    for sample in samples:
        fhs.append(open('%s.methylation.txt' % sample, 'w'))
        fhs[-1].write("#chrom\tpos0\tpos1\tpct\tn_same\tn_converted\n")

    import re
    # regexp to remove the ^, $ stuff
    end_re = re.compile("\^.|\$")
    deletion_re = re.compile("(-\d+[ACTGNatcgn]+)")

    insertion_re = re.compile("\+(\d+)([ACTGNatcgn]+)")
    deletion_re = re.compile("-(\d+)([ACTGNatcgn]+)")

    def del_sub(m):
        """ with a base string like "-1ACTG" just return CTG
        for -3CCCTG return TG
        """
        n = int(m.groups()[0])
        bases = m.groups()[1]
        return bases[n:]

{'fraction_non_alphanumeric': 0.16000522466039707, 'fraction_numerical': 0.013845350052246604, 'mean_word_length': 2.579710144927536, 'pattern_counts': {'":': 12, '<': 3, '<?xml version=': 0, '>': 17, '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': '5193520', 'n_tokens_mistral': 2712, 'n_tokens_neox': 2608, 'n_words': 722}