<commit_before>"""
This module contains various functions and classes for handling variant
information, including utilities for generating combinations of variants, 
haplotypes, and genotypes.
"""
import copy
import logging
import vcf

logger = logging.getLogger('Log')

#####################################################################

class VariantCandidateReader(object):
    """
    A class to read variant information from vcf files, and return
    a batch of stream of variant instances.

    Required PyVCF
    """

    def __init__(self, filenames, options = None):

        """
        Constructor. Takes the vcf files (bgzip format) and return
        vcf stream for the variant instances.
        """

        self.options     = options
        self.vcf_readers = []

        for filename in filenames:

            if '.gz' not in filename:
                logger.error(
                    '\nSource File %s does not look like a bgzipped VCF '
                    '(i.e. name does not end in .gz)\nYou must supply a '
                    'VCF that has been compressed with bgzip, and indexed '
                    'with tabix' % (filename))
                logger.error(
                    'Compress the VCF using bgzip: '
                    'bgzip %s --> %s.gz' % (filename, filename))
                logger.error(
                    'Remember to use the "tabix -p vcf %s.gz" command to '
                    'index the compressed file' % (filename))
                logger.error(
                    'This should create an index file: %s.gz.tbi\n' % (filename))
                raise ValueError('\nInput VCF source file %s was not compressed '
                                 'and indexed' % (filename))
            else:
                # How To: how to close the open vcf files' handle?
                self.vcf_readers.append(vcf.Reader(filename = filename))

    def variants(self, chrom = None, start = None, end = None, nosnp = True):
        """
        *(chrom = None, start = None, end = None, nosnp = True)*
        Generator funtion. Yields variants in order of genomic coordinate.

        fetch variants in a region using 0-based indexing.
        
        The region is specified by :term:`chrom`, *start* and *end*.

        fetch returns an empty string if the region is out of range or
        addresses an unknown *chrom*.

        If *chrom* is given and *start* is None, the sequence from the
        first base is returned. Similarly, if *end* is None, the sequence
        until the last base is returned.

        """
        varlist = []
        # TO DO: All the variant in the same position could make by a net
        for vcf_reader in self.vcf_readers:

            for r in vcf_reader.fetch(chrom, start, end):

                # Continue, if the ALT == '.'
                if r.ALT[0] is None: continue

                # ignore the information that we don't care
                r.INFO    = None
                r.FORMAT  = None
                r.samples = None

                if nosnp and r.is_snp:

                    continue
                elif r.is_snp:

                    varlist.append(r)
                else:
                    # Copy the VCF row, and all the change will just happen 
                    # in this new copy. It would be very useful when we assign
                    # new record to 'varlist'
                    record = copy.copy(r)

                    # TO DO: This trim strategy is greedy algorithm, which is 
                    # not the best method. 
                    # e.g: Assume the turth is [A, ATC] 
                    # If we see [AA, ATCA] then will be [A, TCA] instead of
                    # [A, ATC] after triming 
                    for alt in r.ALT:
                        # Non snp variants may leading and/or trailing bases 
                        # trimming.
                        # For indel the first base will always be the same with
                        # REF. That will not be trim!

                        pos, ref = r.POS, r.REF
                        
                        # Trim the leading bases, should keep at lest 1 base 
                        while (alt.sequence[0].upper() == ref[0].upper() and 
                               len(ref) > 1 and len(alt) > 1):
                            alt.sequence = alt.sequence[1:]
                            ref  = ref[1:]
                            pos += 1

                        # Trim the trailing bases, should keep at lest 1 base
                        while (alt.sequence[-1].upper() == ref[-1].upper() and
                               len(ref) > 1 and len(alt) > 1):
                            alt.sequence = alt.sequence[:-1]
                            ref = ref[:-1]

                        if len(ref) > 0 and len(alt) > 0:
                            record.POS = pos
                            record.REF = ref
                            record.ALT = [alt]
                            # After the 'for loop', there may contain some 
                            # duplication positions in 'varlist'
                            varlist.append(record)

        varlist = self._dedup(varlist)
        logger.debug('Found %s variants in region %s in source file' 
                     % (len(varlist), '%s:%s-%s' % (chrom, start, end)))

        # It's a list of '_Record' which type is defined by 'PyVCF'
        return varlist

    def _dedup(self, varlist):
        """
        Delete conflict variants at the same position.  

        And keep the small one
        """
        if len(varlist) == 0: 
            return [] 

        vdict = {}
        # Firstly, we should the duplicate positions
        for i, v in enumerate(varlist):
            k = v.CHROM + ':' + str(v.POS)
            vdict[k] = vdict.get(k, []) + [i]

        variants = []
        for k, idx in vdict.items():

            if len(idx) == 1:
                variants.append(varlist[idx[0]])
            else:
                """
                delete conflict here
                """
                prevar = varlist[idx[0]]
                for i in idx[1:]:

                    if prevar.REF == varlist[i].REF:

                        prevar.ALT += varlist[i].ALT

                    elif varlist[i].is_snp:

                        continue
                    else:
                        # Keep the smaller one
                        pre_max_vlen  = max([abs(len(prevar.REF) - len(v)) 
                                             for v in prevar.AL])
                        this_max_vlen = max([abs(len(varlist[i].REF) - len(v)) 
                                             for v in varlist[i].ALT])
                        if this_max_vlen < pre_max_vlen:
                            prevar = varlist[i]
                        
                prevar.ALT = list(set(prevar.ALT)) # Unique the sequence
                variants.append(prevar)

        return sorted(list(set(variants))) # Sorted by reference pos order 
        

def get_sequence_context(fa_stream, variant):
    """
    Return the sequence surrounding this variant's position.
    """

    start = max(0, variant.POS - 10)
    return fa_stream.fetch(variant.CHROM, start, variant.POS + 10)

def homoRunForOneVariant(fa_stream, variant):

    """
    Calculate and return the length of the largest homopolymer<commit_msg>a bug<commit_after>"""
This module contains various functions and classes for handling variant
information, including utilities for generating combinations of variants, 
haplotypes, and genotypes.
"""
import copy
import logging
import vcf

logger = logging.getLogger('Log')

#####################################################################

class VariantCandidateReader(object):
    """
    A class to read variant information from vcf files, and return
    a batch of stream of variant instances.

    Required PyVCF
    """

    def __init__(self, filenames, options = None):

        """
        Constructor. Takes the vcf files (bgzip format) and return
        vcf stream for the variant instances.
        """

        self.options     = options
        self.vcf_readers = []

        for filename in filenames:

            if '.gz' not in filename:
                logger.error(
                    '\nSource File %s does not look like a bgzipped VCF '
                    '(i.e. name does not end in .gz)\nYou must supply a '
                    'VCF that has been compressed with bgzip, and indexed '
                    'with tabix' % (filename))
                logger.error(
                    'Compress the VCF using bgzip: '
                    'bgzip %s --> %s.gz' % (filename, filename))
                logger.error(
                    'Remember to use the "tabix -p vcf %s.gz" command to '
                    'index the compressed file' % (filename))
                logger.error(
                    'This should create an index file: %s.gz.tbi\n' % (filename))
                raise ValueError('\nInput VCF source file %s was not compressed '
                                 'and indexed' % (filename))
            else:
                # How To: how to close the open vcf files' handle?
                self.vcf_readers.append(vcf.Reader(filename = filename))

    def variants(self, chrom = None, start = None, end = None, nosnp = True):
        """
        *(chrom = None, start = None, end = None, nosnp = True)*
        Generator funtion. Yields variants in order of genomic coordinate.

        fetch variants in a region using 0-based indexing.
        
        The region is specified by :term:`chrom`, *start* and *end*.

        fetch returns an empty string if the region is out of range or
        addresses an unknown *chrom*.

        If *chrom* is given and *start* is None, the sequence from the
        first base is returned. Similarly, if *end* is None, the sequence
        until the last base is returned.

        """
        varlist = []
        # TO DO: All the variant in the same position could make by a net
        for vcf_reader in self.vcf_readers:

            for r in vcf_reader.fetch(chrom, start, end):

                # Continue, if the ALT == '.'
                if r.ALT[0] is None: continue

                # ignore the information that we don't care
                r.INFO    = None
                r.FORMAT  = None
                r.samples = None

                if nosnp and r.is_snp:

                    continue
                elif r.is_snp:

                    varlist.append(r)
                else:
                    # Copy the VCF row, and all the change will just happen 
                    # in this new copy. It would be very useful when we assign
                    # new record to 'varlist'
                    record = copy.copy(r)

                    # TO DO: This trim strategy is greedy algorithm, which is 
                    # not the best method. 
                    # e.g: Assume the turth is [A, ATC] 
                    # If we see [AA, ATCA] then will be [A, TCA] instead of
                    # [A, ATC] after triming 
                    for alt in r.ALT:
                        # Non snp variants may leading and/or trailing bases 
                        # trimming.
                        # For indel the first base will always be the same with
                        # REF. That will not be trim!

                        pos, ref = r.POS, r.REF
                        
                        # Trim the leading bases, should keep at lest 1 base 
                        while (alt.sequence[0].upper() == ref[0].upper() and 
                               len(ref) > 1 and len(alt) > 1):
                            alt.sequence = alt.sequence[1:]
                            ref  = ref[1:]
                            pos += 1

                        # Trim the trailing bases, should keep at lest 1 base
                        while (alt.sequence[-1].upper() == ref[-1].upper() and
                               len(ref) > 1 and len(alt) > 1):
                            alt.sequence = alt.sequence[:-1]
                            ref = ref[:-1]

                        if len(ref) > 0 and len(alt) > 0:
                            record.POS = pos
                            record.REF = ref
                            record.ALT = [alt]
                            # After the 'for loop', there may contain some 
                            # duplication positions in 'varlist'
                            varlist.append(record)

        varlist = self._dedup(varlist)
        logger.debug('Found %s variants in region %s in source file' 
                     % (len(varlist), '%s:%s-%s' % (chrom, start, end)))

        # It's a list of '_Record' which type is defined by 'PyVCF'
        return varlist

    def _dedup(self, varlist):
        """
        Delete conflict variants at the same position.  

        And keep the small one
        """
        if len(varlist) == 0: 
            return [] 

        vdict = {}
        # Firstly, we should the duplicate positions
        for i, v in enumerate(varlist):
            k = v.CHROM + ':' + str(v.POS)
            vdict[k] = vdict.get(k, []) + [i]

        variants = []
        for k, idx in vdict.items():

            if len(idx) == 1:
                variants.append(varlist[idx[0]])
            else:
                """
                delete conflict here
                """
                prevar = varlist[idx[0]]
                for i in idx[1:]:

                    if prevar.is_snp and (not varlist[i].is_snp):
                        # Ignore SNP,if we find other variant type  
                        prevar = varlist[i]
                    elif prevar.REF == varlist[i].REF:
                        # Merge variants if they have the same REF, SNP
                        # will be merge here, too
                        prevar.ALT += varlist[i].ALT
                    elif varlist[i].is_snp:
                        # Ignore SNP,if we find other variant type
                        continue
                    else:
                        # Keep the smaller one
                        pre_max_vlen  = max([abs(len(prevar.REF) - len(v)) 
                                             for v in prevar.AL])
                        this_max_vlen = max([abs(len(varlist[i].REF) - len(v)) 
                                             for v in varlist[i].ALT])
                        if this_max_vlen < pre_max_vlen:
                            prevar = varlist[i]
                        
                prevar.ALT = list(set(prevar.ALT)) # Unique the sequence
                variants.append(prevar)

        return sorted(list(set(variants))) # Sorted by reference pos order 
        

def get_sequence_context(fa_stream, variant):
    """
    Return the sequence surrounding this variant's position.
    """

    start = max(0, variant.POS - 10)
    return fa_stream.fetch(variant.CHROM, start, variant.POS + 10)

def homoRunForOneVariant(fa_stream, variant):

    """
    Calculate and return the length of the largest homopolymer
{'fraction_non_alphanumeric': 0.09227311785312975, 'fraction_numerical': 0.0038336968735541016, 'mean_word_length': 1.2619225594259231, 'pattern_counts': {'":': 0, '<': 5, '<?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': 14, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3909215', 'n_tokens_mistral': 4069, 'n_tokens_neox': 3691, 'n_words': 1353}