Full extraction of seq flanking the sites.
514412
def extract_full(rec, sites, flank, fw):
    
    for s in sites:
        newid = "{0}:{1}".format(rec.name, s)
        left = max(s - flank, 0)
        right = min(s + flank, len(rec))
        frag = rec.seq[left:right].strip("Nn")
        newrec = SeqRecord(frag, id=newid, description="")
        SeqIO.write([newrec], fw, "fasta")
632753
def get_start_stops(transcript_sequence, start_codons=None, stop_codons=None):
        
        transcript_sequence = transcript_sequence.upper()  # for comparison with codons below
        if not start_codons:
            start_codons = ['ATG']
        if not stop_codons:
            stop_codons = ['TAA', 'TAG', 'TGA']

        seq_frames = {1: {'starts': [], 'stops': []},
                      2: {'starts': [], 'stops': []},
                      3: {'starts': [], 'stops': []}}

        for codons, positions in ((start_codons, 'starts'),
                                  (stop_codons, 'stops')):
            if len(codons) > 1:
                pat = re.compile('|'.join(codons))
            else:
                pat = re.compile(codons[0])

            for m in re.finditer(pat, transcript_sequence):
                # Increment position by 1, Frame 1 starts at position 1 not 0
                start = m.start() + 1
                rem = start % 3
                if rem == 1:  # frame 1
                    seq_frames[1][positions].append(start)
                elif rem == 2:  # frame 2
                    seq_frames[2][positions].append(start)
                elif rem == 0:  # frame 3
                    seq_frames[3][positions].append(start)
        return seq_frames
218926
def get_gene_seqs(database_path, gene):
    
    gene_path = database_path + "/" + gene + ".fsa"
    gene_seq = ""
    # Open fasta file
    with open(gene_path) as gene_file:
        header = gene_file.readline()
        for line in gene_file:
            seq = line.strip()
            gene_seq += seq
    return gene_seq
1364850
private static Object[] extractExtendedFeatures(Configuration configuration, int length)
    {
        Object[] featureMap = new Object[length];

        State state = configuration.state;
        Sentence sentence = configuration.sentence;

        int b0Position = 0;
        int b1Position = 0;
        int b2Position = 0;
        int s0Position = 0;

        long svr = 0; // stack right valency
        long svl = 0; // stack left valency
        long bvl = 0; // buffer left valency

        long b0w = 0;
        long b0p = 0;

        long b1w = 0;
        long b1p = 0;

        long b2w = 0;
        long b2p = 0;

        long s0w = 0;
        long s0p = 0;
        long s0l = 0;

        long bl0p = 0;
        long bl0w = 0;
        long bl0l = 0;

        long bl1w = 0;
        long bl1p = 0;
        long bl1l = 0;

        long sr0p = 0;
        long sr0w = 0;
        long sr0l = 0;

        long sh0w = 0;
        long sh0p = 0;
        long sh0l = 0;

        long sl0p = 0;
        long sl0w = 0;
        long sl0l = 0;

        long sr1w = 0;
        long sr1p = 0;
        long sr1l = 0;

        long sh1w = 0;
        long sh1p = 0;

        long sl1w = 0;
        long sl1p = 0;
        long sl1l = 0;

        long sdl = 0;
        long sdr = 0;
        long bdl = 0;

        int[] words = sentence.getWords();
        int[] tags = sentence.getTags();

        if (0 < state.bufferSize())
        {
            b0Position = state.bufferHead();
            b0w = b0Position == 0 ? 0 : words[b0Position - 1];
            b0w += 2;
            b0p = b0Position == 0 ? 0 : tags[b0Position - 1];
            b0p += 2;
            bvl = state.leftValency(b0Position);

            int leftMost = state.leftMostModifier(state.getBufferItem(0));
            if (leftMost >= 0)
            {
                bl0p = leftMost == 0 ? 0 : tags[leftMost - 1];
                bl0p += 2;
                bl0w = leftMost == 0 ? 0 : words[leftMost - 1];
                bl0w += 2;
                bl0l = state.getDependent(leftMost);
                bl0l += 2;

                int l2 = state.leftMostModifier(leftMost);
                if (l2 >= 0)
                {
                    bl1w = l2 == 0 ? 0 : words[l2 - 1];
                    bl1w += 2;
                    bl1p = l2 == 0 ? 0 : tags[l2 - 1];
                    bl1p += 2;
                    bl1l = state.getDependent(l2);
                    bl1l += 2;
                }
            }

            if (1 < state.bufferSize())
            {
                b1Position = state.getBufferItem(1);
                b1w = b1Position == 0 ? 0 : words[b1Position - 1];
                b1w += 2;
                b1p = b1Position == 0 ? 0 : tags[b1Position - 1];
                b1p += 2;

                if (2 < state.bufferSize())
                {
                    b2Position = state.getBufferItem(2);

                    b2w = b2Position == 0 ? 0 : words[b2Position - 1];
                    b2w += 2;
                    b2p = b2Position == 0 ? 0 : tags[b2Position - 1];
                    b2p += 2;
                }
            }
        }

        if (0 < state.stackSize())
        {
            s0Position = state.stackTop();
            s0w = s0Position == 0 ? 0 : words[s0Position - 1];
            s0w += 2;
            s0p = s0Position == 0 ? 0 : tags[s0Position - 1];
            s0p += 2;
            s0l = state.getDependent(s0Position);
            s0l += 2;

            svl = state.leftValency(s0Position);
            svr = state.rightValency(s0Position);

            int leftMost = state.leftMostModifier(s0Position);
            if (leftMost >= 0)
            {
                sl0p = leftMost == 0 ? 0 : tags[leftMost - 1];
                sl0p += 2;
                sl0w = leftMost == 0 ? 0 : words[leftMost - 1];
                sl0w += 2;
                sl0l = state.getDependent(leftMost);
                sl0l += 2;
            }

            int rightMost = state.rightMostModifier(s0Position);
            if (rightMost >= 0)
            {
                sr0p = rightMost == 0 ? 0 : tags[rightMost - 1];
                sr0p += 2;
                sr0w = rightMost == 0 ? 0 : words[rightMost - 1];
                sr0w += 2;
                sr0l = state.getDependent(rightMost);
                sr0l += 2;
            }

            int headIndex = state.getHead(s0Position);
            if (headIndex >= 0)
            {
                sh0w = headIndex == 0 ? 0 : words[headIndex - 1];
                sh0w += 2;
                sh0p = headIndex == 0 ? 0 : tags[headIndex - 1];
                sh0p += 2;
                sh0l = state.getDependent(headIndex);
                sh0l += 2;
            }

            if (leftMost >= 0)
            {
                int l2 = state.leftMostModifier(leftMost);
                if (l2 >= 0)
                {
                    sl1w = l2 == 0 ? 0 : words[l2 - 1];
                    sl1w += 2;
                    sl1p = l2 == 0 ? 0 : tags[l2 - 1];
                    sl1p += 2;
                    sl1l = state.getDependent(l2);
                    sl1l += 2;
                }
            }
            if (headIndex >= 0)
            {
                if (state.hasHead(headIndex))
                {
                    int h2 = state.getHead(headIndex);
                    sh1w = h2 == 0 ? 0 : words[h2 - 1];
                    sh1w += 2;
                    sh1p = h2 == 0 ? 0 : tags[h2 - 1];
                    sh1p += 2;
                }
            }
            if (rightMost >= 0)
            {
                int r2 = state.rightMostModifier(rightMost);
                if (r2 >= 0)
                {
                    sr1w = r2 == 0 ? 0 : words[r2 - 1];
                    sr1w += 2;
                    sr1p = r2 == 0 ? 0 : tags[r2 - 1];
                    sr1p += 2;
                    sr1l = state.getDependent(r2);
                    sr1l += 2;
                }
            }
        }
        int index = 0;

        long b0wp = b0p;
        b0wp |= (b0w << 8);
        long b1wp = b1p;
        b1wp |= (b1w << 8);
        long s0wp = s0p;
        s0wp |= (s0w << 8);
        long b2wp = b2p;
        b2wp |= (b2w << 8);

        /**
         * From single words
         */
        if (s0w != 1)
        {
            featureMap[index++] = s0wp;
            featureMap[index++] = s0w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p;

        if (b0w != 1)
        {
            featureMap[index++] = b0wp;
            featureMap[index++] = b0w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = b0p;

        if (b1w != 1)
        {
            featureMap[index++] = b1wp;
            featureMap[index++] = b1w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = b1p;

        if (b2w != 1)
        {
            featureMap[index++] = b2wp;
            featureMap[index++] = b2w;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
        }
        featureMap[index++] = b2p;

        /**
         * from word pairs
         */
        if (s0w != 1 && b0w != 1)
        {
            featureMap[index++] = (s0wp << 28) | b0wp;
            featureMap[index++] = (s0wp << 20) | b0w;
            featureMap[index++] = (s0w << 28) | b0wp;
        }
        else
        {
            featureMap[index++] = null;
            featureMap[index++] = null;
            featureMap[index++] = null;
        }

        if (s0w != 1)
        {
            featureMap[index++] = (s0wp << 8) | b0p;
        }
        else
        {
            featureMap[index++] = null;
        }

        if (b0w != 1)
        {
            featureMap[index++] = (s0p << 28) | b0wp;
        }
        else
        {
            featureMap[index++] = null;
        }

        if (s0w != 1 && b0w != 1)
        {
            featureMap[index++] = (s0w << 20) | b0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (s0p << 8) | b0p;
        featureMap[index++] = (b0p << 8) | b1p;

        /**
         * from three words
         */
        featureMap[index++] = (b0p << 16) | (b1p << 8) | b2p;
        featureMap[index++] = (s0p << 16) | (b0p << 8) | b1p;
        featureMap[index++] = (sh0p << 16) | (s0p << 8) | b0p;
        featureMap[index++] = (s0p << 16) | (sl0p << 8) | b0p;
        featureMap[index++] = (s0p << 16) | (sr0p << 8) | b0p;
        featureMap[index++] = (s0p << 16) | (b0p << 8) | bl0p;

        /**
         * distance
         */
        long distance = 0;
        if (s0Position > 0 && b0Position > 0)
            distance = Math.abs(b0Position - s0Position);
        if (s0w != 1)
        {
            featureMap[index++] = s0w | (distance << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (distance << 8);
        if (b0w != 1)
        {
            featureMap[index++] = b0w | (distance << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = b0p | (distance << 8);
        if (s0w != 1 && b0w != 1)
        {
            featureMap[index++] = s0w | (b0w << 20) | (distance << 40);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (b0p << 8) | (distance << 28);

        /**
         * Valency information
         */
        if (s0w != 1)
        {
            featureMap[index++] = s0w | (svr << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (svr << 8);
        if (s0w != 1)
        {
            featureMap[index++] = s0w | (svl << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = s0p | (svl << 8);
        if (b0w != 1)
        {
            featureMap[index++] = b0w | (bvl << 20);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = b0p | (bvl << 8);

        /**
         * Unigrams
         */
        if (sh0w != 1)
        {
            featureMap[index++] = sh0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sh0p;
        featureMap[index++] = s0l;
        if (sl0w != 1)
        {
            featureMap[index++] = sl0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sl0p;
        featureMap[index++] = sl0l;
        if (sr0w != 1)
        {
            featureMap[index++] = sr0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sr0p;
        featureMap[index++] = sr0l;
        if (bl0w != 1)
        {
            featureMap[index++] = bl0w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = bl0p;
        featureMap[index++] = bl0l;

        /**
         * From third order features
         */
        if (sh1w != 1)
        {
            featureMap[index++] = sh1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sh1p;
        featureMap[index++] = sh0l;
        if (sl1w != 1)
        {
            featureMap[index++] = sl1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sl1p;
        featureMap[index++] = sl1l;
        if (sr1w != 1)
        {
            featureMap[index++] = sr1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = sr1p;
        featureMap[index++] = sr1l;
        if (bl1w != 1)
        {
            featureMap[index++] = bl1w;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = bl1p;
        featureMap[index++] = bl1l;
        featureMap[index++] = s0p | (sl0p << 8) | (sl1p << 16);
        featureMap[index++] = s0p | (sr0p << 8) | (sr1p << 16);
        featureMap[index++] = s0p | (sh0p << 8) | (sh1p << 16);
        featureMap[index++] = b0p | (bl0p << 8) | (bl1p << 16);

        /**
         * label set
         */
        if (s0Position >= 0)
        {
            sdl = state.leftDependentLabels(s0Position);
            sdr = state.rightDependentLabels(s0Position);
        }

        if (b0Position >= 0)
        {
            bdl = state.leftDependentLabels(b0Position);
        }

        if (s0w != 1)
        {
            featureMap[index++] = (s0w + "|" + sdr);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (s0p + "|" + sdr);
        if (s0w != 1)
        {
            featureMap[index++] = s0w + "|" + sdl;
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (s0p + "|" + sdl);
        if (b0w != 1)
        {
            featureMap[index++] = (b0w + "|" + bdl);
        }
        else
        {
            featureMap[index++] = null;
        }
        featureMap[index++] = (b0p + "|" + bdl);
        return featureMap;
    }
79767
def get_transcript(self,exon_bounds='max'):
    """"""
    out = Transcript()
    out.junctions = [x.get_junction() for x in self.junction_groups]
    # check for single exon transcript
    if len(out.junctions) == 0:
      leftcoord = min([x.exons[0].range.start for x in self.transcripts])
      rightcoord = max([x.exons[-1].range.end for x in self.transcripts])
      e = Exon(GenomicRange(x.exons[0].range.chr,leftcoord,rightcoord))
      e.set_is_leftmost()
      e.set_is_rightmost()
      out.exons.append(e)
      return out
    # get internal exons
    self.exons = []
    for i in range(0,len(self.junction_groups)-1):
      j1 = self.junction_groups[i].get_junction()
      j2 = self.junction_groups[i+1].get_junction()
      e = Exon(GenomicRange(j1.right.chr,j1.right.end,j2.left.start))
      e.set_left_junc(j1)
      e.set_right_junc(j2)
      #print str(i)+" to "+str(i+1)
      out.exons.append(e)
    # get left exon
    left_exons = [y for y in [self.transcripts[e[0]].junctions[e[1]].get_left_exon() for e in self.junction_groups[0].evidence] if y]
    if len(left_exons) == 0:
      sys.stderr.write("ERROR no left exon\n")
      sys.exit()
    e_left = Exon(GenomicRange(out.junctions[0].left.chr,\
                               min([x.range.start for x in left_exons]),
                               out.junctions[0].left.start))
    e_left.set_right_junc(out.junctions[0])
    out.exons.insert(0,e_left)
    # get right exon
    right_exons = [y for y in [self.transcripts[e[0]].junctions[e[1]].get_right_exon() for e in self.junction_groups[-1].evidence] if y]
    if len(right_exons) == 0:
      sys.stderr.write("ERROR no right exon\n")
      sys.exit()
    e_right = Exon(GenomicRange(out.junctions[-1].right.chr,\
                               out.junctions[-1].right.end,\
                               max([x.range.end for x in right_exons])))
    e_right.set_left_junc(out.junctions[-1])
    out.exons.append(e_right)
    return out
1281403
def _get_purecn_files(paired, work_dir, require_exist=False):
    
    out_base = os.path.join(work_dir, "%s-purecn" % (dd.get_sample_name(paired.tumor_data)))
    out = {"plot": {}}
    all_files = []
    for plot in ["chromosomes", "local_optima", "segmentation", "summary"]:
        if plot == "summary":
            cur_file = "%s.pdf" % out_base
        else:
            cur_file = "%s_%s.pdf" % (out_base, plot)
        if not require_exist or os.path.exists(cur_file):
            out["plot"][plot] = cur_file
            all_files.append(os.path.basename(cur_file))
    for key, ext in [["hetsummary", ".csv"], ["dnacopy", "_dnacopy.seg"], ["genes", "_genes.csv"],
                     ["log", ".log"], ["loh", "_loh.csv"], ["rds", ".rds"],
                     ["variants", "_variants.csv"]]:
        cur_file = "%s%s" % (out_base, ext)
        if not require_exist or os.path.exists(cur_file):
            out[key] = cur_file
            all_files.append(os.path.basename(cur_file))
    return out_base, out, all_files
670565
def create_dexseq_annotation(gff, count_file):
    
    out_file = count_file + ".ann"
    if file_exists(out_file):
        return out_file
    with file_transaction(out_file) as tx_out:
        with open(tx_out, 'w') as out_handle:
            with open(gff) as in_handle:
                for line in in_handle:
                    cols = line.strip().split("\t")
                    if cols[2] == "exonic_part":
                        exon = [f for f in cols[8].split(";") if f.strip().startswith("exonic_part_number")]
                        gene = [f for f in cols[8].split(";") if f.strip().startswith("gene_id")]
                        exon = exon[0].replace("\"", "").split()[1]
                        gene = gene[0].replace("\"", "").split()[1]
                        length = int(cols[4]) - int(cols[3]) + 1
                        line = "%s:%s\t%s\t%s\t%s\t%s\t%s\t%s\n" % (gene, exon, gene,
                                                                    cols[0], cols[3],
                                                                    cols[4],
                                                                    length, cols[6])
                        out_handle.write(line)
947891
def read_alignment(out_sam, loci, seqs, out_file):
    """"""
    hits = defaultdict(list)
    with open(out_file, "w") as out_handle:
        samfile = pysam.Samfile(out_sam, "r")
        for a in samfile.fetch():
            if not a.is_unmapped:
                nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                a = makeBED(a)
                if not a:
                    continue
                ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
        for hit in hits.values():
            nm = hit[0][0]
            for l in hit:
                if nm == l[0]:
                    print(l[1], file=out_handle)
    return out_file
356759
def get_genomic_seq_for_region(self, chrom, start_pos, end_pos):
        
        
        headers = {"content-type": "text/plain"}
        
        self.attempt = 0
        ext = "/sequence/region/human/{}:{}..{}:1".format(chrom, start_pos, end_pos)
        
        return self.ensembl_request(ext, headers)
178254
def seqs_from_fastacmd(acc_list, blast_db,is_protein=True):
    """"""
    fasta_cmd_res = fasta_cmd_get_seqs(acc_list, blast_db=blast_db, \
        is_protein=is_protein)
    recs = FastaCmdFinder(fasta_cmd_res['StdOut'])
    result = {}
    for rec in recs:
        try:
            result[rec[0][1:].strip()] = ''.join(map(strip, rec[1:]))
        except IndexError:  #maybe we didn't get a sequence?
            pass
    fasta_cmd_res.cleanUp()
    return result
1068667
def gene_panel(self, ax, feature):
        
        from gffutils.contrib.plotting import Gene
        extent = [feature.start, feature.stop]
        nearby_genes = self.db.region(
            (feature.chrom, feature.start, feature.stop), featuretype='gene')
        ybase = 0
        ngenes = 0
        for nearby_gene in nearby_genes:
            ngenes += 1
            extent.extend([nearby_gene.start, nearby_gene.stop])
            gene_collection = Gene(
                self.db,
                nearby_gene,
                transcripts=['mRNA'],
                cds=['CDS'],
                utrs=['exon'],
                ybase=ybase,
                color="0.5", picker=5)
            gene_collection.name = nearby_gene.id
            gene_collection.add_to_ax(ax)
            ybase += gene_collection.max_y

        xmin = min(extent)
        xmax = max(extent)
        ymax = ngenes

        # 1% padding seems to work well
        padding = (xmax - xmin) * 0.01
        ax.axis('tight')

        # add lines indicating extent of current feature
        vline_kwargs = dict(color='k', linestyle='--')
        ax.axvline(feature.start, **vline_kwargs)
        ax.axvline(feature.stop, **vline_kwargs)

        # Make a new feature to represent the region plus surrounding genes
        interval = pybedtools.create_interval_from_list(feature.fields)
        interval.start = xmin - padding
        interval.stop = xmax + padding
        interval.strand = '.'
        return interval
947816
def rnaseq_prep_samples(config, run_info_yaml, parallel, dirs, samples):
    
    pipeline = dd.get_in_samples(samples, dd.get_analysis)
    trim_reads_set = any([tz.get_in(["algorithm", "trim_reads"], d) for d in dd.sample_data_iterator(samples)])
    resources = ["picard"]
    needs_trimming = (_is_smallrnaseq(pipeline) or trim_reads_set)
    if needs_trimming:
        resources.append("atropos")
    with prun.start(_wres(parallel, resources),
                    samples, config, dirs, "trimming",
                    max_multicore=1 if not needs_trimming else None) as run_parallel:
        with profile.report("organize samples", dirs):
            samples = run_parallel("organize_samples", [[dirs, config, run_info_yaml,
                                                            [x[0]["description"] for x in samples]]])
            samples = run_parallel("prepare_sample", samples)
        if needs_trimming:
            with profile.report("adapter trimming", dirs):
                if _is_smallrnaseq(pipeline):
                    samples = run_parallel("trim_srna_sample", samples)
                else:
                    samples = run_parallel("trim_sample", samples)
    return samples
715814
def _rnaseq_checkpoints(samples):
    
    checkpoints = {}
    checkpoints["rnaseq"] = True
    checkpoints["vc"] = any([dd.get_variantcaller(d) for d in samples])
    return checkpoints
20023
def truncate_seq_pair(tokens_a, tokens_b, max_num_tokens, rng):
    """"""
    while True:
        total_length = len(tokens_a) + len(tokens_b)
        if total_length <= max_num_tokens:
            break

        trunc_tokens = tokens_a if len(tokens_a) > len(tokens_b) else tokens_b
        assert len(trunc_tokens) >= 1

        # We want to sometimes truncate from the front and sometimes from the
        # back to add more randomness and avoid biases.
        if rng.random() < 0.5:
            del trunc_tokens[0]
        else:
            trunc_tokens.pop()
705163
def _retranslate(seq):
        '''
        
        '''
        if len(seq['vdj_nt']) % 3 != 0:
            trunc = len(seq['vdj_nt']) % 3
            seq['vdj_nt'] = seq['vdj_nt'][:-trunc]
        seq['vdj_aa'] = Seq(seq['vdj_nt'], generic_dna).translate()
241545
def subset_to_genome(in_file, out_file, data):
    
    if not utils.file_uptodate(out_file, in_file):
        contigs = set([x.name for x in ref.file_contigs(dd.get_ref_file(data))])
        with utils.open_gzipsafe(in_file) as in_handle:
            with file_transaction(data, out_file) as tx_out_file:
                with open(tx_out_file, "w") as out_handle:
                    for line in in_handle:
                        parts = line.split()
                        if parts and parts[0] in contigs:
                            out_handle.write(line)
    return out_file
707707
def map_to_precursors_on_fly(seqs, names, loci, args):
    """"""
    precursor = precursor_sequence(loci, args.ref).upper()
    dat = dict()
    for s, n in itertools.izip(seqs, names):
        res = pyMatch.Match(precursor, str(s), 1, 3)
        if res > -1:
            dat[n] = [res, res + len(s)]
    logger.debug("mapped in %s: %s out of %s" % (loci, len(dat), len(seqs)))
    return dat