Guess species abundance file path produced by Mothur
509214
def _derive_species_abundance_path(self):
        """"""
        base, ext = path.splitext(self._input_filename)
        return '%s.unique.%s.sabund' % (base, self.__get_method_abbrev())
909589
def find_enriched(sample_entities=None,
                  background_entities=None,
                  object_category=None,
                  **kwargs):

    
    if sample_entities is None:
        sample_entites = []

    (sample_counts, sample_results) = get_counts(entities=sample_entities,
                                                 object_category=object_category,
                                                 min_count=2,
                                                 **kwargs)
    print(str(sample_counts))

    sample_fcs = sample_results['facet_counts']
    taxon_count_dict = sample_fcs[M.SUBJECT_TAXON]

    taxon=None
    for (t,tc) in taxon_count_dict.items():
        # TODO - throw error if multiple taxa
        taxon = t

    if background_entities is None:
        objects = list(sample_counts.keys())
        print("OBJECTS="+str(objects))
        background_entities = get_background(objects, taxon, object_category)

    # TODO: consider caching
    (bg_counts,_) = get_counts(entities=background_entities,
                               object_category=object_category,
                               **kwargs)

    sample_n = len(sample_entities) # TODO - annotated only?
    pop_n = len(background_entities)
    # adapted from goatools
    for (sample_termid,sample_count) in sample_counts.items():
        pop_count = bg_counts[sample_termid]

        # https://en.wikipedia.org/wiki/Fisher's_exact_test
        #              Cls  NotCls
        # study/sample [a,      b]
        # rest of ref  [c,      d]
        #              
        a = sample_count
        b = sample_n - sample_count
        c = pop_count - sample_count
        d = pop_n - pop_count - b
        print("ABCD="+str((sample_termid,a,b,c,d,sample_n)))
        _, p_uncorrected = sp.stats.fisher_exact( [[a, b], [c, d]])
        print("P="+str(p_uncorrected))
774691
private static int countWildcardPhraseFrequencies(File dir, String pattern)
        throws Exception {
        
        File[] files = dir.listFiles();

        int total = 0;

        for (int i = 0; i < files.length; i++) {
            File f = files[i];
            if (f.isDirectory()) {
                total += countWildcardPhraseFrequencies(f, pattern);
            } else if (f.getName().endsWith(".txt")) {
                Scanner sc = new Scanner(f);
                while (sc.hasNext()) {
                    String line = sc.nextLine();
                    if (line.matches(pattern)) {
                        total++;
                    }
                }
            }
        }
        return total;
    }
215964
def handle_program_options():
    """"""
    parser = argparse.ArgumentParser(description="Calculate the alpha diversity\
                                     of a set of samples using one or more \
                                     metrics and output a kernal density \
                                     estimator-smoothed histogram of the \
                                     results.")
    parser.add_argument("-m", "--map_file",
                        help="QIIME mapping file.")
    parser.add_argument("-i", "--biom_fp",
                        help="Path to the BIOM table")
    parser.add_argument("-c", "--category",
                        help="Specific category from the mapping file.")
    parser.add_argument("-d", "--diversity", default=["shannon"], nargs="+",
                        help="The alpha diversity metric. Default \
                             value is 'shannon', which will calculate the Shannon\
                             entropy. Multiple metrics can be specified (space separated).\
                             The full list of metrics is available at:\
                             http://scikit-bio.org/docs/latest/generated/skbio.diversity.alpha.html.\
                             Beta diversity metrics will be supported in the future.")
    parser.add_argument("--x_label", default=[None], nargs="+",
                        help="The name of the diversity metric to be displayed on the\
                        plot as the X-axis label. If multiple metrics are specified,\
                        then multiple entries for the X-axis label should be given.")
    parser.add_argument("--color_by",
                        help="A column name in the mapping file containing\
                              hexadecimal (#FF0000) color values that will\
                              be used to color the groups. Each sample ID must\
                              have a color entry.")
    parser.add_argument("--plot_title", default="",
                        help="A descriptive title that will appear at the top \
                        of the output plot. Surround with quotes if there are\
                        spaces in the title.")
    parser.add_argument("-o", "--output_dir", default=".",
                        help="The directory plots will be saved to.")
    parser.add_argument("--image_type", default="png",
                        help="The type of image to save: png, svg, pdf, eps, etc...")
    parser.add_argument("--save_calculations",
                        help="Path and name of text file to store the calculated "
                        "diversity metrics.")
    parser.add_argument("--suppress_stats", action="store_true", help="Do not display "
                        "significance testing results which are shown by default.")
    parser.add_argument("--show_available_metrics", action="store_true",
                        help="Supply this parameter to see which alpha diversity metrics "
                             " are available for usage. No calculations will be performed"
                             " if this parameter is provided.")
    return parser.parse_args()
1121658
def _derive_rank_abundance_path(self):
        """"""
        base, ext = path.splitext(self._input_filename)
        return '%s.unique.%s.rabund' % (base, self.__get_method_abbrev())
152015
function(options, filePaths, results, numResults) {
		var content = "{\n\t\"numResults\": " + numResults + ",\n" + 
			"\t\"creationDate\": \"" + _getISODateString() + "\",\n" + 
			"\t\"results\": {\n";

		var group = [];
		for (var file in results) {
			var groupStr = "\t\t\"" + file + "\": [\n";

			var matchGroup = [];
			for (var i=0; i<results[file].length; i++) {
				matchGroup.push("\t\t\t{\n" + 
					"\t\t\t\t\"line\": " + results[file][i].line + ",\n" + 
					"\t\t\t\t\"match\": " + "\"" + _cleanStr(results[file][i].match) + "\"" + 
					"\n\t\t\t}");
			}
			groupStr += matchGroup.join(",\n") + "\n";
			groupStr += "\t\t]"
			group.push(groupStr);
		}
		content += group.join(",\n");
		content += "\n\t}"

		if (options.outputExaminedFiles) {
			content += ",\n\t\"examinedFiles\": [\n";
			var files = [];
			for (var i=0; i<filePaths.length; i++) {
				files.push("\t\t\"" + _cleanStr(filePaths[i]) + "\"");
			}
			content += files.join(",\n");
			content += "\n\t]";
		}

		content += "\n}";

		return content;
	}
1167286
def dnld_gafs(species_list, prt=sys.stdout, loading_bar=True):
    """"""
    # Example GAF files in  http://current.geneontology.org/annotations/:
    #   http://current.geneontology.org/annotations/mgi.gaf.gz
    #   http://current.geneontology.org/annotations/fb.gaf.gz
    #   http://current.geneontology.org/annotations/goa_human.gaf.gz
    http = "http://current.geneontology.org/annotations"
    # There are two filename patterns for gene associations on geneontology.org
    fin_gafs = []
    cwd = os.getcwd()
    for species_txt in species_list: # e.g., goa_human mgi fb
        gaf_base = '{ABC}.gaf'.format(ABC=species_txt) # goa_human.gaf
        gaf_cwd = os.path.join(cwd, gaf_base) # {CWD}/goa_human.gaf
        wget_cmd = "{HTTP}/{GAF}.gz".format(HTTP=http, GAF=gaf_base)
        dnld_file(wget_cmd, gaf_cwd, prt, loading_bar)
        fin_gafs.append(gaf_cwd)
    return fin_gafs
1226685
def find_clades(self):
        ""
        # index names from the first tree
        ndict = {j: i for i, j in enumerate(self.names)}
        namedict = {i: j for i, j in enumerate(self.names)}

        # store counts
        clade_counts = {}
        for tidx, ncopies in self.treedict.items():
            
            # testing on unrooted trees is easiest but for some reason slow
            ttree = self.treelist[tidx].unroot()

            # traverse over tree
            for node in ttree.treenode.traverse('preorder'):
                bits = np.zeros(len(ttree), dtype=np.bool_)
                for child in node.iter_leaf_names():
                    bits[ndict[child]] = True

                # get bit string and its reverse
                bitstring = bits.tobytes()
                revstring = np.invert(bits).tobytes()

                # add to clades first time, then check for inverse next hits
                if bitstring in clade_counts:
                    clade_counts[bitstring] += ncopies
                else:
                    if revstring not in clade_counts:
                        clade_counts[bitstring] = ncopies
                    else:
                        clade_counts[revstring] += ncopies

        # convert to freq
        for key, val in clade_counts.items():
            clade_counts[key] = val / float(len(self.treelist))

        ## return in sorted order
        self.namedict = namedict
        self.clade_counts = sorted(
            clade_counts.items(),
            key=lambda x: x[1],
            reverse=True)
96210
def _get_loftee(data):
    
    ancestral_file = tz.get_in(("genome_resources", "variation", "ancestral"), data)
    if not ancestral_file or not os.path.exists(ancestral_file):
        ancestral_file = "false"
    vep = config_utils.get_program("vep", data["config"])
    args = ["--plugin", "LoF,human_ancestor_fa:%s,loftee_path:%s" %
            (ancestral_file, os.path.dirname(os.path.realpath(vep)))]
    return args
489321
def orthologs_iterator(fo, version):
    """"""

    species_list = config["bel_resources"].get("species_list", [])

    fo.seek(0)
    with gzip.open(fo, "rt") as f:
        for line in f:
            edge = json.loads(line)
            if "metadata" in edge:
                source = edge["metadata"]["source"]
                continue

            if "ortholog" in edge:
                edge = edge["ortholog"]
                subj_tax_id = edge["subject"]["tax_id"]
                obj_tax_id = edge["object"]["tax_id"]

                # Skip if species not listed in species_list
                if species_list and subj_tax_id and subj_tax_id not in species_list:
                    continue
                if species_list and obj_tax_id and obj_tax_id not in species_list:
                    continue

                # Converted to ArangoDB legal chars for _key
                subj_key = arangodb.arango_id_to_key(edge["subject"]["id"])
                subj_id = edge["subject"]["id"]

                # Converted to ArangoDB legal chars for _key
                obj_key = arangodb.arango_id_to_key(edge["object"]["id"])
                obj_id = edge["object"]["id"]

                # Subject node
                yield (
                    arangodb.ortholog_nodes_name,
                    {
                        "_key": subj_key,
                        "name": subj_id,
                        "tax_id": edge["subject"]["tax_id"],
                        "source": source,
                        "version": version,
                    },
                )
                # Object node
                yield (
                    arangodb.ortholog_nodes_name,
                    {
                        "_key": obj_key,
                        "name": obj_id,
                        "tax_id": edge["object"]["tax_id"],
                        "source": source,
                        "version": version,
                    },
                )

                arango_edge = {
                    "_from": f"{arangodb.ortholog_nodes_name}/{subj_key}",
                    "_to": f"{arangodb.ortholog_nodes_name}/{obj_key}",
                    "_key": bel.utils._create_hash(f"{subj_id}>>{obj_id}"),
                    "type": "ortholog_to",
                    "source": source,
                    "version": version,
                }

                yield (arangodb.ortholog_edges_name, arango_edge)
514098
def _contextualise(self):
        ''''''
        # loop through hierarchy identifying unique lineages
        # TODO: gain other contextual information, not just ident
        deja_vues = []
        for rank in reversed(self.taxonomy):
            # return named clades -- '' are ignored
            clades = [e for e in self.hierarchy[rank] if e[1]]
            # print 'Rank: {0} - {1}'.format(rank, len(clades))
            # get unique lineages at this level
            uniques = [e for e in clades if len(e[0]) == 1]
            # removed those already seen
            uniques = [e for e in uniques if e[0][0].ident not in deja_vues]
            # add each to self[ident]['cident']
            for e in uniques:
                ident = e[0][0].ident
                self[ident]['cident'] = e[1]
                deja_vues.append(ident)
1335122
def _lstree(files, dirs):
    """"""
    for f, sha1 in files:
        yield "100644 blob {}\t{}\0".format(sha1, f)

    for d, sha1 in dirs:
        yield "040000 tree {}\t{}\0".format(sha1, d)
1023486
def _prepare_sample(data, run_folder):
    
    want = set(["description", "files", "genome_build", "name", "analysis", "upload", "algorithm"])
    out = {}
    for k, v in data.items():
        if k in want:
            out[k] = _relative_paths(v, run_folder)
    if "algorithm" not in out:
        analysis, algorithm = _select_default_algorithm(out.get("analysis"))
        out["algorithm"] = algorithm
        out["analysis"] = analysis
    description = "%s-%s" % (out["name"], clean_name(out["description"]))
    out["name"] = [out["name"], description]
    out["description"] = description
    return out
934135
def main(argv):
    
    import argparse
    import codecs
    # have to be ready to deal with utf-8 names
    out = codecs.getwriter('utf-8')(sys.stdout)
    description = '''Takes a series of at least 2 OTT ids and reports the OTT of their least inclusive taxonomic ancestor and that taxon's ancestors.'''
    parser = argparse.ArgumentParser(prog='ot-taxo-mrca-to-root', description=description)
    parser.add_argument('ids', nargs='+', type=int, help='OTT IDs')
    args = parser.parse_args(argv)
    id_list = args.ids
    last_id = id_list.pop()
    anc_list = get_taxonomic_ancestor_ids(last_id)
    common_anc = set(anc_list)
    for curr_id in id_list:
        curr_anc_set = set(get_taxonomic_ancestor_ids(curr_id))
        common_anc &= curr_anc_set
        if not common_anc:
            break
    for anc_id in anc_list:
        if anc_id in common_anc:
            out.write('{}\n'.format(anc_id))
517458
def num_species(self):
        
        all_headers = reduce(lambda x, y: set(x) | set(y),
                             (rec.get_names() for rec in self.records))
        return len(all_headers)
1044955
private static Pair<Map<Integer, TermParameter>, Map<Integer, TermParameter>> inferOrthologs(
            Kam kam, KAMStore kAMStore, SpeciesDialect dialect,
            EdgeFilter inferf, Set<KamNode> species,
            Map<KamNode, KamNode> ortho) {

        final List<org.openbel.framework.common.model.Namespace> spl = dialect
                .getSpeciesNamespaces();
        final Set<String> rlocs = constrainedHashSet(spl.size());
        for (final org.openbel.framework.common.model.Namespace n : spl) {
            rlocs.add(n.getResourceLocation());
        }

        Map<Integer, TermParameter> ntp = sizedHashMap(species.size());
        Map<Integer, TermParameter> etp = sizedHashMap(species.size());
        for (final KamNode snode : species) {
            if (snode != null) {
                // XXX term parameter looked up 2x; may impact perf/determinism
                // TODO redesign orthologousNodes / inferOrthologs
                TermParameter p = findParameter(kam, kAMStore, snode, rlocs);

                // recurse incoming connections from species node
                recurseConnections(kam, snode, p, inferf, REVERSE, ortho, ntp, etp);
                // recurse outgoing connections from species node
                recurseConnections(kam, snode, p, inferf, FORWARD, ortho, ntp, etp);
            }
        }

        return new Pair<Map<Integer, TermParameter>, Map<Integer, TermParameter>>(
                ntp, etp);
    }
1288602
def augustus(args):
    
    p = OptionParser(augustus.__doc__)
    p.add_option("--species", default="maize",
                 help="Use species model for prediction")
    p.add_option("--hintsfile", help="Hint-guided AUGUSTUS")
    p.add_option("--nogff3", default=False, action="store_true",
                 help="Turn --gff3=off")
    p.set_home("augustus")
    p.set_cpus()
    opts, args = p.parse_args(args)

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

    fastafile, = args
    cpus = opts.cpus
    mhome = opts.augustus_home
    gff3 = not opts.nogff3
    suffix = ".gff3" if gff3 else ".out"
    cfgfile = op.join(mhome, "config/extrinsic/extrinsic.M.RM.E.W.cfg")

    outdir = mkdtemp(dir=".")
    fs = split([fastafile, outdir, str(cpus)])

    augustuswrap_params = partial(augustuswrap, species=opts.species,
                            gff3=gff3, cfgfile=cfgfile,
                            hintsfile=opts.hintsfile)
    g = Jobs(augustuswrap_params, fs.names)
    g.run()

    gff3files = [x.rsplit(".", 1)[0] + suffix for x in fs.names]
    outfile = fastafile.rsplit(".", 1)[0] + suffix
    FileMerger(gff3files, outfile=outfile).merge()
    shutil.rmtree(outdir)

    if gff3:
        from jcvi.annotation.reformat import augustus as reformat_augustus
        reformat_outfile = outfile.replace(".gff3", ".reformat.gff3")
        reformat_augustus([outfile, "--outfile={0}".format(reformat_outfile)])