%prog enrich omgfile groups ntaxa > enriched.omg

    Enrich OMG output by pulling genes misses by OMG.
964867
def enrich(args):
    
    p = OptionParser(enrich.__doc__)
    p.add_option("--ghost", default=False, action="store_true",
                 help="Add ghost homologs already used [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) != 3:
        sys.exit(not p.print_help())

    omgfile, groupsfile, ntaxa = args
    ntaxa = int(ntaxa)
    ghost = opts.ghost

    # Get gene pair => weight mapping
    weights = get_edges()
    info = get_info()
    # Get gene => taxon mapping
    info = dict((k, v.split()[5]) for k, v in info.items())

    groups = Grouper()

    fp = open(groupsfile)
    for row in fp:
        members = row.strip().split(",")
        groups.join(*members)

    logging.debug("Imported {0} families with {1} members.".\
                    format(len(groups), groups.num_members))

    seen = set()
    omggroups = Grouper()
    fp = open(omgfile)
    for row in fp:
        genes, idxs = row.split()
        genes = genes.split(",")
        seen.update(genes)
        omggroups.join(*genes)

    nmembers = omggroups.num_members
    logging.debug("Imported {0} OMG families with {1} members.".\
                    format(len(omggroups), nmembers))
    assert nmembers == len(seen)

    alltaxa = set(str(x) for x in range(ntaxa))
    recruited = []
    fp = open(omgfile)
    for row in fp:
        genes, idxs = row.split()
        genes = genes.split(",")
        a = genes[0]

        idxs = set(idxs.split(","))
        missing_taxa = alltaxa - idxs
        if not missing_taxa:
            print(row.rstrip())
            continue

        leftover = groups[a]
        if not ghost:
            leftover = set(leftover) - seen

        if not leftover:
            print(row.rstrip())
            continue

        leftover_sorted_by_taxa = dict((k, \
                             [x for x in leftover if info[x] == k]) \
                                for k in missing_taxa)

        #print genes, leftover
        #print leftover_sorted_by_taxa
        solutions = []
        for solution in product(*leftover_sorted_by_taxa.values()):
            score = sum(weights.get((a, b), 0) for a in solution for b in genes)
            if score == 0:
                continue
            score += sum(weights.get((a, b), 0) for a, b in combinations(solution, 2))
            solutions.append((score, solution))
            #print solution, score

        best_solution = max(solutions) if solutions else None
        if best_solution is None:
            print(row.rstrip())
            continue

        #print "best ==>", best_solution
        best_score, best_addition = best_solution
        genes.extend(best_addition)
        recruited.extend(best_addition)

        genes = sorted([(info[x], x) for x in genes])
        idxs, genes = zip(*genes)

        if ghost:  # decorate additions so it's clear that they were added
            pgenes = []
            for g in genes:
                if g in recruited and g in seen:
                    pgenes.append("|{0}|".format(g))
                else:
                    pgenes.append(g)
            genes = pgenes

        print("\t".join((",".join(genes), ",".join(idxs))))
        if not ghost:
            seen.update(best_addition)

    logging.debug("Recruited {0} new genes.".format(len(recruited)))
234146
def example_generator(self, encoder, tmp_dir, task_id):
    
    filepaths = self.text_filepaths_for_task(tmp_dir, task_id)
    if task_id >= self.num_train_shards:
      # this is dev data - limit the total length.
      max_chars_per_file = self.max_dev_chars // (
          self.num_dev_shards * len(filepaths))
    else:
      max_chars_per_file = None
    tokens = []
    for ftext in self.file_generator(
        filepaths, max_chars_per_file=max_chars_per_file):
      tokens.extend(encoder.encode(ftext))
      pos = 0
      while pos + self.sequence_length <= len(tokens):
        yield {"targets": tokens[pos:pos + self.sequence_length]}
        pos += self.sequence_length
      if pos > 0:
        tokens = tokens[pos:]
    if self.remainder_policy == "pad":
      if tokens:
        targets = tokens + [0] * (self.sequence_length - len(tokens))
        yield {"targets": targets}
    else:
      assert self.remainder_policy == "drop"
759381
def _group_by_cnv_method(batches):
    
    CnvGroup = collections.namedtuple("CnvGroup", "items, work_dir, access_file, region_file")
    out = []
    groups = collections.defaultdict(list)
    for batch, items in batches.items():
        for data in items:
            work_dir = utils.safe_makedir(os.path.join(dd.get_work_dir(data), "structural", "bins", batch))
            cnv_file = get_base_cnv_regions(data, work_dir, "transcripts100", include_gene_names=False)
            if cnv_file:
                break
        assert cnv_file, ("Did not find coverage regions for batch %s: %s" %
                          (batch, " ".join([dd.get_sample_name(d) for d in items])))
        groups[(cnv_file, dd.get_prep_method(data))].append((items, data, work_dir))
    for (cnv_file, _), cur_group in groups.items():
        group_items = reduce(operator.add, [xs[0] for xs in cur_group])
        access_file = tz.get_in(["config", "algorithm", "callable_regions"], cur_group[0][1])
        out.append(CnvGroup(group_items, cur_group[0][2], access_file, cnv_file))
    return out
265222
def summary(args):
    
    from jcvi.formats.base import DictFile
    from jcvi.utils.cbook import percentage, Registry

    p = OptionParser(summary.__doc__)
    p.add_option("--extra", help="Cross with extra tsv file [default: %default]")
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    frfile, statusfile = args
    status = DictFile(statusfile)
    fp = open(frfile)
    registry = Registry()  # keeps all the tags for any given gene
    for row in fp:
        seqid, gene, tag = row.split()
        if tag == '.':
            registry[gene].append("outside")
        else:
            registry[gene].append("inside")
            if tag[0] == '[':
                registry[gene].append("no_syntenic_model")
                if tag.startswith("[S]"):
                    registry[gene].append("[S]")
                    gstatus = status.get(gene, None)
                    if gstatus == 'complete':
                        registry[gene].append("complete")
                    elif gstatus == 'pseudogene':
                        registry[gene].append("pseudogene")
                    elif gstatus == 'partial':
                        registry[gene].append("partial")
                    else:
                        registry[gene].append("gmap_fail")
                elif tag.startswith("[NS]"):
                    registry[gene].append("[NS]")
                    if "random" in tag or "Scaffold" in tag:
                        registry[gene].append("random")
                    else:
                        registry[gene].append("real_ns")
                elif tag.startswith("[NF]"):
                    registry[gene].append("[NF]")
            else:
                registry[gene].append("syntenic_model")

    inside = registry.count("inside")
    outside = registry.count("outside")
    syntenic = registry.count("syntenic_model")
    non_syntenic = registry.count("no_syntenic_model")
    s = registry.count("[S]")
    ns = registry.count("[NS]")
    nf = registry.count("[NF]")
    complete = registry.count("complete")
    pseudogene = registry.count("pseudogene")
    partial = registry.count("partial")
    gmap_fail = registry.count("gmap_fail")
    random = registry.count("random")
    real_ns = registry.count("real_ns")

    complete_models = registry.get_tag("complete")
    pseudogenes = registry.get_tag("pseudogene")
    partial_deletions = registry.get_tag("partial")

    m = "{0} inside synteny blocks\n".format(inside)
    m += "{0} outside synteny blocks\n".format(outside)
    m += "{0} has syntenic gene\n".format(syntenic)
    m += "{0} lack syntenic gene\n".format(non_syntenic)
    m += "{0} has sequence match in syntenic location\n".format(s)
    m += "{0} has sequence match in non-syntenic location\n".format(ns)
    m += "{0} has sequence match in un-ordered scaffolds\n".format(random)
    m += "{0} has sequence match in real non-syntenic location\n".format(real_ns)
    m += "{0} has no sequence match\n".format(nf)
    m += "{0} syntenic sequence - complete model\n".format(percentage(complete, s))
    m += "{0} syntenic sequence - partial model\n".format(percentage(partial, s))
    m += "{0} syntenic sequence - pseudogene\n".format(percentage(pseudogene, s))
    m += "{0} syntenic sequence - gmap fail\n".format(percentage(gmap_fail, s))
    print(m, file=sys.stderr)

    aa = ["complete_models", "partial_deletions", "pseudogenes"]
    bb = [complete_models, partial_deletions, pseudogenes]
    for a, b in zip(aa, bb):
        fw = open(a, "w")
        print("\n".join(b), file=fw)
        fw.close()

    extra = opts.extra
    if extra:
        registry.update_from(extra)

    fp.seek(0)
    fw = open("registry", "w")
    for row in fp:
        seqid, gene, tag = row.split()
        ts = registry[gene]
        print("\t".join((seqid, gene, tag, "-".join(ts))), file=fw)
    fw.close()

    logging.debug("Registry written.")
659823
private NodeMetadata enrich(NodeMetadata nodeMetadata, Template template) {
    final NodeMetadataBuilder nodeMetadataBuilder = NodeMetadataBuilder
        .fromNodeMetadata(nodeMetadata);

    if (nodeMetadata.getHardware() == null) {
      nodeMetadataBuilder.hardware(template.getHardware());
    }

    if (nodeMetadata.getImageId() == null) {
      nodeMetadataBuilder.imageId(template.getImage().getId());
    }

    if (nodeMetadata.getLocation() == null) {
      nodeMetadataBuilder.location(template.getLocation());
    }

    return nodeMetadataBuilder.build();
  }
1323543
def run_all2(protgroup, memornot, subsequences, base_outdir,
            protgroup_dict, protein_feathers_dir, date, errfile, impute_counts=True,
            cutoff_num_proteins=0, core_only_genes=None,
            length_filter_pid=.8, remove_correlated_feats=True,
            force_rerun_counts=False, force_rerun_percentages=False, force_rerun_pca=False):
    """"""
    import ssbio.utils

    # Need to set multiprocessing limit for scipy/numpy stuff if parallelizing anything
    import os
    os.environ['OMP_NUM_THREADS'] = '1'

    # First, filter down the protein group to the membrane/nonmembrane definition
    prots_filtered_feathers = get_protein_feather_paths(protgroup=protgroup, memornot=memornot,
                                                        protgroup_dict=protgroup_dict,
                                                        protein_feathers_dir=protein_feathers_dir,
                                                        core_only_genes=core_only_genes)
    num_proteins = len(prots_filtered_feathers)
    if num_proteins <= cutoff_num_proteins:
        return

    # Make output directories
    protscale = 'proteome_unscaled'
    outdir_d0 = ssbio.utils.make_dir(op.join(base_outdir, protscale))
    outdir_d1 = ssbio.utils.make_dir(op.join(outdir_d0, '-'.join(memornot)))
    outdir_final = ssbio.utils.make_dir(op.join(outdir_d1, '-'.join(protgroup)))

    if impute_counts:
        big_strain_counts_df = get_proteome_counts_impute_missing(prots_filtered_feathers=prots_filtered_feathers,
                                                                  outpath=op.join(outdir_final,
                                                                                  '{}-subsequence_proteome_IMP.fthr'.format(
                                                                                      date)),
                                                                  length_filter_pid=length_filter_pid,
                                                                  force_rerun=force_rerun_counts)

        big_strain_percents_df = get_proteome_percentages(counts_df=big_strain_counts_df,
                                                          outpath=op.join(outdir_final,
                                                                          '{}-subsequence_proteome_perc_IMP.fthr'.format(
                                                                                  date)),
                                                          force_rerun=force_rerun_percentages)
        pca_pickle = op.join(outdir_final, '{}-subsequence_pca.pckl'.format(date))
    # Divide by totals to get percentages in a new dataframe
    else:
        try:
            big_strain_percents_df = get_proteome_correct_percentages(prots_filtered_feathers=prots_filtered_feathers,
                                                                      outpath=op.join(outdir_final,
                                                                                      '{}-subsequence_proteome_perc_AVG.fthr'.format(
                                                                                              date)),
                                                                      length_filter_pid=length_filter_pid,
                                                                      force_rerun=force_rerun_percentages)
            pca_pickle = op.join(outdir_final, '{}-subsequence_pca_AVG.pckl'.format(date))
        except:
            with open(errfile, "a") as myfile:
                myfile.write('PERCENTAGES ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return



    if ssbio.utils.force_rerun(flag=force_rerun_pca, outfile=pca_pickle):

        # Then, get filters for rows of the loaded feathers for interested subsequences
        keep_subsequences = get_interested_subsequences(subsequences=subsequences)

        # Some numbers: number of features
        num_feats = len(big_strain_percents_df)

        # Make an unwieldy title
        big_title = 'LOC={0}; PROTGROUP={1};\n' \
                    'NUMPROTS={2}; NUMFEATS={3}'.format('-'.join(memornot),
                                                                        '-'.join(protgroup),
                                                                        num_proteins,
                                                                        num_feats)

        # Run PCA and make plots
        runner = PCAMultiROS(features_df=big_strain_percents_df, observations_df=pd.DataFrame(), plot_title=big_title)
        try:
            runner.clean_data(keep_features=keep_subsequences, remove_correlated_feats=remove_correlated_feats)
        except:
            with open(errfile, "a") as myfile:
                myfile.write(
                    'CLEAN ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
            return
        # try:
        runner.run_pca()
        # except:
        #     with open(errfile, "a") as myfile:
        #         myfile.write(
        #             'PCA ERR: ' + '-'.join(memornot) + '\t' + '-'.join(protgroup) + "\n")
        #     return
        with open(pca_pickle, 'wb') as f:
            pickle.dump(runner, f)
    else:
        with open(pca_pickle, 'rb') as f:
            runner = pickle.load(f)
458327
def map_all(self, prot_alignment, nucl_sequences):
        
        zipped = itertools.zip_longest(prot_alignment, nucl_sequences)
        for p, n in zipped:
            if p is None:
                raise ValueError("Exhausted protein sequences")
            elif n is None:
                raise ValueError("Exhausted nucleotide sequences")
            yield self.map_alignment(p, n)
999398
def enrich_proteins_with_rnas(graph):
    
    for protein_node in list(graph):
        if not isinstance(protein_node, Protein):
            continue

        if protein_node.variants:
            continue

        rna_node = protein_node.get_rna()
        graph.add_translation(rna_node, protein_node)
223991
def token_generator(self, texts, **kwargs):
        
        for text_idx, text in enumerate(texts):
            if self.lower:
                text = text.lower()
            for char in text:
                yield text_idx, char
505623
def generate_skip_gram_data_set(self, token_list):
        '''
        
        '''
        n_gram_tuple_zip = self.generate_tuple_zip(token_list, 3)
        skip_gram_list = []
        for pre, point, post in n_gram_tuple_zip:
            skip_gram_list.append((point, pre))
            skip_gram_list.append((point, post))
        return zip(skip_gram_list)
989057
def fold_enrichment(self):
        """"""
        return self.k / (self.K*(self.cutoff/float(self.N)))
457691
def liftover(args):
    
    p = OptionParser(liftover.__doc__)
    p.add_option("--checkvalid", default=False, action="store_true",
                help="Check minscore, period and length")
    opts, args = p.parse_args(args)

    if len(args) != 2:
        sys.exit(not p.print_help())

    refbed, fastafile = args
    genome = pyfasta.Fasta(fastafile)
    edits = []
    fp = open(refbed)
    for i, row in enumerate(fp):
        s = STRLine(row)
        seq = genome[s.seqid][s.start - 1: s.end].upper()
        s.motif = get_motif(seq, len(s.motif))
        s.fix_counts(seq)
        if opts.checkvalid and not s.is_valid():
            continue
        edits.append(s)
        if i % 10000 == 0:
            print(i, "lines read", file=sys.stderr)

    edits = natsorted(edits, key=lambda x: (x.seqid, x.start))
    for e in edits:
        print(str(e))
266211
def seqingroups(groups,seq):
  ''
  if not (groups and seq): return None
  bestmatch=None,None,0
  if any(len(g)<2 for g in groups): raise ValueError('some subgroups have length < 2')
  for i,g in filter(lambda x:x[1][0],enumerate(groups)): # i.e. we're only interested in groups with common elements
    # begin starts at 0 so begin+1 starts at 1. (first elt of each group is the bool indicator)
    begin=0
    while 1:
      try: begin=g.index(seq[0],begin+1)
      except ValueError: break
      jmax=min(len(g)-begin,len(seq))
      for j in range(jmax):
        if g[begin+j]!=seq[j]: break
      else: j+=1 # so matchlen works below
      matchlen=min(j,jmax)
      if matchlen<bestmatch[2]: continue
      bestmatch=[i,begin,matchlen] # note: begin is an offset including the initial bool
  return bestmatch if bestmatch[2] else None
155014
static public void addGmodGFF3GeneFeatures(LinkedHashMap<String, ChromosomeSequence> chromosomeSequenceList, FeatureList listGenes) throws Exception {


		// key off mRNA as being a known feature that may or may not have a parent gene


		FeatureList mRNAFeatures = listGenes.selectByType("mRNA");
		LinkedHashMap<String,FeatureList> featureIDHashMap = FeatureHelper.buildFeatureAtrributeIndex("ID", listGenes);
		LinkedHashMap<String,FeatureList> featureParentHashMap = FeatureHelper.buildFeatureAtrributeIndex("Parent", listGenes);

		for (FeatureI f : mRNAFeatures) {
			String geneID;
			String geneNote = null;
			String geneSource = null;
			String sequenceName = null;
			ChromosomeSequence seq = null;
			GeneSequence geneSequence = null;

			Feature mRNAFeature = (Feature) f;
			String mRNAID = mRNAFeature.getAttribute("ID");
			String mRNAsource = mRNAFeature.source();
			String mRNANote = mRNAFeature.getAttribute("Note");
			String mRNAParent = mRNAFeature.getAttribute("Parent");
			if (mRNAParent != null && mRNAParent.length() > 0) {
			   // FeatureList geneFeatureList = listGenes.selectByAttribute("ID", mRNAParent);
				FeatureList geneFeatureList = featureIDHashMap.get(mRNAParent);
				Feature geneFeature = (Feature) geneFeatureList.get(0);
				geneID = geneFeature.getAttribute("ID");
				geneNote = geneFeature.getAttribute("Note");
				geneSource = geneFeature.source();
				sequenceName = geneFeature.seqname();

				//
			} else {
				//deal with cases where no parent gene is given
				geneID = mRNAID;
				geneSource = mRNAsource;
				sequenceName = mRNAFeature.seqname();
			}

			seq = chromosomeSequenceList.get(sequenceName);

			AccessionID geneAccessionID = new AccessionID(geneID);

		  //  FeatureList mRNAChildren = listGenes.selectByAttribute("Parent", mRNAID);
			FeatureList mRNAChildren = featureParentHashMap.get(mRNAID);
			FeatureList cdsFeatures = mRNAChildren.selectByType("CDS");
			FeatureI feature = cdsFeatures.get(0);
			Strand strand = Strand.POSITIVE;

			if (feature.location().isNegative()) {
				strand = Strand.NEGATIVE;
			}
			cdsFeatures = cdsFeatures.sortByStart();







			//String seqName = feature.seqname();
			FeatureI startCodon = null;
			FeatureI stopCodon = null;
			Integer startCodonBegin = null;
			Integer stopCodonEnd = null;
			String startCodonName = "";
			String stopCodonName = "";
			FeatureList startCodonList = mRNAChildren.selectByType("five_prime_UTR");
			if (startCodonList != null && startCodonList.size() > 0) {
				startCodon = startCodonList.get(0);
				if (strand == Strand.NEGATIVE) {
					startCodonBegin = startCodon.location().bioEnd();
				} else {
					startCodonBegin = startCodon.location().bioStart();
				}
				startCodonName = startCodon.getAttribute("ID");
			}

			FeatureList stopCodonList = mRNAChildren.selectByType("three_prime_UTR");

			if (stopCodonList != null && stopCodonList.size() > 0) {
				stopCodon = stopCodonList.get(0);
				if (strand == Strand.NEGATIVE) {
					stopCodonEnd = stopCodon.location().bioStart();
				} else {
					stopCodonEnd = stopCodon.location().bioEnd();
				}
				stopCodonName = stopCodon.getAttribute("ID");

			}




			if (startCodonBegin == null) {
				if (strand == Strand.NEGATIVE) {
					FeatureI firstFeature = cdsFeatures.get(0);

					startCodonBegin = firstFeature.location().bioEnd();
				} else {
					FeatureI firstFeature = cdsFeatures.get(0);

					startCodonBegin = firstFeature.location().bioStart();
				}
			}

			if (stopCodonEnd == null) {
				if (strand == Strand.NEGATIVE) {
					FeatureI lastFeature = cdsFeatures.get(cdsFeatures.size() - 1);
					stopCodonEnd = lastFeature.location().bioStart();
				} else {
					FeatureI lastFeature = cdsFeatures.get(cdsFeatures.size() - 1);
					stopCodonEnd = lastFeature.location().bioEnd();
				}
			}
			//for gtf ordering can be strand based so first is last and last is first
			if (startCodonBegin > stopCodonEnd) {
				int temp = startCodonBegin;
				startCodonBegin = stopCodonEnd;
				stopCodonEnd = temp;
			}



			AccessionID transcriptAccessionID = new AccessionID(mRNAID);
			geneSequence = seq.getGene(geneID);
			if (geneSequence == null) {
				geneSequence = seq.addGene(geneAccessionID, startCodonBegin, stopCodonEnd, strand);
				geneSequence.setSource(geneSource);
				if (geneNote != null && geneNote.length() > 0) {
					geneSequence.addNote(geneNote);
				}
			} else {

				if (startCodonBegin < geneSequence.getBioBegin()) {
					geneSequence.setBioBegin(startCodonBegin);
				}
				if (stopCodonEnd > geneSequence.getBioBegin()) {
					geneSequence.setBioEnd(stopCodonEnd);
				}

			}
			TranscriptSequence transcriptSequence = geneSequence.addTranscript(transcriptAccessionID, startCodonBegin, stopCodonEnd);
			transcriptSequence.setSource(mRNAsource);
			if (mRNANote != null && mRNANote.length() > 0) {
				transcriptSequence.addNote(mRNANote);

			}
			if (startCodon != null) {
				if (startCodonName == null || startCodonName.length() == 0) {
					startCodonName = mRNAID + "-start_codon-" + startCodon.location().bioStart() + "-" + startCodon.location().bioEnd();
				}
				transcriptSequence.addStartCodonSequence(new AccessionID(startCodonName), startCodon.location().bioStart(), startCodon.location().bioEnd());
			}
			if (stopCodon != null) {
				if (stopCodonName == null || stopCodonName.length() == 0) {
					stopCodonName = mRNAID + "-stop_codon-" + stopCodon.location().bioStart() + "-" + stopCodon.location().bioEnd();
				}
				transcriptSequence.addStopCodonSequence(new AccessionID(stopCodonName), stopCodon.location().bioStart(), stopCodon.location().bioEnd());
			}

			for (FeatureI cdsFeature : cdsFeatures) {
				Feature cds = (Feature) cdsFeature;
				String cdsNote = cdsFeature.getAttribute("Note");
				String cdsSource = cds.source();
				String cdsName = cds.getAttribute("ID");
				if (cdsName == null || cdsName.length() == 0) {
					cdsName = mRNAID + "-cds-" + cds.location().bioStart() + "-" + cds.location().bioEnd();
				}
				AccessionID cdsAccessionID = new AccessionID(cdsName);
				ExonSequence exonSequence = geneSequence.addExon(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd());
				exonSequence.setSource(cdsSource);
				if (cdsNote != null && cdsNote.length() > 0) {
					exonSequence.addNote(cdsNote);
				}
				transcriptSequence.addCDS(cdsAccessionID, cdsFeature.location().bioStart(), cdsFeature.location().bioEnd(), cds.frame());
			}
			geneSequence.addIntronsUsingExons();

		}

	}
603006
def _balanced_strand_gen(base_aligns, limit):
    ''''''
    predicate=lambda align: align.is_reverse
    gen1, gen2 = itertools.tee((predicate(align), align) for align in base_aligns)
    for_strand_gen = (align for pred, align in gen1 if not pred)
    rev_strand_gen = (align for pred, align in gen2 if pred)
    combined_gen = iter_zip(for_strand_gen, rev_strand_gen)
    interleaved_gen = (align for for_rev_aligns in combined_gen for align in for_rev_aligns)
    return itertools.islice(interleaved_gen, limit)
557515
def stack(args):
    
    p = OptionParser(stack.__doc__)
    p.add_option("--top", default=10, type="int",
                 help="Draw the first N chromosomes [default: %default]")
    p.add_option("--stacks",
                 default="Exons,Introns,DNA_transposons,Retrotransposons",
                 help="Features to plot in stackplot [default: %default]")
    p.add_option("--switch",
                 help="Change chr names based on two-column file [default: %default]")
    add_window_options(p)
    opts, args, iopts = p.set_image_options(args, figsize="8x8")

    if len(args) != 1:
        sys.exit(not p.print_help())

    fastafile, = args
    top = opts.top
    window, shift, subtract, merge = check_window_options(opts)
    switch = opts.switch
    if switch:
        switch = DictFile(opts.switch)

    stacks = opts.stacks.split(",")
    bedfiles = get_beds(stacks)
    binfiles = get_binfiles(bedfiles, fastafile, shift,
                            subtract=subtract, merge=merge)

    sizes = Sizes(fastafile)
    s = list(sizes.iter_sizes())[:top]
    maxl = max(x[1] for x in s)
    margin = .08
    inner = .02   # y distance between tracks

    pf = fastafile.rsplit(".", 1)[0]
    fig = plt.figure(1, (iopts.w, iopts.h))
    root = fig.add_axes([0, 0, 1, 1])

    # Gauge
    ratio = draw_gauge(root, margin, maxl)

    # Per chromosome
    yinterval = (1 - 2 * margin) / (top + 1)
    xx = margin
    yy = 1 - margin
    for chr, clen in s:
        yy -= yinterval
        xlen = clen / ratio
        cc = chr
        if "_" in chr:
            ca, cb = chr.split("_")
            cc = ca[0].upper() + cb

        if switch and cc in switch:
            cc = "\n".join((cc, "({0})".format(switch[cc])))

        root.add_patch(Rectangle((xx, yy), xlen, yinterval - inner, color=gray))
        ax = fig.add_axes([xx, yy, xlen, yinterval - inner])

        nbins = clen / shift
        if clen % shift:
            nbins += 1

        stackplot(ax, binfiles, nbins, palette, chr, window, shift)
        root.text(xx - .04, yy + .5 * (yinterval - inner), cc, ha="center", va="center")

        ax.set_xlim(0, nbins)
        ax.set_ylim(0, 1)
        ax.set_axis_off()

    # Legends
    yy -= yinterval
    xx = margin
    for b, p in zip(bedfiles, palette):
        b = b.rsplit(".", 1)[0].replace("_", " ")
        b = Registration.get(b, b)

        root.add_patch(Rectangle((xx, yy), inner, inner, color=p, lw=0))
        xx += 2 * inner
        root.text(xx, yy, b, size=13)
        xx += len(b) * .012 + inner

    root.set_xlim(0, 1)
    root.set_ylim(0, 1)
    root.set_axis_off()

    image_name = pf + "." + iopts.format
    savefig(image_name, dpi=iopts.dpi, iopts=iopts)
286870
protected void enrichWithGenerics(Set<BioPAXElement> seed, Set<BioPAXElement> all)
	{
		Set addition;
		if (type == Type.TO_GENERAL) addition = access(upperGenAcc, seed, all);
		else addition = access(lowerGenAcc, seed, all);
		
		all.addAll(addition);
		seed.addAll(addition);

		enrichWithCM(seed, all);
	}