Join together alignment blocks to create a semi human projected local 
    alignment (small reference sequence deletions are kept as supported by 
    the local alignment).
180606
def do_interval( sources, index, out, ref_src, start, end, seq_db, missing_data, strand ):
    
    ref_src_size = None
    # Make sure the reference component is also the first in the source list
    assert sources[0].split('.')[0] == ref_src.split('.')[0], "%s != %s" \
        % ( sources[0].split('.')[0], ref_src.split('.')[0] )
    # Determine the overall length of the interval
    base_len = end - start
    # Counter for the last reference species base we have processed
    last_stop = start
    # Rows in maf blocks come in in arbitrary order, we'll convert things
    # to the destred order of the tiled block
    source_to_index = dict( ( name, i ) for ( i, name ) in enumerate( sources ) )
    # This gets all the maf blocks overlapping our interval of interest
    # NOTE: Unlike maf_tile we're expecting 
    # things to be single coverage in the reference species, so we won't 
    # sort by score and lay down.
    blocks = index.get( ref_src, start, end )
    # The last component seen for each species onto which we are tiling
    last_components = [ None ] * len( sources )
    last_status = [ None ] * len( sources )
    cols_needing_fill = [ 0 ] * len( sources )
    # The list of strings in which we build up the tiled alignment
    tiled_rows = [ "" for i in range( len( sources ) ) ]
    # Enumerate the (ordered) list of blocks
    for i, block in enumerate( blocks ):
        # Check for overlap in reference species
        ref = block.get_component_by_src_start( ref_src )
        if ref.start < last_stop:
            if ref.end < last_stop: 
                continue
            block = block.slice_by_component( ref, last_stop, min( end, ref.end ) )
            ref = block.get_component_by_src_start( ref_src )
        block = block.slice_by_component( ref, max( start, ref.start ), min( end, ref.end ) )
        ref = block.get_component_by_src_start( ref_src )
        # print block
        assert last_components[0] is None or ref.start >= last_components[0].end, \
            "MAF must be sorted and single coverage in reference species!"
        assert ref.strand == "+", \
            "MAF must have all reference species blocks on the plus strand"
        # Store the size of the reference sequence for building fake block   
        if ref_src_size is None:
            ref_src_size = ref.src_size
        # Handle the reference component seperately, it has no synteny status
        # but we will try to fill in missing sequence
        if ref.start > last_stop:
            # Need to fill in some reference sequence
            chunk_len = ref.start - last_stop
            text = bx.seq.nib.NibFile( open( seq_db[ ref_src ] ) ).get( last_stop, chunk_len ) 
            tiled_rows[0] += text
            for source in sources[1:]:
                cols_needing_fill[ source_to_index[ source ] ] += chunk_len
        # Do reference component
        chunk_len = len( ref.text )
        tiled_rows[0] += ref.text
        # Do each other component
        for source in sources[1:]:
            source_index = source_to_index[ source ]
            comp = block.get_component_by_src_start( source )
            if comp:
                if comp.synteny_left is None:
                    left_status, left_length = None, -1
                else:
                    left_status, left_length = comp.synteny_left
                if comp.synteny_right is None:
                    right_status, right_length = None, -1
                else:
                    right_status, right_length = comp.synteny_right
                # We have a component, do we need to do some filling?
                cols_to_fill = cols_needing_fill[ source_index ]
                if cols_to_fill > 0:
                    # Adjacent components should have matching status
                    ## assert last_status[ source_index ] is None or last_status[ source_index ] == left_status, \
                    ##     "left status (%s) does not match right status (%s) of last component for %s" \
                    ##         % ( left_status, last_status[ source_index ], source )
                    if left_status is None:
                        fill_char = guess_fill_char( last_components[source_index], comp )
                    else:
                        fill_char = get_fill_char( left_status )
                    tiled_rows[ source_index ] += ( fill_char * cols_to_fill )
                    cols_needing_fill[ source_index ] = 0
                # Okay, filled up to current position, now append the text
                tiled_rows[ source_index ] += comp.text
                assert len( tiled_rows[ source_index ] ) == len( tiled_rows[ 0 ] ), \
                    "length of tiled row should match reference row"
                last_components[ source_index ] = comp
                last_status[ source_index ] = right_status
            else:
                # No component, we'll have to fill this region when we know
                # the status
                cols_needing_fill[ source_index ] += chunk_len
        last_stop = ref.end
    # No more components, clean up the ends
    if last_stop < end:
        # Need to fill in some reference sequence
        chunk_len = end - last_stop
        tiled_rows[0] += bx.seq.nib.NibFile( open( seq_db[ ref_src ] ) ).get( last_stop, chunk_len ) 
        for source in sources[1:]:
            cols_needing_fill[ source_to_index[ source ] ] += chunk_len
    # Any final filling that needs to be done?
    for source in sources[1:]:
        source_index = source_to_index[ source ]
        fill_needed = cols_needing_fill[ source_index ]
        if fill_needed > 0:
            if last_components[ source_index ] is None:
                # print >>sys.stderr, "Never saw any components for %s, filling with @" % source
                fill_char = '@'
            else:
                if last_status[ source_index ] is None:
                    fill_char = '*'
                else:
                    fill_char = get_fill_char( last_status[ source_index ] )
            tiled_rows[ source_index ] += fill_char * fill_needed
        assert len( tiled_rows[ source_index ] ) == len( tiled_rows[ 0 ] ), \
            "length of tiled row should match reference row"
    # Okay, now make up the fake alignment from the tiled rows.
    tiled_rows = remove_all_gap_columns( tiled_rows )
    a = align.Alignment()
    for i, name in enumerate( sources ):
        text = "".join( tiled_rows[i] )
        size = len( text ) - text.count( "-" )
        if i == 0:
            if ref_src_size is None: ref_src_size = bx.seq.nib.NibFile( open( seq_db[ ref_src ] ) ).length
            c = align.Component( ref_src, start, end-start, "+", ref_src_size, text )
        else:
            c = align.Component( name + ".fake", 0, size, "?", size, text )
        a.add_component( c )
    if strand == '-':
        a = a.reverse_complement()
    out.write( a )
1053237
def align(s1,s2,test=False):
    
    import operator
    from Bio import pairwise2
    alignments = pairwise2.align.localms(s1.upper(),s2.upper(),1,-1,-5,-5)
    if test:
        print(alignments)
    alignsymb=np.nan
    score=np.nan
    sorted_alignments = sorted(alignments, key=operator.itemgetter(2))
    for a in alignments:
        alignstr=pairwise2.format_alignment(*a)
        alignsymb=alignstr.split('\n')[1]
        score=a[2]
        if test:
            print(alignstr)
        break
    return alignsymb,score
604776
def create_features(bam_in, loci_file, reference, out_dir):
    
    lenvec_plus = op.join(out_dir, 'genomic_lenvec.plus')
    lenvec_minus = op.join(out_dir, 'genomic_lenvec.minus')
    compute_genomic_cmd = ("compute_genomic_lenvectors "
                           "{bam_in} {lenvec_plus} "
                           "{lenvec_minus} "
                           "{min_len} "
                           "{max_len} ")
    index_genomic_cmd = ("index_genomic_lenvectors "
                         "{lenvec} ")
    genomic_lenvec = op.join(out_dir, 'genomic_lenvec')
    feat_len_file = op.join(out_dir, 'feat_lengths.txt')
    compute_locus_cmd = ("compute_locus_lenvectors "
                         "{loci_file} "
                         "{genomic_lenvec} "
                         "{min_len} "
                         "{max_len} "
                         "> {feat_len_file}")
    cov_S_file = op.join(out_dir, 'loci.cov_anti')
    coverage_anti_cmd = ("coverageBed -S -counts -b "
                         "{bam_in} -a {loci_file} "
                         "> {cov_S_file}")
    feat_posentropy = op.join(out_dir, 'feat_posentropy.txt')
    entropy_cmd = ("compute_locus_entropy.rb "
                   "{counts_reads} "
                   "> {feat_posentropy}")
    with utils.chdir(out_dir):
        run(compute_genomic_cmd.format(min_len=min_trimmed_read_len, max_len=max_trimmed_read_len, **locals()), "Run compute_genomic")
        run(index_genomic_cmd.format(lenvec=lenvec_plus), "Run index in plus")
        run(index_genomic_cmd.format(lenvec=lenvec_minus), "Run index in minus")
        run(compute_locus_cmd.format(min_len=min_trimmed_read_len, max_len=max_trimmed_read_len, **locals()), "Run compute locus")
        run(coverage_anti_cmd.format(**locals()), "Run coverage antisense")
        feat_antisense = _order_antisense_column(cov_S_file, min_trimmed_read_len)

        counts_reads = _reads_per_position(bam_in, loci_file, out_dir)
        run(entropy_cmd.format(**locals()), "Run entropy")

        rnafold = calculate_structure(loci_file, reference)
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
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
1272196
def get_conservation(block):
    
    consensus = block['sequences'][0]['seq']
    assert all(c.isupper() for c in consensus), \
            "So-called consensus contains indels!"
    # remove all non-consensus positions -- now alignment is easy
    cleaned = [[c for c in s['seq'] if not c.islower()]
               for s in block['sequences'][1:]]
    height = float(len(cleaned))
    # validation
    for row in cleaned:
        if len(row) != len(consensus):
            raise ValueError("Aligned sequence length (%s) doesn't match "
                             "consensus (%s)"
                             % (len(row), len(consensus)))
    # transpose & go
    columns = zip(*cleaned)
    return dict((idx + 1, columns[idx].count(cons_char) / height)
                for idx, cons_char in enumerate(consensus))
737567
def merge(assembled_gtfs, ref_file, gtf_file, num_cores, data):
    
    assembled_file = tempfile.NamedTemporaryFile(delete=False).name
    with open(assembled_file, "w") as temp_handle:
        for assembled in assembled_gtfs:
            temp_handle.write(assembled + "\n")
    out_dir = os.path.join("assembly", "cuffmerge")
    merged_file = os.path.join(out_dir, "merged.gtf")
    out_file = os.path.join(out_dir, "assembled.gtf")
    if file_exists(out_file):
        return out_file
    if not file_exists(merged_file):
        with file_transaction(data, out_dir) as tmp_out_dir:
            cmd = ("cuffmerge -o {tmp_out_dir} --ref-gtf {gtf_file} "
                   "--num-threads {num_cores} --ref-sequence {ref_file} "
                   "{assembled_file}")
            cmd = cmd.format(**locals())
            message = ("Merging the following transcript assemblies with "
                       "Cuffmerge: %s" % ", ".join(assembled_gtfs))
            do.run(cmd, message)
    clean, _ = clean_assembly(merged_file)
    fixed = fix_cufflinks_attributes(gtf_file, clean, data)
    classified = annotate_gtf.annotate_novel_coding(fixed, gtf_file, ref_file,
                                                    data)
    filtered = annotate_gtf.cleanup_transcripts(classified, gtf_file, ref_file)
    shutil.move(filtered, out_file)
    return out_file
965669
def _load_data_alignment(self, chain1, chain2):
        
        parser = PDB.PDBParser(QUIET=True)
        ppb = PDB.PPBuilder()

        structure1 = parser.get_structure(chain1, self.pdb1)
        structure2 = parser.get_structure(chain2, self.pdb2)

        seq1 = str(ppb.build_peptides(structure1)[0].get_sequence())
        seq2 = str(ppb.build_peptides(structure2)[0].get_sequence())

        # Alignment parameters taken from PconsFold renumbering script.
        align = pairwise2.align.globalms(seq1, seq2, 2, -1, -0.5, -0.1)[0]
        indexes = set(i for i, (s1, s2) in enumerate(zip(align[0], align[1]))
                      if s1 != '-' and s2 != '-')
        coord1 = np.hstack([np.concatenate((r['CA'].get_coord(), (1,)))[:, None]
                            for i, r in enumerate(structure1.get_residues())
                            if i in indexes and 'CA' in r]).astype(DTYPE,
                                                                   copy=False)
        coord2 = np.hstack([np.concatenate((r['CA'].get_coord(), (1,)))[:, None]
                            for i, r in enumerate(structure2.get_residues())
                            if i in indexes and 'CA' in r]).astype(DTYPE,
                                                                   copy=False)

        self.coord1 = coord1
        self.coord2 = coord2
        self.N = len(seq1)
750060
def run(align_bams, items, ref_file, assoc_files, region=None, out_file=None):
    
    paired = vcfutils.get_paired_bams(align_bams, items)
    assert paired and not paired.normal_bam, ("Pisces supports tumor-only variant calling: %s" %
                                              (",".join([dd.get_sample_name(d) for d in items])))
    vrs = bedutils.population_variant_regions(items)
    target = shared.subset_variant_regions(vrs, region,
                                            out_file, items=items, do_merge=True)
    min_af = float(dd.get_min_allele_fraction(paired.tumor_data)) / 100.0
    if not utils.file_exists(out_file):
        base_out_name = utils.splitext_plus(os.path.basename(paired.tumor_bam))[0]
        raw_file = "%s.vcf" % utils.splitext_plus(out_file)[0]
        with file_transaction(paired.tumor_data, raw_file) as tx_out_file:
            ref_dir = _prep_genome(os.path.dirname(tx_out_file), paired.tumor_data)
            out_dir = os.path.dirname(tx_out_file)
            cores = dd.get_num_cores(paired.tumor_data)
            emit_min_af = min_af / 10.0
            cmd = ("pisces --bampaths {paired.tumor_bam} --genomepaths {ref_dir} --intervalpaths {target} "
                   "--maxthreads {cores} --minvf {emit_min_af} --vffilter {min_af} "
                   "--ploidy somatic --gvcf false -o {out_dir}")
            # Recommended filtering for low frequency indels
            # https://github.com/bcbio/bcbio-nextgen/commit/49d0cbb1f6dcbea629c63749e2f9813bd06dcee3#commitcomment-29765373
            cmd += " -RMxNFilter 5,9,0.35"
            # For low frequency UMI tagged variants, set higher variant thresholds
            # https://github.com/Illumina/Pisces/issues/14#issuecomment-399756862
            if min_af < (1.0 / 100.0):
                cmd += " --minbasecallquality 30"
            do.run(cmd.format(**locals()), "Pisces tumor-only somatic calling")
            shutil.move(os.path.join(out_dir, "%s.vcf" % base_out_name),
                        tx_out_file)
        vcfutils.bgzip_and_index(raw_file, paired.tumor_data["config"],
                                 prep_cmd="sed 's#%s.bam#%s#' | %s" %
                                 (base_out_name, dd.get_sample_name(paired.tumor_data),
                                  vcfutils.add_contig_to_header_cl(dd.get_ref_file(paired.tumor_data), out_file)))
    return vcfutils.bgzip_and_index(out_file, paired.tumor_data["config"])
575844
def variant_context_w_alignment(am, var, margin=20, tx_ac=None):
    """"""
    from uta_align.align.algorithms import align, cigar_alignment

    fh = full_house(am, var, tx_ac=tx_ac)
    tm = am._fetch_AlignmentMapper(fh['n'].ac, fh['g'].ac, am.alt_aln_method)
    strand = tm.strand
    span_g = _ival_to_span(fh['g'].posedit.pos)
    span_g = (span_g[0] - margin, span_g[1] + margin)
    ival_g = Interval(SimplePosition(span_g[0]), SimplePosition(span_g[1]))
    ival_n = tm.g_to_n(ival_g)
    assert ival_n.start.offset == 0 and ival_n.end.offset == 0, "limited to coding variants"
    span_n = _ival_to_span(ival_n)
    ival_c = tm.g_to_c(ival_g)
    span_c = _ival_to_span(ival_c)
    seq_gt = am.hdp.get_seq(fh['g'].ac, span_g[0] - 1, span_g[1])
    seq_gb = complement(seq_gt)
    seq_n = am.hdp.get_seq(fh['n'].ac, span_n[0] - 1, span_n[1])
    if strand == 1:
        a = align(bytes(seq_gt), bytes(seq_n), b'global', extended_cigar=True)
    else:
        seq_n = ''.join(reversed(seq_n))
        a = align(bytes(seq_gb), bytes(seq_n), b'global', extended_cigar=True)

    aseq_gt, _ = cigar_alignment(seq_gt, a.query, a.cigar, hide_match=False)
    aseq_gb, aseq_n = cigar_alignment(seq_gb, a.query, a.cigar, hide_match=False)
    aln_str = _reformat_aln_str(cigar_alignment(a.ref, a.query, a.cigar, hide_match=True)[1])
    s_dir = '>' if strand == 1 else '<'

    lines = [
        [
            1,
            0,
            seq_line_fmt(
                var=fh['c'],
                span=span_c if strand == 1 else list(reversed(span_c)),
                content='',
                dir=s_dir),
        ],
        [
            2,
            0,
            seq_line_fmt(
                var=fh['n'],
                span=span_n if strand == 1 else list(reversed(span_n)),
                content=aseq_n,
                dir=s_dir),
        ],
        [
            3,
            0,
            _line_fmt.format(pre='', content=aln_str, post=a.cigar.to_string(), comment=''),
        ],
        [
            4,
            1,
            seq_line_fmt(var=fh['g'], span=span_g, content=aseq_gt, dir='>'),
        ],
        [
            4,
            2,
            seq_line_fmt(var=fh['g'], span=span_g, content=aseq_gb, dir='<'),
        ],
        [
            5,
            0,
            pointer_line(var=fh['g'], span=span_g),
        ],
    ]
    if strand == -1:
        lines.sort(key=lambda e: (-e[0], e[1]))
    return '\n'.join(r[2] for r in lines)
921357
def align(args):
    
    from jcvi.formats.fastq import guessoffset

    p = OptionParser(align.__doc__)
    p.add_option("--rnaseq", default=False, action="store_true",
                 help="Input is RNA-seq reads, turn splicing on")
    p.add_option("--native", default=False, action="store_true",
                 help="Convert GSNAP output to NATIVE format")
    p.set_home("eddyyeh")
    p.set_outdir()
    p.set_cpus()
    opts, args = p.parse_args(args)

    if len(args) == 2:
        logging.debug("Single-end alignment")
    elif len(args) == 3:
        logging.debug("Paired-end alignment")
    else:
        sys.exit(not p.print_help())

    dbfile, readfile = args[:2]
    outdir = opts.outdir
    assert op.exists(dbfile) and op.exists(readfile)
    prefix = get_prefix(readfile, dbfile)
    logfile = op.join(outdir, prefix + ".log")
    gsnapfile = op.join(outdir, prefix + ".gsnap")
    nativefile = gsnapfile.rsplit(".", 1)[0] + ".unique.native"

    if not need_update((dbfile, readfile), gsnapfile):
        logging.error("`{0}` exists. `gsnap` already run.".format(gsnapfile))
    else:
        dbdir, dbname = check_index(dbfile)
        cmd = "gsnap -D {0} -d {1}".format(dbdir, dbname)
        cmd += " -B 5 -m 0.1 -i 2 -n 3"  # memory, mismatch, indel penalty, nhits
        if opts.rnaseq:
            cmd += " -N 1"
        cmd += " -t {0}".format(opts.cpus)
        cmd += " --gmap-mode none --nofails"
        if readfile.endswith(".gz"):
            cmd += " --gunzip"
        try:
            offset = "sanger" if guessoffset([readfile]) == 33 else "illumina"
            cmd += " --quality-protocol {0}".format(offset)
        except AssertionError:
            pass
        cmd += " " + " ".join(args[1:])
        sh(cmd, outfile=gsnapfile, errfile=logfile)

    if opts.native:
        EYHOME = opts.eddyyeh_home
        if need_update(gsnapfile, nativefile):
            cmd = op.join(EYHOME, "convert2native.pl")
            cmd += " --gsnap {0} -o {1}".format(gsnapfile, nativefile)
            cmd += " -proc {0}".format(opts.cpus)
            sh(cmd)

    return gsnapfile, logfile
1081527
def pair_hmm_align_unaligned_seqs(seqs, moltype=DNA_cogent, params={}):
    

    seqs = LoadSeqs(data=seqs, moltype=moltype, aligned=False)
    try:
        s1, s2 = seqs.values()
    except ValueError:
        raise ValueError(
            "Pairwise aligning of seqs requires exactly two seqs.")

    try:
        gap_open = params['gap_open']
    except KeyError:
        gap_open = 5
    try:
        gap_extend = params['gap_extend']
    except KeyError:
        gap_extend = 2
    try:
        score_matrix = params['score_matrix']
    except KeyError:
        score_matrix = make_dna_scoring_dict(
            match=1, transition=-1, transversion=-1)

    return local_pairwise(s1, s2, score_matrix, gap_open, gap_extend)
1159976
def liftover_to_genome(pass_pos, gtf):
    """"""

    fixed_pos = []
    for pos in pass_pos:
        if pos["chrom"] not in gtf:
            continue
        db_pos = gtf[pos["chrom"]][0]
        mut = _parse_mut(pos["sv"])
        print([db_pos, pos])
        if db_pos[3] == "+":
            pos['pre_pos'] = db_pos[1] + pos["pre_pos"] + 1
        else:
            pos['pre_pos'] = db_pos[2] - (pos["pre_pos"] - 1)
        pos['chrom'] = db_pos[0]
        pos['nt'] = list(mut[0])
        fixed_pos.append(pos)

    _print_header(fixed_pos)
    for pos in fixed_pos:
        print_vcf(pos)
1025370
public static void fillAlignedAtomArrays(AFPChain afpChain, Atom[] ca1,
			Atom[] ca2, Atom[] ca1aligned, Atom[] ca2aligned) {

		int pos=0;
		int[] blockLens = afpChain.getOptLen();
		int[][][] optAln = afpChain.getOptAln();
		assert(afpChain.getBlockNum() <= optAln.length);

		for (int block=0; block < afpChain.getBlockNum(); block++) {
			for(int i=0;i<blockLens[block];i++) {
				int pos1 = optAln[block][0][i];
				int pos2 = optAln[block][1][i];
				Atom a1 = ca1[pos1];
				Atom a2 = (Atom) ca2[pos2].clone();
				ca1aligned[pos] = a1;
				ca2aligned[pos] = a2;
				pos++;
			}
		}

		// this can happen when we load an old XML serialization which did not support modern ChemComp representation of modified residues.
		if (pos != afpChain.getOptLength()){
			logger.warn("AFPChainScorer getTMScore: Problems reconstructing alignment! nr of loaded atoms is " + pos + " but should be " + afpChain.getOptLength());
			// we need to resize the array, because we allocated too many atoms earlier on.
			ca1aligned = (Atom[]) resizeArray(ca1aligned, pos);
			ca2aligned = (Atom[]) resizeArray(ca2aligned, pos);
		}

	}
12299
def map_alignment(self, prot_seq, nucl_seq):
        
        if prot_seq.id != nucl_seq.id:
            logging.warning(
                'ID mismatch: %s != %s. Are the sequences in the same order?',
                prot_seq.id, nucl_seq.id)

        # Ungap nucleotides
        codons = batch(str(nucl_seq.seq.ungap('-')), 3)
        codons = [''.join(i) for i in codons]
        codon_iter = iter(codons)

        ungapped_prot = str(prot_seq.seq).replace('-', '')

        if len(ungapped_prot) != len(codons):
            table = self.translation_table.forward_table
            prot_str = ' '.join(' ' + p + ' ' for p in ungapped_prot)
            codon_str = ' '.join(codons)
            trans_str = ' '.join(' ' + table.get(codon, 'X') + ' '
                                 for codon in codons)
            raise ValueError("""Length of codon sequence ({0}) does not match \
length of protein sequence ({1}) for {2}
Protein:       {3}
Codons:        {4}
Trans. Codons: {5}""".format(len(codons), len(ungapped_prot), nucl_seq.id, prot_str,
                            codon_str, trans_str))

        try:
            nucl_align = ['---' if p == '-' else next(codon_iter)
                          for p in str(prot_seq.seq)]
        except StopIteration:
            assert False  # Should be checked above

        result = SeqRecord(Seq(''.join(nucl_align)), id=nucl_seq.id,
                           description=nucl_seq.description)

        # Validate
        self._validate_translation(prot_seq.seq.upper(), result.seq.upper())

        return result
835536
def alignment_to_partials(alignment, missing_data=None):
    
    partials_dict = {}
    for (name, sequence) in alignment.get_sequences():
        datatype = 'dna' if alignment.is_dna() else 'protein'
        partials_dict[name] = seq_to_partials(sequence, datatype)

    if missing_data is not None:
        l = len(alignment)
        for name in missing_data:
            if name not in partials_dict:
                partials_dict[name] = seq_to_partials('-'*l, datatype)
    return partials_dict
1072342
def make_V_and_J_mask_mapping(self, genV, genJ):
        
        #construct mapping between allele/gene names and index for custom V_usage_masks
        V_allele_names = [V[0] for V in genV]
        V_mask_mapping = {}
        for v in set([x.split('*')[0] for x in V_allele_names]):
            V_mask_mapping[v] = []
        for v in set(['V'.join((x.split('*')[0]).split('V')[1:]) for x in V_allele_names]):
            V_mask_mapping[v] = []
        for i, v in enumerate(V_allele_names):
            V_mask_mapping[v] = [i]
            V_mask_mapping['V'.join((v.split('*')[0]).split('V')[1:])].append(i)
            V_mask_mapping[v.split('*')[0]].append(i)      
        
        #construct mapping between allele/gene names and index for custom J_usage_masks
        J_allele_names = [J[0] for J in genJ]
        J_mask_mapping = {}
        for j in set([x.split('*')[0] for x in J_allele_names]):
            J_mask_mapping[j] = []
        for j in set(['J'.join((x.split('*')[0]).split('J')[1:]) for x in J_allele_names]):
            J_mask_mapping[j] = []
        for i, j in enumerate(J_allele_names):
            J_mask_mapping[j] = [i]
            J_mask_mapping['J'.join((j.split('*')[0]).split('J')[1:])].append(i)
            J_mask_mapping[j.split('*')[0]].append(i)
            
        self.V_allele_names = V_allele_names
        self.V_mask_mapping = V_mask_mapping
        self.J_allele_names = J_allele_names
        self.J_mask_mapping = J_mask_mapping