find 16S rRNA gene sequence coordinates
515491
def find_coordinates(hmms, bit_thresh):
    
    # get coordinates from cmsearch output
    seq2hmm = parse_hmm(hmms, bit_thresh)
    seq2hmm = best_model(seq2hmm)
    group2hmm = {} # group2hmm[seq][group] = [model, strand, coordinates, matches, gaps]
    for seq, info in list(seq2hmm.items()):
        group2hmm[seq] = {}
        # info = [model, [[hit1], [hit2], ...]]
        for group_num, group in enumerate(hit_groups(info[1])):
            # group is a group of hits to a single 16S gene
            # determine matching strand based on best hit
            best = sorted(group, reverse = True, key = itemgetter(-1))[0]
            strand = best[5]
            coordinates = [i[0] for i in group] + [i[1] for i in group]
            coordinates = [min(coordinates), max(coordinates), strand]
            # make sure all hits are to the same strand
            matches = [i for i in group if i[5] == strand]
            # gaps = [[gstart, gend], [gstart2, gend2]]
            gaps = check_gaps(matches)
            group2hmm[seq][group_num] = [info[0], strand, coordinates, matches, gaps]
    return group2hmm
85137
def load_just_in_time_genome_alignment(path, ref_spec, extensions=None,
                                       index_exts=None, fail_no_index=True,
                                       verbose=False):
    """"""
    if index_exts is None and fail_no_index:
      raise ValueError("Failure on no index specified for loading genome " +
                       "alignment, but no index extensions specified")

    extensions = __trim_extensions_dot(extensions)
    index_exts = __trim_extensions_dot(index_exts)

    partial_chrom_files = {}
    whole_chrom_files = {}
    for fn in os.listdir(path):
      pth = os.path.join(path, fn)
      if os.path.isfile(pth):
        base, ext = os.path.splitext(pth)
        ext = __trim_extension_dot(ext)
        if extensions is None or ext in extensions:
          idx_path = __find_index(pth, index_exts)
          if idx_path is None and fail_no_index:
            raise PyokitIOError("No index file for " + fn)
          chrom, start, end = __split_genomic_interval_filename(fn)
          assert((start is None and end is None) or
                 (start is not None and end is not None))
          if start is None:
            if chrom in whole_chrom_files:
              raise PyokitIOError("multiple files for chrom " + chrom)
            whole_chrom_files[chrom] = (pth, idx_path)
          else:
            k = (chrom, start, end)
            if k in partial_chrom_files:
              other_pth = partial_chrom_files[k][0]
              raise PyokitIOError("multiple files for " + str(k) + " --> " +
                                  pth + " and " + other_pth)
            partial_chrom_files[k] = (pth, idx_path)

    def factory(k):
      pth, idx = k
      return build_genome_alignment_from_file(pth, ref_spec, idx, verbose)

    return JustInTimeGenomeAlignment(whole_chrom_files, partial_chrom_files,
                                     factory)
1013855
def get_sequences(options, fasta_file_contents):
    ''' 
    '''
    errors = []
    fasta_files_str = ", ".join(fasta_file_contents.keys())
    fasta_records = None
    reverse_mapping = {}

    try:
        fasta_records, reverse_mapping = parse_FASTA_files(options, fasta_file_contents)
        if not fasta_records:
            errors.append("No protein sequences found in the FASTA file(s) %s." % fasta_files_str)
    except Exception, e:
        e = '\n'.join([l for l in traceback.format_exc(), str('e') if l.strip()])
        errors.append("Error parsing FASTA file(s) %s:\n%s" % (fasta_files_str, str(e)))

    if not fasta_records:
        return None, {}, errors

    colorprinter.message('Found %d protein sequence(s).' % len(fasta_records))
    return fasta_records, reverse_mapping, errors
191635
def find_databases(databases):
    
    # 16 ribosomal proteins in their expected order
    proteins = ['L15', 'L18', 'L6', 'S8', 'L5', 'L24', 'L14',
            'S17', 'L16', 'S3', 'L22', 'S19', 'L2', 'L4', 'L3', 'S10']
    # curated databases
    protein_databases = {
                'L14': 'rpL14_JGI_MDM.filtered.faa',
                'L15': 'rpL15_JGI_MDM.filtered.faa',
                'L16': 'rpL16_JGI_MDM.filtered.faa',
                'L18': 'rpL18_JGI_MDM.filtered.faa',
                'L22': 'rpL22_JGI_MDM.filtered.faa',
                'L24': 'rpL24_JGI_MDM.filtered.faa',
                'L2': 'rpL2_JGI_MDM.filtered.faa',
                'L3': 'rpL3_JGI_MDM.filtered.faa',
                'L4': 'rpL4_JGI_MDM.filtered.faa',
                'L5': 'rpL5_JGI_MDM.filtered.faa',
                'L6': 'rpL6_JGI_MDM.filtered.faa',
                'S10': 'rpS10_JGI_MDM.filtered.faa',
                'S17': 'rpS17_JGI_MDM.filtered.faa',
                'S19': 'rpS19_JGI_MDM.filtered.faa',
                'S3': 'rpS3_JGI_MDM.filtered.faa',
                'S8': 'rpS8_JGI_MDM.filtered.faa'}
    protein_databases = {key: '%s/%s' % (databases, database) \
            for key, database in list(protein_databases.items())}
    return proteins, protein_databases
218593
def get_column(self, chrom, position,
                 missing_seqs=MissingSequenceHandler.TREAT_AS_ALL_GAPS,
                 species=None):
    """"""
    blocks = self.get_blocks(chrom, position, position + 1)
    if len(blocks) == 0:
      raise NoSuchAlignmentColumnError("Request for column on chrom " +
                                       chrom + " at position " +
                                       str(position) + " not possible; " +
                                       "genome alignment not defined at " +
                                       "that locus.")
    if len(blocks) > 1:
      raise NoUniqueColumnError("Request for column on chrom " + chrom +
                                " at position " + str(position) + "not " +
                                "possible; ambiguous alignment of that locus.")

    return blocks[0].get_column_absolute(position, missing_seqs, species)
73950
def _seqfeature_to_coral(feature):
    '''

    '''
    # Some genomic sequences don't have a label attribute
    # TODO: handle genomic cases differently than others. Some features lack
    # a label but should still be incorporated somehow.
    qualifiers = feature.qualifiers
    if 'label' in qualifiers:
        feature_name = qualifiers['label'][0]
    elif 'locus_tag' in qualifiers:
        feature_name = qualifiers['locus_tag'][0]
    else:
        raise FeatureNameError('Unrecognized feature name')
    # Features with gaps are special, require looking at subfeatures
    # Assumption: subfeatures are never more than one level deep
    if feature.location_operator == 'join':
        # Feature has gaps. Have to figure out start/stop from subfeatures,
        # calculate gap indices. A nested feature model may be required
        # eventually.
        # Reorder the sub_feature list by start location
        # Assumption: none of the subfeatures overlap so the last entry in
        # the reordered list also has the final stop point of the feature.
        # FIXME: Getting a deprecation warning about using sub_features
        # instead of feature.location being a CompoundFeatureLocation
        reordered = sorted(feature.location.parts,
                           key=lambda location: location.start)
        starts = [int(location.start) for location in reordered]
        stops = [int(location.end) for location in reordered]
        feature_start = starts.pop(0)
        feature_stop = stops.pop(-1)
        starts = [start - feature_start for start in starts]
        stops = [stop - feature_start for stop in stops]
        feature_gaps = list(zip(stops, starts))
    else:
        # Feature doesn't have gaps. Ignore subfeatures.
        feature_start = int(feature.location.start)
        feature_stop = int(feature.location.end)
        feature_gaps = []
    feature_type = _process_feature_type(feature.type)
    if feature.location.strand == -1:
        feature_strand = 1
    else:
        feature_strand = 0
    if 'gene' in qualifiers:
        gene = qualifiers['gene']
    else:
        gene = []
    if 'locus_tag' in qualifiers:
        locus_tag = qualifiers['locus_tag']
    else:
        locus_tag = []
    coral_feature = coral.Feature(feature_name, feature_start,
                                  feature_stop, feature_type,
                                  gene=gene, locus_tag=locus_tag,
                                  qualifiers=qualifiers,
                                  strand=feature_strand,
                                  gaps=feature_gaps)
    return coral_feature
248485
def strainer(sequencepath):
    
    metadata_list = list()
    assert os.path.isdir(sequencepath), 'Cannot locate sequence path as specified: {}' \
        .format(sequencepath)
    # Get the sequences in the sequences folder into a list. Note that they must have a file extension that
    # begins with .fa
    strains = sorted(glob.glob(os.path.join(sequencepath, '*.fa*')))
    # Populate the metadata object. This object will be populated to mirror the objects created in the
    # genome assembly pipeline. This way this script will be able to be used as a stand-alone, or as part
    # of a pipeline
    assert strains, 'Could not find any files with an extension starting with "fa" in {}. Please check ' \
                    'to ensure that your sequence path is correct'.format(sequencepath)
    for sample in strains:
        # Create the object
        metadata = MetadataObject()
        # Set the base file name of the sequence. Just remove the file extension
        filename = os.path.splitext(os.path.split(sample)[1])[0]
        # Set the .name attribute to be the file name
        metadata.name = filename
        # Create the .general attribute
        metadata.general = GenObject()
        metadata.commands = GenObject()
        metadata.general.outputdirectory = os.path.join(sequencepath, filename)
        # Set the .general.bestassembly file to be the name and path of the sequence file
        metadata.general.bestassemblyfile = os.path.join(metadata.general.outputdirectory, '{sn}.fasta'
                                                         .format(sn=filename))
        make_path(metadata.general.outputdirectory)
        # Create a symlink to the directory
        relative_symlink(sample,
                         metadata.general.outputdirectory)
        metadata.general.logout = os.path.join(metadata.general.outputdirectory, 'out')
        metadata.general.logerr = os.path.join(metadata.general.outputdirectory, 'err')
        # Append the metadata for each sample to the list of samples
        metadata_list.append(metadata)
    return strains, metadata_list
1097873
def sort_length(source_file, source_file_type, direction=1):
    
    direction_text = 'ascending' if direction == 1 else 'descending'

    logging.info('Indexing sequences by length: %s', direction_text)

    # Adapted from the Biopython tutorial example.

    # Get the lengths and ids, and sort on length
    len_and_ids = sorted((len(rec), rec.id)
                         for rec in SeqIO.parse(source_file, source_file_type))

    if direction == 0:
        ids = reversed([seq_id for (length, seq_id) in len_and_ids])
    else:
        ids = [seq_id for (length, seq_id) in len_and_ids]
    del len_and_ids  # free this memory

    # SeqIO.index does not handle gzip instances
    if isinstance(source_file, gzip.GzipFile):
        tmpfile = tempfile.NamedTemporaryFile()
        source_file.seek(0)
        tmpfile.write(source_file.read())
        tmpfile.seek(0)
        source_file = tmpfile

    record_index = SeqIO.index(source_file.name, source_file_type)

    for seq_id in ids:
        yield record_index[seq_id]
4955
def copies(mapping, s2bins, rna, min_rna = 800, mismatches = 0):
    
    cov = {} # cov[scaffold] = [bases, length]
    s2bins, bins2s = parse_s2bins(s2bins)
    rna_cov = parse_rna(rna, s2bins, min_rna)
    s2bins, bins2s = filter_missing_rna(s2bins, bins2s, rna_cov)
    # count bases mapped to scaffolds and rRNA gene regions
    for line in mapping:
        line = line.strip().split()
        # get scaffold lengths
        if line[0].startswith('@'):
            if line[0].startswith('@SQ') is False:
                continue
            s = line[1].split(':')[1]
            l = int(line[2].split(':')[1])
            # check if scaffold is binned
            if s not in s2bins:
                continue
            if s not in cov:
                cov[s] = [0, l]
        # check mismatch threshold
        mm = count_mismatches(line)
        if mm is False or mm > mismatches:
            continue
        # check that scaffold is in bin
        s, bases = line[2], len(line[9])
        if s not in cov:
            continue
        cov[s][0] += bases 
        rna_cov = rna_bases(rna_cov, s, bases, line) 
    print('# mismatches threshold: %s' % (mismatches))
    header = ['#rRNA scaffold', 'rRNA genes >=%sbp on scaffold' % (min_rna), \
            'rRNA coverage', \
            'bin', 'bin info', 'bin coverage', \
            'rRNAs >=%sbp in bin' % (min_rna), \
            'rRNA coverage/bin coverage', \
            'estimated number of copies']
    print('\t'.join(header))
    for bin, scaffolds in list(bins2s.items()):
        rna_count = sum([len(rna_cov[s][2]) for s in scaffolds if s in rna_cov])
        for s in scaffolds:
            if s not in rna_cov:
                continue
            out = []
            counts = rna_cov[s]
            bin_cov = calc_bin_cov(bins2s[bin], cov)
            num_genes = len(counts[2])
            rna_coverage = float(float(counts[0])/float(counts[1]))
            if bin_cov == 0:
                rna_div_bin = 0
            else:
                rna_div_bin = float(rna_coverage/bin_cov)
            est = int(max([rna_count, counts, rna_div_bin]))
            out = [s, num_genes, rna_coverage, bin, bin_cov, rna_count, rna_div_bin, est]
            print('\t'.join([str(i) for i in out]))
438180
def get_coords(variant):
    
    coordinates = {
        'chrom': None,
        'end_chrom': None,
        'sv_length': None,
        'sv_type': None,
        'pos': None,
        'end': None,
    }
    chrom = variant.CHROM
    if chrom.startswith(('chr', 'CHR', 'Chr')):
        chrom = chrom[3:]
    coordinates['chrom'] = chrom
    end_chrom = chrom
    
    pos = int(variant.POS)
    alt = variant.ALT[0]

    # Get the end position
    # This will be None for non-svs
    end_pos = variant.INFO.get('END')
    if end_pos:
        end = int(end_pos)
    else:
        end = int(variant.end)
    coordinates['end'] = end
    
    sv_type = variant.INFO.get('SVTYPE')
    length = variant.INFO.get('SVLEN')
    if length:
        sv_len = abs(length)
    else:
        sv_len = end - pos

    # Translocations will sometimes have a end chrom that differs from chrom
    if sv_type == 'BND':
        other_coordinates = alt.strip('ACGTN[]').split(':')
        end_chrom = other_coordinates[0]
        if end_chrom.startswith(('chr', 'CHR', 'Chr')):
            end_chrom = end_chrom[3:]

        end = int(other_coordinates[1])

        #Set 'infinity' to length if translocation
        sv_len = float('inf')

    # Insertions often have length 0 in VCF
    if (sv_len == 0 and alt != '<INS>'):
        sv_len = len(alt)

    if (pos == end) and (sv_len > 0):
        end = pos + sv_len

    position = Position(chrom, pos)
    end_position = Position(end_chrom, end)
    
    # If 'start' is greater than 'end', switch positions
    if is_greater(position, end_position):
        end_chrom = position.chrom
        end = position.pos
        
        chrom = end_position.chrom
        pos = end_position.pos
    
    coordinates['end_chrom'] = end_chrom
    coordinates['pos'] = pos
    coordinates['end'] = end
    coordinates['sv_length'] = sv_len
    coordinates['sv_type'] = sv_type
    
    return coordinates
599145
def get_sequences_from_cluster(c1, c2, data):
    """"""
    seqs1 = data[c1]['seqs']
    seqs2 = data[c2]['seqs']
    seqs = list(set(seqs1 + seqs2))
    names = []
    for s in seqs:
        if s in seqs1 and s in seqs2:
            names.append("both")
        elif s in seqs1:
            names.append(c1)
        else:
            names.append(c2)
    return seqs, names
293558
def parse_rna(rna, s2bins, min_rna):
    
    rna_cov = {}
    for seq in parse_fasta(rna):
        # check that length passes threshold
        length = len(seq[1])
        if length < min_rna:
            continue
        # check if sequence is binnned
        s = seq[0].split('>')[1].split()[0]
        if s not in s2bins:
            continue
        if s not in rna_cov:
            rna_cov[s] = [0, 0, []]
        position = [int(i) for i in seq[0].rsplit('pos=', 1)[1].split()[0].split('-')]
        rna_cov[s][2].append(position)
        rna_cov[s][1] += length
    return rna_cov
163510
def extract_ngram_filter(pos_seq, regex=SimpleNP, minlen=1, maxlen=8):
	""""""
	ss = coarse_tag_str(pos_seq)

	def gen():
		for s in xrange(len(ss)):
			for n in xrange(minlen, 1 + min(maxlen, len(ss) - s)):
				e = s + n
				substr = ss[s:e]
				if re.match(regex + "$", substr):
					yield (s, e)

	return list(gen())
1000213
def cdsparse(self, record):
        
        try:
            # Find genes that are present in all strains of interest - the number of times the gene is found is
            # equal to the number of strains. Earlier parsing ensures that the same gene is not present in a strain
            # more than once
            if self.genes[self.genenames[record.id]] == len(self.runmetadata.samples):
                # Add the gene names and sequences to the appropriate dictionaries
                try:
                    self.genesequence[self.genenames[record.id]].add(str(record.seq))
                # Initialise the dictionary as required, then populate as above
                except KeyError:
                    self.genesequence[self.genenames[record.id]] = set()
                    self.genesequence[self.genenames[record.id]].add(str(record.seq))
                try:
                    self.coresequence[str(record.seq)].add(record.id)
                except KeyError:
                    self.coresequence[str(record.seq)] = set()
                    self.coresequence[str(record.seq)].add(record.id)
        except KeyError:
            pass
346145
def _get_read_depths(cls, read_depths_file, sequence_name, position):
        ''''''
        assert os.path.exists(read_depths_file)
        assert os.path.exists(read_depths_file + '.tbi')
        tbx = pysam.TabixFile(read_depths_file)
        try:
            rows = [x for x in tbx.fetch(sequence_name, position, position + 1)]
        except:
            return None

        if len(rows) > 1: # which happens with indels, mutiple lines for same base of reference
            test_rows = [x for x in rows if x.rstrip().split()[3] != '.']
            if len(test_rows) != 1:
                rows = [rows[-1]]
            else:
                rows = test_rows

        if len(rows) == 1:
            r, p, ref_base, alt_base, ref_counts, alt_counts = rows[0].rstrip().split()
            bases = ref_base if alt_base == '.' else ref_base + ',' + alt_base
            return bases, int(ref_counts), alt_counts
        else:
            return None
384834
def dna_to_re(seq):
    

    seq = seq.replace('K', '[GT]')
    seq = seq.replace('M', '[AC]')
    seq = seq.replace('R', '[AG]')
    seq = seq.replace('Y', '[CT]')
    seq = seq.replace('S', '[CG]')
    seq = seq.replace('W', '[AT]')
    seq = seq.replace('B', '[CGT]')
    seq = seq.replace('V', '[ACG]')
    seq = seq.replace('H', '[ACT]')
    seq = seq.replace('D', '[AGT]')
    seq = seq.replace('X', '[GATC]')
    seq = seq.replace('N', '[GATC]')

    return re.compile(seq)
1026293
def load_sequences_to_strains(self, joblib=False, cores=1, force_rerun=False):
        """"""
        log.info('Loading sequences to strain GEM-PROs...')
        if joblib:
            result = DictList(Parallel(n_jobs=cores)(delayed(self._load_sequences_to_strain)(s, force_rerun) for s in self.strain_ids))
        else:
            result = []
            for s in tqdm(self.strain_ids):
                result.append(self._load_sequences_to_strain(s, force_rerun))

        for strain_id, gp_seqs_path in result:
            self.strain_infodict[strain_id]['gp_seqs_path'] = gp_seqs_path