Annotate miRNAs using miRBase database with seqbuster tool
17183
def sample_annotation(data):
    
    names = data["rgnames"]['sample']
    tools = dd.get_expression_caller(data)
    work_dir = os.path.join(dd.get_work_dir(data), "mirbase")
    out_dir = os.path.join(work_dir, names)
    utils.safe_makedir(out_dir)
    out_file = op.join(out_dir, names)
    if dd.get_mirbase_hairpin(data):
        mirbase = op.abspath(op.dirname(dd.get_mirbase_hairpin(data)))
        if utils.file_exists(data["collapse"]):
            data['transcriptome_bam'] = _align(data["collapse"],
                                               dd.get_mirbase_hairpin(data),
                                               out_file,
                                               data)
            data['seqbuster'] = _miraligner(data["collapse"], out_file,
                                            dd.get_species(data),
                                            mirbase,
                                            data['config'])
            data["mirtop"] = _mirtop(data['seqbuster'],
                                     dd.get_species(data),
                                     mirbase,
                                     out_dir,
                                     data['config'])
        else:
            logger.debug("Trimmed collapsed file is empty for %s." % names)
    else:
        logger.debug("No annotation file from miRBase.")

    sps = dd.get_species(data) if dd.get_species(data) else "None"
    logger.debug("Looking for mirdeep2 database for %s" % names)
    if file_exists(op.join(dd.get_work_dir(data), "mirdeep2", "novel", "hairpin.fa")):
        data['seqbuster_novel'] = _miraligner(data["collapse"], "%s_novel" % out_file, sps,
                                              op.join(dd.get_work_dir(data),
                                                      "mirdeep2", "novel"),
                                              data['config'])

    if "trna" in tools:
        data['trna'] = _mint_trna_annotation(data)

    data = spikein.counts_spikein(data)
    return [[data]]
636978
def _cluster(bam_file, ma_file, out_dir, reference, annotation_file=None):
    
    seqcluster = op.join(get_bcbio_bin(), "seqcluster")
    # cl = ["cluster", "-o", out_dir, "-m", ma_file, "-a", bam_file, "-r", reference]
    if annotation_file:
        annotation_file = "-g " + annotation_file
    else:
        annotation_file = ""

    if not file_exists(op.join(out_dir, "counts.tsv")):
        cmd = ("{seqcluster} cluster -o {out_dir} -m {ma_file} -a {bam_file} -r {reference} {annotation_file}")
        do.run(cmd.format(**locals()), "Running seqcluster.")
    counts = op.join(out_dir, "counts.tsv")
    stats = op.join(out_dir, "read_stats.tsv")
    json = op.join(out_dir, "seqcluster.json")
    return {'out_dir': out_dir, 'count_file': counts, 'stat_file': stats, 'json': json}
1113790
def _coord(sequence, start, mirna, precursor, iso):
    
    dif = abs(mirna[0] - start)
    if start < mirna[0]:
        iso.t5 = sequence[:dif].upper()
    elif start > mirna[0]:
        iso.t5 = precursor[mirna[0] - 1:mirna[0] - 1 + dif].lower()
    elif start == mirna[0]:
        iso.t5 = "NA"
    if dif > 4:
        logger.debug("start > 3 %s %s %s %s %s" % (start, len(sequence), dif, mirna, iso.format()))
        return None

    end = start + (len(sequence) - len(iso.add)) - 1
    dif = abs(mirna[1] - end)
    if iso.add:
        sequence = sequence[:-len(iso.add)]
    # if dif > 3:
    #    return None
    if end > mirna[1]:
        iso.t3 = sequence[-dif:].upper()
    elif end < mirna[1]:
        iso.t3 = precursor[mirna[1] - dif:mirna[1]].lower()
    elif end == mirna[1]:
        iso.t3 = "NA"
    if dif > 4:
        logger.debug("end > 3 %s %s %s %s %s" % (len(sequence), end, dif, mirna, iso.format()))
        return None
    logger.debug("%s %s %s %s %s %s" % (start, len(sequence), end, dif, mirna, iso.format()))
    return True
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
                )
399374
def read_annotations(**kws):
    """"""
    if 'gaf' not in kws and 'gene2go' not in kws:
        return
    gene2gos = None
    if 'gaf' in kws:
        gene2gos = read_gaf(kws['gaf'], prt=sys.stdout)
        if not gene2gos:
            raise RuntimeError("NO ASSOCIATIONS LOADED FROM {F}".format(F=kws['gaf']))
    elif 'gene2go' in kws:
        assert 'taxid' in kws, 'taxid IS REQUIRED WHEN READING gene2go'
        gene2gos = read_ncbi_gene2go(kws['gene2go'], taxids=[kws['taxid']])
        if not gene2gos:
            raise RuntimeError("NO ASSOCIATIONS LOADED FROM {F} FOR TAXID({T})".format(
                F=kws['gene2go'], T=kws['taxid']))
    return gene2gos
124144
def seq_str(self):
        """"""

        if not self.seq:
            return None

        return ssbio.protein.sequence.utils.cast_to_str(self.seq)
706195
def seq_info(names, id2names, insertions, sequences):
    
    seqs = {} # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns]], ...]]
    for name in names:
        id = id2names[name]
        gene = name.split('fromHMM::', 1)[0].rsplit(' ', 1)[1]
        model = name.split('fromHMM::', 1)[1].split('=', 1)[1].split()[0]
        i_gene_pos = insertions[id] # coordinates of each insertion wrt gene
        i_model_pos = name.split('fromHMM::', 1)[1].split('model-pos(ins-len)=')[1].split()[0].split(';') # model overlap
        i_info = []
        for i, ins in enumerate(i_gene_pos):
            model_pos = i_model_pos[i].split('-')[1].split('(')[0]
            length = i_model_pos[i].split('(')[1].split(')')[0]
            iheader = '>%s_%s insertion::seq=%s type=insertion strand=n/a gene-pos=%s-%s model-pos=%s'\
                    % (id, (i + 1), (i + 1), ins[0], ins[1], model_pos)
            iseq = sequences[id][1][ins[0]:(ins[1] + 1)]
            iseq = [iheader, iseq]
            info = [ins, model_pos, length, iseq, [], []]
            i_info.append(info)
        seqs[id] = [gene, model, i_info]
    return seqs
31071
def from_seqfeature(s, **kwargs):
    
    source = s.qualifiers.get('source', '.')[0]
    score = s.qualifiers.get('score', '.')[0]
    seqid = s.qualifiers.get('seqid', '.')[0]
    frame = s.qualifiers.get('frame', '.')[0]
    strand = _feature_strand[s.strand]

    # BioPython parses 1-based GenBank positions into 0-based for use within
    # Python.  We need to convert back to 1-based GFF format here.
    start = s.location.start.position + 1
    stop = s.location.end.position
    featuretype = s.type
    id = s.id
    attributes = dict(s.qualifiers)
    attributes.pop('source', '.')
    attributes.pop('score', '.')
    attributes.pop('seqid', '.')
    attributes.pop('frame', '.')
    return Feature(seqid, source, featuretype, start, stop, score, strand,
                   frame, attributes, id=id, **kwargs)
1313512
def _setup_genome_annotations(g, args, ann_groups):
    
    available_anns = g.get("annotations", []) + g.pop("annotations_available", [])
    anns = []
    for orig_target in args.datatarget:
        if orig_target in ann_groups:
            targets = ann_groups[orig_target]
        else:
            targets = [orig_target]
        for target in targets:
            if target in available_anns:
                anns.append(target)
    g["annotations"] = anns
    if "variation" not in args.datatarget and "validation" in g:
        del g["validation"]
    return g
188179
def create_psm_lookup(fn, fastafn, mapfn, header, pgdb, unroll=False,
                      specfncol=None, decoy=False,
                      fastadelim=None, genefield=None):
    
    proteins = store_proteins_descriptions(pgdb, fastafn, fn, mapfn, header,
                                           decoy, fastadelim, genefield)
    mzmlmap = pgdb.get_mzmlfile_map()
    sequences = {}
    for psm in tsvreader.generate_tsv_psms(fn, header):
        seq = tsvreader.get_psm_sequence(psm, unroll)
        sequences[seq] = 1
    pgdb.store_pepseqs(((seq,) for seq in sequences))
    pepseqmap = pgdb.get_peptide_seq_map()
    psms = []
    for row, psm in enumerate(tsvreader.generate_tsv_psms(fn, header)):
        specfn, psm_id, scan, seq, score = tsvreader.get_psm(psm, unroll,
                                                             specfncol)
        if len(psms) % DB_STORE_CHUNK == 0:
            pgdb.store_psms(psms)
            psms = []
        psms.append({'rownr': row,
                     'psm_id': psm_id,
                     'seq': pepseqmap[seq],
                     'score': score,
                     'specfn': mzmlmap[specfn],
                     'scannr': scan,
                     'spec_id': '{}_{}'.format(mzmlmap[specfn], scan),
                     })
    pgdb.store_psms(psms)
    pgdb.index_psms()
    store_psm_protein_relations(fn, header, pgdb, proteins)
451196
def print_seqs(seqs, id2name, name):
    
    orfs = open('%s.orfs.faa' % (name.rsplit('.', 1)[0]), 'w')
    introns = open('%s.introns.fa' % (name.rsplit('.', 1)[0]), 'w')
    insertions = open('%s.insertions.fa' % (name.rsplit('.', 1)[0]), 'w')
    for seq in seqs:
        for i, ins in enumerate(seqs[seq][2], 1):
            model_pos = ins[1]
            if ins[6] is True: # orf(s) in ins[4]
                for orf in ins[4]:
                    orf_info = orf[0].split('>')[1].split()
                    id = orf_info[0].split('_', 1)[0]
                    name = id2name[id]
                    annotation = orf_info[1]
                    strand = orf_info[7]
                    if strand == '1':
                        strand = '+'
                    else:
                        strand = '-'
                    start, stop = sorted([int(orf_info[3]), int(orf_info[5])])
                    header = '>%s insertion::seq=%s type=%s strand=%s gene-pos=%s-%s model-pos=%s'\
                            % (name, i, annotation, strand, start, stop, model_pos)
                    print('\n'.join([header, orf[1]]), file=orfs)
            if ins[7] is True: # intron(s) in ins[5]
                for intron in ins[5]:
                    id, type, strand, pos = intron[0].split('>', 1)[1].split()
                    name = id2name[id.split('_')[0]]
                    header = '>%s insertion::seq=%s type=%s strand=%s gene-pos=%s model-pos=%s'\
                            % (name, i, type, strand, pos, model_pos)
                    print('\n'.join([header, intron[1]]), file=introns)
            insertion = ins[3]
            id, info = insertion[0].split('>')[1].split(' ', 1)
            name = id2name[id.split('_')[0]]
            header = '>%s %s' % (name, info)
            print('\n'.join([header, insertion[1]]), file=insertions)
    insertions.close()        
    orfs.close()
    introns.close()
1216751
def seqrecord(self, allele, locus):
        
        try:
            db = self.server[self.dbversion + "_" + locus]
        except:
            self.logger.error("The database " + self.dbversion + "_"
                              + locus + " does not exist!")
            return ''

        seqrecord = db.lookup(name=allele)
        return seqrecord
1000013
def seqannotation(self, seqrecord, allele, loc):
        
        #seqrecord = self.seqrecord(allele, loc)
        complete_annotation = get_features(seqrecord)
        annotation = Annotation(annotation=complete_annotation,
                                method='match',
                                complete_annotation=True)

        if self.alignments:
            alignment = {f: self.annoated_alignments[loc][allele][f]['Seq']
                         for f in self.annoated_alignments[loc][allele].keys()}
            annotation.aligned = alignment

        return annotation
309871
def get_annotated_chain_sequence_string(self, chain_id, use_seqres_sequences_if_possible, raise_Exception_if_not_found = True):
        ''''''
        if use_seqres_sequences_if_possible and self.seqres_sequences and self.seqres_sequences.get(chain_id):
            return ('SEQRES', self.seqres_sequences[chain_id])
        elif self.atom_sequences.get(chain_id):
            return ('ATOM', self.atom_sequences[chain_id])
        elif raise_Exception_if_not_found:
            raise Exception('Error: Chain %s expected but not found.' % (str(chain_id)))
        else:
            return None
749780
def dannots2dalignbed2dannotsagg(cfg):
    
    datatmpd=cfg['datatmpd']
    
    daannotp=f'{datatmpd}/08_dannot.tsv'
    cfg['daannotp']=daannotp
    dannotsaggp=cfg['dannotsaggp']
    logging.info(basename(daannotp))
    if ((not exists(daannotp)) and (not exists(dannotsaggp))) or cfg['force']:
        gff_renamed_cols=[c+' annotation' if c in set(bed_colns).intersection(gff_colns) else c for c in gff_colns]
        dannots=pd.read_csv(cfg['annotationsbedp'],sep='\t',
                   names=bed_colns+gff_renamed_cols,
                           low_memory=False)
        dannots=del_Unnamed(dannots)

        dannots=dannots.set_index('id')
        dannots['annotations count']=1
        # separate ids from attribute columns
        dannots=lambda2cols(dannots,lambdaf=gffatributes2ids,
                    in_coln='attributes',
                to_colns=['gene name','gene id','transcript id','protein id','exon id'])
        dannots=dannots.drop(['attributes']+[c for c in gff_renamed_cols if 'annotation' in c],axis=1)
        logging.debug('or this step takes more time?')
        to_table(dannots,daannotp)
#         to_table_pqt(dannots,daannotp)
    else:
#         dannots=read_table_pqt(daannotp)
        dannots=read_table(daannotp)
        dannots=del_Unnamed(dannots)
        
    logging.info(basename(dannotsaggp))
    if not exists(dannotsaggp) or cfg['force']:
        if not 'dannots' in locals():
#             dannots=read_table_pqt(daannotp)
            dannots=pd.read_table(daannotp,low_memory=False)
        dannots=del_Unnamed(dannots)
        dannots=dannots.reset_index()
        logging.debug('aggregating the annotations')
        from rohan.dandage.io_sets import unique 
        cols2aggf={'annotations count':np.sum,
                  'type': unique,
                  'gene name': unique,
                  'gene id': unique,
                  'transcript id': unique,
                  'protein id': unique,
                  'exon id': unique}
        dannotsagg=dannots.groupby('id').agg(cols2aggf)
        dannotsagg['annotations count']=dannotsagg['annotations count']-1
        dannotsagg.loc[dannotsagg['annotations count']==0,'region']='intergenic'
        dannotsagg.loc[dannotsagg['annotations count']!=0,'region']='genic'
        logging.debug('end of the slowest step')            
        del dannots    
        dannotsagg=dannotsagg.reset_index()
#         to_table_pqt(dannotsagg,dannotsaggp)
        dannotsagg.to_csv(dannotsaggp,sep='\t')        
    return cfg
475557
def get_gene_info(ensembl_ids=None, hgnc_symbols=None):
    
    uniq_ensembl_ids = set(ensembl_id for ensembl_id in (ensembl_ids or []))
    uniq_hgnc_symbols = set(hgnc_symbol for hgnc_symbol in (hgnc_symbols or []))
    genes = []
    gene_data = []
    
    if uniq_ensembl_ids:
        for ensembl_id in uniq_ensembl_ids:
            for res in query_gene(ensembl_id=ensembl_id):
                gene_data.append(res) 

    elif uniq_hgnc_symbols:
        for hgnc_symbol in uniq_hgnc_symbols:
            query_res = query_gene(hgnc_symbol=hgnc_symbol)
            if query_res:
                for res in query_res:
                    gene_data.append(res)
            else:
                # If no result we add just the symbol
                gene_data.append({
                    'hgnc_symbol': hgnc_symbol,
                    'hgnc_id': None,
                    'ensembl_id': None,
                    'description': None,
                    'chrom': 'unknown',
                    'start': 0,
                    'stop': 0,
                    'hi_score': None,
                    'constraint_score': None,
                })
    for gene in gene_data:
        genes.append(Gene(
            symbol=gene ['hgnc_symbol'],
            hgnc_id=gene['hgnc_id'],
            ensembl_id=gene['ensembl_id'],
            description=gene['description'],
            chrom=gene['chrom'],
            start=gene['start'],
            stop=gene['stop'],
            location=get_cytoband_coord(gene['chrom'], gene['start']),
            hi_score=gene['hi_score'],
            constraint_score=gene['constraint_score'],
            omim_number=get_omim_number(gene['hgnc_symbol'])
            ))

    return genes
693187
def find_seq_homologues(self, return_raw=False):
        
        
        
        # Get sequence from structure/chain
        # We could use Bio.PDB.Polypeptide?
        
        from Bio.SCOP.Raf import to_one_letter_code
        
        s = self
        seq_iter = s.get_residues()
        seq_str = ''

        for aa in seq_iter:
            if aa.resname in to_one_letter_code:
                seq_str += to_one_letter_code[aa.resname]
        
        # Use BLAST to find homologous sequences with associated
        # structures in the PDB database.
        # Perhaps include local BLAST?
        
        from Bio.Blast.NCBIWWW import qblast

        # Adapt for short query sequences if needed
        # From http://www.ncbi.nlm.nih.gov/blast/producttable.shtml#shortp
        
        if len(seq_str) < 15:
            word_size = 2
            expect = 20000
            matrix_name = 'PAM30'
            filter = None
        else:
            word_size = 3
            expect = 10.0
            matrix_name = 'BLOSUM62'
            filter = 'SEG'
             
        query_result = qblast(  "blastp", "pdb", seq_str, 
                                word_size, expect, matrix_name, filter)
        
        if return_raw:
           return query_result 
            
        # Parse BLAST result to yield results
        # PDBID : (E-Value, Identity, Positives, Gaps, Alignment)
        
        from Bio.Blast import NCBIXML
        
        blast_records = NCBIXML.read(query_result)
        
        results = []
        
        for alignment in blast_records.alignments:
            for hsps in alignment.hsps:
                id_perc = "%s/%s" %(hsps.identities, alignment.length)
                pos_perc = "%s/%s" %(hsps.positives, alignment.length)
                gaps_perc = "%s/%s" %(hsps.gaps, alignment.length)

                pdb_id = alignment.title.split('|')[3]
                e_value = hsps.expect
                
                results.append( ( pdb_id,
                                "%2.5e" %e_value,
                                id_perc, 
                                pos_perc, 
                                gaps_perc,
                                '\n'.join([hsps.query, hsps.match, hsps.sbjct]) ) )

        return results