<commit_before>"""
File: intervals.py
Author: <NAME>
Purpose: Treat the genome as a series of intervals, some of which have a known subspecific origin.
        Over subsets of all samples, find the frequency of all pairwise combinations of origins for intervals.
        E.g. find the number of samples that are domesticus at interval A and musculus at interval B
"""

import numpy as np
import csv
import string

# integer representations of subspecies
NO_SPECIES = 0
SPECIES_TO_INT = {'None': NO_SPECIES, 'Dom': 0o001, 'Mus': 0o010, 'Cast': 0o100}
MAX_POSITION = 200000000  # 200m, bigger than the largest mouse chromosome
CHROMO_TO_INT = {str(x): x for x in xrange(20)}
CHROMO_TO_INT['X'] = 21
CHROMO_TO_INT['Y'] = 22
CHROMO_TO_INT['MT'] = 23
CHROMO_SIZES = [200000000] * len(CHROMO_TO_INT)  # upper bound on size of each chromosome


def genome_index(chromosome, position):
    """ Converts chromosome and position to a single number
    :param chromosome: integer representation of chromosome
    :param position: position on chromosome
    :return: integer denoting chromosome and position
    """
    return ((chromosome - 1) * MAX_POSITION) + position


def make_source_matrix(strain_name, chromosomes):
    """ Constructs and saves the pairwise combinations of subspecific origins of intervals
    :param strain_name: string, name of the strain/sample
    :param chromosomes: dictionary with the intervals for each chromosome
    """
    num_intervals = sum([len(ints) for ints in chromosomes.itervalues()])
    intervals = np.empty(num_intervals, dtype=np.uint32)
    source_array = np.empty(num_intervals, dtype=np.uint8)
    sources = np.empty([num_intervals, num_intervals], dtype=np.uint8)
    interval_num = 0
    for chromosome, interval_list in sorted(chromosomes.iteritems(), key=lambda x: x[0]):
        for species, end in interval_list:
            intervals[interval_num] = genome_index(chromosome, end)
            source_array[interval_num] = species
            interval_num += 1
    for interval_num in xrange(num_intervals):
        sources[interval_num, :] = source_array
    for interval_num in xrange(num_intervals):
        sources[:, interval_num] = (source_array + sources[:, interval_num])
    np.save(strain_name + '_intervals.npy', intervals)
    np.save(strain_name + '_sources.npy', sources)


def rgb_to_subspecies(file_list, suffix='_ss_names'):
    """ replaces color values with subspecies names
    :param file_list: list of files to convert
    :param suffix: the converted files are saved as oldfilename_suffix.extension
    :return: list of filenames of converted files
    """
    new_filenames = []
    for file_name in file_list:
        dot_pos = file_name.find('.')
        with open(file_name) as in_file, open(file_name[:dot_pos] + suffix + file_name[dot_pos:], 'w+') as out_file:
            for row in in_file:
                row = row.replace('255 0 0', 'Mus')
                row = row.replace('0 255 0', 'Cast')
                row = row.replace('0 0 255', 'Dom')
                out_file.write(row)
    return new_filenames


def parse_csvs(file_list):
    """ parses downloaded haplotypes from the Mouse Phylogeny Viewer
    :param file_list: list of filenames
    :return: dictionary of {strain name: {chromosome number: list of (subspecies id, interval end)}}
    """
    # csv headers
    STRAIN = 'strain'
    CHROMOSOME = 'chromosome'
    START = 'start_position'
    END = 'end_position'
    SUBSPECIES = 'subspecies'
    strains = {}
    for filename in file_list:
        with open(filename) as csvfile:
            reader = csv.DictReader(csvfile)
            for row in reader:
                chromosomes = strains.setdefault(row[STRAIN], {})
                intervals = chromosomes.setdefault(CHROMO_TO_INT[row[CHROMOSOME]], [])
                lastEnd = 0 if not intervals else intervals[-1]
                # add null interval if there is a gap between intervals with assigned subspecies
                if not int(row[START]) - 1 <= lastEnd <= int(row[START]) + 1:
                    intervals.append((NO_SPECIES, int(row[START])))
                intervals.append((SPECIES_TO_INT[row[SUBSPECIES]], int(row[END])))
    # add null interval to end of each chromosome
    for chromosomes in strains.itervalues():
        for chromosome, intervals in chromosomes.iteritems():
            intervals.append((NO_SPECIES, CHROMO_SIZES[chromosome]))
    return strains


def preprocess(file_list):
    """ creates the preprocessed matrix for each sample. Run once
    :param file_list: list of file names of csv files with subspecific origins
    :return: list of the strain names that were processed
    """
    strain_names = []
    # remove chars from strain names that can't be in filenames (e.g. /)
    valid_chars = "-_.() %s%s" % (string.ascii_letters, string.digits)
    for strain_name, chromosomes in parse_csvs(file_list).iteritems():
        strain_name = ''.join(c for c in strain_name if c in valid_chars)
        make_source_matrix(strain_name, chromosomes)
        strain_names.append(strain_name)
    return strain_names


def source_frequencies(strain_names):
    """ creates matrices of the pairwise subspecific origin frequencies
    Andrew: this is the only function you should need to edit
    :param strain_names: list of strain names to analyze (must be a subset of the output from preprocess())
    """
    interval_lists = []
    for strain_name in strain_names:
        interval_lists.append(list(np.load(strain_name + '_intervals.npy')))
    # compute elementary intervals
    elem_intervals = []
    while interval_lists:
        elem_intervals.append(min(intervals[0] for intervals in interval_lists))
        for i, intervals in enumerate(interval_lists):
            if intervals[0] == elem_intervals[-1]:
                intervals.pop(0)
            if not intervals:
                interval_lists.pop(i)
    # dictionary of matrices, one for each combination of two sources. i.e. 1 matrix has the counts for the dom-mus combo
    # elementary intervals are along the axes
    # Each element contains the number of samples that have that matrix's combo at that pair of intervals
    source_counts = {}
    for source1 in SPECIES_TO_INT.itervalues():
        if source1 != NO_SPECIES:  # don't count combinations for which one source is null
            for source2 in SPECIES_TO_INT.itervalues():
                if source2 != NO_SPECIES:
                    source_counts[source1 + source2] = np.zeros([len(elem_intervals), len(elem_intervals)],
                                                                dtype=np.uint16)
    # fill in matrices
    for strain_name in strain_names:
        intervals = np.load(strain_name + '_intervals.npy')
        sources = np.load(strain_name + '_sources.npy')
        elem_row = 0
        for row, row_end in enumerate(intervals):
            while elem_row < len(elem_intervals) and elem_intervals[elem_row] <= row_end:
                elem_col = 0
                for source, col_end in zip(sources[row], intervals):
                    while elem_col < len(elem_intervals) and elem_intervals[elem_col] <= col_end:
                        if source in source_counts:
                            source_counts[source][elem_row, elem_col] += 1
                        elem_col += 1
                elem_row += 1
    print sum([source for source in source_counts.itervalues()])
    print elem_intervals
    print source_counts


def main():
    # you should only need to run these once
    converted_file_names = rgb_to_subspecies(['put', 'file', 'names', 'here'])
    file_names = converted_file_names + ['put files', 'you do not need to convert', 'here']
    strain_names = preprocess(file_names)
    print strain_names
    # the following you'll need to run multiple times
    # for any subset you want of strain names, run something like:
    source_frequencies(strain_names)


if __name__ == '__main__':
    main()<commit_msg>use real chromosome lengths<commit_after>"""
File: intervals.py
Authors: <NAME>, <NAME>
Purpose: Treat the genome as a series of intervals, some of which have a known subspecific origin.
        Over subsets of all samples, find the frequency of all pairwise combinations of origins for intervals.
        E.g. find the number of samples that are domesticus at interval A and musculus at interval B
"""

import os
import sys
import numpy as np
import csv
import string
import glob

from time import clock
from collections import OrderedDict

# integer representations of subspecies
NO_SPECIES = 0
SPECIES_TO_INT = OrderedDict({'None': NO_SPECIES, 'dom': 0o001, 'mus': 0o010, 'cas': 0o100})
INTRA_SUBSPECIFIC = [ v+v for k,v in SPECIES_TO_INT.iteritems() ] + [ -1*(v+v) for k,v in SPECIES_TO_INT.iteritems() ]
MAX_POSITION = 200000000  # 200m, bigger than the largest mouse chromosome
CHROMO_TO_INT = OrderedDict()
for x in range(1,21):
    CHROMO_TO_INT[ str(x) ] = x
#CHROMO_TO_INT['X'] = 20
CHROMO_TO_INT['Y'] = 21
CHROMO_TO_INT['MT'] = 22
CHROMO_SIZES = [    197195432, 181748087, 159599783, 155630120, 152537259,
                    149517037, 152524553, 131738871, 124076172, 129993255,
                    121843856, 121257530, 120284312, 125194864, 103494974,
                    98319150,   95272651,  90772031,  61342430, 166650296,
                    91744698,      16299 ] # these are mm9 cooords
CHROMO_OFFSETS = np.cumsum([0] + CHROMO_SIZES)

def genome_index(chromosome, position):
    """ Converts chromosome and position to a single number
    :param chromosome: integer representation of chromosome
    :param position: position on chromosome
    :return: integer denoting chromosome and position
    """
    return CHROMO_OFFSETS[chromosome-1] + position


def make_source_matrix(strain_name, chromosomes):
    """ Constructs and saves the pairwise combinations of subspecific origins of intervals
    :param strain_name: string, name of the strain/sample
    :param chromosomes: dictionary with the intervals for each chromosome
    """
    num_intervals = sum([len(ints) for ints in chromosomes.itervalues()])
    intervals = np.empty(num_intervals, dtype=np.uint32)
    source_array = np.empty(num_intervals, dtype=np.int8)
    sources = np.empty([num_intervals, num_intervals], dtype=np.int8)
    interval_num = 0
    for chromosome, interval_list in sorted(chromosomes.iteritems(), key=lambda x: x[0]):
        for species, end in interval_list:
            intervals[interval_num] = genome_index(chromosome, end)
            source_array[interval_num] = species
            interval_num += 1
    for interval_num in xrange(num_intervals):
        sources[interval_num, :] = source_array
    for interval_num in xrange(num_intervals):
        sources[:, interval_num] = (source_array + sources[:, interval_num])
    ## prox::dist is upper triangle; dist::prox is lower triangle
    sources_new = np.empty([num_intervals, num_intervals], dtype=np.int8)
    sources_new = np.triu(sources) + -1*np.tril(sources, k = -1)
    np.save(os.path.join("intervals", strain_name + '_intervals.npy'), intervals)
    np.save(os.path.join("intervals", strain_name + '_sources.npy'), sources_new)


def rgb_to_subspecies(file_list, suffix='_ss_names'):
    """ replaces color values with subspecies names
    :param file_list: list of files to convert
    :param suffix: the converted files are saved as oldfilename_suffix.extension
    :return: list of filenames of converted files
    """
    new_filenames = []
    for file_name in file_list:
        dot_pos = file_name.find('.')
        this_file = file_name[:dot_pos] + suffix + file_name[dot_pos:]
        new_filenames.append(this_file)
        with open(file_name) as in_file, open(this_file, 'w+') as out_file:
            first = 0
            for row in in_file:
                if not first:
                    out_file.write(row)
                    first += 1
                else:
                    row = row.strip()
                    # skip regions with ambiguous origin (mostly applies to wild-caught mice)
                    if not (row.endswith(",255 0 0") or row.endswith(",0 255 0") or row.endswith(",0 0 255")):
                        continue
                    else:
                        row = row.replace('255 0 0', 'mus')
                        row = row.replace('0 255 0', 'cas')
                        row = row.replace('0 0 255', 'dom')
                        out_file.write(row + "\n")
                        first += 1

    return new_filenames


def parse_csvs(file_list):
    """ parses downloaded haplotypes from the Mouse Phylogeny Viewer
    :param file_list: list of filenames
    :return: dictionary of {strain name: {chromosome number: list of (subspecies id, interval end)}}
    """
    # csv headers
    STRAIN = 'strain'
    CHROMOSOME = 'chrom'
    START = 'start'
    END = 'end'
    SUBSPECIES = 'subspecies'
    strains = {}
    for filename in file_list:
        with open(filename) as csvfile:
            reader = csv.DictReader(csvfile)
            for row in reader:
                chromosomes = strains.setdefault(row[STRAIN], OrderedDict())
                intervals = chromosomes.setdefault(CHROMO_TO_INT[row[CHROMOSOME]], [])
                lastEnd = 0 if not intervals else intervals[-1][1]
                last_species = 0 if not intervals else intervals[-1][0]
                print row
                if last_species != SPECIES_TO_INT[row[SUBSPECIES]]:
                # add null interval if there is a gap between intervals with assigned subspecies
                    if not int(row[START]) - 1 <= lastEnd <= int(row[START]) + 1:
                        intervals.append((NO_SPECIES, int(row[START])))
                    intervals.append((SPECIES_TO_INT[row[SUBSPECIES]], int(row[END])))
                lastEnd = row[END]
                last_species = SPECIES_TO_INT[row[SUBSPECIES]]
    # add null interval to end of each chromosome
    for chromosomes in strains.itervalues():
        for chromosome, intervals in chromosomes.iteritems():
            intervals.append((NO_SPECIES, CHROMO_SIZES[chromosome-1]))
    print strains
    return strains


def preprocess(file_list):
    """ creates the preprocessed matrix for each sample. Run once
    :param file_list: list of file names of csv files with subspecific origins
    :return: list of the strain names that were processed
    """
    strain_names = []
    # remove chars from strain names that can't be in filenames (e.g. /)
    valid_chars = "-_.()%s%s" % (string.ascii_letters, string.digits)
    for strain_name, chromosomes in parse_csvs(file_list).iteritems():
        strain_name = ''.join(c for c in strain_name if c in valid_chars)
        make_source_matrix(strain_name, chromosomes)
        strain_names.append(strain_name)
    return strain_names


def source_frequencies(strain_names, in_path = ""):
    """ creates matrices of the pairwise subspecific origin frequencies
    Andrew: this is the only function you should need to edit
    :param strain_names: list of strain names to analyze (must be a subset of the output from preprocess())
    :param in_path: where to look for pre-computed intervals from preprocess()
    """
    sys.stderr.write("Loading intervals...\n")
    interval_lists = []
    interval_dict = {}
    for strain_name in strain_names:
        interval_lists.append( list(np.load(os.path.join(in_path, strain_name + '_intervals.npy'))) )
        interval_dict[strain_name] = np.copy(interval_lists[-1])

    # compute elementary intervals
    sys.stderr.write("Computing elementary intervals...")
    elem_intervals = []
    while interval_lists:
        elem_intervals.append(min(intervals[0] for intervals in interval_lists))
        for i, intervals in enumerate(interval_lists):
            #print intervals
            if intervals[0] == elem_intervals[-1]:
                intervals.pop(0)
            if not intervals:
                interval_lists.pop(i)
    sys.stderr.write(" {} total intervals\n".format(len(elem_intervals)))

    # dictionary of matrices, one for each combination of two sources. i.e. 1 matrix has the counts for the dom-mus combo
    # elementary intervals are along the axes
    # Each element contains the number of samples that have that matrix's combo at that pair of intervals
    source_counts = {}
    for species1, source1 in SPECIES_TO_INT.iteritems():
        for species2, source2 in SPECIES_TO_INT.iteritems():
            source_counts[ source1+source2 ] = np.zeros([len(elem_intervals), len(elem_intervals)],
                                                        dtype=np.int16)
            source_counts[ -1*(source1+source2) ] = np.zeros([len(elem_intervals), len(elem_intervals)],
                                                        dtype=np.int16)
    ## fill in matrices
    sys.stderr.write("Counting incidence of subspecies pairs...\n")
    for strain_name in strain_names:
        sys.stderr.write("\t-- {}\n".format(strain_name))
        start = clock()
        intervals = interval_dict[strain_name]
        sources = np.load(os.path.join(in_path, strain_name + "_sources.npy"))
        elem_row = 0
        # map this strain's intervals onto the elementary intervals
        breaks = np.searchsorted(elem_intervals, intervals)
        for row, row_end in enumerate(intervals):
            row_ind = range(breaks[row], breaks[row]+1)
            for col, col_end in enumerate(intervals):
                source = sources[row,col]
                col_ind = range(breaks[col], breaks[col]+1)
                if source in source_counts:
                    for i in col_ind:
                        source_counts[source][ row_ind,i ] += 1
        end = clock()
        sys.stderr.write("\t\t... in {} seconds\n".format(end-start))

    return ( source_counts, elem_intervals )

def calc_genomic_area(counts, intervals):
    '''
    Compute the total genomic 'area' occupied by each combination of subspecies.
    :param counts: dictionary of incidence matrices, one per subspecies combo
    :param intervals: the 'elementary intervals' over which the counts were computed
    '''
    # compute area of each cell in the interval grid
    intervals = np.array([0] + intervals, dtype = np.float32)/1.0e6
    areas = np.zeros([ len(intervals)-1, len(intervals)-1 ], dtype = np.float32)
    for row in range(1, len(intervals)):
        for col in range(row, len(intervals)):
            areas[row-1,col-1] = (intervals[row]-intervals[row-1]) * (intervals[col]-intervals[col-1])
            if col > row:
                areas[col-1,row-1] = areas[row-1,col-1]

    print intervals[:10]
    print (areas[ range(0,10),: ])[ :,range(0,10) ]

    areas_masked = OrderedDict()
    # get denominator for genomic area
    genome_len = sum(float(x)/1.0e6 for x in CHROMO_SIZES)-(3.0*20.0)
    denom = genome_len**2
    for combo, vals in counts.iteritems():
        factor = 2.0 if combo in INTRA_SUBSPECIFIC else 1.0
        areas_masked.update({ combo: np.sum((vals > 0)*areas*factor)/denom })
    return areas_masked

def list_available_strains(globber):
    '''
    Given path to a databse of pre-computed subspecies origin intervals,
    list the strains that are available.
    :param globber: file pattern for finding strains, possibly including wildcards
    '''
    sources = []
    for s in glob.glob(globber + "_sources.npy"):
        sources.append( os.path.basename(s).replace("_sources.npy","") )

    intervals = []
    for s in glob.glob(globber + "_intervals.npy"):
        intervals.append( os.path.basename(s).replace("_intervals.npy","") )

    avail = set(intervals) & set(sources)
    return list(avail)

def main():

    # processed_files = rgb_to_subspecies(["subspecific_origins.csv"])
    # strain_names = preprocess(processed_files)
    # print strain_names
    # strain_names = preprocess(["cc.subspecies.csv"])
    #strain_names = preprocess(["test.csv"])
    #print strain_names
    #sys.exit(1)

    # the following you'll need to run multiple times
    # for any subset you want of strain names, run something like:
    strain_names = list_available_strains("intervals/STU*")[0:4]

    sys.stderr.write("Computing pairwise subspecies origins for these strains:\n")
    sys.stderr.write("\t-- {}\n".format("\n\t-- ".join(strain_names)))

    pairs = source_frequencies(strain_names, in_path = "intervals/")
    combos = OrderedDict()
    for s1, i1 in SPECIES_TO_INT.iteritems():
        for s2, i2 in SPECIES_TO_INT.iteritems():
            combos[ s1 + "_" + s2 ] = i1+i2
            combos[ s2 + "_" + s1 ] = -1*(i1+i2)

    sys.stderr.write("\nSumming up...\n")
    areas = calc_genomic_area(pairs[0], pairs[1])
    for combo, code in combos.iteritems():
        if code in areas:
            print "\t{} ({}): {: 1.5f}".format(combo, code, areas[code])

if __name__ == '__main__':
    main()
{'fraction_non_alphanumeric': 0.09931047692013024, 'fraction_numerical': 0.0215475962459299, 'mean_word_length': 2.835629017447199, 'pattern_counts': {'":': 0, '<': 14, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2519169', 'n_tokens_mistral': 6705, 'n_tokens_neox': 6125, 'n_words': 1834}