intersect transcription position with annotation files
4588
def anncluster(c, clus_obj, db, type_ann, feature_id="name"):
    """"""
    id_sa, id_ea, id_id, id_idl, id_sta = 1, 2, 3, 4, 5
    if type_ann == "bed":
        id_sb = 7
        id_eb = 8
        id_stb = 11
        id_tag = 9
    ida = 0
    clus_id = clus_obj.clus
    loci_id = clus_obj.loci
    db = os.path.splitext(db)[0]
    logger.debug("Type:%s\n" % type_ann)
    for cols in c.features():
        if type_ann == "gtf":
            cb, sb, eb, stb, db, tag = read_gtf_line(cols[6:], feature_id)
        else:
            sb = int(cols[id_sb])
            eb = int(cols[id_eb])
            stb = cols[id_stb]
            tag = cols[id_tag]
        id = int(cols[id_id])
        idl = int(cols[id_idl])
        if (id in clus_id):
            clus = clus_id[id]
            sa = int(cols[id_sa])
            ea = int(cols[id_ea])
            ida += 1
            lento5, lento3, strd = _position_in_feature([sa, ea, cols[id_sta]], [sb, eb, stb])
            if db in loci_id[idl].db_ann:
                ann = annotation(db, tag, strd, lento5, lento3)
                tdb = loci_id[idl].db_ann[db]
                tdb.add_db_ann(ida, ann)
                loci_id[idl].add_db(db, tdb)
            else:
                ann = annotation(db, tag, strd, lento5, lento3)
                tdb = dbannotation(1)
                tdb.add_db_ann(ida, ann)
                loci_id[idl].add_db(db, tdb)
            clus_id[id] = clus
    clus_obj.clus = clus_id
    clus_obj.loci = loci_id
    return clus_obj
699016
def download_annotation(genomebuild, gene_file): 
    
    pred_bin = "genePredToBed"
    pred = find_executable(pred_bin)
    if not pred:
        sys.stderr.write("{} not found in path!\n".format(pred_bin))
        sys.exit(1)

    tmp = NamedTemporaryFile(delete=False, suffix=".gz")

    anno = []
    f = urlopen(UCSC_GENE_URL.format(genomebuild))
    p = re.compile(r'\w+.Gene.txt.gz')
    for line in f.readlines():
        m = p.search(line.decode())
        if m:
            anno.append(m.group(0))

    sys.stderr.write("Retrieving gene annotation for {}\n".format(genomebuild))
    url = ""
    for a in ANNOS:
        if a in anno:
            url = UCSC_GENE_URL.format(genomebuild) + a
            break
    if url:
        sys.stderr.write("Using {}\n".format(url))
        urlretrieve(
                url,
                tmp.name
                )
         
        with gzip.open(tmp.name) as f:
            cols = f.readline().decode(errors='ignore').split("\t")

        start_col = 1
        for i,col in enumerate(cols):
            if col == "+" or col == "-":
                start_col = i - 1
                break
        end_col = start_col + 10
       
        cmd = "zcat {} | cut -f{}-{} | {} /dev/stdin {}"
        print(cmd.format(tmp.name, start_col, end_col, pred, gene_file))
        sp.call(cmd.format(
            tmp.name, start_col, end_col, pred, gene_file), 
            shell=True)

    else:
        sys.stderr.write("No annotation found!")
1056032
def coordinate_documents(docs)
      regex = document_url_regex
      approved = {}
      docs.each do |doc|
        lang = doc.data['lang'] || @default_lang
        url = doc.url.gsub(regex, '/')
        doc.data['permalink'] = url
        next if @file_langs[url] == @active_lang
        next if @file_langs[url] == @default_lang && lang != @active_lang
        approved[url] = doc
        @file_langs[url] = lang
      end
      approved.values
    end
1015902
def cli(ctx, exon_a, exon_b, organism="", sequence=""):
    
    return ctx.gi.annotations.merge_exons(exon_a, exon_b, organism=organism, sequence=sequence)
641892
def _get_region(self, contig, start, end, fout, min_length=250):
        ''''''
        sam_reader = pysam.Samfile(self.bam, "rb")
        trimming_end = (start == 0)
        for read in sam_reader.fetch(contig, start, end):
            read_interval = pyfastaq.intervals.Interval(read.pos, read.reference_end - 1)
            seq = mapping.aligned_read_to_read(read, ignore_quality=not self.fastq_out, revcomp=False)

            if trimming_end:
                bases_off_start = 0
                bases_off_end = max(0, read.reference_end - 1 - end)
                #seq.seq = seq.seq[:read.query_alignment_end - bases_off_end]
                seq = seq.subseq(0, read.query_alignment_end - bases_off_end)
            else:
                bases_off_start = max(0, start - read.pos + 1)
                #seq.seq = seq.seq[bases_off_start  + read.query_alignment_start:]
                seq = seq.subseq(bases_off_start  + read.query_alignment_start, len(seq))

            if read.is_reverse:
                seq.revcomp()

            if len(seq) >= min_length:
                print(seq, file=fout)
991145
def bamtag(sam):
    ''' 
    '''
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    if is_python3():
        queryalignment = next(track)
    else:
        queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            aln.tags += [('XC', match.group('CB'))]
        if "molecular" in annotations:
            aln.tags += [('RX', match.group('MB'))]
        if "sample" in annotations:
            aln.tags += [('XS', match.group('SB'))]

        out_file.write(aln)

    total_time = time.time() - start_time
    logger.info('BAM tag conversion done - {:.3}s, {:,} alns/min'.format(total_time, int(60. * count / total_time)))
    logger.info("Processed %d alignments." % count)
500016
def construct_transformed_regex(annotations):
    
    re_string = '.*'
    if "cellular" in annotations:
        re_string += ":CELL_(?P<CB>.*)"
    if "molecular" in annotations:
        re_string += ":UMI_(?P<MB>\w*)"
    if "sample" in annotations:
        re_string += ":SAMPLE_(?P<SB>\w*)"
    if re_string == ".*":
        logger.error("No annotation present on this file, aborting.")
        sys.exit(1)
    return re_string
690464
def _reverse_annotations(old_record, new_record):
    
    # Copy the annotations over
    for k, v in list(old_record.annotations.items()):
        # Trim if appropriate
        if isinstance(v, (tuple, list)) and len(v) == len(old_record):
            assert len(v) == len(old_record)
            v = v[::-1]
        new_record.annotations[k] = v

    # Letter annotations must be lists / tuples / strings of the same
    # length as the sequence
    for k, v in list(old_record.letter_annotations.items()):
        assert len(v) == len(old_record)
        new_record.letter_annotations[k] = v[::-1]
991154
def bam2es(
        bam_fn,
        es_fo,
        allowed_delta,
    ):
        

        es_fo.write("# RN:   read name" + os.linesep)
        es_fo.write("# Q:    is mapped with quality" + os.linesep)
        es_fo.write("# Chr:  chr id" + os.linesep)
        es_fo.write("# D:    direction" + os.linesep)
        es_fo.write("# L:    leftmost nucleotide" + os.linesep)
        es_fo.write("# R:    rightmost nucleotide" + os.linesep)
        es_fo.write("# Cat:  category of alignment assigned by LAVEnder" + os.linesep)
        es_fo.write("#         M_i    i-th segment is correctly mapped" + os.linesep)
        es_fo.write("#         m      segment should be unmapped but it is mapped" + os.linesep)
        es_fo.write("#         w      segment is mapped to a wrong location" + os.linesep)
        es_fo.write("#         U      segment is unmapped and should be unmapped" + os.linesep)
        es_fo.write("#         u      segment is unmapped and should be mapped" + os.linesep)
        es_fo.write("# Segs: number of segments" + os.linesep)
        es_fo.write("# " + os.linesep)
        es_fo.write("# RN\tQ\tChr\tD\tL\tR\tCat\tSegs" + os.linesep)

        with pysam.AlignmentFile(bam_fn, "rb") as sam:
            references_dict = {}

            for i in range(len(sam.references)):
                references_dict[sam.references[i]] = i + 1

            for read in sam:
                rnf_read_tuple = rnftools.rnfformat.ReadTuple()
                rnf_read_tuple.destringize(read.query_name)

                left = read.reference_start + 1
                right = read.reference_end
                chrom_id = references_dict[sam.references[read.reference_id]]

                nb_of_segments = len(rnf_read_tuple.segments)

                if rnf_read_tuple.segments[0].genome_id == 1:
                    should_be_mapped = True
                else:
                    should_be_mapped = False

                # read is unmapped
                if read.is_unmapped:
                    # read should be mapped
                    if should_be_mapped:
                        category = "u"
                    # read should be unmapped
                    else:
                        category = "U"
                # read is mapped
                else:
                    # read should be mapped
                    if should_be_mapped:
                        exists_corresponding_segment = False

                        for j in range(len(rnf_read_tuple.segments)):
                            segment = rnf_read_tuple.segments[j]
                            if (
                                (segment.left == 0 or abs(segment.left - left) <= allowed_delta)
                                and (segment.right == 0 or abs(segment.right - right) <= allowed_delta)
                                and (segment.left != 0 or segment.right == 0)
                                and (chrom_id == 0 or chrom_id == segment.chr_id)
                            ):
                                exists_corresponding_segment = True
                                segment = str(j + 1)
                                break

                        # read was mapped to correct location
                        if exists_corresponding_segment:  # exists ok location?
                            category = "M_" + segment
                        # read was mapped to incorrect location
                        else:
                            category = "w"
                    # read should be unmapped
                    else:
                        category = "m"

                es_fo.write(
                    "\t".join(
                        map(
                            str,
                            [
                                # read name
                                read.query_name,
                                # aligned?
                                "unmapped" if read.is_unmapped else "mapped_" + str(read.mapping_quality),
                                # reference id
                                chrom_id,
                                # direction
                                "R" if read.is_reverse else "F",
                                # left
                                left,
                                # right
                                right,
                                # assigned category
                                category,
                                # count of segments
                                nb_of_segments
                            ]
                        )
                    ) + os.linesep
                )
1052733
def generate_intersect_subparser(subparsers):
    """"""
    parser = subparsers.add_parser(
        'intersect', description=constants.INTERSECT_DESCRIPTION,
        epilog=constants.INTERSECT_EPILOG, formatter_class=ParagraphFormatter,
        help=constants.INTERSECT_HELP)
    parser.set_defaults(func=ngram_intersection)
    utils.add_common_arguments(parser)
    utils.add_db_arguments(parser)
    utils.add_corpus_arguments(parser)
    utils.add_query_arguments(parser)
26442
def get_annotations(fname, prefix=None):
    ""
    annot_dict = json.load(open(fname))
    id2images, id2bboxes, id2cats = {}, collections.defaultdict(list), collections.defaultdict(list)
    classes = {}
    for o in annot_dict['categories']:
        classes[o['id']] = o['name']
    for o in annot_dict['annotations']:
        bb = o['bbox']
        id2bboxes[o['image_id']].append([bb[1],bb[0], bb[3]+bb[1], bb[2]+bb[0]])
        id2cats[o['image_id']].append(classes[o['category_id']])
    for o in annot_dict['images']:
        if o['id'] in id2bboxes:
            id2images[o['id']] = ifnone(prefix, '') + o['file_name']
    ids = list(id2images.keys())
    return [id2images[k] for k in ids], [[id2bboxes[k], id2cats[k]] for k in ids]
247026
def _codon_positions(self, feature):
        
        results = self._transcript_feature_positions(feature)
        if len(results) != 3:
            raise ValueError(
                "Expected 3 positions for %s of %s but got %d" % (
                    feature,
                    self.id,
                    len(results)))
        return results
517588
function transferAnnotations (doc, path, offset, newPath, newOffset) {
  var index = doc.getIndex('annotations')
  var annotations = index.get(path, offset)
  for (let i = 0; i < annotations.length; i++) {
    let a = annotations[i]
    var isInside = (offset > a.start.offset && offset < a.end.offset)
    var start = a.start.offset
    var end = a.end.offset
    // 1. if the cursor is inside an annotation it gets either split or truncated
    if (isInside) {
      // create a new annotation if the annotation is splittable
      if (a.canSplit()) {
        let newAnno = a.toJSON()
        newAnno.id = uuid(a.type + '_')
        newAnno.start.path = newPath
        newAnno.start.offset = newOffset
        newAnno.end.path = newPath
        newAnno.end.offset = newOffset + a.end.offset - offset
        doc.create(newAnno)
      }
      // in either cases truncate the first part
      let newStartOffset = a.start.offset
      let newEndOffset = offset
      // if after truncate the anno is empty, delete it
      if (newEndOffset === newStartOffset) {
        doc.delete(a.id)
      // ... otherwise update the range
      } else {
        // TODO: Use coordintate ops!
        if (newStartOffset !== start) {
          doc.set([a.id, 'start', 'offset'], newStartOffset)
        }
        if (newEndOffset !== end) {
          doc.set([a.id, 'end', 'offset'], newEndOffset)
        }
      }
    // 2. if the cursor is before an annotation then simply transfer the annotation to the new node
    } else if (a.start.offset >= offset) {
      // TODO: Use coordintate ops!
      // Note: we are preserving the annotation so that anything which is connected to the annotation
      // remains valid.
      doc.set([a.id, 'start', 'path'], newPath)
      doc.set([a.id, 'start', 'offset'], newOffset + a.start.offset - offset)
      doc.set([a.id, 'end', 'path'], newPath)
      doc.set([a.id, 'end', 'offset'], newOffset + a.end.offset - offset)
    }
  }

  // TODO: fix support for container annotations
  // // same for container annotation anchors
  // index = doc.getIndex('container-annotation-anchors');
  // var anchors = index.get(path);
  // var containerAnnoIds = [];
  // forEach(anchors, function(anchor) {
  //   containerAnnoIds.push(anchor.id);
  //   var start = anchor.offset;
  //   if (offset <= start) {
  //     // TODO: Use coordintate ops!
  //     let coor = anchor.isStart?'start':'end'
  //     doc.set([anchor.id, coor, 'path'], newPath);
  //     doc.set([anchor.id, coor, 'offset'], newOffset + anchor.offset - offset);
  //   }
  // });
  // // check all anchors after that if they have collapsed and remove the annotation in that case
  // forEach(uniq(containerAnnoIds), function(id) {
  //   var anno = doc.get(id);
  //   var annoSel = anno.getSelection();
  //   if(annoSel.isCollapsed()) {
  //     // console.log("...deleting container annotation because it has collapsed" + id);
  //     doc.delete(id);
  //   }
  // });
}
213923
def dump_annotations(type_info, files):
    
    with open(type_info) as f:
        data = json.load(f)
    for item in data:
        path, line, func_name = item['path'], item['line'], item['func_name']
        if files and path not in files:
            for f in files:
                if path.startswith(os.path.join(f, '')):
                    break
            else:
                continue  # Outer loop
        print("%s:%d: in %s:" % (path, line, func_name))
        type_comments = item['type_comments']
        signature = unify_type_comments(type_comments)
        arg_types = signature['arg_types']
        return_type = signature['return_type']
        print("    # type: (%s) -> %s" % (", ".join(arg_types), return_type))
790320
def subset_bamfile(sam, barcodes):
    
    from pysam import AlignmentFile

    start_time = time.time()

    sam_file = open_bamfile(sam)
    out_file = AlignmentFile("-", "wh", template=sam_file)
    track = sam_file.fetch(until_eof=True)

    # peek at first alignment to determine the annotations
    queryalignment = track.next()
    annotations = detect_alignment_annotations(queryalignment)
    track = itertools.chain([queryalignment], track)

    re_string = construct_transformed_regex(annotations)
    parser_re = re.compile(re_string)
    barcodes = set(barcode.strip() for barcode in barcodes)

    for count, aln in enumerate(track, start=1):
        if count and not count % 1000000:
            logger.info("Processed %d alignments." % count)

        match = parser_re.match(aln.qname)
        tags = aln.tags

        if "cellular" in annotations:
            cb = match.group('CB')
            if cb in barcodes:
                out_file.write(aln)
1062372
def parse_fasta_annotations(fastas, annot_tables, trans_table):
    
    if annot_tables is not False:
        annots = {}
        for table in annot_tables:
            for cds in open(table):
                ID, start, end, strand = cds.strip().split()
                annots[ID] = [start, end, int(strand)]
    for fasta in fastas:
        for seq in parse_fasta(fasta):
            if ('# ;gc_cont' not in seq[0] and '# ID=' not in seq[0]) and annot_tables is False:
                print('# specify fasta from Prodigal or annotations table (-t)', file=sys.stderr)
                exit()
            if 'ID=' in seq[0]:
                ID = seq[0].rsplit('ID=', 1)[1].split(';', 1)[0]
                contig = seq[0].split()[0].split('>')[1].rsplit('_%s' % (ID), 1)[0]
            else:
                contig = seq[0].split()[0].split('>')[1].rsplit('_', 1)[0]
            locus = seq[0].split()[0].split('>')[1]
            # annotation info from Prodigal
            if ('# ;gc_cont' in seq[0] or '# ID=' in seq[0]):
                info = seq[0].split(' # ')
                start, end, strand = int(info[1]), int(info[2]), info[3]
                if strand == '1':
                    strand = 1
                else:
                    strand = -1
                product = [''.join(info[4].split()[1:])]
            # annotation info from table
            else:
                start, end, strand = annots[locus]
                product = seq[0].split(' ', 1)[1]
            info = {'transl_table':[trans_table], \
                    'translation':[seq[1]], \
                    'product':product}
            yield contig, [locus, [start, end, strand], info]
141796
def align_two_alignments(aln1, aln2, moltype, params=None):
    
    #create SequenceCollection object from seqs
    aln1 = Alignment(aln1,MolType=moltype)
    #Create mapping between abbreviated IDs and full IDs
    aln1_int_map, aln1_int_keys = aln1.getIntMap()
    #Create SequenceCollection from int_map.
    aln1_int_map = Alignment(aln1_int_map,MolType=moltype)

    #create Alignment object from aln
    aln2 = Alignment(aln2,MolType=moltype)
    #Create mapping between abbreviated IDs and full IDs
    aln2_int_map, aln2_int_keys = aln2.getIntMap(prefix='seqn_')
    #Create SequenceCollection from int_map.
    aln2_int_map = Alignment(aln2_int_map,MolType=moltype)

    #Update aln1_int_keys with aln2_int_keys
    aln1_int_keys.update(aln2_int_keys)

    #Create Mafft app.
    app = Clustalw(InputHandler='_input_as_multiline_string',\
        params=params,
        SuppressStderr=True)
    app.Parameters['-align'].off()
    app.Parameters['-infile'].off()
    app.Parameters['-profile'].on()

    #Add aln_int_map as profile1
    app.Parameters['-profile1'].on(\
        app._tempfile_as_multiline_string(aln1_int_map.toFasta()))

    #Add seq_int_map as profile2
    app.Parameters['-profile2'].on(\
        app._tempfile_as_multiline_string(aln2_int_map.toFasta()))
    #Get results using int_map as input to app
    res = app()

    #Get alignment as dict out of results
    alignment = dict(ClustalParser(res['Align'].readlines()))

    #Make new dict mapping original IDs
    new_alignment = {}
    for k,v in alignment.items():
        new_alignment[aln1_int_keys[k]]=v
    #Create an Alignment object from alignment dict
    new_alignment = Alignment(new_alignment,MolType=moltype)
    #Clean up
    res.cleanUp()
    remove(app.Parameters['-profile1'].Value)
    remove(app.Parameters['-profile2'].Value)
    del(aln1,aln1_int_map,aln1_int_keys,\
        aln2,aln2_int_map,aln2_int_keys,app,res,alignment)

    return new_alignment