genobstacles command parser
495826
def cmd_genobstacles(self, args):
        ''''''
        usage = "usage: genobstacles <start|stop|restart|clearall|status|set>"
        if len(args) == 0:
            print(usage)
            return
        if args[0] == "set":
            gen_settings.command(args[1:])
        elif args[0] == "start":
            if self.have_home:
                self.start()
            else:
                self.pending_start = True
        elif args[0] == "stop":
            self.stop()
            self.pending_start = False
        elif args[0] == "restart":
            self.stop()
            self.start()
        elif args[0] == "status":
            print(self.status())
        elif args[0] == "remove":
            latlon = self.module('map').click_position
            if self.last_click is not None and self.last_click == latlon:
                return
            self.last_click = latlon
            if latlon is not None:
                closest = None
                closest_distance = 1000
                for a in self.aircraft:
                    dist = a.distance_from(latlon[0], latlon[1])
                    if dist < closest_distance:
                        closest_distance = dist
                        closest = a
                if closest is not None:
                    self.aircraft.remove(closest)
                else:
                    print("No obstacle found at click point")
                    
        elif args[0] == "dropcloud":
            self.cmd_dropobject(Weather())
        elif args[0] == "dropeagle":
            self.cmd_dropobject(BirdOfPrey())
        elif args[0] == "dropbird":
            self.cmd_dropobject(BirdMigrating())
        elif args[0] == "dropplane":
            self.cmd_dropobject(Aircraft())
        elif args[0] == "clearall":
            self.clearall()
        else:
            print(usage)
708835
def main():
    """"""

    parser = OptionParser(conflict_handler="resolve")

    parser.add_option('--humanTRA', '--human_T_alpha', action='store_true', dest='humanTRA', default=False, help='use default human TRA model (T cell alpha chain)')
    parser.add_option('--humanTRB', '--human_T_beta', action='store_true', dest='humanTRB', default=False, help='use default human TRB model (T cell beta chain)')
    parser.add_option('--mouseTRB', '--mouse_T_beta', action='store_true', dest='mouseTRB', default=False, help='use default mouse TRB model (T cell beta chain)')
    parser.add_option('--humanIGH', '--human_B_heavy', action='store_true', dest='humanIGH', default=False, help='use default human IGH model (B cell heavy chain)')
    parser.add_option('--set_custom_model_VDJ', dest='vdj_model_folder', metavar='PATH/TO/FOLDER/', help='specify PATH/TO/FOLDER/ for a custom VDJ generative model')
    parser.add_option('--set_custom_model_VJ', dest='vj_model_folder', metavar='PATH/TO/FOLDER/', help='specify PATH/TO/FOLDER/ for a custom VJ generative model')

    parser.add_option('-i', '--infile', dest = 'infile_name',metavar='PATH/TO/FILE', help='read in CDR3 sequences (and optionally V/J masks) from PATH/TO/FILE')
    parser.add_option('-o', '--outfile', dest = 'outfile_name', metavar='PATH/TO/FILE', help='write CDR3 sequences and pgens to PATH/TO/FILE')
    parser.add_option('--seq_in', '--seq_index', type='int', metavar='INDEX', dest='seq_in_index', default = 0, help='specifies sequences to be read in are in column INDEX. Default is index 0 (the first column).')

    parser.add_option('--v_in', '--v_mask_index', type='int', metavar='INDEX', dest='V_mask_index', help='specifies V_masks are found in column INDEX in the input file. Default is no V mask.')
    parser.add_option('--j_in', '--j_mask_index', type='int', metavar='INDEX', dest='J_mask_index', help='specifies J_masks are found in column INDEX in the input file. Default is no J mask.')

    parser.add_option('--v_mask', type='string', dest='V_mask', help='specify V usage to condition Pgen on for seqs read in as arguments.')
    parser.add_option('--j_mask', type='string', dest='J_mask', help='specify J usage to condition Pgen on for seqs read in as arguments.')

    parser.add_option('-m', '--max_number_of_seqs', type='int',metavar='N', dest='max_number_of_seqs', help='compute Pgens for at most N sequences.')
    parser.add_option('--lines_to_skip', type='int',metavar='N', dest='lines_to_skip', default = 0, help='skip the first N lines of the file. Default is 0.')
    parser.add_option('-a', '--alphabet_filename', dest='alphabet_filename', metavar='PATH/TO/FILE', help="specify PATH/TO/FILE defining a custom 'amino acid' alphabet. Default is no custom alphabet.")
    parser.add_option('--seq_type_out', type='choice',metavar='SEQ_TYPE', dest='seq_type_out',  choices=['all', 'ntseq', 'nucleotide', 'aaseq', 'amino_acid'], help="if read in sequences are ntseqs, declare what type of sequence to compute pgen for. Default is all. Choices: 'all', 'ntseq', 'nucleotide', 'aaseq', 'amino_acid'")
    parser.add_option('--skip_off','--skip_empty_off', action='store_true', dest = 'skip_empty', default=True, help='stop skipping empty or blank sequences/lines (if for example you want to keep line index fidelity between the infile and outfile).')

    parser.add_option('--display_off', action='store_false', dest='display_seqs', default=True, help='turn the sequence display off (only applies in write-to-file mode). Default is on.')
    parser.add_option('--num_lines_for_display', type='int', metavar='N', default = 50, dest='num_lines_for_display', help='N lines of the output file are displayed when sequence display is on. Also used to determine the number of sequences to average over for speed and time estimates.')
    parser.add_option('--time_updates_off', action='store_false', dest='time_updates', default=True, help='turn time updates off (only applies when sequence display is disabled).')
    parser.add_option('--seqs_per_time_update', type='float', metavar='N', default = 100, dest='seqs_per_time_update', help='specify the number of sequences between time updates. Default is 1e5.')

    parser.add_option('-d', '--delimiter', type='choice', dest='delimiter',  choices=['tab', 'space', ',', ';', ':'], help="declare infile delimiter. Default is tab for .tsv input files, comma for .csv files, and any whitespace for all others. Choices: 'tab', 'space', ',', ';', ':'")
    parser.add_option('--raw_delimiter', type='str', dest='delimiter', help="declare infile delimiter as a raw string.")
    parser.add_option('--delimiter_out', type='choice', dest='delimiter_out',  choices=['tab', 'space', ',', ';', ':'], help="declare outfile delimiter. Default is tab for .tsv output files, comma for .csv files, and the infile delimiter for all others. Choices: 'tab', 'space', ',', ';', ':'")
    parser.add_option('--raw_delimiter_out', type='str', dest='delimiter_out', help="declare for the delimiter outfile as a raw string.")
    parser.add_option('--gene_mask_delimiter', type='choice', dest='gene_mask_delimiter',  choices=['tab', 'space', ',', ';', ':'], help="declare gene mask delimiter. Default comma unless infile delimiter is comma, then default is a semicolon. Choices: 'tab', 'space', ',', ';', ':'")
    parser.add_option('--raw_gene_mask_delimiter', type='str', dest='gene_mask_delimiter', help="declare delimiter of gene masks as a raw string.")
    parser.add_option('--comment_delimiter', type='str', dest='comment_delimiter', help="character or string to indicate comment or header lines to skip.")


    (options, args) = parser.parse_args()

    #Check that the model is specified properly
    main_folder = os.path.dirname(__file__)

    default_models = {}
    default_models['humanTRA'] = [os.path.join(main_folder, 'default_models', 'human_T_alpha'),  'VJ']
    default_models['humanTRB'] = [os.path.join(main_folder, 'default_models', 'human_T_beta'), 'VDJ']
    default_models['mouseTRB'] = [os.path.join(main_folder, 'default_models', 'mouse_T_beta'), 'VDJ']
    default_models['humanIGH'] = [os.path.join(main_folder, 'default_models', 'human_B_heavy'), 'VDJ']

    num_models_specified = sum([1 for x in default_models.keys() + ['vj_model_folder', 'vdj_model_folder'] if getattr(options, x)])

    if num_models_specified == 1: #exactly one model specified
        try:
            d_model = [x for x in default_models.keys() if getattr(options, x)][0]
            model_folder = default_models[d_model][0]
            recomb_type = default_models[d_model][1]
        except IndexError:
            if options.vdj_model_folder: #custom VDJ model specified
                model_folder = options.vdj_model_folder
                recomb_type = 'VDJ'
            elif options.vj_model_folder: #custom VJ model specified
                model_folder = options.vj_model_folder
                recomb_type = 'VJ'
    elif num_models_specified == 0:
        print 'Need to indicate generative model.'
        print 'Exiting...'
        return -1
    elif num_models_specified > 1:
        print 'Only specify one model'
        print 'Exiting...'
        return -1

    #Check that all model and genomic files exist in the indicated model folder
    if not os.path.isdir(model_folder):
        print 'Check pathing... cannot find the model folder: ' + model_folder
        print 'Exiting...'
        return -1

    params_file_name = os.path.join(model_folder,'model_params.txt')
    marginals_file_name = os.path.join(model_folder,'model_marginals.txt')
    V_anchor_pos_file = os.path.join(model_folder,'V_gene_CDR3_anchors.csv')
    J_anchor_pos_file = os.path.join(model_folder,'J_gene_CDR3_anchors.csv')

    for x in [params_file_name, marginals_file_name, V_anchor_pos_file, J_anchor_pos_file]:
        if not os.path.isfile(x):
            print 'Cannot find: ' + x
            print 'Please check the files (and naming conventions) in the model folder ' + model_folder
            print 'Exiting...'
            return -1

    alphabet_filename = options.alphabet_filename #used if a custom alphabet is to be specified
    if alphabet_filename is not None:
        if not os.path.isfile(alphabet_filename):
            print 'Cannot find custom alphabet file: ' + infile_name
            print 'Exiting...'
            return -1

    #Load up model based on recomb_type
    #VDJ recomb case --- used for TCRB and IGH
    if recomb_type == 'VDJ':
        genomic_data = load_model.GenomicDataVDJ()
        genomic_data.load_igor_genomic_data(params_file_name, V_anchor_pos_file, J_anchor_pos_file)
        generative_model = load_model.GenerativeModelVDJ()
        generative_model.load_and_process_igor_model(marginals_file_name)
        pgen_model = generation_probability.GenerationProbabilityVDJ(generative_model, genomic_data, alphabet_filename)
    #VJ recomb case --- used for TCRA and light chain
    elif recomb_type == 'VJ':
        genomic_data = load_model.GenomicDataVJ()
        genomic_data.load_igor_genomic_data(params_file_name, V_anchor_pos_file, J_anchor_pos_file)
        generative_model = load_model.GenerativeModelVJ()
        generative_model.load_and_process_igor_model(marginals_file_name)
        pgen_model = generation_probability.GenerationProbabilityVJ(generative_model, genomic_data, alphabet_filename)

    aa_alphabet = ''.join(pgen_model.codons_dict.keys())

    if options.infile_name is not None:
        infile_name = options.infile_name

        if not os.path.isfile(infile_name):
            print 'Cannot find input file: ' + infile_name
            print 'Exiting...'
            return -1

    if options.outfile_name is not None:
        outfile_name = options.outfile_name
        if os.path.isfile(outfile_name):
            if not raw_input(outfile_name + ' already exists. Overwrite (y/n)? ').strip().lower() in ['y', 'yes']:
                print 'Exiting...'
                return -1

    #Parse delimiter
    delimiter = options.delimiter
    if delimiter is None: #Default case
        if options.infile_name is None:
            delimiter = '\t'
        elif infile_name.endswith('.tsv'): #parse TAB separated value file
            delimiter = '\t'
        elif infile_name.endswith('.csv'): #parse COMMA separated value file
            delimiter = ','
    else:
        try:
            delimiter = {'tab': '\t', 'space': ' ', ',': ',', ';': ';', ':': ':'}[delimiter]
        except KeyError:
            pass #Other string passed as the delimiter.

    #Parse delimiter_out
    delimiter_out = options.delimiter_out
    if delimiter_out is None: #Default case
        if delimiter is None:
            delimiter_out = '\t'
        else:
            delimiter_out = delimiter
        if options.outfile_name is None:
            pass
        elif outfile_name.endswith('.tsv'): #output TAB separated value file
            delimiter_out = '\t'
        elif outfile_name.endswith('.csv'): #output COMMA separated value file
            delimiter_out = ','
    else:
        try:
            delimiter_out = {'tab': '\t', 'space': ' ', ',': ',', ';': ';', ':': ':'}[delimiter_out]
        except KeyError:
            pass #Other string passed as the delimiter.

    #Parse gene_delimiter
    gene_mask_delimiter = options.gene_mask_delimiter
    if gene_mask_delimiter is None: #Default case
        gene_mask_delimiter = ','
        if delimiter == ',':
            gene_mask_delimiter = ';'
    else:
        try:
            gene_mask_delimiter = {'tab': '\t', 'space': ' ', ',': ',', ';': ';', ':': ':'}[gene_mask_delimiter]
        except KeyError:
            pass #Other string passed as the delimiter.


    #More options
    time_updates = options.time_updates
    display_seqs = options.display_seqs
    num_lines_for_display = options.num_lines_for_display
    seq_in_index = options.seq_in_index #where in the line the sequence is after line.split(delimiter)
    lines_to_skip = options.lines_to_skip #one method of skipping header
    comment_delimiter = options.comment_delimiter #another method of skipping header
    seqs_per_time_update = options.seqs_per_time_update
    max_number_of_seqs = options.max_number_of_seqs
    V_mask_index = options.V_mask_index #Default is not conditioning on V identity
    J_mask_index = options.J_mask_index #Default is not conditioning on J identity
    skip_empty = options.skip_empty

    seq_type_out = options.seq_type_out #type of pgens to be computed. Can be ntseq, aaseq, or both
    if seq_type_out is not None:
        seq_type_out = {'all': None, 'ntseq': 'ntseq', 'nucleotide': 'ntseq', 'aaseq': 'aaseq', 'amino_acid': 'aaseq'}[seq_type_out]

    if options.infile_name is None: #No infile specified -- args should be the input seqs
        print_warnings = True
        seqs = args
        seq_types = [determine_seq_type(seq, aa_alphabet) for seq in seqs]
        unrecognized_seqs = [seq for i, seq in enumerate(seqs) if seq_types[i] is None]
        if len(unrecognized_seqs) > 0 and print_warnings:
            print 'The following sequences/arguments were not recognized: ' + ', '.join(unrecognized_seqs)
        seqs = [seq for i, seq in enumerate(seqs) if seq_types[i] is not None]
        seq_types = [seq_type for seq_type in seq_types if seq_type is not None]


        #Format V and J masks -- uniform for all argument input sequences
        try:
            V_mask = options.V_mask.split(',')
            unrecognized_v_genes = [v for v in V_mask if v not in pgen_model.V_mask_mapping.keys()]
            V_mask = [v for v in V_mask if v in pgen_model.V_mask_mapping.keys()]
            if len(unrecognized_v_genes) > 0:
                print 'These V genes/alleles are not recognized: ' + ', '.join(unrecognized_v_genes)
            if len(V_mask) == 0:
                print 'No recognized V genes/alleles in the provided V_mask. Continuing without conditioning on V usage.'
                V_mask = None
        except AttributeError:
            V_mask = options.V_mask #Default is None, i.e. not conditioning on V identity

        try:
            J_mask = options.J_mask.split(',')
            unrecognized_j_genes = [j for j in J_mask if j not in pgen_model.J_mask_mapping.keys()]
            J_mask = [j for j in J_mask if j in pgen_model.J_mask_mapping.keys()]
            if len(unrecognized_j_genes) > 0:
                print 'These J genes/alleles are not recognized: ' + ', '.join(unrecognized_j_genes)
            if len(J_mask) == 0:
                print 'No recognized J genes/alleles in the provided J_mask. Continuing without conditioning on J usage.'
                J_mask = None
        except AttributeError:
            J_mask = options.J_mask #Default is None, i.e. not conditioning on J identity

        print ''
        start_time = time.time()
        for seq, seq_type in zip(seqs, seq_types):
            if seq_type == 'aaseq':
                c_pgen = pgen_model.compute_aa_CDR3_pgen(seq, V_mask, J_mask, print_warnings)
                print 'Pgen of the amino acid sequence ' + seq + ': ' + str(c_pgen)
                print ''
            elif seq_type == 'regex':
                c_pgen = pgen_model.compute_regex_CDR3_template_pgen(seq, V_mask, J_mask, print_warnings)
                print 'Pgen of the regular expression sequence ' + seq + ': ' + str(c_pgen)
                print ''
            elif seq_type == 'ntseq':
                if seq_type_out is None or seq_type_out == 'ntseq':
                    c_pgen_nt = pgen_model.compute_nt_CDR3_pgen(seq, V_mask, J_mask, print_warnings)
                    print 'Pgen of the nucleotide sequence ' + seq + ': ' + str(c_pgen_nt)
                if seq_type_out is None or seq_type_out == 'aaseq':
                    c_pgen_aa = pgen_model.compute_aa_CDR3_pgen(nt2aa(seq), V_mask, J_mask, print_warnings)
                    print 'Pgen of the amino acid sequence nt2aa(' + seq + ') = ' + nt2aa(seq) + ': ' + str(c_pgen_aa)
                print ''

        c_time = time.time() - start_time
        if c_time > 86400: #more than a day
            c_time_str = '%d days, %d hours, %d minutes, and %.2f seconds.'%(int(c_time)/86400, (int(c_time)/3600)%24, (int(c_time)/60)%60, c_time%60)
        elif c_time > 3600: #more than an hr
            c_time_str = '%d hours, %d minutes, and %.2f seconds.'%((int(c_time)/3600)%24, (int(c_time)/60)%60, c_time%60)
        elif c_time > 60: #more than a min
            c_time_str = '%d minutes and %.2f seconds.'%((int(c_time)/60)%60, c_time%60)
        else:
            c_time_str = '%.2f seconds.'%(c_time)

        print 'Completed pgen computation in: ' + c_time_str

    else: #Read sequences in from file
        print_warnings = False #Most cases of reading in from file should have warnings disabled
        seqs = []
        seq_types = []
        V_usage_masks = []
        J_usage_masks = []

        infile = open(infile_name, 'r')

        for i, line in enumerate(infile):
            if comment_delimiter is not None: #Default case -- no comments/header delimiter
                if line.startswith(comment_delimiter): #allow comments
                    continue
            if i < lines_to_skip:
                continue

            if delimiter is None: #Default delimiter is any whitespace
                split_line = line.split()
            else:
                split_line = line.split(delimiter)

            #Find the seq
            try:
                seq = split_line[seq_in_index].strip()
                if len(seq.strip()) == 0:
                    if skip_empty:
                        continue
                    else:
                        seqs.append(seq) #keep the blank seq as a placeholder
                        seq_types.append('aaseq')
                else:
                    seqs.append(seq)
                    seq_types.append(determine_seq_type(seq, aa_alphabet))
            except IndexError: #no index match for seq
                if skip_empty and len(line.strip()) == 0:
                    continue
                print 'seq_in_index is out of range'
                print 'Exiting...'
                infile.close()
                return -1

            #Find and format V_usage_mask
            if V_mask_index is None:
                V_usage_masks.append(None) #default mask
            else:
                try:
                    V_usage_mask = split_line[V_mask_index].strip().split(gene_mask_delimiter)
                    #check that all V gene/allele names are recognized
                    if all([v in pgen_model.V_mask_mapping for v in V_usage_mask]):
                        V_usage_masks.append(V_usage_mask)
                    else:
                        print str(V_usage_mask) + " is not a usable V_usage_mask composed exclusively of recognized V gene/allele names"
                        print 'Unrecognized V gene/allele names: ' + ', '.join([v for v in V_usage_mask if not v in pgen_model.V_mask_mapping.keys()])
                        print 'Exiting...'
                        infile.close()
                        return -1
                except IndexError: #no index match for V_mask_index
                    print 'V_mask_index is out of range'
                    print 'Exiting...'
                    infile.close()
                    return -1

            #Find and format J_usage_mask
            if J_mask_index is None:
                J_usage_masks.append(None) #default mask
            else:
                try:
                    J_usage_mask = split_line[J_mask_index].strip().split(gene_mask_delimiter)
                    #check that all V gene/allele names are recognized
                    if all([j in pgen_model.J_mask_mapping for j in J_usage_mask]):
                        J_usage_masks.append(J_usage_mask)
                    else:
                        print str(J_usage_mask) + " is not a usable J_usage_mask composed exclusively of recognized J gene/allele names"
                        print 'Unrecognized J gene/allele names: ' + ', '.join([j for j in J_usage_mask if not j in pgen_model.J_mask_mapping.keys()])
                        print 'Exiting...'
                        infile.close()
                        return -1
                except IndexError: #no index match for J_mask_index
                    print 'J_mask_index is out of range'
                    print 'Exiting...'
                    infile.close()
                    return -1

            if max_number_of_seqs is not None:
                if len(seqs) >= max_number_of_seqs:
                    break


        unrecognized_seqs = [seq for i, seq in enumerate(seqs) if seq_types[i] is None]
        if len(unrecognized_seqs) > 0 and len(unrecognized_seqs) < len(seqs):
            if print_warnings or options.outfile_name is not None:
                print 'Some strings read in were not parsed as sequences -- they will be omitted.'
                print 'Examples of improperly read strings: '
                for unrecognized_seq in unrecognized_seqs[:10]:
                    print unrecognized_seq
            seqs = [seq for i, seq in enumerate(seqs) if seq_types[i] is not None]
            V_usage_masks = [V_usage_mask for i, V_usage_mask in enumerate(V_usage_masks) if seq_types[i] is not None]
            seq_types = [seq_type for seq_type in seq_types if seq_type is not None]
        elif len(unrecognized_seqs) > 0 and len(unrecognized_seqs) == len(seqs):
            print 'None of the read in strings were parsed as sequences. Check input file.'
            print 'Examples of improperly read strings:'
            for unrecognized_seq in unrecognized_seqs[:10]:
                print unrecognized_seq
            print 'Exiting...'
            return -1

        infile.close()


        if options.outfile_name is not None: #OUTFILE SPECIFIED, allow printed info/display

            print 'Successfully read in and formatted ' + str(len(seqs)) + ' sequences and any V or J usages.'
            if display_seqs:
                sys.stdout.write('\r'+'Continuing to Pgen computation in 3... ')
                sys.stdout.flush()
                time.sleep(0.4)
                sys.stdout.write('\r'+'Continuing to Pgen computation in 2... ')
                sys.stdout.flush()
                time.sleep(0.4)
                sys.stdout.write('\r'+'Continuing to Pgen computation in 1... ')
                sys.stdout.flush()
                time.sleep(0.4)
            else:
                print 'Continuing to Pgen computation.'
                print_warnings = True #Display is off, can print warnings

            if display_seqs:
                lines_for_display = []
                times_for_speed_calc = [time.time()]

            outfile = open(outfile_name, 'w')
            start_time = time.time()
            for i, seq in enumerate(seqs):
                if seq_types[i] == 'aaseq':
                    #Compute Pgen and print out
                    c_pgen_line = seq + delimiter_out + str(pgen_model.compute_aa_CDR3_pgen(seq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                if seq_types[i] == 'regex':
                    #Compute Pgen and print out
                    c_pgen_line = seq + delimiter_out + str(pgen_model.compute_regex_CDR3_template_pgen(seq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                elif seq_types[i] == 'ntseq':
                    ntseq = seq
                    if len(ntseq) % 3 == 0: #inframe sequence
                        aaseq = nt2aa(ntseq)
                        #Compute Pgen and print out based on recomb_type and seq_type_out
                        if seq_type_out is None:
                            c_pgen_line = ntseq + delimiter_out + str(pgen_model.compute_nt_CDR3_pgen(ntseq, V_usage_masks[i], J_usage_masks[i], print_warnings)) + delimiter_out + aaseq + delimiter_out +  str(pgen_model.compute_aa_CDR3_pgen(aaseq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                        elif seq_type_out == 'ntseq':
                            c_pgen_line = ntseq + delimiter_out + str(pgen_model.compute_nt_CDR3_pgen(ntseq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                        elif seq_type_out == 'aaseq':
                            c_pgen_line = aaseq + delimiter_out + str(pgen_model.compute_aa_CDR3_pgen(aaseq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                    else: #out of frame sequence -- Pgens are 0 and use 'out_of_frame' for aaseq
                        if seq_type_out is None:
                            c_pgen_line = ntseq + delimiter_out + '0' + delimiter_out + 'out_of_frame' + delimiter_out + '0'
                        elif seq_type_out == 'ntseq':
                            c_pgen_line = ntseq + delimiter_out + '0'
                        elif seq_type_out == 'aaseq':
                            c_pgen_line = 'out_of_frame' + delimiter_out + '0'

                outfile.write(c_pgen_line + '\n')

                #Print time update
                if display_seqs:
                    cc_time = time.time()
                    c_time = cc_time - start_time
                    times_for_speed_calc = [cc_time] + times_for_speed_calc[:num_lines_for_display]
                    c_avg_speed = (len(times_for_speed_calc)-1)/float(times_for_speed_calc[0] - times_for_speed_calc[-1])

                    #eta = ((len(seqs) - (i+1))/float(i+1))*c_time

                    eta = (len(seqs) - (i+1))/c_avg_speed

                    lines_for_display = [c_pgen_line] + lines_for_display[:num_lines_for_display]


                    c_time_str = '%s hours, %s minutes, and %s seconds.'%(repr(int(c_time)/3600).rjust(3), repr((int(c_time)/60)%60).rjust(2), repr(int(c_time)%60).rjust(2))
                    eta_str = '%s hours, %s minutes, and %s seconds.'%(repr(int(eta)/3600).rjust(3), repr((int(eta)/60)%60).rjust(2), repr(int(eta)%60).rjust(2))
                    time_str = 'Time to compute Pgen on %s seqs: %s \nEst. time for remaining %s seqs: %s'%(repr(i+1).rjust(9), c_time_str, repr(len(seqs) - (i + 1)).rjust(9), eta_str)
                    speed_str = 'Current Pgen computation speed: %s seqs/min'%(repr(round((len(times_for_speed_calc)-1)*60/float(times_for_speed_calc[0] - times_for_speed_calc[-1]), 2)).rjust(8))
                    display_str = '\n'.join(lines_for_display[::-1]) + '\n' + '-'*80 + '\n' + time_str + '\n' + speed_str + '\n' + '-'*80
                    print '\033[2J' + display_str
                elif (i+1)%seqs_per_time_update == 0 and time_updates:
                    c_time = time.time() - start_time
                    eta = ((len(seqs) - (i+1))/float(i+1))*c_time
                    if c_time > 86400: #more than a day
                        c_time_str = '%d days, %d hours, %d minutes, and %.2f seconds.'%(int(c_time)/86400, (int(c_time)/3600)%24, (int(c_time)/60)%60, c_time%60)
                    elif c_time > 3600: #more than an hr
                        c_time_str = '%d hours, %d minutes, and %.2f seconds.'%((int(c_time)/3600)%24, (int(c_time)/60)%60, c_time%60)
                    elif c_time > 60: #more than a min
                        c_time_str = '%d minutes and %.2f seconds.'%((int(c_time)/60)%60, c_time%60)
                    else:
                        c_time_str = '%.2f seconds.'%(c_time)

                    if eta > 86400: #more than a day
                        eta_str = '%d days, %d hours, %d minutes, and %.2f seconds.'%(int(eta)/86400, (int(eta)/3600)%24, (int(eta)/60)%60, eta%60)
                    elif eta > 3600: #more than an hr
                        eta_str = '%d hours, %d minutes, and %.2f seconds.'%((int(eta)/3600)%24, (int(eta)/60)%60, eta%60)
                    elif eta > 60: #more than a min
                        eta_str = '%d minutes and %.2f seconds.'%((int(eta)/60)%60, eta%60)
                    else:
                        eta_str = '%.2f seconds.'%(eta)

                    print 'Pgen computed for %d sequences in: %s Estimated time remaining: %s'%(i+1, c_time_str, eta_str)

            c_time = time.time() - start_time
            if c_time > 86400: #more than a day
                c_time_str = '%d days, %d hours, %d minutes, and %.2f seconds.'%(int(c_time)/86400, (int(c_time)/3600)%24, (int(c_time)/60)%60, c_time%60)
            elif c_time > 3600: #more than an hr
                c_time_str = '%d hours, %d minutes, and %.2f seconds.'%((int(c_time)/3600)%24, (int(c_time)/60)%60, c_time%60)
            elif c_time > 60: #more than a min
                c_time_str = '%d minutes and %.2f seconds.'%((int(c_time)/60)%60, c_time%60)
            else:
                c_time_str = '%.2f seconds.'%(c_time)
            print 'Completed Pgen computation for %d sequences: in %s'%(len(seqs), c_time_str)

            outfile.close()

        else: #NO OUTFILE -- print directly to stdout
            start_time = time.time()
            for i, seq in enumerate(seqs):
                if seq_types[i] == 'aaseq':
                    #Compute Pgen and print out
                    c_pgen_line = seq + delimiter_out + str(pgen_model.compute_aa_CDR3_pgen(seq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                if seq_types[i] == 'regex':
                    #Compute Pgen and print out
                    c_pgen_line = seq + delimiter_out + str(pgen_model.compute_regex_CDR3_template_pgen(seq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                elif seq_types[i] == 'ntseq':
                    ntseq = seq
                    if len(ntseq) % 3 == 0: #inframe sequence
                        aaseq = nt2aa(ntseq)
                        #Compute Pgen and print out based on recomb_type and seq_type_out
                        if seq_type_out is None:
                            c_pgen_line = ntseq + delimiter_out + str(pgen_model.compute_nt_CDR3_pgen(ntseq, V_usage_masks[i], J_usage_masks[i], print_warnings)) + delimiter_out + aaseq + delimiter_out +  str(pgen_model.compute_aa_CDR3_pgen(aaseq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                        elif seq_type_out == 'ntseq':
                            c_pgen_line = ntseq + delimiter_out + str(pgen_model.compute_nt_CDR3_pgen(ntseq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                        elif seq_type_out == 'aaseq':
                            c_pgen_line = aaseq + delimiter_out + str(pgen_model.compute_aa_CDR3_pgen(aaseq, V_usage_masks[i], J_usage_masks[i], print_warnings))
                    else: #out of frame sequence -- Pgens are 0 and use 'out_of_frame' for aaseq
                        if seq_type_out is None:
                            c_pgen_line = ntseq + delimiter_out + '0' + delimiter_out + 'out_of_frame' + delimiter_out + '0'
                        elif seq_type_out == 'ntseq':
                            c_pgen_line = ntseq + delimiter_out + '0'
                        elif seq_type_out == 'aaseq':
                            c_pgen_line = 'out_of_frame' + delimiter_out + '0'

                print c_pgen_line
689777
def main(argv=None):
    

    if argv is None:
        argv = sys.argv

    # setup command line parser
    parser = U.OptionParser(version="%prog version: $Id$",
                            usage=usage,
                            description=globals()["__doc__"])

    group = U.OptionGroup(parser, "count-specific options")

    parser.add_option("--wide-format-cell-counts", dest="wide_format_cell_counts",
                      action="store_true",
                      default=False,
                      help=("output the cell counts in a wide format "
                            "(rows=genes, columns=cells)"))

    parser.add_option_group(group)

    # add common options (-h/--help, ...) and parse command line
    (options, args) = U.Start(parser, argv=argv, add_group_dedup_options=False)

    options.per_gene = True  # hardcodes counting to per-gene only

    U.validateSamOptions(options, group=False)

    if options.random_seed:
        np.random.seed(options.random_seed)

    if options.stdin != sys.stdin:
        in_name = options.stdin.name
        options.stdin.close()
    else:
        raise ValueError("Input on standard in not currently supported")

    if options.in_sam:
        in_mode = "r"
    else:
        in_mode = "rb"

    infile = pysam.Samfile(in_name, in_mode)

    # write out to tempfile and then sort to stdout
    tmpfilename = U.getTempFilename(dir=options.tmpdir)
    tmpfile = U.openFile(tmpfilename, mode="w")

    nInput, nOutput, input_reads = 0, 0, 0

    gene_tag = options.gene_tag
    metacontig2contig = None

    if options.chrom:
        inreads = infile.fetch(reference=options.chrom)
    else:
        if options.gene_transcript_map:
            metacontig2contig = sam_methods.getMetaContig2contig(
                infile, options.gene_transcript_map)
            metatag = "MC"
            inreads = sam_methods.metafetcher(infile, metacontig2contig, metatag)
            gene_tag = metatag
        else:
            inreads = infile.fetch()

    bundle_iterator = sam_methods.get_bundles(
        options,
        only_count_reads=True,
        metacontig_contig=metacontig2contig)

    # set up UMIClusterer functor with methods specific to
    # specified options.method

    processor = network.UMIClusterer(options.method)

    for bundle, key, status in bundle_iterator(inreads):
        if status == "single_read":
            continue

        gene, cell = key

        umis = bundle.keys()
        counts = {umi: bundle[umi]["count"] for umi in umis}

        nInput += sum(counts.values())

        while nInput >= input_reads + 1000000:
            input_reads += 1000000
            U.info("Parsed %i input reads" % input_reads)

        # group the umis
        groups = processor(
            counts,
            threshold=options.threshold)

        gene_count = len(groups)

        if options.per_cell:
            tmpfile.write("%s\n" % "\t".join((gene, cell.decode(), str(gene_count))))
        else:
            tmpfile.write("%s\n" % "\t".join((gene, str(gene_count))))

        nOutput += gene_count

    tmpfile.close()

    if options.per_cell:

        gene_counts_dict = {}

        with U.openFile(tmpfilename, mode="r") as inf:
            genes = set()
            cells = set()
            for line in inf:
                gene, cell, gene_count = line.strip().split("\t")
                genes.add(gene)
                cells.add(cell)

                if gene not in gene_counts_dict:
                    gene_counts_dict[gene] = {}

                gene_counts_dict[gene][cell] = gene_count

        if options.wide_format_cell_counts:  # write out in wide format

            options.stdout.write(
                "%s\t%s\n" % ("gene", "\t".join(sorted(cells))))

            for gene in sorted(genes):
                counts = []
                for cell in sorted(cells):
                    if cell in gene_counts_dict[gene]:
                        counts.append(gene_counts_dict[gene][cell])
                    else:
                        counts.append(0)
                options.stdout.write(
                    "%s\t%s\n" % (gene, "\t".join(map(str, counts))))

        else:  # write out in long format
            options.stdout.write("%s\t%s\t%s\n" % ("gene", "cell", "count"))
            for gene in sorted(genes):
                for cell in sorted(list(gene_counts_dict[gene].keys())):
                    options.stdout.write("%s\t%s\t%s\n" % (
                        gene, cell, gene_counts_dict[gene][cell]))
    else:
        options.stdout.write("%s\t%s\n" % ("gene", "count"))

        with U.openFile(tmpfilename, mode="r") as inf:
            for line in inf:
                options.stdout.write(line)

    os.unlink(tmpfilename)

    # output reads events and benchmark information.
    for event in bundle_iterator.read_events.most_common():
        U.info("%s: %s" % (event[0], event[1]))

    U.info("Number of (post deduplication) reads counted: %i" % nOutput)

    U.Stop()
594565
def genstis(outname):
    
    pattern="""class stisS%d(countrateCase):
    def setUp(self):
        self.obsmode="%s"
        self.spectrum="%s"
        self.setglobal(__file__)
        self.runpy()\n"""

    speclist=['/grp/hst/cdbs/calspec/gd71_mod_005.fits',
              '/grp/hst/cdbs/calspec/gd153_mod_004.fits',
              '/grp/hst/cdbs/calspec/g191b2b_mod_004.fits']

    glist={'g140l':'fuvmama','g230l':'nuvmama','g430l':'ccd','g750l':'ccd',
    'g230lb':'ccd'}

    out=open(outname,'a')
    out.write("""from pytools import testutil
import sys
from basecase import calcphotCase, calcspecCase, countrateCase,SpecSourcerateSpecCase\n
""")
    count=0
    for g in glist:
        for sp in speclist:
            obsmode='stis,%s,fuvmama,s52x2'%g
            defn=pattern%(count,obsmode,sp)
            out.write(defn)
            count+=1

    out.write("""\n\n
if __name__ == '__main__':
    if 'debug' in sys.argv:
        testutil.debug(__name__)
    else:
        testutil.testall(__name__,2)
""")

    out.close()
992478
def main():

    gtselect_keys = ['tmin', 'tmax', 'emin', 'emax', 'zmax', 'evtype', 'evclass',
                     'phasemin', 'phasemax', 'convtype', 'rad', 'ra', 'dec']

    gtmktime_keys = ['roicut', 'filter']

    usage = "usage: %(prog)s [options] "
    description = "Run gtselect and gtmktime on one or more FT1 files.  "
    ""
    parser = argparse.ArgumentParser(usage=usage, description=description)

    add_lsf_args(parser)

    for k in gtselect_keys:

        if k in ['evtype', 'evclass', 'convtype']:
            parser.add_argument('--%s' % k, default=None, type=int, help='')
        else:
            parser.add_argument('--%s' % k, default=None, type=float, help='')

    for k in gtmktime_keys:
        parser.add_argument('--%s' % k, default=None, type=str, help='')

    parser.add_argument('--rock_angle', default=None, type=float, help='')

    parser.add_argument('--outdir', default=None, type=str,
                        help='Path to output directory used when merge=False.')
    parser.add_argument('--output', default=None, type=str,
                        help='Path to output file used when merge=True.')
    parser.add_argument('--scfile', default=None, type=str, help='')

    parser.add_argument('--dry_run', default=False, action='store_true')
    parser.add_argument('--overwrite', default=False, action='store_true')
    parser.add_argument('--merge', default=False, action='store_true',
                        help='Merge input FT1 files into N files where N is determined '
                        'by files_per_split.')

    parser.add_argument('--files_per_split', default=100,
                        type=int, help='Set the number of files to combine in each '
                        'split of the input file list.')

    parser.add_argument('--file_idx_min', default=None,
                        type=int, help='Set the number of files to assign to '
                        'each batch job.')

    parser.add_argument('--file_idx_max', default=None,
                        type=int, help='Set the number of files to assign to '
                        'each batch job.')

    parser.add_argument('files', nargs='+', default=None,
                        help='List of files.')

    args = parser.parse_args()

    batch = vars(args).pop('batch')
    files = vars(args).pop('files')
    args.outdir = os.path.abspath(args.outdir)
    files = [os.path.abspath(f) for f in files]

    ft1_files = get_files(files, ['.fit', '.fits'])
    for i, f in enumerate(ft1_files):
        if re.search('^root\:\/\/', f) is None:
            ft1_files[i] = os.path.abspath(f)

    input_files = []
    output_files = []
    files_idx_min = []
    files_idx_max = []
    opts = []

    if args.file_idx_min is not None and args.file_idx_max is not None:

        files_idx_min = [args.file_idx_min]
        files_idx_max = [args.file_idx_max]
        input_files = [files]
        output_files = [args.output]

    elif args.merge:
        if not args.output:
            raise Exception('No output file defined.')

        nfiles = len(ft1_files)
        njob = int(np.ceil(nfiles / float(args.files_per_split)))
        for ijob, i in enumerate(range(0, nfiles, args.files_per_split)):

            if args.outdir is not None:
                mkdir(args.outdir)
                outdir = os.path.abspath(args.outdir)
            else:
                outdir = os.path.dirname(os.path.dirname(args.output))

            outfile = os.path.splitext(os.path.basename(args.output))[0]
            outfile += '_%03i.fits' % (ijob)
            outfile = os.path.join(outdir, outfile)
            input_files += [files]
            output_files += [outfile]
            files_idx_min += [i]
            files_idx_max += [i + args.files_per_split]
            opts += [vars(args).copy()]
            opts[-1]['output'] = outfile
            opts[-1]['file_idx_min'] = i
            opts[-1]['file_idx_max'] = i + args.files_per_split

    else:
        input_files = ft1_files
        files_idx_min = [i for i in range(len(ft1_files))]
        files_idx_max = [i + 1 for i in range(len(ft1_files))]
        output_files = [os.path.join(
            args.outdir, os.path.basename(x)) for x in ft1_files]
        opts = [vars(args).copy() for x in ft1_files]

    if batch:
        submit_jobs('fermipy-select',
                    input_files, opts, output_files, overwrite=args.overwrite,
                    dry_run=args.dry_run)
        sys.exit(0)

    logger = Logger.configure(os.path.basename(__file__), None, logging.INFO)
    logger.info('Starting.')

    if args.scfile is not None:
        args.scfile = os.path.abspath(args.scfile)

    cwd = os.getcwd()
    user = os.environ['USER']
    tmpdir = tempfile.mkdtemp(prefix=user + '.', dir='/scratch')
    os.chdir(tmpdir)

    logger.info('tmpdir %s', tmpdir)
    logger.info('outdir %s', args.outdir)
    logger.info('output %s', args.output)

    for infiles, outfile, idx_min, idx_max in zip(input_files, output_files,
                                                  files_idx_min, files_idx_max):

        logger.info('infiles %s', pprint.pformat(infiles))
        logger.info('outfile %s', outfile)
        infiles = get_files(infiles, ['.fit', '.fits'])
        if idx_min is not None:
            infiles = infiles[idx_min:idx_max]

        for i, f in enumerate(infiles):

            if re.search('^root\:\/\/', f) is None:
                continue
            os.system('xrdcp %s %s' % (f, f.split('/')[-1]))
            infiles[i] = os.path.join(tmpdir, f.split('/')[-1])

        kw = {k: args.__dict__[k] for k in gtselect_keys}
        if kw['emax'] is None:
            kw['emax'] = 1E6

        create_filelist(infiles, 'list.txt')
        kw['infile'] = 'list.txt'
        kw['outfile'] = 'out.fits'
        staged_outfile = kw['outfile']
        run_gtapp('gtselect', logger, kw)

        kw = {k: args.__dict__[k] for k in gtmktime_keys}
        if kw['roicut'] is None:
            kw['roicut'] = 'no'

        if kw['filter'] is None:
            kw['filter'] = 'DATA_QUAL==1 && LAT_CONFIG==1'
            if args.rock_angle is not None:
                kw['filter'] += ' && ABS(ROCK_ANGLE)<%(rock)s ' % dict(
                    rock=args.rock_angle)
        kw['evfile'] = 'out.fits'
        kw['outfile'] = 'out_filtered.fits'
        if args.scfile is not None:
            kw['scfile'] = args.scfile
            staged_outfile = kw['outfile']
            run_gtapp('gtmktime', logger, kw)

        logger.info('cp %s %s', staged_outfile, outfile)
        shutil.copy(staged_outfile, outfile)

    os.chdir(cwd)
    logger.info('Deleting %s', tmpdir)
    shutil.rmtree(tmpdir)
    logger.info('Done.')
1032106
def main():
    
    # Define Parser object and add to toil
    url_prefix = 'https://s3-us-west-2.amazonaws.com/cgl-pipeline-inputs/'
    parser = argparse.ArgumentParser(description=main.__doc__, formatter_class=argparse.RawTextHelpFormatter)
    parser.add_argument('--config', required=True,
                        help='Path to configuration file for samples, one per line. UUID,URL_to_bamfile. '
                             'The URL may be a standard "http://", a "file://<abs_path>", or "s3://<bucket>/<key>"')
    parser.add_argument('--gtf', help='URL to annotation GTF file',
                        default=url_prefix + 'rnaseq_cgl/gencode.v23.annotation.gtf')
    parser.add_argument('--gtf-pickle', help='Pickled GTF file',
                        default=url_prefix + 'spladder/gencode.v23.annotation.gtf.pickle')
    parser.add_argument('--gtf-m53', help='M53 preprocessing annotation table',
                        default=url_prefix + 'spladder/gencode.v23.annotation.gtf.m53')
    parser.add_argument('--positions', help='URL to SNP positions over genes file (TSV)',
                        default=url_prefix + 'spladder/positions_fixed.tsv')
    parser.add_argument('--genome', help='URL to Genome fasta',
                        default=url_prefix + 'rnaseq_cgl/hg38_no_alt.fa')
    parser.add_argument('--genome-index', help='Index file (fai) of genome',
                        default=url_prefix + 'spladder/hg38_no_alt.fa.fai')
    parser.add_argument('--ssec', default=None, help='Path to master key used for downloading encrypted files.')
    parser.add_argument('--output-s3-dir', default=None, help='S3 Directory of the form: s3://bucket/directory')
    parser.add_argument('--output-dir', default=None, help='full path where final results will be output')
    parser.add_argument('--sudo', action='store_true', default=False,
                        help='Set flag if sudo is required to run Docker.')
    parser.add_argument('--star-index', help='URL to download STAR Index built from HG38/gencodev23 annotation.',
                        default=url_prefix + 'rnaseq_cgl/starIndex_hg38_no_alt.tar.gz')
    parser.add_argument('--fwd-3pr-adapter', help="Sequence for the FWD 3' Read Adapter.", default='AGATCGGAAGAG')
    parser.add_argument('--rev-3pr-adapter', help="Sequence for the REV 3' Read Adapter.", default='AGATCGGAAGAG')
    Job.Runner.addToilOptions(parser)
    args = parser.parse_args()
    # Sanity Checks
    if args.config:
        assert os.path.isfile(args.config), 'Config not found at: {}'.format(args.config)
    if args.ssec:
        assert os.path.isfile(args.ssec), 'Encryption key not found at: {}'.format(args.config)
    if args.output_s3_dir:
        assert args.output_s3_dir.startswith('s3://'), 'Wrong format for output s3 directory'
    # Program checks
    for program in ['curl', 'docker']:
        assert which(program), 'Program "{}" must be installed on every node.'.format(program)

    Job.Runner.startToil(Job.wrapJobFn(parse_input_samples, args), args)
839679
def tx2genedict(gtf, keep_version=False):
    
    d = {}
    with open_gzipsafe(gtf) as in_handle:
        for line in in_handle:
            if "gene_id" not in line or "transcript_id" not in line:
                continue
            geneid = line.split("gene_id")[1].split(" ")[1]
            geneid = _strip_non_alphanumeric(geneid)
            txid = line.split("transcript_id")[1].split(" ")[1]
            txid = _strip_non_alphanumeric(txid)
            if keep_version and "transcript_version" in line:
                txversion = line.split("transcript_version")[1].split(" ")[1]
                txversion = _strip_non_alphanumeric(txversion)
                txid  += "." + txversion
            if has_transcript_version(line) and not keep_version:
                txid = _strip_feature_version(txid)
                geneid = _strip_feature_version(geneid)
            d[txid] = geneid
    return d
495930
def genargs() -> ArgumentParser:
    
    parser = ArgumentParser()
    parser.add_argument("infile", help="Input ShExC specification")
    parser.add_argument("-nj", "--nojson", help="Do not produce json output", action="store_true")
    parser.add_argument("-nr", "--nordf", help="Do not produce rdf output", action="store_true")
    parser.add_argument("-j", "--jsonfile", help="Output ShExJ file (Default: {infile}.json)")
    parser.add_argument("-r", "--rdffile", help="Output ShExR file (Default: {infile}.{fmt suffix})")
    parser.add_argument("--context", help="Alternative @context")
    parser.add_argument("-f", "--format",
                        choices=list(set(x.name for x in rdflib_plugins(None, rdflib_Serializer)
                                         if '/' not in str(x.name))),
                        help="Output format (Default: turtle)", default="turtle")
    return parser
1150902
def main(args):
   """"""
   txome = Transcriptome()
   #read the reference gpd if one is gven
   if args.reference: 
      rinf = None
      if re.search('\.gz$',args.reference):
         rinf = gzip.open(args.reference)
      else:
         rinf = open(args.reference)
      sys.stderr.write("Reading in reference\n")
      z = 0
      # populate txome with reference transcripts for each chromosome
      for line in rinf:
         z += 1
         gpd = GPD(line)
         gpd.set_payload(z)
         if z%100 == 0:  sys.stderr.write(str(z)+"          \r")
         txome.add_transcript(gpd)
      rinf.close()
      sys.stderr.write(str(z)+"          \r")
      sys.stderr.write("\n")
   txome.sort_transcripts()
   sys.stderr.write("Buffering mappings\n")
   inf = sys.stdin
   if args.input != '-':
      if re.search('\.gz$',args.input):
         inf = gzip.open(args.input)
      else:
         inf = open(args.input)
   tof = gzip.open(args.tempdir+'/reads.gpd.gz','w')
   for line in inf: tof.write(line.rstrip()+"\n")
   tof.close()
  
   sys.stderr.write("1. Process by overlapping locus.\n")
   annotated_reads(txome,gzip.open(args.tempdir+'/reads.gpd.gz'),'initial',args)

   #of.close()

   """Now sequences we should be able to annotate as partial have
      consensus transcripts.  Lets put those to annotation scrutiny again.
   """
   sys.stderr.write("2. Reannotated partially annotated by overlapping locus.\n")
   annotated_reads(txome,gzip.open(args.tempdir+'/initial/partial_annotated_multiexon.sorted.gpd.gz'),'repartial',args)


   txs = {}
   for transcript in txome.transcripts:
      txs[transcript.name] = transcript
   detected = {}
   tinf = gzip.open(args.tempdir+'/initial/annotated.txt.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      detected[f[1]] = True
   tinf.close()
   tinf = gzip.open(args.tempdir+'/repartial/annotated.txt.gz')
   remove_partial = {}
   for line in tinf:
      f = line.rstrip().split("\t")
      detected[f[1]] = True
      remove_partial[f[1]] = True
   tinf.close()
   tof = gzip.open(args.tempdir+'/candidate.gpd.gz','w')
   for name in detected:
      tof.write(txs[name].get_gpd_line()+"\n")
   tinf = gzip.open(args.tempdir+'/initial/partial_annotated_multiexon.sorted.gpd.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      if f[1] not in remove_partial:
         tof.write(line)
   tinf = gzip.open(args.tempdir+'/initial/unannotated_singleexon.sorted.gpd.gz')
   for line in tinf:
      f = line.rstrip().split("\t")
      f[0] = str(uuid.uuid4())
      tof.write("\t".join(f)+"\n")
   tinf.close()
   tof.close()
   sort_gpd(args.tempdir+'/candidate.gpd.gz',args.tempdir+'/candidate.sorted.gpd.gz',args)
   """We can ignore the partial annotations that have been detected"""
   ntxome = Transcriptome()
   tinf = gzip.open(args.tempdir+'/candidate.sorted.gpd.gz')
   for line in tinf: ntxome.add_transcript(GPD(line))     
   annotated_reads(ntxome,gzip.open(args.tempdir+'/initial/unannotated_multiexon.sorted.gpd.gz'),'unannotated',args)
   """now we know which unannotated reads actually have annotations"""
   tinf.close()

   tof = gzip.open(args.tempdir+'/final.gpd.gz','w')
   tinf = gzip.open(args.tempdir+'/candidate.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tinf = gzip.open(args.tempdir+'/unannotated/unannotated_multiexon.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tinf = gzip.open(args.tempdir+'/unannotated/partial_annotated_multiexon.sorted.gpd.gz')   
   for line in tinf: tof.write(line)
   tinf.close()
   tof.close()
   sort_gpd(args.tempdir+'/final.gpd.gz',args.tempdir+'/final.sorted.gpd.gz',args)
   """Prepare outputs"""
   of = sys.stdout
   if args.output:
      if re.search('\.gz$',args.output):
         of = gzip.open(args.output,'w')
      else:
         of = open(args.output,'w')
   tinf = gzip.open(args.tempdir+'/final.sorted.gpd.gz')
   for line in tinf: of.write(line)
   tinf.close()
   of.close()
671554
def _get_input_args(bam_file, data, out_base, background):
    
    if dd.get_genome_build(data) in ["hg19"]:
        return ["--PileupFile", _create_pileup(bam_file, data, out_base, background)]
    else:
        return ["--BamFile", bam_file]
1341348
def csv_writer(parser, keep, extract, args):
    """"""
    # The output
    output = sys.stdout if args.output == "-" else open(args.output, "w")

    try:
        # Getting the samples
        samples = np.array(parser.get_samples(), dtype=str)
        k = _get_sample_select(samples=samples, keep=keep)

        # Writing the CSV header
        print("sample_id", "variant_id", "chromosome", "position", "reference",
              "coded", "dosage", "hard_call", sep=",", file=output)

        # The data generator
        generator = _get_generator(parser=parser, extract=extract, keep=k,
                                   check_maf=args.maf)

        # The number of markers extracted
        nb_extracted = 0

        for data in generator:
            # Keeping only the required genotypes
            genotypes = data.genotypes

            # The hard call mapping
            hard_call_mapping = {
                0: "{ref}/{ref}".format(ref=data.reference),
                1: "{ref}/{alt}".format(ref=data.reference, alt=data.coded),
                2: "{alt}/{alt}".format(alt=data.coded),
            }

            for sample, geno in zip(samples[k], genotypes):
                # Is the genotype missing
                is_missing = np.isnan(geno)

                # Hard coding (NaN values are empty string)
                hard_coded = None
                if is_missing:
                    geno = ""
                    hard_coded = ""
                else:
                    hard_coded = hard_call_mapping[int(round(geno, 0))]

                print(sample, data.variant.name, data.variant.chrom,
                      data.variant.pos, data.reference, data.coded,
                      geno, hard_coded, sep=",", file=output)

            nb_extracted += 1

        if nb_extracted == 0:
            logger.warning("No markers matched the extract list")

    finally:
        output.close()
983376
def entrez(args):
    
    p = OptionParser(entrez.__doc__)

    allowed_databases = {"fasta": ["genome", "nuccore", "nucgss", "protein", "nucest"],
                         "asn.1": ["genome", "nuccore", "nucgss", "protein", "gene"],
                         "xml": ["genome", "nuccore", "nucgss", "nucest", "gene"],
                         "gb": ["genome", "nuccore", "nucgss"],
                         "est": ["nucest"],
                         "gss": ["nucgss"],
                         "acc": ["nuccore"],
                         }

    valid_formats = tuple(allowed_databases.keys())
    valid_databases = ("genome", "nuccore", "nucest",
                       "nucgss", "protein", "gene")

    p.add_option("--noversion", dest="noversion",
                 default=False, action="store_true",
                 help="Remove trailing accession versions")
    p.add_option("--format", default="fasta", choices=valid_formats,
                 help="download format [default: %default]")
    p.add_option("--database", default="nuccore", choices=valid_databases,
                 help="search database [default: %default]")
    p.add_option("--retmax", default=1000000, type="int",
                 help="how many results to return [default: %default]")
    p.add_option("--skipcheck", default=False, action="store_true",
                 help="turn off prompt to check file existence [default: %default]")
    p.add_option("--batchsize", default=500, type="int",
                 help="download the results in batch for speed-up [default: %default]")
    p.set_outdir(outdir=None)
    p.add_option("--outprefix", default="out",
                 help="output file name prefix [default: %default]")
    p.set_email()
    opts, args = p.parse_args(args)

    if len(args) != 1:
        sys.exit(p.print_help())

    filename, = args
    if op.exists(filename):
        pf = filename.rsplit(".", 1)[0]
        list_of_terms = [row.strip() for row in open(filename)]
        if opts.noversion:
            list_of_terms = [x.rsplit(".", 1)[0] for x in list_of_terms]
    else:
        pf = filename
        # the filename is the search term
        list_of_terms = [filename.strip()]

    fmt = opts.format
    database = opts.database
    batchsize = opts.batchsize

    assert database in allowed_databases[fmt], \
        "For output format '{0}', allowed databases are: {1}".\
        format(fmt, allowed_databases[fmt])
    assert batchsize >= 1, "batchsize must >= 1"

    if " " in pf:
        pf = opts.outprefix

    outfile = "{0}.{1}".format(pf, fmt)

    outdir = opts.outdir
    if outdir:
        mkdir(outdir)

    # If noprompt, will not check file existence
    if not outdir:
        fw = must_open(outfile, "w", checkexists=True,
                       skipcheck=opts.skipcheck)
        if fw is None:
            return

    seen = set()
    totalsize = 0
    for id, size, term, handle in batch_entrez(list_of_terms, retmax=opts.retmax,
                                               rettype=fmt, db=database, batchsize=batchsize,
                                               email=opts.email):
        if outdir:
            outfile = urljoin(outdir, "{0}.{1}".format(term, fmt))
            fw = must_open(outfile, "w", checkexists=True,
                           skipcheck=opts.skipcheck)
            if fw is None:
                continue

        rec = handle.read()
        if id in seen:
            logging.error("Duplicate key ({0}) found".format(rec))
            continue

        totalsize += size
        print(rec, file=fw)
        print(file=fw)

        seen.add(id)

    if seen:
        print("A total of {0} {1} records downloaded.".
              format(totalsize, fmt.upper()), file=sys.stderr)

    return outfile
63542
def gatk(args):
    
    p = OptionParser(gatk.__doc__)
    p.add_option("--indelrealign", default=False, action="store_true",
                 help="Perform indel realignment")
    p.set_home("gatk")
    p.set_home("picard")
    p.set_phred()
    p.set_cpus(cpus=24)
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    bamfile, ref = args
    pf = bamfile.rsplit(".", 1)[0]
    mm = MakeManager()
    picard = "java -Xmx32g -jar {0}/picard.jar".format(opts.picard_home)
    tk = "java -Xmx32g -jar {0}/GenomeAnalysisTK.jar".format(opts.gatk_home)
    tk += " -R {0}".format(ref)

    # Step 0 - build reference
    dictfile = ref.rsplit(".", 1)[0] + ".dict"
    cmd1 = picard + " CreateSequenceDictionary"
    cmd1 += " R={0} O={1}".format(ref, dictfile)
    cmd2 = "samtools faidx {0}".format(ref)
    mm.add(ref, dictfile, (cmd1, cmd2))

    # Step 1 - sort bam
    sortedbamfile = pf + ".sorted.bam"
    cmd = picard + " SortSam"
    cmd += " INPUT={0} OUTPUT={1}".format(bamfile, sortedbamfile)
    cmd += " SORT_ORDER=coordinate CREATE_INDEX=true"
    mm.add(bamfile, sortedbamfile, cmd)

    # Step 2 - mark duplicates
    dedupbamfile = pf + ".dedup.bam"
    cmd = picard + " MarkDuplicates"
    cmd += " INPUT={0} OUTPUT={1}".format(sortedbamfile, dedupbamfile)
    cmd += " METRICS_FILE=dedup.log CREATE_INDEX=true"
    mm.add(sortedbamfile, dedupbamfile, cmd)

    if opts.indelrealign:
        # Step 3 - create indel realignment targets
        intervals = pf + ".intervals"
        cmd = tk + " -T RealignerTargetCreator"
        cmd += " -I {0} -o {1}".format(dedupbamfile, intervals)
        mm.add(dedupbamfile, intervals, cmd)

        # Step 4 - indel realignment
        realignedbamfile = pf + ".realigned.bam"
        cmd = tk + " -T IndelRealigner"
        cmd += " -targetIntervals {0}".format(intervals)
        cmd += " -I {0} -o {1}".format(dedupbamfile, realignedbamfile)
        mm.add((dictfile, intervals), realignedbamfile, cmd)
    else:
        realignedbamfile = dedupbamfile

    # Step 5 - SNP calling
    vcf = pf + ".vcf"
    cmd = tk + " -T HaplotypeCaller"
    cmd += " -I {0}".format(realignedbamfile)
    cmd += " --genotyping_mode DISCOVERY"
    cmd += " -stand_emit_conf 10 -stand_call_conf 30"
    cmd += " -nct {0}".format(opts.cpus)
    cmd += " -o {0}".format(vcf)
    if opts.phred == "64":
        cmd += " --fix_misencoded_quality_scores"
    mm.add(realignedbamfile, vcf, cmd)

    # Step 6 - SNP filtering
    filtered_vcf = pf + ".filtered.vcf"
    cmd = tk + " -T VariantFiltration"
    cmd += " -V {0}".format(vcf)
    cmd += ' --filterExpression "DP < 10 || DP > 300 || QD < 2.0 || FS > 60.0 || MQ < 40.0"'
    cmd += ' --filterName "LOWQUAL"'
    cmd += ' --genotypeFilterExpression "isHomVar == 1"'
    cmd += ' --genotypeFilterName "HOMOVAR"'
    cmd += ' --genotypeFilterExpression "isHet == 1"'
    cmd += ' --genotypeFilterName "HET"'
    cmd += " -o {0}".format(filtered_vcf)
    mm.add(vcf, filtered_vcf, cmd)

    mm.write()
1010297
def get_argument_parser():
    
    file_mv = cli.file_mv

    desc = 'Extracts gene-level expression data from StringTie output.'
    parser = cli.get_argument_parser(desc)

    parser.add_argument(
        '-s', '--stringtie-file', type=str, required=True, metavar=file_mv,
        help="""Path of the StringTie output file ."""
    )

    parser.add_argument(
        '-g', '--gene-file', type=str, required=True, metavar=file_mv,
        help="""File containing a list of protein-coding genes."""
    )

    parser.add_argument(
        '--no-novel-transcripts', action='store_true',
        help="""Ignore novel transcripts."""
    )

    # parser.add_argument(
    #     '--ambiguous-transcripts', default = 'ignore',
    #      help='Strategy for counting expression of ambiguous novel '
    #            'transcripts.'
    # )
    # possible strategies for ambiguous transcripts: 'ignore','highest','all'

    parser.add_argument(
        '-o', '--output-file', type=str, required=True, metavar=file_mv,
        help="""Path of output file."""
    )

    cli.add_reporting_args(parser)

    return parser
759632
def parse_ggKbase_tables(tables, id_type):
    
    g2info = {}
    for table in tables:
        for line in open(table):
            line = line.strip().split('\t')
            if line[0].startswith('name'):
                header = line
                header[4] = 'genome size (bp)'
                header[12] = '#SCGs'
                header[13] = '#SCG duplicates'
                continue
            name, code, info = line[0], line[1], line
            info = [to_int(i) for i in info]
            if id_type is False: # try to use name and code ID
                if 'UNK' in code or 'unknown' in code:
                    code = name
                if (name != code) and (name and code in g2info):
                    print('# duplicate name or code in table(s)', file=sys.stderr)
                    print('# %s and/or %s' % (name, code), file=sys.stderr)
                    exit()
                if name not in g2info:
                    g2info[name] = {item:stat for item, stat in zip(header, info)}
                if code not in g2info:
                    g2info[code] = {item:stat for item, stat in zip(header, info)}
            else:
                if id_type == 'name':
                    ID = name
                elif id_type == 'code':
                    ID = code
                else:
                    print('# specify name or code column using -id', file=sys.stderr)
                    exit()
                ID = ID.replace(' ', '')
                g2info[ID] = {item:stat for item, stat in zip(header, info)}
                if g2info[ID]['genome size (bp)'] == '':
                    g2info[ID]['genome size (bp)'] = 0
    return g2info
172990
def GeneReader( fh, format='gff' ):
    

    known_formats = ( 'gff', 'gtf', 'bed')
    if format not in known_formats: 
        print('%s format not in %s' % (format, ",".join( known_formats )), file=sys.stderr)
        raise Exception('?')

    if format == 'bed':
        for line in fh:    
            f = line.strip().split()
            chrom = f[0]
            chrom_start = int(f[1])
            name = f[4]
            strand = f[5]
            cdsStart = int(f[6])
            cdsEnd = int(f[7])
            blockCount = int(f[9])
            blockSizes = [ int(i) for i in f[10].strip(',').split(',') ]
            blockStarts = [ chrom_start + int(i) for i in f[11].strip(',').split(',') ]

            # grab cdsStart - cdsEnd
            gene_exons = []
            for base,offset in zip( blockStarts, blockSizes ):
                exon_start = base
                exon_end = base+offset
                gene_exons.append( (exon_start, exon_end) )
            
    genelist = {}
    grouplist = []
    if format == 'gff' or format == 'gtf':
        for line in fh:
            if line.startswith('#'): continue
            fields = line.strip().split('\t')
            if len( fields ) < 9: continue

            # fields

            chrom = fields[0]
            ex_st = int( fields[3] ) - 1 # make zero-centered
            ex_end = int( fields[4] ) #+ 1 # make exclusive
            strand = fields[6]

            if format == 'gtf':
                group = fields[8].split(';')[0]
            else:
                group = fields[8]

            if group not in grouplist: grouplist.append( group )
            if group not in genelist:
                genelist[group] = (chrom, strand, [])
            exons_i = 2
            genelist[group][exons_i].append( ( ex_st, ex_end ) )

        sp = lambda a,b: cmp( a[0], b[0] )

        #for gene in genelist.values():
        for gene in grouplist:
            chrom, strand, gene_exons = genelist[ gene ]
            gene_exons = bitset_union( gene_exons )
            yield chrom, strand, gene_exons, gene
865809
def pasa(args):
    
    from jcvi.formats.base import SetFile
    from jcvi.formats.gff import Gff

    p = OptionParser(pasa.__doc__)
    p.set_home("pasa")
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    fastafile, gffile = args
    transcodergff = fastafile + ".transdecoder.gff3"
    transcodergenomegff = fastafile + ".transdecoder.genome.gff3"
    if need_update((fastafile, gffile), (transcodergff, transcodergenomegff)):
        cmd = "{0}/scripts/pasa_asmbls_to_training_set.dbi".format(opts.pasa_home)
        cmd += " --pasa_transcripts_fasta {0} --pasa_transcripts_gff3 {1}".\
                format(fastafile, gffile)
        sh(cmd)

    completeids = fastafile.rsplit(".", 1)[0] + ".complete.ids"
    if need_update(transcodergff, completeids):
        cmd = "grep complete {0} | cut -f1 | sort -u".format(transcodergff)
        sh(cmd, outfile=completeids)

    complete = SetFile(completeids)
    seen = set()
    completegff = transcodergenomegff.rsplit(".", 1)[0] + ".complete.gff3"
    fw = open(completegff, "w")
    gff = Gff(transcodergenomegff)
    for g in gff:
        a = g.attributes
        if "Parent" in a:
            id = a["Parent"][0]
        else:
            id = a["ID"][0]
        asmbl_id = id.split("|")[0]
        if asmbl_id not in complete:
            continue
        print(g, file=fw)
        if g.type == "gene":
            seen.add(id)

    fw.close()
    logging.debug("A total of {0} complete models extracted to `{1}`.".\
                    format(len(seen), completegff))