## ## MCMC sampler for Mixture-of-Isoforms (MISO) model ## ## Yarden Katz ## ## The sampler uses a Metropolis-Hastings sampling scheme, combined with ## a Gibbs sampling step. ## import scipy import misopy from misopy.reads_utils import count_aligned_reads, \ count_isoform_assignments from misopy.read_simulator import simulate_reads, print_reads_summary, \ read_counts_to_read_list, \ get_reads_summary import misopy.hypothesis_test as ht from misopy.Gene import Gene, Exon from misopy.py2c_gene import * # C MISO interface import pysplicing from scipy import * from numpy import * import cPickle as pickle from scipy.stats import mode import math import time from numpy import numarray import os import sys from collections import defaultdict import glob import logging import logging.handlers loggers = {} def get_logger(logger_name, log_outdir, level=logging.WARNING, include_stdout=True): """ Return a logging object. """ global loggers # Avoid race-conditions try: os.makedirs(log_outdir) except OSError: pass if loggers.get(logger_name): return loggers.get(logger_name) logger = logging.getLogger(logger_name) formatter = \ logging.Formatter('%(asctime)s - %(name)s - %(levelname)s - %(message)s', datefmt='%m/%d/%Y %I:%M:%S %p') # Do not log to file #if log_outdir is not None: # log_filename = os.path.join(log_outdir, "%s.log" %(logger_name)) # fh = logging.FileHandler(log_filename) # fh.setLevel(level) # fh.setFormatter(formatter) # logger.addHandler(fh) logging.root.setLevel(level) # Optionally add handler that streams all logs # to stdout if include_stdout: ch = logging.StreamHandler(sys.stdout) ch.setLevel(level) ch.setFormatter(formatter) logger.addHandler(ch) logger.info("Created logger %s" %(logger_name)) loggers.update({logger_name: logger}) return logger ## ## Helper statistics/linear algebra functions ## def set_diag(a, v): for i, elt in enumerate(a): a[i, i] = v return a def maxi(l): m = max(l) for i, v in enumerate(l): if m == v: return i def mini(l): m = min(l) for i, v in enumerate(l): if m == v: return i def exp_logsumexp(a): return exp(a - logsumexp(a)) def vect_logsumexp(a, axis=None): if axis is None: # Use the scipy.maxentropy version. return logsumexp(a) a = asarray(a) shp = list(a.shape) shp[axis] = 1 a_max = a.max(axis=axis) s = log(exp(a - a_max.reshape(shp)).sum(axis=axis)) lse = a_max + s return lse def print_assignment_summary(assignments): counts = defaultdict(int) for a in assignments: counts[a] += 1 for k, v in counts.iteritems(): print "Total of %d in isoform %d" %(v, k) def float_array_to_str(array_of_floats): """ Convert a float numpy array to a string for printing purposes. """ str_float_array = '[' + ' '.join(['%.3f' %(val) for val in array_of_floats]) + ']' return str_float_array def get_paired_end_sampler_params(num_isoforms, mean_frag_len, frag_variance, read_len, overhang_len=1): """ Return parameters for MISO sampler, in paired-end mode. """ hyperparameters = ones(num_isoforms) proposal_diag = 0.05 sigma = set_diag(zeros([num_isoforms-1, num_isoforms-1]), proposal_diag) sampler_params = {'read_len': read_len, 'overhang_len': overhang_len, 'uniform_proposal': False, 'sigma_proposal': sigma, 'mean_frag_len': mean_frag_len, 'frag_variance': frag_variance} return sampler_params def get_single_end_sampler_params(num_isoforms, read_len, overhang_len=1): """ Return parameters for MISO sampler, in single-end mode. """ hyperparameters = ones(num_isoforms) proposal_diag = 0.05 sigma = set_diag(zeros([num_isoforms-1, num_isoforms-1]), proposal_diag) sampler_params = {'read_len': read_len, 'overhang_len': overhang_len, 'uniform_proposal': False, 'sigma_proposal': sigma} return sampler_params class MISOSampler: def __init__(self, params, paired_end=False, log_dir=None): """ Make a sampler with the given parameters. """ self.params = params self.paired_end = paired_end # set default fragment length distribution parameters if self.paired_end: if ((not 'mean_frag_len' in self.params) or \ (not 'frag_variance' in self.params)): raise Exception, "Must set mean_frag_len and frag_variance when " \ "running in sampler on paired-end data." self.mean_frag_len = self.params['mean_frag_len'] self.frag_variance = self.params['frag_variance'] if log_dir != None: self.log_dir = os.path.abspath(os.path.expanduser(log_dir)) self.log_dir = os.path.join(log_dir, 'logs') # Avoid race-conditions try: os.makedirs(self.log_dir) except OSError: pass self.miso_logger = get_logger('miso_logger', self.log_dir) self.miso_logger.info("Instantiated sampler.") def run_sampler(self, num_iters, reads, gene, hyperparameters, params, output_file, num_chains=6, burn_in=1000, lag=2, prior_params=None, # By default, use sampler with read classes (collapsed) # to get speed boost for single-end reads # (To revert to old reassigning sampler, use # pysplicing.MISO_ALGO_REASSIGN) algorithm=pysplicing.MISO_ALGO_CLASSES, start_cond=pysplicing.MISO_START_AUTO, stop_cond=pysplicing.MISO_STOP_FIXEDNO, verbose=True): """ Fast version of MISO MCMC sampler. Calls C version and returns results. """ num_isoforms = len(gene.isoforms) self.num_isoforms = num_isoforms if prior_params == None: prior_params = (1.0,) * num_isoforms read_positions = reads[0] read_cigars = reads[1] self.num_reads = len(read_positions) if self.num_reads == 0: print "No reads for gene: %s" %(gene.label) return output_file = output_file + ".miso" # If output filename exists, don't run sampler if os.path.isfile(os.path.normpath(output_file)): print "Output filename %s exists, not running MISO." \ %(output_file) return None self.params['iters'] = num_iters self.params['burn_in'] = burn_in self.params['lag'] = lag # Define local variables related to reads and overhang self.overhang_len = self.params['overhang_len'] self.read_len = self.params['read_len'] t1 = 0 t2 = 0 if verbose: t1 = time.time() #self.miso_logger.info("Running sampler...") #self.miso_logger.info(" - num_iters: " + str(num_iters)) #self.miso_logger.info(" - burn-in: " + str(burn_in)) #self.miso_logger.info(" - lag: " + str(lag)) #self.miso_logger.info(" - paired-end? " + str(self.paired_end)) #self.miso_logger.info(" - gene: " + str(gene)) rejected_proposals = 0 accepted_proposals = 0 psi_vectors = [] all_psi_proposals = [] if params['uniform_proposal']: self.miso_logger.debug("UNIFORM independent proposal being used.") proposal_type = "unif" else: self.miso_logger.debug("Non-uniform proposal being used.") self.miso_logger.debug(" - sigma_proposal: " + str(params['sigma_proposal'])) proposal_type = "drift" init_psi = ones(num_isoforms)/float(num_isoforms) # Do not process genes with one isoform if num_isoforms == 1: one_iso_msg = "Gene %s has only one isoform; skipping..." \ %(gene.label) self.miso_logger.warning(one_iso_msg) return # Convert Python Gene object to C c_gene = py2c_gene(gene) ## ## Run C MISO ## read_positions = tuple([r+1 for r in read_positions]) if self.paired_end: # Number of standard deviations in insert length # distribution to consider when assigning reads # to isoforms num_sds = 4L # Run paired-end miso_results = pysplicing.MISOPaired(c_gene, 0L, read_positions, read_cigars, long(self.read_len), float(self.mean_frag_len), float(self.frag_variance), float(num_sds), long(num_iters), long(burn_in), long(lag), prior_params, long(self.overhang_len), long(num_chains), start_cond, stop_cond) else: # Run single-end miso_results = pysplicing.MISO(c_gene, 0L, read_positions, read_cigars, long(self.read_len), long(num_iters), long(burn_in), long(lag), prior_params, long(self.overhang_len), long(num_chains), start_cond, stop_cond, pysplicing.MISO_ALGO_REASSIGN) # algorithm) # Psi samples psi_vectors = transpose(array(miso_results[0])) # Log scores of accepted samples kept_log_scores = transpose(array(miso_results[1])) # Read classes read_classes = miso_results[2] # Read class statistics read_class_data = miso_results[3] # Assignments of reads to isoforms assignments = miso_results[4] # Statistics and parameters about sampler run run_stats = miso_results[5] # Assignments of reads to classes. # read_classes[n] represents the read class that has # read_assignments[n]-many reads. reads_data = (read_classes, read_class_data) assignments = array(assignments) # Skip events where all reads are incompatible with the annotation; # do not output a file for those. if all(assignments == -1): print "All reads incompatible with annotation, skipping..." return accepted_proposals = run_stats[4] rejected_proposals = run_stats[5] percent_acceptance = (float(accepted_proposals)/(accepted_proposals + \ rejected_proposals)) * 100 #self.miso_logger.info("Percent acceptance (including burn-in): %.4f" %(percent_acceptance)) #self.miso_logger.info("Number of iterations recorded: %d" %(len(psi_vectors))) # Write MISO output to file print "Outputting samples to: %s..." %(output_file) self.miso_logger.info("Outputting samples to: %s" %(output_file)) self.output_miso_results(output_file, gene, reads_data, assignments, psi_vectors, kept_log_scores, num_iters, burn_in, lag, percent_acceptance, proposal_type) if verbose: t2 = time.time() print "Event took %.2f seconds" %(t2 - t1) def output_miso_results(self, output_file, gene, reads_data, assignments, psi_vectors, kept_log_scores, num_iters, burn_in, lag, percent_acceptance, proposal_type): """ Output results of MISO to a file. """ output = open(output_file, 'w') # Get a string representation of the isoforms - use '_' # in the delimiter regardless iso_delim = '_' if type(gene.isoforms[0].desc) == list: str_isoforms = '[' + ",".join(["\'" + iso_delim.join(iso.desc) + "\'" \ for iso in gene.isoforms]) + ']' else: str_isoforms = '[' + ",".join(["\'" + iso.desc + "\'" \ for iso in gene.isoforms]) + ']' num_isoforms = len(gene.isoforms) # And of the exon lengths exon_lens = ",".join(["(\'%s\',%d)" %(p.label, p.len) \ for p in gene.parts]) ## Compile header with information about isoforms and internal parameters used ## by the sampler, and also information about read counts and number of ## reads assigned to each isoform. read_classes, read_class_counts = reads_data read_counts_list = [] for class_num, class_type in enumerate(read_classes): class_counts = read_class_counts[class_num] # Get the read class type in string format class_str = str(tuple([int(c) for c in class_type])).replace(" ", "") # Get the read class counts in string format class_counts_str = "%s" %(int(read_class_counts[class_num])) # Put class and counts together curr_str = "%s:%s" %(class_str, class_counts_str) read_counts_list.append(curr_str) # Get a summary of the raw read counts supporting each isoform read_counts_str = ",".join(read_counts_list) assigned_counts = count_isoform_assignments(assignments) # Get number of reads assigned to each isoform assigned_counts_str = ",".join(["%d:%d" %(c[0], c[1]) \ for c in assigned_counts]) # coordinates where mRNAs start mRNA_starts = [] mRNA_ends = [] for iso in gene.isoforms: mRNA_starts.append(iso.genomic_start) mRNA_ends.append(iso.genomic_end) mRNA_start_coords = ",".join([str(start) for start in mRNA_starts]) mRNA_end_coords = ",".join([str(end) for end in mRNA_ends]) chrom = gene.chrom if chrom == None: chrom = "NA" strand = gene.strand if strand == None: strand = "NA" header = "#isoforms=%s\texon_lens=%s\titers=%d\tburn_in=%d\tlag=%d\t" \ "percent_accept=%.2f\tproposal_type=%s\t" \ "counts=%s\tassigned_counts=%s\tchrom=%s\tstrand=%s\tmRNA_starts=%s\tmRNA_ends=%s\n" \ %(str_isoforms, exon_lens, num_iters, burn_in, lag, percent_acceptance, proposal_type, read_counts_str, assigned_counts_str, # Fields related to gene/event chrom, strand, mRNA_start_coords, mRNA_end_coords) output.write(header) # Output samples and their associated log scores, as well as read counts results_fields = ["sampled_psi", "log_score"] results_header = "%s\n" %("\t".join(results_fields)) output.write(results_header) for psi_sample, curr_log_score in zip(psi_vectors, kept_log_scores): psi_sample_str = ",".join(["%.4f" %(psi) for psi in psi_sample]) output_line = "%s\t%.2f\n" %(psi_sample_str, curr_log_score) output.write(output_line) output.close() print "Completed outputting." # return [percent_acceptance, array(psi_vectors), array(kept_log_scores)] def run_sampler_on_event(gene, ni, ne, nb, read_len, overhang_len, num_iters, output_dir, confidence_level=.95): """ Run sampler on a two-isoform gene event. """ print "Running sampler on a two-isoform event..." print " - Gene label: ", gene.label, gene print " - NI, NE, NB: %d, %d, %d" %(ni, ne, nb) print "Using default sampler parameters." if gene.chrom != None: # Index output by chromosome print "Indexing by chromosome..." output_dir = os.path.join(output_dir, gene.chrom) if not os.path.exists(output_dir): os.makedirs(output_dir) output_filename = os.path.join(output_dir, gene.label) samples = [] cred_interval = [] num_isoforms = len(gene.isoforms) burn_in = 500 lag = 10 hyperparameters = ones(num_isoforms) proposal_diag = 0.05 sigma = set_diag(zeros([num_isoforms-1, num_isoforms-1]), proposal_diag) sampler_params = {'read_len': read_len, 'overhang_len': overhang_len, 'uniform_proposal': False, 'sigma_proposal': sigma} sampler = MISOSampler(sampler_params, log_dir=output_dir) reads = read_counts_to_read_list(ni, ne, nb) t1 = time.time() sampler_results = sampler.run_sampler(num_iters, reads, gene, hyperparameters, sampler_params, output_filename, burn_in=burn_in, lag=lag) if not sampler_results: return (samples, cred_interval) samples = sampler_results[1] # Compute credible intervals cred_interval = ht.compute_credible_intervals(samples, confidence_level=confidence_level) t2 = time.time() print " - Sampler run took %s seconds." %(str(t2-t1)) # return samples and credible intervals return (samples, cred_interval) def profile_miso(): from Gene import make_gene gene = make_gene([150, 100, 150], [[1, 2, 3], [1, 3]]) read_len = 36 overhang_len = 4 output_dir = "profiler-test" for x in range(10): print "x = %d" %(x) a, b = run_sampler_on_event(gene, 500, 50, 40, read_len, overhang_len, 10000, output_dir) def main(): return # import cProfile as profile # import pstats # output_file = "profile" # profile.run('profile_miso()', output_file) # p = pstats.Stats(output_file) # print "name: " # print p.sort_stats('name') # print "all stats: " # p.print_stats() # print "cumulative (top 10): " # p.sort_stats('cumulative').print_stats(20) if __name__ == '__main__': main()