<commit_before>#!/usr/bin/env python
"""
Given a BAM file, this script will only allow reads that meet filtering 
criteria to be written to output. The output is another BAM file with the 
reads not matching the criteria removed.

Currently, the available filters are:
    -mapped                    Keep only mapped reads
    
    -mismatch num              # mismatches or indels
                               indel always counts as 1 regardless of length
                               (requires NM tag in reads)
    -mismatch_ref num ref.fa   # mismatches or indel - looks up mismatches 
                               directly in a ref FASTA file (if NM not 
                               available)
                            
    -nosecondary               Remove reads that have the 0x100 flag set
    -noqcfail                  Remove reads that have the 0x200 flag set
    
    -mask bitmask              Remove reads that match the mask (base 10/hex)

    -exclude ref:start-end     Remove reads in this region (1-based start)
    -excludebed file.bed       Remove reads that are in any of the regions 
                               from the given BED file
                               
    -include ref:start-end     Remove reads NO this region (can only be one)
    -includebed file.bed       Remove reads that are NOT any of the regions 
                               from the given BED file
                               Note: If this is a large dataset, use 
                               "bamutils extract" instead.

    -whitelist fname           Remove reads that aren't on this list (by name)
    -blacklist fname           Remove reads that are on this list (by name)
    
    -eq  tag_name value
    -lt  tag_name value
    -lte tag_name value
    -gt  tag_name value
    -gte tag_name value

    Where tag_name should be the full name, plus the type eg: AS:i
    
Common tags to filter by:
    AS:i    Alignment score
    NM:i    Edit distance (each indel counts as many as it's length)
    IH:i    Number of alignments
    
"""

import os,sys,gzip
from support.eta import ETA
import pysam

def usage():
    print __doc__
    print """
Usage: bamutils filter in.bam out.bam {-failed out.txt} criteria...

Options:

  -failed fname    A text file containing the read names of all reads 
                   that were removed with filtering

Example: 
bamutils filter filename.bam output.bam -mapped -gte AS:i 1000

This will remove all unmapped reads, as well as any reads that have an AS:i 
value less than 1000.
"""
    sys.exit(1)

bam_cigar = ['M','I','D','N','S','H','P']

class Blacklist(object):
    def __init__(self,fname):
        self.fname = fname
        self.notallowed = []
        with open(fname) as f:
            for line in f:
                self.notallowed.append(line.strip())
    def filter(self,bam,read):
        return read.qname not in self.notallowed
    def __repr__(self):
        return 'Blacklist: %s' % (self.fname)

class Whitelist(object):
    def __init__(self,fname):
        self.fname = fname
        self.allowed = []
        with open(fname) as f:
            for line in f:
                self.allowed.append(line.strip())
    def filter(self,bam,read):
        return read.qname in self.allowed
    def __repr__(self):
        return 'Whitelist: %s' % (self.fname)


class IncludeRegion(object):
    _excludes = []
    _last = None
    def __init__(self,region):
        IncludeRegion._excludes.append(ExcludeRegion(region))
        # self.excl = ExcludeRegion(region)
    def filter(self,bam,read):
        if read == IncludeRegion._last:
            return True
        <commit_msg>updated whitelist/blacklist to support tab/space delimited (such as output from bamutils reads)<commit_after>#!/usr/bin/env python
"""
Given a BAM file, this script will only allow reads that meet filtering 
criteria to be written to output. The output is another BAM file with the 
reads not matching the criteria removed.

Currently, the available filters are:
    -mapped                    Keep only mapped reads
    
    -mismatch num              # mismatches or indels
                               indel always counts as 1 regardless of length
                               (requires NM tag in reads)
    -mismatch_ref num ref.fa   # mismatches or indel - looks up mismatches 
                               directly in a ref FASTA file (if NM not 
                               available)
                            
    -nosecondary               Remove reads that have the 0x100 flag set
    -noqcfail                  Remove reads that have the 0x200 flag set
    
    -mask bitmask              Remove reads that match the mask (base 10/hex)

    -exclude ref:start-end     Remove reads in this region (1-based start)
    -excludebed file.bed       Remove reads that are in any of the regions 
                               from the given BED file
                               
    -include ref:start-end     Remove reads NO this region (can only be one)
    -includebed file.bed       Remove reads that are NOT any of the regions 
                               from the given BED file
                               Note: If this is a large dataset, use 
                               "bamutils extract" instead.

    -whitelist fname           Remove reads that aren't on this list (by name)
    -blacklist fname           Remove reads that are on this list (by name)
                                 These lists can be whitespace-delimited with 
                                 the read name as the first column.
    
    -eq  tag_name value
    -lt  tag_name value
    -lte tag_name value
    -gt  tag_name value
    -gte tag_name value

    Where tag_name should be the full name, plus the type eg: AS:i
    
Common tags to filter by:
    AS:i    Alignment score
    NM:i    Edit distance (each indel counts as many as it's length)
    IH:i    Number of alignments
    
"""

import os,sys,gzip
from support.eta import ETA
import pysam

def usage():
    print __doc__
    print """
Usage: bamutils filter in.bam out.bam {-failed out.txt} criteria...

Options:
  -failed fname    A text file containing the read names of all reads 
                   that were removed with filtering

Example: 
bamutils filter filename.bam output.bam -mapped -gte AS:i 1000

This will remove all unmapped reads, as well as any reads that have an AS:i 
value less than 1000.
"""
    sys.exit(1)

bam_cigar = ['M','I','D','N','S','H','P']

class Blacklist(object):
    def __init__(self,fname):
        self.fname = fname
        self.notallowed = []
        with open(fname) as f:
            for line in f:
                self.notallowed.append(line.strip().split()[0])
    def filter(self,bam,read):
        return read.qname not in self.notallowed
    def __repr__(self):
        return 'Blacklist: %s' % (self.fname)

class Whitelist(object):
    def __init__(self,fname):
        self.fname = fname
        self.allowed = []
        with open(fname) as f:
            for line in f:
                self.allowed.append(line.strip().split()[0])
    def filter(self,bam,read):
        return read.qname in self.allowed
    def __repr__(self):
        return 'Whitelist: %s' % (self.fname)


class IncludeRegion(object):
    _excludes = []
    _last = None
    def __init__(self,region):
        IncludeRegion._excludes.append(ExcludeRegion(region))
        # self.excl = ExcludeRegion(region)
    def filter(self,bam,read):
        if read == IncludeRegion._last:
            return True
        
{'fraction_non_alphanumeric': 0.06871035940803383, 'fraction_numerical': 0.005813953488372093, 'mean_word_length': 1.7940199335548173, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1644677', 'n_tokens_mistral': 2184, 'n_tokens_neox': 2063, 'n_words': 802}