' + factor_str + "
"
+
+ return factor_str
+
def default_motifs():
"""Return list of Motif instances from default motif database."""
@@ -1393,14 +1475,8 @@ def parse_motifs(motifs):
motifs : list
List of Motif instances.
"""
- if isinstance(motifs, six.string_types):
- with open(motifs) as f:
- if motifs.endswith("pwm") or motifs.endswith("pfm"):
- motifs = read_motifs(f, fmt="pwm")
- elif motifs.endswith("transfac"):
- motifs = read_motifs(f, fmt="transfac")
- else:
- motifs = read_motifs(f)
+ if isinstance(motifs, str):
+ return read_motifs(motifs)
elif isinstance(motifs, Motif):
motifs = [motifs]
else:
@@ -1426,6 +1502,8 @@ def _add_factors_from_handle(motifs, handle):
m2f_direct = {}
m2f_indirect = {}
for line in open(map_file):
+ if line.startswith("#"):
+ continue
try:
motif, *factor_info = line.strip().split("\t")
if len(factor_info) == 1:
@@ -1438,7 +1516,11 @@ def _add_factors_from_handle(motifs, handle):
except Exception:
pass
+ m2f = pd.read_csv(map_file, sep="\t", comment="#", index_col=0)
+
for motif in motifs:
+ if motif.id in m2f.index:
+ motif.factor_info = m2f.loc[motif.id]
if motif.id in m2f_direct:
motif.factors[DIRECT_NAME] = m2f_direct[motif.id]
if motif.id in m2f_indirect:
@@ -1518,7 +1600,7 @@ def read_motifs(infile=None, fmt="pfm", as_dict=False):
if fmt == "pwm":
fmt = "pfm"
- if infile is None or isinstance(infile, six.string_types):
+ if infile is None or isinstance(infile, str):
infile = pfmfile_location(infile)
with open(infile) as f:
motifs = _read_motifs_from_filehandle(f, fmt)
diff --git a/gimmemotifs/plot.py b/gimmemotifs/plot.py
index d161c138..d27f9d21 100644
--- a/gimmemotifs/plot.py
+++ b/gimmemotifs/plot.py
@@ -4,7 +4,6 @@
# the terms of the MIT License, see the file COPYING included with this
# distribution.
""" Various plotting functions """
-from __future__ import print_function
from PIL import Image
import seaborn as sns
from mpl_toolkits.axes_grid1 import ImageGrid
@@ -33,8 +32,7 @@
def axes_off(ax):
- """Get rid of all axis ticks, lines, etc.
- """
+ """Get rid of all axis ticks, lines, etc."""
ax.set_frame_on(False)
ax.axes.get_yaxis().set_visible(False)
ax.axes.get_xaxis().set_visible(False)
@@ -355,7 +353,7 @@ def _get_motif_tree(tree, data, circle=True, vmin=None, vmax=None):
m = 25 / data.values.max()
for node in t.traverse("levelorder"):
- val = data[[l.name for l in node.get_leaves()]].values.mean()
+ val = data[[leaf.name for leaf in node.get_leaves()]].values.mean()
style = NodeStyle()
style["size"] = 0
diff --git a/gimmemotifs/rank.py b/gimmemotifs/rank.py
index 15076a86..5ae37690 100644
--- a/gimmemotifs/rank.py
+++ b/gimmemotifs/rank.py
@@ -9,6 +9,7 @@
import subprocess as sp
import pandas as pd
import numpy as np
+from scipy.stats import rankdata, norm
try:
from scipy.special import factorial
@@ -77,23 +78,77 @@ def qStuart(r):
return factorial(N) * v[N]
-def rankagg(df, method="stuart"):
- """Return aggregated ranks.
+def _rank_int(series, c=3.0 / 8, stochastic=True):
+ # Based on code by Edward Mountjoy
+ # See: https://github.com/edm1/rank-based-INT
+ """Perform rank-based inverse normal transformation on pandas series.
+ If stochastic is True ties are given rank randomly, otherwise ties will
+ share the same value. NaN values are ignored.
+ Args:
+ param1 (pandas.Series): Series of values to transform
+ param2 (Optional[float]): Constand parameter (Bloms constant)
+ param3 (Optional[bool]): Whether to randomise rank of ties
+
+ Returns:
+ pandas.Series
+ """
+
+ # Check input
+ assert isinstance(series, pd.Series)
+ assert isinstance(c, float)
+ assert isinstance(stochastic, bool)
+
+ # Set seed
+ np.random.seed(123)
+
+ # Take original series indexes
+ orig_idx = series.index
+
+ # Drop NaNs
+ series = series.loc[~pd.isnull(series)]
+
+ # Get ranks
+ if stochastic:
+ # Shuffle by index
+ series = series.loc[np.random.permutation(series.index)]
+ # Get rank, ties are determined by their position in the series (hence
+ # why we randomised the series)
+ rank = rankdata(series, method="ordinal")
+ else:
+ # Get rank, ties are averaged
+ rank = rankdata(series, method="average")
+
+ # Convert numpy array back to series
+ rank = pd.Series(rank, index=series.index)
+
+ # Convert rank to normal distribution
+ transformed = rank.apply(_rank_to_normal, c=c, n=len(rank))
+
+ return transformed[orig_idx]
+
+def _rank_to_normal(rank, c, n):
+ # Standard quantile function
+ x = (rank - c) / (n - 2 * c + 1)
+ return norm.ppf(x)
+
+
+def _rankagg_int(df):
+ # Convert values to ranks
+ df_int = df.apply(_rank_int)
+ # Combine z-score using Stouffer's method
+ df_int = (df_int.sum(1) / np.sqrt(df_int.shape[1])).to_frame()
+ df_int.columns = ["z-score"]
+ return df_int
+
+
+def _rankagg_stuart(df):
+ """
Implementation is ported from the RobustRankAggreg R package
References:
Kolde et al., 2012, DOI: 10.1093/bioinformatics/btr709
Stuart et al., 2003, DOI: 10.1126/science.1087447
-
- Parameters
- ----------
- df : pandas.DataFrame
- DataFrame with values to be ranked and aggregated
-
- Returns
- -------
- pandas.DataFrame with aggregated ranks
"""
rmat = pd.DataFrame(index=df.iloc[:, 0])
@@ -105,3 +160,49 @@ def rankagg(df, method="stuart"):
rmat = rmat.apply(sorted, 1, result_type="expand")
p = rmat.apply(qStuart, 1)
return pd.DataFrame({"score": p}, index=rmat.index)
+
+
+def rankagg(df, method="int_stouffer", include_reverse=True, log_transform=True):
+ """Return aggregated ranks.
+
+ Stuart implementation is ported from the RobustRankAggreg R package
+
+ References:
+ Kolde et al., 2012, DOI: 10.1093/bioinformatics/btr709
+ Stuart et al., 2003, DOI: 10.1126/science.1087447
+
+ Parameters
+ ----------
+ df : pandas.DataFrame
+ DataFrame with values to be ranked and aggregated
+ method : str, optional
+ Either "int_stouffer" or "stuart". The "int_stouffer" method is based on combining z-scores
+ from a inverse normal transform of ranks using Stouffer's method.
+
+ Returns
+ -------
+ pandas.DataFrame with aggregated ranks
+ """
+ method = method.lower()
+ if method not in ["stuart", "int_stouffer"]:
+ raise ValueError("unknown method for rank aggregation")
+
+ if method == "stuart":
+ df_asc = pd.DataFrame()
+ df_desc = pd.DataFrame()
+ for col in df.columns:
+ df_asc[col] = (
+ df.sample(frac=1).sort_values(col, ascending=False).index.values
+ )
+ if include_reverse:
+ df_desc[col] = (
+ df.sample(frac=1).sort_values(col, ascending=True).index.values
+ )
+
+ df_result = -np.log10(_rankagg_stuart(df_asc))
+ if include_reverse:
+ df_result += np.log10(_rankagg_stuart(df_desc))
+
+ return df_result
+ if method == "int_stouffer":
+ return _rankagg_int(df)
diff --git a/gimmemotifs/report.py b/gimmemotifs/report.py
index 239eb441..e41951dc 100644
--- a/gimmemotifs/report.py
+++ b/gimmemotifs/report.py
@@ -16,11 +16,21 @@
import numpy as np
import pandas as pd
from statsmodels.stats.multitest import multipletests
+from pandas.core.indexing import _non_reducing_slice
+from pandas.io.formats.style import Styler
+import seaborn as sns
+from matplotlib import colors
+import matplotlib.pyplot as plt
+
+try:
+ import emoji
+except ImportError:
+ pass
from gimmemotifs.comparison import MotifComparer
from gimmemotifs.fasta import Fasta
from gimmemotifs.motif import read_motifs
-from gimmemotifs.config import MotifConfig, DIRECT_NAME, INDIRECT_NAME
+from gimmemotifs.config import MotifConfig
from gimmemotifs.plot import roc_plot
from gimmemotifs.stats import calc_stats, add_star, write_stats
from gimmemotifs import __version__
@@ -28,6 +38,588 @@
logger = logging.getLogger("gimme.report")
+FACTOR_TOOLTIP = "factors
(direct or predicted)
"
+
+
+def _wrap_html_str(x):
+ if " " not in x:
+ return x
+
+ min_pos, max_pos = 0, len(x)
+ if ">" in x and "[^<>]*<").search(x)
+ min_pos, max_pos = m.start(), m.end()
+
+ positions = [m.start() for m in re.compile(" ").finditer(x)]
+
+ positions = [p for p in positions if min_pos < p < max_pos]
+
+ if len(positions) == 0:
+ return x
+
+ pos = sorted(positions, key=lambda p: abs(p - len(x) / 2))[0]
+ x = x[:pos] + "(direct or predicted)
" + x[pos + 1 :] + return x + + +def relative_luminance(rgba): + """ + Calculate relative luminance of a color. + The calculation adheres to the W3C standards + (https://www.w3.org/WAI/GL/wiki/Relative_luminance) + Parameters + ---------- + color : rgb or rgba tuple + Returns + ------- + float + The relative luminance as a value from 0 to 1 + """ + r, g, b = ( + x / 12.92 if x <= 0.03928 else ((x + 0.055) / 1.055 ** 2.4) for x in rgba[:3] + ) + return 0.2126 * r + 0.7152 * g + 0.0722 * b + + +def contrasting_text_color(color, text_color_threshold=0.408): + dark = relative_luminance(color) < text_color_threshold + text_color = "#f1f1f1" if dark else "#000000" + return text_color + + +class ExtraStyler(Styler): + """ + Extra styles for a DataFrame or Series based on pandas.styler using HTML and CSS. + """ + + loader = jinja2.ChoiceLoader( + [jinja2.FileSystemLoader(MotifConfig().get_template_dir()), Styler.loader] + ) + env = jinja2.Environment(loader=loader) + template = env.get_template("table.tpl") + + def __init__(self, *args, **kwargs): + self._data_todo = [] + self.circle_styles = None + self.palette_styles = None + self.col_heading_style = { + "name": "col_heading", + "props": [("border-bottom", "1px solid #e0e0e0")], + } + super(ExtraStyler, self).__init__(*args, **kwargs) + self.display_data = self.data.copy() + + # self.template = + + self._font = "Nunito Sans" + + @property + def font(self): + return self._font + + @font.setter + def font(self, font_name): + self._font = font_name + + def set_font(self, font_name): + """ + Set the font that will be used. + + Parameters + ---------- + font_name : str + Should be a font name available though the Google Font API. + + Returns + ------- + self : ExtraStyler + + Notes + ----- + ``font_name`` can contain spaces, eg. "Nunito Sans". + + Examples + -------- + >>> df = pd.DataFrame(np.random.randn(4, 2), columns=['a', 'b']) + >>> ExtraStyler(df).font("Roboto) + """ + self.font = font_name + return self + + def _current_index(self, subset): + subset = pd.IndexSlice[:, :] if subset is None else subset + subset = _non_reducing_slice(subset) + selected = self.data.loc[subset] + idx_slice = pd.IndexSlice[ + self.data.index.get_indexer(selected.index), + self.data.columns.get_indexer(selected.columns), + ] + return idx_slice + + def _translate(self): + self._compute_data() + d = super()._translate() + circle_styles = self.circle_styles or [] + palette_styles = self.palette_styles or [] + col_heading_style = self.col_heading_style or [] + d.update( + { + "font": self.font, + "circle_styles": circle_styles, + "palette_styles": palette_styles, + "col_heading_style": col_heading_style, + } + ) + return d + + def _compute_data(self): + r = self + for func, args, kwargs in self._data_todo: + r = func(self)(*args, **kwargs) + r.data = r.display_data + return r + + def _tooltip(self, tip, subset=None, part=None): + subset = pd.IndexSlice[:, :] if subset is None else subset + subset = _non_reducing_slice(subset) + + if part is None: + part = "data" + + if part == "data": + self.display_data.loc[subset] = ( + "
"
+ + self.display_data.loc[subset].astype(str)
+ + "
"
+ )
+ elif part == "columns":
+ idx = self._current_index(subset)[1]
+ rename = dict(
+ zip(
+ self.display_data.columns[idx],
+ ""
+ + self.display_data.columns[idx].astype(str)
+ + "
",
+ )
+ )
+ self.display_data.rename(columns=rename, inplace=True)
+ elif part == "index":
+ idx = self._current_index(subset)[0]
+ rename = dict(
+ zip(
+ self.display_data.index[idx],
+ ""
+ + self.display_data.index[idx].astype(str)
+ + "
",
+ )
+ )
+ self.display_data.rename(index=rename, inplace=True)
+ else:
+ raise ValueError(f"unknown value for part: {part}")
+ return self
+
+ def _wrap_iterable(self, it):
+ return [_wrap_html_str(val) for val in it]
+
+ def _wrap(self, subset=None, axis=0):
+ subset = pd.IndexSlice[:, :] if subset is None else subset
+ subset = _non_reducing_slice(subset)
+
+ if axis in [0, "columns"]:
+ idx = self._current_index(subset)[1]
+ rename = dict(
+ zip(
+ self.display_data.columns[idx],
+ self._wrap_iterable(self.display_data.columns[idx]),
+ )
+ )
+ self.display_data.rename(columns=rename, inplace=True)
+ elif axis in [1, "index"]:
+ idx = self._current_index(subset)[0]
+ rename = dict(
+ zip(
+ self.display_data.index[idx],
+ self._wrap_iterable(self.display_data.index[idx]),
+ )
+ )
+ self.display_data.rename(index=rename, inplace=True)
+ else:
+ raise ValueError(f"unknown value for axis: {axis}")
+ return self
+
+ def _convert_to_image(self, subset=None, height=30):
+ subset = pd.IndexSlice[:, :] if subset is None else subset
+ subset = _non_reducing_slice(subset)
+
+ self.display_data.loc[subset] = (
+ f'"
+ + text
+ + "
"
+ )
+ else:
+ text = self.display_data.iloc[idx].astype(str) if show_text else ""
+ self.display_data.iloc[idx] = (
+ f"" + text + "
"
+ )
+
+ return self
+
+ def add_circle(self, **kwargs):
+ self._data_todo.append((lambda instance: instance._circle, (), kwargs))
+ return self
+
+ def wrap(self, **kwargs):
+ self._data_todo.append((lambda instance: instance._wrap, (), kwargs))
+ return self
+
+ def add_tooltip(self, tip, **kwargs):
+ self._data_todo.append((lambda instance: instance._tooltip, (tip,), kwargs))
+ return self
+
+ def convert_to_image(self, **kwargs):
+ self._data_todo.append(
+ (lambda instance: instance._convert_to_image, (), kwargs)
+ )
+ return self
+
+ def rename(self, columns=None, index=None):
+ self.display_data = self.display_data.rename(columns=columns, index=index)
+ return self
+
+ def _emoji_score(self, series, emoji_str=None, bins=None):
+ if emoji_str is None:
+ emoji_str = ":star:"
+ if bins is None:
+ bins = 3
+
+ if isinstance(bins, int):
+ labels = range(1, bins + 1)
+ else:
+ labels = range(1, len(bins))
+
+ return [
+ emoji.emojize(emoji_str * val, use_aliases=True)
+ for val in pd.cut(series, bins=bins, labels=labels)
+ ]
+
+ def _emoji_scale(self, series, emojis=None, bins=None):
+ emoji_dict = {
+ "thumbs": [":thumbsdown:", ":thumbsup:"],
+ "check": [":cross_mark:", ":white_check_mark:"],
+ "smiley": [
+ ":crying_face:",
+ ":slightly_frowning_face:",
+ ":neutral_face:",
+ ":slightly_smiling_face:",
+ ":grin:",
+ ],
+ "black_square": [
+ ":black_small_square:",
+ ":black_medium_small_square:",
+ ":black_medium_square:",
+ ":black_large_square:",
+ ],
+ "white_square": [
+ ":white_small_square:",
+ ":white_medium_small_square:",
+ ":white_medium_square:",
+ ":white_large_square:",
+ ],
+ }
+
+ if emojis is None:
+ emojis = "smiley"
+
+ if emojis in emoji_dict:
+ labels = emoji_dict[emojis]
+ if bins is None:
+ bins = len(labels)
+
+ return [
+ emoji.emojize(val, use_aliases=True)
+ for val in pd.cut(series, bins=bins, labels=labels)
+ ]
+
+ def emoji_scale(self, subset=None, emojis=None, bins=None, axis=0):
+ subset = pd.IndexSlice[:, :] if subset is None else subset
+ subset = _non_reducing_slice(subset)
+
+ idx = self._current_index(subset=subset)
+
+ result = self.display_data.iloc[idx].apply(
+ self._emoji_scale, axis=axis, result_type="expand", args=(emojis, bins)
+ )
+ self.display_data.iloc[idx] = result.values
+
+ return self.align(subset=subset, location="center", axis=axis)
+
+ def emoji_score(self, subset=None, emoji_str=None, bins=None, axis=0):
+ subset = pd.IndexSlice[:, :] if subset is None else subset
+ subset = _non_reducing_slice(subset)
+
+ idx = self._current_index(subset=subset)
+ result = self.display_data.iloc[idx].apply(
+ self._emoji_score, axis=axis, result_type="expand", args=(emoji_str, bins)
+ )
+ self.display_data.iloc[idx] = result.values
+
+ return self.align(subset=subset, location="left", axis=axis)
+
+ def emojify(self, subset=None):
+ subset = pd.IndexSlice[:, :] if subset is None else subset
+ subset = _non_reducing_slice(subset)
+
+ idx = self._current_index(subset=subset)
+ result = self.display_data.iloc[idx].applymap(emoji.emojize)
+ self.display_data.iloc[idx] = result.values
+
+ return self
+
+ def scaled_background_gradient(
+ self,
+ subset=None,
+ cmap="RdBu_r",
+ low=0,
+ high=0,
+ center_zero=False,
+ vmin=None,
+ vmax=None,
+ ):
+ subset = pd.IndexSlice[:, :] if subset is None else subset
+ subset = _non_reducing_slice(subset)
+
+ vmax = (
+ self.data.loc[subset]
+ .replace({np.inf: np.nan, -np.inf: np.nan})
+ .max(skipna=True)
+ .max()
+ if vmax is None
+ else vmax
+ )
+ vmin = (
+ self.data.loc[subset]
+ .replace({np.inf: np.nan, -np.inf: np.nan})
+ .min(skipna=True)
+ .min()
+ if vmin is None
+ else vmin
+ )
+
+ if center_zero:
+ vmax = max(abs(vmax), abs(vmin))
+ vmin = -vmax
+
+ r = self
+ for col in self.data.loc[subset].columns:
+ r = r.apply(
+ self._background_gradient,
+ subset=pd.IndexSlice[subset[0], col],
+ cmap=cmap,
+ m=vmin,
+ M=vmax,
+ low=low,
+ high=high,
+ )
+
+ return r
+
def get_roc_values(motif, fg_file, bg_file, genome):
"""Calculate ROC AUC values for ROC plots."""
@@ -167,8 +759,8 @@ class ReportMotif(object):
+ os.path.basename(roc_img_file % (motif.id, bg))
+ ".png"
}
- rm.bg[bg][u"roc_img_link"] = {
- u"href": "images/"
+ rm.bg[bg]["roc_img_link"] = {
+ "href": "images/"
+ os.path.basename(roc_img_file % (motif.id, bg))
+ ".png"
}
@@ -253,93 +845,114 @@ def create_denovo_motif_report(
)
-def maelstrom_html_report(outdir, infile, pfmfile=None, threshold=4):
- df = pd.read_table(infile, index_col=0)
- df = df[np.any(abs(df) >= threshold, 1)]
+def motif_to_factor_series(series, pfmfile=None, motifs=None):
+ if motifs is None:
+ motifs = read_motifs(pfmfile, as_dict=True)
- motifs = read_motifs(pfmfile)
+ if isinstance(series, pd.Index):
+ index = series
+ else:
+ index = series.index
- df.rename_axis(None, inplace=True)
- cols = df.columns
+ factors = [motifs[motif].format_factors(html=True) for motif in series]
+ return pd.Series(data=factors, index=index)
+
+
+def motif_to_img_series(series, pfmfile=None, motifs=None, outdir=".", subdir="logos"):
+ if motifs is None:
+ motifs = read_motifs(pfmfile, as_dict=True)
+
+ if not os.path.exists(outdir):
+ os.makedirs(outdir)
+ if not os.path.exists(os.path.join(outdir, subdir)):
+ os.makedirs(os.path.join(outdir, subdir))
+
+ img_series = []
+ for motif in series:
+ if motif not in motifs:
+ raise ValueError(f"Motif {motif} does not occur in motif database")
+ fname = subdir + "/{}.png".format(re.sub(r"[^a-zA-Z0-9\-]+", "_", motif))
+ if not os.path.exists(fname):
+ motifs[motif].plot_logo(fname=os.path.join(outdir, fname))
+ img_series.append(fname)
+
+ if isinstance(series, pd.Index):
+ index = series
+ else:
+ index = series.index
+ return pd.Series(data=img_series, index=index)
- motifs = read_motifs(pfmfile)
- idx = [motif.id for motif in motifs]
- direct = [
- ",".join(sorted(set([x.upper() for x in motif.factors[DIRECT_NAME]])))
- for motif in motifs
- ]
- indirect = [
- ",".join(sorted(set([x.upper() for x in motif.factors[INDIRECT_NAME]])))
- for motif in motifs
- ]
- m2f = pd.DataFrame({DIRECT_NAME: direct, INDIRECT_NAME: indirect}, index=idx)
-
- factor_cols = [DIRECT_NAME, INDIRECT_NAME]
- if True:
- for factor_col in factor_cols:
- f = m2f[factor_col].str.len() > 30
- m2f[factor_col] = (
- ''
- + m2f[factor_col].str.slice(0, 30)
- )
- m2f.loc[f, factor_col] += "(...)"
- m2f[factor_col] += "
"
- df = df.join(m2f)
- df["logo"] = [
- '
'.format(re.sub("[()/]", "_", x))
- for x in list(df.index)
+def maelstrom_html_report(outdir, infile, pfmfile=None, threshold=3):
+
+ # Read the maelstrom text report
+ df = pd.read_table(infile, index_col=0)
+
+ # Columns with maelstrom rank aggregation value
+ value_cols = df.columns[
+ ~df.columns.str.contains("corr") & ~df.columns.str.contains("% with motif")
]
+ # Columns with correlation values
+ corr_cols = df.columns[df.columns.str.contains("corr")]
- if not os.path.exists(outdir + "/logos"):
- os.makedirs(outdir + "/logos")
- for motif in motifs:
- if motif.id in df.index:
- motif.plot_logo(
- fname=outdir + "/logos/{}.png".format(re.sub("[()/]", "_", motif.id))
- )
+ df = df[np.any(abs(df[value_cols]) >= threshold, 1)]
- template_dir = MotifConfig().get_template_dir()
- js = open(
- os.path.join(template_dir, "sortable/sortable.min.js"), encoding="utf-8"
- ).read()
- css = open(
- os.path.join(template_dir, "sortable/sortable-theme-slick.css"),
- encoding="utf-8",
- ).read()
- df = df[factor_cols + ["logo"] + list(cols)]
-
- df_styled = df.style
- absmax = np.max((abs(df[cols].max().max()), abs(df[cols].min().min())))
- target = absmax * 1.75
-
- for col in cols:
- smin = df[col].min()
- smax = df[col].max()
- diff = smax - smin
- low = abs((-target - smin) / diff)
- high = (target - smax) / diff
- df_styled = df_styled.background_gradient(
- cmap="RdBu_r", low=low, high=high, subset=[col]
- )
+ # Add motif logo's
+ df.insert(
+ 0,
+ "logo",
+ motif_to_img_series(df.index, pfmfile=pfmfile, outdir=outdir, subdir="logos"),
+ )
+ # Add factors that can bind to the motif
+ df.insert(0, "factors", motif_to_factor_series(df.index, pfmfile=pfmfile))
- df_styled = df_styled.set_precision(3)
- df_styled = df_styled.set_table_attributes("data-sortable")
- df_styled = df_styled.render()
- df_styled = df_styled.replace(
- "data-sortable", 'class="sortable-theme-slick" data-sortable'
+ rename_columns = {"factors": FACTOR_TOOLTIP}
+
+ df_styled = (
+ ExtraStyler(df)
+ .set_precision(2)
+ .convert_to_image(
+ subset=["logo"],
+ height=30,
+ )
+ .scaled_background_gradient(
+ subset=value_cols, center_zero=True, low=1 / 1.75, high=1 / 1.75
+ )
+ .border(subset=list(value_cols[:1]), location="left")
+ .border(part="columns", location="bottom")
+ .set_table_attributes('class="sortable-theme-slick" data-sortable')
+ .align(subset=list(value_cols), location="center")
+ .set_font("Nunito Sans")
+ .rename(columns=rename_columns)
)
+ if len(corr_cols) > 0:
+ df_styled = (
+ df_styled.wrap(subset=list(corr_cols))
+ .align(subset=list(corr_cols), location="center")
+ .scaled_background_gradient(
+ subset=corr_cols,
+ cmap="PuOr_r",
+ center_zero=True,
+ low=1 / 1.75,
+ high=1 / 1.75,
+ )
+ )
+
+ for col in df.columns:
+ if "% with motif" in col:
+ df_styled = (
+ df_styled.add_circle(subset=[col], cmap="Purples", vmax=100, size=40)
+ .wrap(subset=[col])
+ .align(subset=[col], location="center")
+ .border(subset=[col], location="left")
+ .to_precision_str(subset=[col])
+ )
+
+ df_styled = df_styled.wrap().render()
+
with open(outdir + "/gimme.maelstrom.report.html", "w", encoding="utf-8") as f:
- f.write("\n")
- f.write("\n".format(css))
- f.write("\n")
- f.write("\n")
f.write(df_styled)
- f.write("\n".format(js))
- f.write("\n")
def roc_html_report(
@@ -353,102 +966,104 @@ def roc_html_report(
):
df = pd.read_table(infile, index_col=0)
df.rename_axis(None, inplace=True)
- df["corrected P-value"] = multipletests(df["P-value"], method="fdr_bh")[1]
+
+ motifs = read_motifs(pfmfile, as_dict=True)
+ if use_motifs is not None:
+ motifs = {k: v for k, v in motifs.items() if k in use_motifs}
+ idx = list(motifs.keys())
+ df = df.loc[idx]
+
+ df.insert(2, "corrected P-value", multipletests(df["P-value"], method="fdr_bh")[1])
+ df.insert(3, "-log10 P-value", -np.log10(df["corrected P-value"]))
+ df = df[df["corrected P-value"] <= threshold]
cols = [
+ "factors",
"logo",
- "# matches",
- "# matches background",
- "P-value",
- "log10 P-value",
- "corrected P-value",
+ "% matches input",
+ "%matches background",
+ "-log10 P-value",
"ROC AUC",
"PR AUC",
"Enr. at 1% FPR",
"Recall at 10% FDR",
]
- motifs = read_motifs(pfmfile)
- if use_motifs is not None:
- motifs = [m for m in motifs if m.id in use_motifs]
-
- idx = [motif.id for motif in motifs]
- df = df.loc[idx]
- direct = [",".join(motif.factors[DIRECT_NAME]) for motif in motifs]
- indirect = [",".join(motif.factors[INDIRECT_NAME]) for motif in motifs]
- m2f = pd.DataFrame({DIRECT_NAME: direct, INDIRECT_NAME: indirect}, index=idx)
-
- factor_cols = [DIRECT_NAME, INDIRECT_NAME]
- if True:
- for factor_col in factor_cols:
- f = m2f[factor_col].str.len() > 30
- m2f[factor_col] = (
- ''
- + m2f[factor_col].str.slice(0, 30)
- )
- m2f.loc[f, factor_col] += "(...)"
- m2f[factor_col] += "
"
- df = df.join(m2f)
- cols = factor_cols + cols
-
- df = df[df["corrected P-value"] <= threshold]
-
if link_matches:
df["# matches"] = (
""
+ df["# matches"].astype(str)
+ ""
)
- df["logo"] = [
- ''.format(re.sub(r"[^-_\w]+", "_", x))
- for x in list(df.index)
- ]
+ # Add motif logo's
+ df.insert(
+ 0,
+ "logo",
+ motif_to_img_series(
+ df.index, pfmfile=pfmfile, motifs=motifs, outdir=outdir, subdir="logos"
+ ),
+ )
+ # Add factors that can bind to the motif
+ df.insert(
+ 0, "factors", motif_to_factor_series(df.index, pfmfile=pfmfile, motifs=motifs)
+ )
+
+ rename_columns = {"factors": FACTOR_TOOLTIP}
df = df[cols]
- if not os.path.exists(outdir + "/logos"):
- os.makedirs(outdir + "/logos")
- for motif in motifs:
- if motif.id in df.index:
- motif.plot_logo(
- fname=outdir
- + "/logos/{}.png".format(re.sub(r"[^-_\w]+", "_", motif.id))
- )
bar_cols = [
- "log10 P-value",
+ "% matches input",
+ "%matches background",
+ "-log10 P-value",
"ROC AUC",
"PR AUC",
"Enr. at 1% FPR",
"Recall at 10% FDR",
]
- template_dir = MotifConfig().get_template_dir()
- js = open(
- os.path.join(template_dir, "sortable/sortable.min.js"), encoding="utf-8"
- ).read()
- css = open(
- os.path.join(template_dir, "sortable/sortable-theme-slick.css"),
- encoding="utf-8",
- ).read()
+
+ df["% matches input"] = df["% matches input"].astype(int)
+ df["%matches background"] = df["%matches background"].astype(int)
+ rename_columns = {"factors": FACTOR_TOOLTIP}
+ df = df.sort_values("ROC AUC", ascending=False)
with open(os.path.join(outdir, outname), "w", encoding="utf-8") as f:
- f.write("\n")
- f.write("\n".format(css))
- f.write("\n")
- f.write("\n")
if df.shape[0] > 0:
f.write(
- df.sort_values("ROC AUC", ascending=False)
- .style.bar(bar_cols)
- .set_precision(3)
- .set_table_attributes("data-sortable")
+ ExtraStyler(df)
+ .convert_to_image(
+ subset=["logo"],
+ height=30,
+ )
+ .add_circle(
+ subset=["% matches input", "%matches background"],
+ vmax=100,
+ cmap="Purples",
+ )
+ .scaled_background_gradient(
+ "-log10 P-value", vmin=0, high=0.3, cmap="Reds"
+ )
+ .scaled_background_gradient(
+ "ROC AUC", vmin=0.5, vmax=1, high=0.3, cmap="Reds"
+ )
+ .scaled_background_gradient(
+ "PR AUC", vmin=0, vmax=1, high=0.3, cmap="Reds"
+ )
+ .scaled_background_gradient(
+ "Enr. at 1% FPR", vmin=1, high=0.3, cmap="Reds"
+ )
+ .scaled_background_gradient(
+ "Recall at 10% FDR", vmin=0, vmax=1, high=0.7, cmap="Reds"
+ )
+ .set_precision(2)
+ .set_table_attributes('class="sortable-theme-slick" data-sortable')
+ .wrap(subset=cols)
+ .align(subset=bar_cols, location="center")
+ .rename(columns=rename_columns)
+ .to_precision_str(subset=["% matches input", "%matches background"])
.render()
- .replace("data-sortable", 'class="sortable-theme-slick" data-sortable')
)
else:
- f.write("No enriched motifs found.")
- f.write("\n".format(js))
- f.write("\n")
+ f.write("No enriched motifs found.")
diff --git a/gimmemotifs/scanner.py b/gimmemotifs/scanner.py
index d45d3906..8397fb8f 100644
--- a/gimmemotifs/scanner.py
+++ b/gimmemotifs/scanner.py
@@ -1,11 +1,11 @@
import os
import re
import sys
+from collections import Counter
from functools import partial
from tempfile import mkdtemp, NamedTemporaryFile
import logging
import multiprocessing as mp
-import six
# "hidden" features, in development
try:
@@ -18,7 +18,8 @@
from genomepy import Genome
from diskcache import Cache
import numpy as np
-from scipy.stats import scoreatpercentile
+from sklearn.preprocessing import scale
+import pandas as pd
from gimmemotifs import __version__
from gimmemotifs.background import RandomGenomicFasta, gc_bin_bedfile
@@ -46,7 +47,6 @@ def _pickle_method(m):
# only used when using cache, should not be a requirement
try:
from dogpile.cache import make_region
- from dogpile.cache.api import NO_VALUE
import xxhash
except ImportError:
pass
@@ -250,8 +250,7 @@ def scan_to_file(
zscore=True,
gcnorm=True,
):
- """Scan an inputfile with motifs.
- """
+ """Scan an inputfile with motifs."""
should_close = False
if filepath_or_buffer is None:
fo = sys.stdout
@@ -341,7 +340,7 @@ def scan_to_best_match(
if genome:
s.set_genome(genome)
- if isinstance(motifs, six.string_types):
+ if isinstance(motifs, str):
motifs = read_motifs(motifs)
logger.debug("scanning %s...", fname)
@@ -370,7 +369,9 @@ def parse_threshold_values(motif_file, cutoff):
return threshold
-def scan_sequence(seq, motifs, nreport, scan_rc):
+def scan_sequence(
+ seq, seq_gc_bin, motifs, nreport, scan_rc, motifs_meanstd=None, zscore=False
+):
ret = []
# scan for motifs
@@ -378,36 +379,37 @@ def scan_sequence(seq, motifs, nreport, scan_rc):
if cutoff is None:
ret.append([])
else:
- result = pwmscan(seq, motif.logodds, cutoff, nreport, scan_rc)
+ if zscore:
+ m_mean, m_std = motifs_meanstd[seq_gc_bin][motif.id]
+ result = pwmscan(
+ seq, motif.logodds, motif.pwm_min_score(), nreport, scan_rc
+ )
+ result = [[(row[0] - m_mean) / m_std, row[1], row[2]] for row in result]
+ result = [row for row in result if row[0] >= cutoff]
+ else:
+ result = pwmscan(seq, motif.logodds, cutoff, nreport, scan_rc)
if cutoff <= motif.pwm_min_score() and len(result) == 0:
result = [[motif.pwm_min_score(), 0, 1]] * nreport
- ret.append(result)
-
- # return results
- return ret
-
-
-def scan_region(region, genome, motifs, nreport, scan_rc):
-
- # retrieve sequence
- chrom, start, end = re.split(r"[:-]", region)
- seq = genome[chrom][int(start) : int(end)].seq.upper()
-
- return scan_sequence(seq, motifs, nreport, scan_rc)
+ ret.append(result)
-def scan_seq_mult(seqs, motifs, nreport, scan_rc):
- ret = []
- for seq in seqs:
- result = scan_sequence(seq.upper(), motifs, nreport, scan_rc)
- ret.append(result)
return ret
-def scan_region_mult(regions, genome, motifs, nreport, scan_rc):
+def scan_seq_mult(
+ seqs, seq_gc_bins, motifs, nreport, scan_rc, motifs_meanstd=None, zscore=False
+):
ret = []
- for region in regions:
- result = scan_region(region, genome, motifs, nreport, scan_rc)
+ for seq, seq_gc_bin in zip(seqs, seq_gc_bins):
+ result = scan_sequence(
+ seq.upper(),
+ seq_gc_bin,
+ motifs,
+ nreport,
+ scan_rc,
+ motifs_meanstd=motifs_meanstd,
+ zscore=zscore,
+ )
ret.append(result)
return ret
@@ -529,7 +531,7 @@ class Scanner(object):
def __init__(self, ncpus=None):
self.config = MotifConfig()
- self.threshold = None
+ self._threshold = None
self.genome = None
self.background = None
self.meanstd = {}
@@ -604,19 +606,15 @@ def _meanstd_from_seqs(self, motifs, seqs):
def _threshold_from_seqs(self, motifs, seqs, fpr):
scan_motifs = [(m, m.pwm_min_score()) for m in motifs]
-
table = []
- for x in self._scan_sequences_with_motif(scan_motifs, seqs, 1, True):
- table.append([row[0][0] for row in x])
+ seq_gc_bins = [self.get_seq_bin(seq) for seq in seqs]
+ for gc_bin, result in zip(
+ seq_gc_bins, self._scan_sequences_with_motif(scan_motifs, seqs, 1, True)
+ ):
+ table.append([gc_bin] + [row[0][0] for row in result])
- for (motif, _), scores in zip(scan_motifs, np.array(table).transpose()):
- if len(scores) > 0:
- opt_score = scoreatpercentile(scores, 100 - (100 * fpr))
- yield motif, opt_score # cutoff
- else:
- raise ValueError(
- "Could not determine threshold for motif {}".format(motif)
- )
+ df = pd.DataFrame(table, columns=["gc_bin"] + [m.id for m in motifs])
+ return df
def set_meanstd(self, gc=False):
if not self.background:
@@ -677,8 +675,22 @@ def set_meanstd(self, gc=False):
lock.release()
+ for gc_bin in self.gc_bins:
+ gc_bin = "{:.2f}-{:.2f}".format(*gc_bin)
+ if gc_bin not in self.meanstd:
+ valid_bins = []
+ for b in self.gc_bins:
+ bstr = "{:.2f}-{:.2f}".format(b[0], b[1])
+ if bstr in self.meanstd:
+ valid_bins.append(((b[0] + b[1]) / 2, bstr))
+
+ v = float(gc_bin.split("-")[1])
+ _, bstr = sorted(valid_bins, key=lambda x: abs(x[0] - v))[0]
+ logger.warn(f"Using {bstr}")
+ self.meanstd[gc_bin] = self.meanstd[bstr]
+
def set_background(
- self, fname=None, genome=None, size=200, nseq=10000, gc=False, gc_bins=None
+ self, fname=None, genome=None, size=200, nseq=None, gc=False, gc_bins=None
):
"""Set the background to use for FPR and z-score calculations.
@@ -705,6 +717,16 @@ def set_background(
size = int(size)
+ if gc_bins is None:
+ if gc:
+ gc_bins = [(0.0, 0.2), (0.8, 1)]
+ for b in np.arange(0.2, 0.799, 0.05):
+ gc_bins.append((b, b + 0.05))
+ else:
+ gc_bins = [(0, 1)]
+ if nseq is None:
+ nseq = max(10000, len(gc_bins) * 1000)
+
if genome and fname:
raise ValueError("Need either genome or filename for background.")
@@ -736,12 +758,6 @@ def set_background(
if not fa:
if gc:
-
- if gc_bins is None:
- gc_bins = [(0.0, 0.2), (0.8, 1)]
- for b in np.arange(0.2, 0.799, 0.05):
- gc_bins.append((b, b + 0.05))
-
with NamedTemporaryFile() as tmp:
logger.info("using {} sequences".format(nseq))
gc_bin_bedfile(
@@ -757,6 +773,12 @@ def set_background(
if gc_bins:
self.gc_bins = gc_bins
+ @property
+ def threshold(self):
+ if self._threshold is None:
+ self.set_threshold()
+ return self._threshold
+
def set_threshold(self, fpr=None, threshold=None, gc=False):
"""Set motif scanning threshold based on background sequences.
@@ -775,6 +797,17 @@ def set_threshold(self, fpr=None, threshold=None, gc=False):
if threshold and fpr:
raise ValueError("Need either fpr or threshold.")
+ if threshold is None and fpr is None:
+ if self.genome:
+ fpr = 0.01
+ logger.info(f"Using default FPR of {fpr}")
+ else:
+ threshold = 0.95
+ logger.info(
+ f"Genome not specified, using default threshold of {threshold}."
+ )
+ logger.info("This is likely not ideal.")
+
if fpr:
fpr = float(fpr)
if not (0.0 < fpr < 1.0):
@@ -783,11 +816,16 @@ def set_threshold(self, fpr=None, threshold=None, gc=False):
if not self.motifs:
raise ValueError("please run set_motifs() first")
- thresholds = {}
motifs = read_motifs(self.motifs)
+ gc_bins = ["{:.2f}-{:.2f}".format(*gc_bin) for gc_bin in self.gc_bins]
if threshold is not None:
- self.threshold = parse_threshold_values(self.motifs, threshold)
+ d = parse_threshold_values(self.motifs, threshold)
+ self._threshold = pd.DataFrame(d, index=[0])
+ self._threshold = self._threshold.join(
+ pd.DataFrame(gc_bins, index=[0] * len(gc_bins), columns=["gc_bin"])
+ )
+ self._threshold = self._threshold.set_index("gc_bin")
return
if not self.background:
@@ -801,36 +839,41 @@ def set_threshold(self, fpr=None, threshold=None, gc=False):
lock.acquire()
with Cache(CACHE_DIR) as cache:
scan_motifs = []
+ self._threshold = None
for motif in motifs:
- k = "{}|{}|{:.4f}".format(motif.hash(), self.background_hash, fpr)
-
- threshold = cache.get(k)
- if threshold is None:
+ k = "{}|{}|{:.4f}|{}".format(
+ motif.hash(), self.background_hash, fpr, ",".join(sorted(gc_bins))
+ )
+ vals = cache.get(k)
+ if vals is None:
scan_motifs.append(motif)
else:
- if np.isclose(threshold, motif.pwm_max_score()):
- thresholds[motif.id] = None
- elif np.isclose(threshold, motif.pwm_min_score()):
- thresholds[motif.id] = 0.0
+ if self._threshold is None:
+ self._threshold = vals.to_frame()
else:
- thresholds[motif.id] = threshold
+ self._threshold[motif.id] = vals
if len(scan_motifs) > 0:
logger.info("determining FPR-based threshold")
- for motif, threshold in self._threshold_from_seqs(
- scan_motifs, seqs, fpr
- ):
- k = "{}|{}|{:.4f}".format(motif.hash(), self.background_hash, fpr)
- cache.set(k, threshold)
- if np.isclose(threshold, motif.pwm_max_score()):
- thresholds[motif.id] = None
- elif np.isclose(threshold, motif.pwm_min_score()):
- thresholds[motif.id] = 0.0
- else:
- thresholds[motif.id] = threshold
+ df = self._threshold_from_seqs(scan_motifs, seqs, fpr).set_index(
+ "gc_bin"
+ )
+ if self._threshold is None:
+ self._threshold = df
+ else:
+ self._threshold = pd.concat((self._threshold, df), axis=1)
+ for motif in scan_motifs:
+ k = "{}|{}|{:.4f}|{}".format(
+ motif.hash(),
+ self.background_hash,
+ fpr,
+ ",".join(sorted(gc_bins)),
+ )
+ cache.set(k, df[motif.id])
lock.release()
- self.threshold_str = "{}_{}_{}".format(fpr, threshold, self.background_hash)
- self.threshold = thresholds
+ self.threshold_str = "{}_{}_{}_{}".format(
+ fpr, threshold, self.background_hash, ",".join(sorted(gc_bins))
+ )
def set_genome(self, genome):
"""
@@ -887,8 +930,11 @@ def best_match(self, seqs, scan_rc=True, zscore=False, gc=False):
yield [m[0] for m in matches]
def get_seq_bin(self, seq):
- useq = seq.upper()
- gc = round((useq.count("G") + useq.count("C")) / len(useq), 2)
+ if len(str(seq)) == 0:
+ gc = 0
+ else:
+ useq = seq.upper()
+ gc = round((useq.count("G") + useq.count("C")) / len(useq), 2)
if gc == 0:
gc = 0.01
for b_start, b_end in self.gc_bins:
@@ -931,16 +977,9 @@ def scan(self, seqs, nreport=100, scan_rc=True, zscore=False, gc=False):
"""
Scan a set of regions or sequences.
"""
-
- if not self.threshold:
- logger.info(
- "Using default threshold of 0.95. " "This is likely not optimal!"
- )
- self.set_threshold(threshold=0.95)
-
seqs = as_fasta(seqs, genome=self.genome)
- it = self._scan_sequences(seqs.seqs, nreport, scan_rc)
+ it = self._scan_sequences(seqs.seqs, nreport, scan_rc, zscore=zscore)
if zscore:
if gc:
@@ -950,165 +989,90 @@ def scan(self, seqs, nreport=100, scan_rc=True, zscore=False, gc=False):
if len(self.meanstd) != 1:
self.set_meanstd(gc=gc)
- gc_seqs = [self.get_seq_bin(seq) for seq in seqs.seqs]
-
logger.debug("Scanning")
- for result, gc_seq in zip(it, gc_seqs):
- if zscore:
- zresult = []
- for i, mrow in enumerate(result):
- try:
- m_mean, m_std = self.get_motif_mean_std(
- gc_seq, self.motif_ids[i]
- )
- except Exception:
- print(self.meanstd)
- print(gc_seq, self.motif_ids[i])
- raise
- mrow = [((x[0] - m_mean) / m_std, x[1], x[2]) for x in mrow]
- zresult.append(mrow)
- yield zresult
- else:
- yield result
+ for result in it:
+ yield result
- def _scan_regions(self, regions, nreport, scan_rc):
- genome = self.genome
- motif_file = self.motifs
- motif_digest = self.checksum.get(motif_file, None)
+ def get_gc_thresholds(self, seqs, motifs=None, zscore=False):
+ # Simple case, only one threshold
+ if np.all(self.threshold.nunique(axis=0) == 1):
+ return self.threshold.iloc[0].to_dict()
- # determine which regions are not in the cache
- scan_regions = regions
- if self.use_cache:
- scan_regions = []
- for region in regions:
- key = str((region, genome, motif_digest, nreport, scan_rc))
- ret = self.cache.get(key)
- if ret == NO_VALUE:
- scan_regions.append(region)
-
- # scan the regions that are not in the cache
- if len(scan_regions) > 0:
-
- g = Genome(genome)
-
- motifs = [(m, self.threshold[m.id]) for m in read_motifs(self.motifs)]
- scan_func = partial(
- scan_region_mult,
- genome=g,
- motifs=motifs,
- nreport=nreport,
- scan_rc=scan_rc,
- )
+ if motifs is None:
+ motifs = read_motifs(self.motifs)
+ seq_gc_bins = [self.get_seq_bin(seq) for seq in seqs]
- for region, ret in self._scan_jobs(scan_func, scan_regions):
- # return values or store values in cache
- if self.use_cache:
- # store values in cache
- key = str(
- (
- region,
- genome,
- motif_digest,
- nreport,
- scan_rc,
- self.threshold_str,
- )
- )
- self.cache.set(key, ret)
- else:
- # return values
- yield ret
+ gc_bin_count = Counter(seq_gc_bins)
- if self.use_cache:
- # return results from cache
- for region in regions:
- key = str(
- (region, genome, motif_digest, nreport, scan_rc, self.threshold_str)
- )
- ret = self.cache.get(key)
- if ret == NO_VALUE or ret is None:
- raise Exception(
- "cache is not big enough to hold all "
- "results, try increasing the cache size "
- "or disable cache"
- )
- yield ret
+ _threshold = self.threshold
+ if zscore:
+ grouped = _threshold.groupby(_threshold.index).apply(scale, axis=0)
+ _threshold = pd.DataFrame(
+ np.vstack(grouped.values),
+ index=_threshold.index,
+ columns=_threshold.columns,
+ )
+
+ nseqs = int(20000 / np.sum(list(gc_bin_count.values())))
+ t = {}
+ maxt = pd.Series([m.pwm_max_score() for m in motifs], index=_threshold.columns)
+ # We do this in a loop as the DataFrame will get too big to fit in memory
+ # when the difference between the number of sequences per gc_bin is very
+ # high.
+ _threshold = _threshold.reset_index()
+ idx = np.hstack(
+ [
+ _threshold[_threshold[_threshold.columns[0]] == gc_bin]
+ .sample(nseqs * count, replace=True, random_state=42)
+ .index.values
+ for gc_bin, count in gc_bin_count.items()
+ ]
+ )
+ for motif in _threshold.columns[1:]:
+ val = _threshold.loc[idx, motif].quantile(0.99, interpolation="higher")
+ if val < maxt.loc[motif]:
+ t[motif] = val
+ else:
+ t[motif] = None
+ return t
def _scan_sequences_with_motif(self, motifs, seqs, nreport, scan_rc):
scan_func = partial(
scan_seq_mult, motifs=motifs, nreport=nreport, scan_rc=scan_rc
)
-
for ret in self._scan_jobs(scan_func, seqs):
yield ret[1]
- def _scan_sequences(self, seqs, nreport, scan_rc):
-
- motif_file = self.motifs
- motif_digest = self.checksum.get(motif_file, None)
-
- scan_seqs = seqs
- if self.use_cache:
- # determine which sequences are not in the cache
- hashes = dict([(s.upper(), xxhash.xxh64(s.upper()).digest()) for s in seqs])
- scan_seqs = []
-
- for seq, seq_hash in hashes.items():
- key = str(
- (seq_hash, motif_digest, nreport, scan_rc, self.threshold_str)
- )
- ret = self.cache.get(key)
- if ret == NO_VALUE or ret is None:
- scan_seqs.append(seq.upper())
-
- # scan the sequences that are not in the cache
- if len(scan_seqs) > 0:
- motifs = [(m, self.threshold[m.id]) for m in read_motifs(self.motifs)]
- scan_func = partial(
- scan_seq_mult, motifs=motifs, nreport=nreport, scan_rc=scan_rc
- )
-
- for seq, ret in self._scan_jobs(scan_func, scan_seqs):
- if self.use_cache:
- h = hashes[seq]
- key = str((h, motif_digest, nreport, scan_rc, self.threshold_str))
- self.cache.set(key, ret)
- else:
- yield ret
-
- if self.use_cache:
- # return results from cache
- for seq in seqs:
- key = str(
- (
- hashes[seq.upper()],
- motif_digest,
- nreport,
- scan_rc,
- self.threshold_str,
- )
- )
- ret = self.cache.get(key)
- if ret == NO_VALUE or ret is None:
- raise Exception(
- "cache is not big enough to hold all "
- "results, try increasing the cache size "
- "or disable cache"
- )
+ def _scan_sequences(self, seqs, nreport, scan_rc, zscore=False):
+ thresholds = self.get_gc_thresholds(seqs, zscore=zscore)
+ motifs = [(m, thresholds[m.id]) for m in read_motifs(self.motifs)]
+ motifs_meanstd = None
+ if zscore:
+ motifs_meanstd = self.meanstd
- yield ret
+ scan_func = partial(
+ scan_seq_mult,
+ motifs=motifs,
+ nreport=nreport,
+ scan_rc=scan_rc,
+ motifs_meanstd=motifs_meanstd,
+ zscore=zscore,
+ )
+ for _, ret in self._scan_jobs(scan_func, seqs):
+ yield ret
def _scan_jobs(self, scan_func, scan_seqs):
batchsize = 1000
+
if self.ncpus > 1:
for i in range((len(scan_seqs) - 1) // batchsize + 1):
batch = scan_seqs[i * batchsize : (i + 1) * batchsize]
chunksize = len(batch) // self.ncpus + 1
jobs = []
for j in range((len(batch) - 1) // chunksize + 1):
- job = self.pool.apply_async(
- scan_func, (batch[j * chunksize : (j + 1) * chunksize],)
- )
+ batch_seqs = batch[j * chunksize : (j + 1) * chunksize]
+ seq_gc_bins = [self.get_seq_bin(seq) for seq in batch_seqs]
+ job = self.pool.apply_async(scan_func, (batch_seqs, seq_gc_bins))
jobs.append(job)
for k, job in enumerate(jobs):
@@ -1117,7 +1081,8 @@ def _scan_jobs(self, scan_func, scan_seqs):
yield region, ret
else:
for i in range((len(scan_seqs) - 1) // batchsize + 1):
- for _j, ret in enumerate(
- scan_func(scan_seqs[i * batchsize : (i + 1) * batchsize])
- ):
+ batch_seqs = scan_seqs[i * batchsize : (i + 1) * batchsize]
+ seq_gc_bins = [self.get_seq_bin(seq) for seq in batch_seqs]
+
+ for _j, ret in enumerate(scan_func(batch_seqs, seq_gc_bins)):
yield scan_seqs[i], ret
diff --git a/gimmemotifs/tools/meme.py b/gimmemotifs/tools/meme.py
index 0319e8ad..3bb69d60 100644
--- a/gimmemotifs/tools/meme.py
+++ b/gimmemotifs/tools/meme.py
@@ -1,5 +1,6 @@
from .motifprogram import MotifProgram
import io
+import os
import re
from subprocess import Popen, PIPE
from tempfile import NamedTemporaryFile
@@ -75,8 +76,11 @@ def _run_program(self, bin, fastafile, params=None):
if not default_params["single"]:
cmd.append(strand)
- # sys.stderr.write(" ".join(cmd) + "\n")
- p = Popen(cmd, bufsize=1, stderr=PIPE, stdout=PIPE)
+ # Fix to run in Docker
+ env = os.environ.copy()
+ env["OMPI_MCA_plm_rsh_agent"] = "sh"
+
+ p = Popen(cmd, bufsize=1, stderr=PIPE, stdout=PIPE, env=env)
stdout, stderr = p.communicate()
motifs = []
diff --git a/gimmemotifs/tools/memew.py b/gimmemotifs/tools/memew.py
index c8d17e45..5f39ddef 100644
--- a/gimmemotifs/tools/memew.py
+++ b/gimmemotifs/tools/memew.py
@@ -1,5 +1,6 @@
from .motifprogram import MotifProgram
import io
+import os
import re
from subprocess import Popen, PIPE
from tempfile import NamedTemporaryFile
@@ -76,8 +77,12 @@ def _run_program(self, bin, fastafile, params=None):
if not default_params["single"]:
cmd.append(strand)
+ # Fix to run in Docker
+ env = os.environ.copy()
+ env["OMPI_MCA_plm_rsh_agent"] = "sh"
+
# sys.stderr.write(" ".join(cmd) + "\n")
- p = Popen(cmd, bufsize=1, stderr=PIPE, stdout=PIPE)
+ p = Popen(cmd, bufsize=1, stderr=PIPE, stdout=PIPE, env=env)
stdout, stderr = p.communicate()
motifs = []
diff --git a/gimmemotifs/tools/weeder.py b/gimmemotifs/tools/weeder.py
index 0361c2d2..0ad4334c 100644
--- a/gimmemotifs/tools/weeder.py
+++ b/gimmemotifs/tools/weeder.py
@@ -76,7 +76,7 @@ def _run_program(self, bin, fastafile, params=None):
shutil.copy(fastafile, name)
fastafile = name
- cmd = "{} -f {} -O".format(self.cmd, fastafile, weeder_organism)
+ cmd = "{} -f {} -O {}".format(self.cmd, fastafile, weeder_organism)
if params["single"]:
cmd += " -ss"
diff --git a/gimmemotifs/utils.py b/gimmemotifs/utils.py
index 17bd6fb9..a1b3c71b 100644
--- a/gimmemotifs/utils.py
+++ b/gimmemotifs/utils.py
@@ -5,8 +5,6 @@
# distribution.
""" Odds and ends that for which I didn't (yet) find another place """
-from __future__ import print_function
-
# Python imports
import os
import re
@@ -15,18 +13,21 @@
import logging
import mmap
import random
-import six
import tempfile
import requests
+from io import TextIOWrapper
+from functools import singledispatch
from subprocess import Popen
from tempfile import NamedTemporaryFile
from shutil import copyfile
# External imports
+import pyfaidx
from scipy import special
import numpy as np
import pybedtools
from genomepy import Genome
+from Bio.SeqIO.FastaIO import SimpleFastaParser
# gimme imports
@@ -49,8 +50,7 @@ def rc(seq):
def narrowpeak_to_bed(inputfile, bedfile, size=0):
- """Convert narrowPeak file to BED file.
- """
+ """Convert narrowPeak file to BED file."""
p = re.compile(r"^(#|track|browser)")
warn_no_summit = True
with open(bedfile, "w") as f_out:
@@ -88,7 +88,7 @@ def pfmfile_location(infile):
"database specified in the config file."
)
- if isinstance(infile, six.string_types):
+ if isinstance(infile, str):
if not os.path.exists(infile):
motif_dir = config.get_motif_dir()
checkfile = os.path.join(motif_dir, infile)
@@ -132,7 +132,7 @@ def phyper_single(k, good, bad, N):
def phyper(k, good, bad, N):
- """ Current hypergeometric implementation in scipy is broken,
+ """Current hypergeometric implementation in scipy is broken,
so here's the correct version.
"""
pvalues = [phyper_single(x, good, bad, N) for x in range(k + 1, N + 1)]
@@ -293,8 +293,8 @@ def motif_localization(fastafile, motif, size, outfile, cutoff=0.9):
def parse_cutoff(motifs, cutoff, default=0.9):
- """ Provide either a file with one cutoff per motif or a single cutoff
- returns a hash with motif id as key and cutoff as value
+ """Provide either a file with one cutoff per motif or a single cutoff
+ returns a hash with motif id as key and cutoff as value
"""
cutoffs = {}
@@ -473,7 +473,7 @@ def get_seqs_type(seqs):
- region file
- BED file
"""
- region_p = re.compile(r"^(.+):(\d+)-(\d+)$")
+ region_p = re.compile(r"^([^\s:]+\@)?(.+):(\d+)-(\d+)$")
if isinstance(seqs, Fasta):
return "fasta"
elif isinstance(seqs, list) or isinstance(seqs, np.ndarray):
@@ -499,24 +499,197 @@ def get_seqs_type(seqs):
raise ValueError("unknown type {}".format(type(seqs).__name__))
-def as_fasta(seqs, genome=None):
- ftype = get_seqs_type(seqs)
- if ftype == "fasta":
- return seqs
- elif ftype == "fastafile":
- return Fasta(seqs)
+# Regular expression to check for region (chr:start-end or genome@chr:start-end)
+region_p = re.compile(r"^[^@]+@([^\s]+):(\d+)-(\d+)$")
+
+
+def _check_minsize(fa, minsize):
+ """
+ Raise ValueError if there is any sequence that is shorter than minsize.
+ If minsize is None the size will not be checked.
+ """
+ if minsize is None:
+ return fa
+
+ for name, seq in fa.items():
+ if len(seq) < minsize:
+ raise ValueError(f"sequence {name} is shorter than {minsize}")
+
+ return fa
+
+
+def _genomepy_convert(to_convert, genome, minsize=None):
+ """
+ Convert a variety of inputs using track2fasta().
+ """
+ if genome is None:
+ raise ValueError("input file is not a FASTA file, need a genome!")
+
+ if isinstance(genome, Genome):
+ g = genome
else:
- if genome is None:
- raise ValueError("need genome to convert to FASTA")
+ g = Genome(genome)
+
+ tmpfile = NamedTemporaryFile()
+ g.track2fasta(to_convert, tmpfile.name)
+
+ fa = as_seqdict(tmpfile.name)
+ return _check_minsize(fa, minsize)
+
+
+def _as_seqdict_genome_regions(regions, minsize=None):
+ """
+ Accepts list of regions where the genome is encoded in the region,
+ using the genome@chrom:start-end format.
+ """
+ genomic_regions = {}
+ for region in regions:
+ genome, region = region.split("@")
+ if genome not in genomic_regions:
+ Genome(genome)
+ genomic_regions[genome] = []
+ genomic_regions[genome].append(region)
+
+ tmpfa = NamedTemporaryFile(mode="w", delete=False)
+ for genome, g_regions in genomic_regions.items():
+ g = Genome(genome)
+
+ fa = g.track2fasta(g_regions)
+
+ for seq in fa:
+ seq.name = f"{genome}@{seq.name}"
+ print(seq.__repr__(), file=tmpfa)
+
+ tmpfa.flush()
+
+ # Open tempfile and restore original sequence order
+ fa = as_seqdict(tmpfa.name)
+ fa = {region: fa[region] for region in regions}
+ return _check_minsize(fa, minsize)
+
+
+@singledispatch
+def as_seqdict(to_convert, genome=None, minsize=None):
+ """
+ Convert input to a dictionary with name as key and sequence as value.
+
+ If the input contains genomic coordinates, the genome needs to be
+ specified. If minsize is specified all sequences will be checked if they
+ are not shorter than minsize. If regions (or a region file) are used as
+ the input, the genome can optionally be specified in the region using the
+ following format: genome@chrom:start-end.
+
+ Current supported input types include:
+ * FASTA, BED and region files.
+ * List or numpy.ndarray of regions.
+ * pyfaidx.Fasta object.
+ * pybedtools.BedTool object.
+
+ Parameters
+ ----------
+ to_convert : list, str, pyfaidx.Fasta or pybedtools.BedTool
+ Input to convert to FASTA-like dictionary
+
+ genome : str, optional
+ Genomepy genome name.
+
+ minsize : int or None, optional
+ If specified, check if all sequences have at least size minsize.
+
+ Returns
+ -------
+ dict with sequence names as key and sequences as value.
+ """
+ raise NotImplementedError(f"Not implement for {type(to_convert)}")
+
+
+@as_seqdict.register(list)
+def _as_seqdict_list(to_convert, genome=None, minsize=None):
+ """
+ Accepts list of regions as input.
+ """
+ if region_p.match(to_convert[0]):
+ return _as_seqdict_genome_regions(to_convert, minsize)
+
+ return _genomepy_convert(to_convert, genome, minsize)
+
+
+@as_seqdict.register(TextIOWrapper)
+def _as_seqdict_file_object(to_convert, genome=None, minsize=None):
+ """
+ Accepts file object as input, should be a FASTA file.
+ """
+ fa = {x: y for x, y in SimpleFastaParser(to_convert)}
+ return _check_minsize(fa, minsize)
+
+
+@as_seqdict.register(str)
+def _as_seqdict_filename(to_convert, genome=None, minsize=None):
+ """
+ Accepts filename as input.
+ """
+ if not os.path.exists(to_convert):
+ raise ValueError("Assuming filename, but it does not exist")
+
+ f = open(to_convert)
+ fa = as_seqdict(f)
+
+ if any(fa):
+ return _check_minsize(fa, minsize)
+
+ with open(to_convert) as f:
+ line = ""
+ while True:
+ line = f.readline()
+ if line == "":
+ break
+ if not line.startswith("#"):
+ break
+
+ if line == "":
+ raise IOError(f"empty file {to_convert}")
+
+ if region_p.match(line.strip()):
+ regions = [myline.strip() for myline in [line] + f.readlines()]
+ return _as_seqdict_genome_regions(regions, minsize=None)
+
+ # Biopython parser resulted in empty dict
+ # Assuming it's a BED or region file
+ return _genomepy_convert(to_convert, genome, minsize)
+
+
+@as_seqdict.register(pyfaidx.Fasta)
+def _as_seqdict_pyfaidx(to_convert, genome=None, minsize=None):
+ """
+ Accepts pyfaidx.Fasta object as input.
+ """
+ fa = {k: str(v) for k, v in to_convert.items()}
+ return _check_minsize(fa, minsize)
+
+
+@as_seqdict.register(pybedtools.BedTool)
+def _as_seqdict_bedtool(to_convert, genome=None, minsize=None):
+ """
+ Accepts pybedtools.BedTool as input.
+ """
+ return _genomepy_convert(
+ ["{}:{}-{}".format(*f[:3]) for f in to_convert], genome, minsize
+ )
+
+
+@as_seqdict.register(np.ndarray)
+def _as_seqdict_array(to_convert, genome=None, minsize=None):
+ """
+ Accepts numpy.ndarray with regions as input.
+ """
+ return as_seqdict(list(to_convert), genome, minsize)
+
- tmpfa = NamedTemporaryFile()
- if isinstance(genome, str):
- genome = Genome(genome)
+def as_fasta(to_convert, genome=None, minsize=None):
+ if isinstance(to_convert, Fasta):
+ return to_convert
- if isinstance(seqs, np.ndarray):
- seqs = list(seqs)
- genome.track2fasta(seqs, tmpfa.name)
- return Fasta(tmpfa.name)
+ return Fasta(fdict=as_seqdict(to_convert, genome, minsize))
def file_checksum(fname):
@@ -539,11 +712,11 @@ def file_checksum(fname):
return checksum
-def join_max(a, l, sep="", suffix=""):
+def join_max(a, length, sep="", suffix=""):
lengths = [len(x) for x in a]
total = 0
for i, size in enumerate(lengths + [0]):
- if total > (l - len(suffix)):
+ if total > (length - len(suffix)):
return sep.join(a[: i - 1]) + suffix
if i > 0:
total += 1
diff --git a/readthedocs.yml b/readthedocs.yml
index 9a6a856f..e458c84e 100644
--- a/readthedocs.yml
+++ b/readthedocs.yml
@@ -1,18 +1,18 @@
# Required
version: 2
-# Build documentation in the docs/ directory with Sphinx
+conda:
+ environment: .rtd-environment.yml
+
+ # Build documentation in the docs/ directory with Sphinx
sphinx:
configuration: docs/conf.py
+build:
+ image: latest
+
python:
- version: 3.7
- install:
- - requirements: docs/requirements.txt
- - method: pip
- path: .
- extra_requirements:
- - docs
- - method: setuptools
- path: another/package
- system_packages: true
+ version: 3.7
+ install:
+ - method: setuptools
+ path: .
diff --git a/requirements.txt b/requirements.txt
index 604878ff..ec991b68 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -10,12 +10,9 @@ pyyaml >= 3.10
pybedtools
statsmodels
scikit-learn
-sklearn-contrib-lightning
seaborn
pysam
xgboost
diskcache
xxhash
-six
-future
pillow
diff --git a/scripts/combine_peaks b/scripts/combine_peaks
index 155d6e83..da2870c0 100644
--- a/scripts/combine_peaks
+++ b/scripts/combine_peaks
@@ -46,7 +46,7 @@ def read_peak_file_to_df(fname):
"qval",
"peak",
]
- df = pd.read_table(fname, names=header)
+ df = pd.read_table(fname, names=header, dtype={"chrom": "str"})
df["chrom"] = df["chrom"].astype(str)
# get the summit
@@ -57,7 +57,7 @@ def read_peak_file_to_df(fname):
df["value"] = df["qval"]
df = df[summit_header]
elif ftype == "bed":
- df = pd.read_table(fname, names=summit_header)
+ df = pd.read_table(fname, names=summit_header, dtype={"chrom": "str"})
if ((df["end"] - df["start"]) != 1).sum() != 0:
raise ValueError(f"{fname} does not contain summits.")
else:
@@ -116,7 +116,7 @@ def combine_peaks(peaks, genome, window, scale_value):
# store summit location + associated value in col4
df_all["col4"] = (
- df_all["chrom"]
+ df_all["chrom"].astype(str)
+ ";"
+ df_all["start"].astype(str)
+ ";"
diff --git a/scripts/coverage_table b/scripts/coverage_table
index 64385b65..14f3128f 100644
--- a/scripts/coverage_table
+++ b/scripts/coverage_table
@@ -9,6 +9,7 @@ import pysam
import numpy as np
import pandas as pd
from sklearn.preprocessing import scale
+import qnorm
from gimmemotifs import __version__
@@ -27,11 +28,12 @@ def make_table(
datafiles,
window,
log_transform=True,
- scale_table=True,
+ normalization="none",
top=0,
topmethod="var",
rmdup=True,
rmrepeats=True,
+ ncpus=12
):
for x in datafiles:
if not os.path.isfile(x):
@@ -49,7 +51,7 @@ def make_table(
data = {}
try:
# Load data in parallel
- pool = multiprocessing.Pool(processes=12)
+ pool = multiprocessing.Pool(processes=ncpus)
jobs = []
for datafile in datafiles:
jobs.append(
@@ -99,9 +101,14 @@ def make_table(
if log_transform:
print("Log transform", file=sys.stderr)
df = np.log1p(df)
- if scale_table:
- print("Scale", file=sys.stderr)
+ if normalization == "scale":
+ print("Normalization by scaling", file=sys.stderr)
df[:] = scale(df, axis=0)
+ if normalization == "quantile":
+ print("Normalization by quantile normalization", file=sys.stderr)
+ df = qnorm.quantile_normalize(df)
+ else:
+ print("No normalization", file=sys.stderr)
if top > 0:
if topmethod == "var":
@@ -165,12 +172,12 @@ if __name__ == "__main__":
action="store_true",
)
parser.add_argument(
- "-s",
- "--scale",
- dest="scale_table",
- help="Scale per datafile",
- default=False,
- action="store_true",
+ "-n",
+ "--normalization",
+ dest="normalization",
+ help="Normalization: none, quantile or scale",
+ default="none",
+ metavar="METHOD",
)
parser.add_argument(
"-t", "--top", dest="top", help="Select regions.", default=0, type=int
@@ -196,6 +203,14 @@ if __name__ == "__main__":
action="store_false",
default=True,
)
+ parser.add_argument(
+ "--nthreads",
+ dest="ncpus",
+ help="Number of threads",
+ metavar="INT",
+ type=int,
+ default=12,
+ )
args = parser.parse_args()
peakfile = args.peakfile
@@ -205,11 +220,12 @@ if __name__ == "__main__":
datafiles,
args.window,
log_transform=args.log_transform,
- scale_table=args.scale_table,
+ normalization=args.normalization,
top=args.top,
topmethod=args.topmethod,
rmdup=args.rmdup,
rmrepeats=args.rmrepeats,
+ ncpus=args.ncpus
)
yesno = {True: "yes", False: "no"}
@@ -222,7 +238,7 @@ output.write("# Window: {}\n".format(args.window))
output.write("# Duplicates removed: {}\n".format(yesno[args.rmdup]))
output.write("# MAPQ 0 removed: {}\n".format(yesno[args.rmrepeats]))
output.write("# Log transformed: {}\n".format(yesno[args.log_transform]))
-output.write("# Scaled: {}\n".format(yesno[args.scale_table]))
+output.write("# Normalization: {}\n".format(args.normalization))
if args.top > 0:
output.write("# Top {} regions selected by {}\n".format(args.top, args.topmethod))
df.to_csv(output, sep="\t", float_format="%0.5f")
diff --git a/setup.py b/setup.py
index 5f9fd60f..b5dc2217 100644
--- a/setup.py
+++ b/setup.py
@@ -124,7 +124,6 @@ def run(self):
"pybedtools",
"statsmodels",
"scikit-learn",
- "sklearn-contrib-lightning",
"seaborn",
"pysam",
"xgboost >= 0.71",
@@ -132,11 +131,10 @@ def run(self):
"diskcache",
"xxhash",
"configparser",
- "six",
- "future",
- "genomepy >= 0.7.2",
+ "genomepy >= 0.8.3",
"tqdm",
"pillow",
"logomaker",
+ "qnorm",
],
)
diff --git a/test/data/rank/ranked.txt b/test/data/rank/ranked.txt
index f5f36f40..bb52cd79 100644
--- a/test/data/rank/ranked.txt
+++ b/test/data/rank/ranked.txt
@@ -1,6 +1,6 @@
a b c d e
-1 bZIP bZIP SOX T-box POU
-2 AP2 AP2 AP2 AP2 AP2
-3 T-box SOX T-box SOX SOX
-4 SOX T-box POU POU T-box
-5 POU POU bZIP bZIP bZIP
+bZIP 5 5 1 1 1
+AP2 4 4 4 4 4
+T-box 3 2 3 5 2
+SOX 2 3 5 3 3
+POU 1 1 2 2 4
diff --git a/test/test_maelstrom.py b/test/test_maelstrom.py
index 269e5b92..53bffc77 100644
--- a/test/test_maelstrom.py
+++ b/test/test_maelstrom.py
@@ -26,13 +26,30 @@ def test1_maelstrom(self):
self.clusters,
"mm10",
self.outdir,
+ filter_redundant=False,
score_table=self.score_table,
count_table=self.count_table,
plot=False,
)
df = pd.read_table(self.outfile, index_col=0, comment="#")
print(df.shape)
- self.assertEquals((623, 4), df.shape)
+
+ self.assertEquals((623, 8), df.shape)
+
+ # Filter redundant motifs
+ run_maelstrom(
+ self.clusters,
+ "mm10",
+ self.outdir,
+ filter_redundant=True,
+ score_table=self.score_table,
+ count_table=self.count_table,
+ plot=False,
+ )
+ df = pd.read_table(self.outfile, index_col=0, comment="#")
+ print(df.shape)
+ self.assertEquals((156, 8), df.shape)
+
for fname in glob(os.path.join(self.outdir, "activity*")):
os.unlink(fname)
diff --git a/test/test_moap.py b/test/test_moap.py
index 775b3534..c6c13be3 100644
--- a/test/test_moap.py
+++ b/test/test_moap.py
@@ -20,7 +20,7 @@ def setUp(self):
def test1_moap(self):
""" Test motif activity prediction """
- for method in ["mwu", "rf", "lightningclassification"]:
+ for method in ["mwu", "rf"]:
df = moap(
self.clusters,
method=method,
@@ -41,7 +41,7 @@ def test1_moap(self):
def test2_moap(self):
""" Test motif activity prediction for two clusters """
- for method in ["mwu", "rf", "lightningclassification"]:
+ for method in ["mwu", "rf"]:
df = moap(
self.clusters2,
method=method,
diff --git a/test/test_rank.py b/test/test_rank.py
index 5cb533e9..cca52132 100644
--- a/test/test_rank.py
+++ b/test/test_rank.py
@@ -2,7 +2,7 @@
import tempfile
import os
import pandas as pd
-from gimmemotifs.rank import rankagg
+from gimmemotifs.rank import rankagg, _rankagg_stuart
class TestRank(unittest.TestCase):
@@ -17,13 +17,15 @@ def setUp(self):
def test1_rankagg(self):
""" Test rank aggregation """
df = pd.read_csv(self.fname, index_col=0, sep="\t")
- result = rankagg(df)
- self.assertEqual("AP2", result.sort_values("score").index[0])
+ result = rankagg(df, method="stuart")
+ self.assertEqual("AP2", result.sort_values("score").index[-1])
+ result = rankagg(df, method="int_stouffer")
+ self.assertEqual("AP2", result.sort_values("z-score").index[-1])
def test2_rankagg(self):
""" Test Python implementation of rank aggregation """
df = pd.read_csv(self.rank_in, index_col=0, sep="\t")
- result = rankagg(df)["score"].values
+ result = _rankagg_stuart(df)["score"].values
ref = pd.read_csv(self.rank_out, index_col=0, sep="\t")["score"].values
for v1, v2 in zip(ref, result):
self.assertAlmostEqual(v1, v2)
diff --git a/test/test_tools.py b/test/test_tools.py
index d6dba2db..9903ce54 100644
--- a/test/test_tools.py
+++ b/test/test_tools.py
@@ -35,6 +35,7 @@ def test_tool(tool_name):
"trawler", # unpredictable, sometimes doesn't find the motif
"weeder", # doesn't work at the moment
"posmo", # motif doesn't predictably look like AP1
+ "dreme", # current dreme in bioconda is broken
]:
return
diff --git a/test/test_utils.py b/test/test_utils.py
index 87f7b9da..c4f72c2f 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -13,7 +13,7 @@ class TestUtils(unittest.TestCase):
""" A test class to test utils functions """
def setUp(self):
- self.genome_dir = "test/data/genome_index"
+ self.genomes_dir = "test/data/genome_index"
self.datadir = "test/data/utils"
def test1_phyper(self):
@@ -31,7 +31,7 @@ def test2_as_fasta(self):
""" convert bed, regions, etc to Fasta """
tmpdir = mkdtemp()
- g = Genome("genome", genome_dir=self.genome_dir)
+ g = Genome("genome", genomes_dir=self.genomes_dir)
fafile = os.path.join(self.datadir, "test.fa")
fa = Fasta(fafile)
diff --git a/versioneer.py b/versioneer.py
index 64fea1c8..cce201c7 100644
--- a/versioneer.py
+++ b/versioneer.py
@@ -1,4 +1,3 @@
-
# Version: 0.18
"""The Versioneer - like a rocketeer, but for versions.
@@ -276,7 +275,6 @@
"""
-from __future__ import print_function
try:
import configparser
except ImportError:
@@ -308,11 +306,13 @@ def get_root():
setup_py = os.path.join(root, "setup.py")
versioneer_py = os.path.join(root, "versioneer.py")
if not (os.path.exists(setup_py) or os.path.exists(versioneer_py)):
- err = ("Versioneer was unable to run the project root directory. "
- "Versioneer requires setup.py to be executed from "
- "its immediate directory (like 'python setup.py COMMAND'), "
- "or in a way that lets it use sys.argv[0] to find the root "
- "(like 'python path/to/setup.py COMMAND').")
+ err = (
+ "Versioneer was unable to run the project root directory. "
+ "Versioneer requires setup.py to be executed from "
+ "its immediate directory (like 'python setup.py COMMAND'), "
+ "or in a way that lets it use sys.argv[0] to find the root "
+ "(like 'python path/to/setup.py COMMAND')."
+ )
raise VersioneerBadRootError(err)
try:
# Certain runtime workflows (setup.py install/develop in a setuptools
@@ -325,8 +325,10 @@ def get_root():
me_dir = os.path.normcase(os.path.splitext(me)[0])
vsr_dir = os.path.normcase(os.path.splitext(versioneer_py)[0])
if me_dir != vsr_dir:
- print("Warning: build in %s is using versioneer.py from %s"
- % (os.path.dirname(me), versioneer_py))
+ print(
+ "Warning: build in %s is using versioneer.py from %s"
+ % (os.path.dirname(me), versioneer_py)
+ )
except NameError:
pass
return root
@@ -348,6 +350,7 @@ def get(parser, name):
if parser.has_option("versioneer", name):
return parser.get("versioneer", name)
return None
+
cfg = VersioneerConfig()
cfg.VCS = VCS
cfg.style = get(parser, "style") or ""
@@ -372,17 +375,18 @@ class NotThisMethod(Exception):
def register_vcs_handler(vcs, method): # decorator
"""Decorator to mark a method as the handler for a particular VCS."""
+
def decorate(f):
"""Store f in HANDLERS[vcs][method]."""
if vcs not in HANDLERS:
HANDLERS[vcs] = {}
HANDLERS[vcs][method] = f
return f
+
return decorate
-def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
- env=None):
+def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False, env=None):
"""Call the given command(s)."""
assert isinstance(commands, list)
p = None
@@ -390,10 +394,13 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
try:
dispcmd = str([c] + args)
# remember shell=False, so use git.cmd on windows, not just git
- p = subprocess.Popen([c] + args, cwd=cwd, env=env,
- stdout=subprocess.PIPE,
- stderr=(subprocess.PIPE if hide_stderr
- else None))
+ p = subprocess.Popen(
+ [c] + args,
+ cwd=cwd,
+ env=env,
+ stdout=subprocess.PIPE,
+ stderr=(subprocess.PIPE if hide_stderr else None),
+ )
break
except EnvironmentError:
e = sys.exc_info()[1]
@@ -418,7 +425,9 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
return stdout, p.returncode
-LONG_VERSION_PY['git'] = '''
+LONG_VERSION_PY[
+ "git"
+] = '''
# This file helps to compute a version number in source trees obtained from
# git-archive tarball (such as those provided by githubs download-from-tag
# feature). Distribution tarballs (built by setup.py sdist) and build
@@ -993,7 +1002,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
# starting in git-1.8.3, tags are listed as "tag: foo-1.0" instead of
# just "foo-1.0". If we see a "tag: " prefix, prefer those.
TAG = "tag: "
- tags = set([r[len(TAG):] for r in refs if r.startswith(TAG)])
+ tags = set([r[len(TAG) :] for r in refs if r.startswith(TAG)])
if not tags:
# Either we're using git < 1.8.3, or there really are no tags. We use
# a heuristic: assume all version tags have a digit. The old git %d
@@ -1002,7 +1011,7 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
# between branches and tags. By ignoring refnames without digits, we
# filter out many common branch names like "release" and
# "stabilization", as well as "HEAD" and "master".
- tags = set([r for r in refs if re.search(r'\d', r)])
+ tags = set([r for r in refs if re.search(r"\d", r)])
if verbose:
print("discarding '%s', no digits" % ",".join(refs - tags))
if verbose:
@@ -1010,19 +1019,26 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
for ref in sorted(tags):
# sorting will prefer e.g. "2.0" over "2.0rc1"
if ref.startswith(tag_prefix):
- r = ref[len(tag_prefix):]
+ r = ref[len(tag_prefix) :]
if verbose:
print("picking %s" % r)
- return {"version": r,
- "full-revisionid": keywords["full"].strip(),
- "dirty": False, "error": None,
- "date": date}
+ return {
+ "version": r,
+ "full-revisionid": keywords["full"].strip(),
+ "dirty": False,
+ "error": None,
+ "date": date,
+ }
# no suitable tags, so version is "0+unknown", but full hex is still there
if verbose:
print("no suitable tags, using unknown + full revision id")
- return {"version": "0+unknown",
- "full-revisionid": keywords["full"].strip(),
- "dirty": False, "error": "no suitable tags", "date": None}
+ return {
+ "version": "0+unknown",
+ "full-revisionid": keywords["full"].strip(),
+ "dirty": False,
+ "error": "no suitable tags",
+ "date": None,
+ }
@register_vcs_handler("git", "pieces_from_vcs")
@@ -1037,8 +1053,7 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
if sys.platform == "win32":
GITS = ["git.cmd", "git.exe"]
- out, rc = run_command(GITS, ["rev-parse", "--git-dir"], cwd=root,
- hide_stderr=True)
+ out, rc = run_command(GITS, ["rev-parse", "--git-dir"], cwd=root, hide_stderr=True)
if rc != 0:
if verbose:
print("Directory %s not under git control" % root)
@@ -1046,10 +1061,19 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
# if there is a tag matching tag_prefix, this yields TAG-NUM-gHEX[-dirty]
# if there isn't one, this yields HEX[-dirty] (no NUM)
- describe_out, rc = run_command(GITS, ["describe", "--tags", "--dirty",
- "--always", "--long",
- "--match", "%s*" % tag_prefix],
- cwd=root)
+ describe_out, rc = run_command(
+ GITS,
+ [
+ "describe",
+ "--tags",
+ "--dirty",
+ "--always",
+ "--long",
+ "--match",
+ "%s*" % tag_prefix,
+ ],
+ cwd=root,
+ )
# --long was added in git-1.5.5
if describe_out is None:
raise NotThisMethod("'git describe' failed")
@@ -1072,17 +1096,16 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
dirty = git_describe.endswith("-dirty")
pieces["dirty"] = dirty
if dirty:
- git_describe = git_describe[:git_describe.rindex("-dirty")]
+ git_describe = git_describe[: git_describe.rindex("-dirty")]
# now we have TAG-NUM-gHEX or HEX
if "-" in git_describe:
# TAG-NUM-gHEX
- mo = re.search(r'^(.+)-(\d+)-g([0-9a-f]+)$', git_describe)
+ mo = re.search(r"^(.+)-(\d+)-g([0-9a-f]+)$", git_describe)
if not mo:
# unparseable. Maybe git-describe is misbehaving?
- pieces["error"] = ("unable to parse git-describe output: '%s'"
- % describe_out)
+ pieces["error"] = "unable to parse git-describe output: '%s'" % describe_out
return pieces
# tag
@@ -1091,10 +1114,12 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
if verbose:
fmt = "tag '%s' doesn't start with prefix '%s'"
print(fmt % (full_tag, tag_prefix))
- pieces["error"] = ("tag '%s' doesn't start with prefix '%s'"
- % (full_tag, tag_prefix))
+ pieces["error"] = "tag '%s' doesn't start with prefix '%s'" % (
+ full_tag,
+ tag_prefix,
+ )
return pieces
- pieces["closest-tag"] = full_tag[len(tag_prefix):]
+ pieces["closest-tag"] = full_tag[len(tag_prefix) :]
# distance: number of commits since tag
pieces["distance"] = int(mo.group(2))
@@ -1105,13 +1130,13 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
else:
# HEX: no tags
pieces["closest-tag"] = None
- count_out, rc = run_command(GITS, ["rev-list", "HEAD", "--count"],
- cwd=root)
+ count_out, rc = run_command(GITS, ["rev-list", "HEAD", "--count"], cwd=root)
pieces["distance"] = int(count_out) # total number of commits
# commit date: see ISO-8601 comment in git_versions_from_keywords()
- date = run_command(GITS, ["show", "-s", "--format=%ci", "HEAD"],
- cwd=root)[0].strip()
+ date = run_command(GITS, ["show", "-s", "--format=%ci", "HEAD"], cwd=root)[
+ 0
+ ].strip()
pieces["date"] = date.strip().replace(" ", "T", 1).replace(" ", "", 1)
return pieces
@@ -1167,16 +1192,22 @@ def versions_from_parentdir(parentdir_prefix, root, verbose):
for i in range(3):
dirname = os.path.basename(root)
if dirname.startswith(parentdir_prefix):
- return {"version": dirname[len(parentdir_prefix):],
- "full-revisionid": None,
- "dirty": False, "error": None, "date": None}
+ return {
+ "version": dirname[len(parentdir_prefix) :],
+ "full-revisionid": None,
+ "dirty": False,
+ "error": None,
+ "date": None,
+ }
else:
rootdirs.append(root)
root = os.path.dirname(root) # up a level
if verbose:
- print("Tried directories %s but none started with prefix %s" %
- (str(rootdirs), parentdir_prefix))
+ print(
+ "Tried directories %s but none started with prefix %s"
+ % (str(rootdirs), parentdir_prefix)
+ )
raise NotThisMethod("rootdir doesn't start with parentdir_prefix")
@@ -1205,11 +1236,13 @@ def versions_from_file(filename):
contents = f.read()
except EnvironmentError:
raise NotThisMethod("unable to read _version.py")
- mo = re.search(r"version_json = '''\n(.*)''' # END VERSION_JSON",
- contents, re.M | re.S)
+ mo = re.search(
+ r"version_json = '''\n(.*)''' # END VERSION_JSON", contents, re.M | re.S
+ )
if not mo:
- mo = re.search(r"version_json = '''\r\n(.*)''' # END VERSION_JSON",
- contents, re.M | re.S)
+ mo = re.search(
+ r"version_json = '''\r\n(.*)''' # END VERSION_JSON", contents, re.M | re.S
+ )
if not mo:
raise NotThisMethod("no version_json in _version.py")
return json.loads(mo.group(1))
@@ -1218,8 +1251,7 @@ def versions_from_file(filename):
def write_to_version_file(filename, versions):
"""Write the given version number to the given _version.py file."""
os.unlink(filename)
- contents = json.dumps(versions, sort_keys=True,
- indent=1, separators=(",", ": "))
+ contents = json.dumps(versions, sort_keys=True, indent=1, separators=(",", ": "))
with open(filename, "w") as f:
f.write(SHORT_VERSION_PY % contents)
@@ -1251,8 +1283,7 @@ def render_pep440(pieces):
rendered += ".dirty"
else:
# exception #1
- rendered = "0+untagged.%d.g%s" % (pieces["distance"],
- pieces["short"])
+ rendered = "0+untagged.%d.g%s" % (pieces["distance"], pieces["short"])
if pieces["dirty"]:
rendered += ".dirty"
return rendered
@@ -1366,11 +1397,13 @@ def render_git_describe_long(pieces):
def render(pieces, style):
"""Render the given version pieces into the requested style."""
if pieces["error"]:
- return {"version": "unknown",
- "full-revisionid": pieces.get("long"),
- "dirty": None,
- "error": pieces["error"],
- "date": None}
+ return {
+ "version": "unknown",
+ "full-revisionid": pieces.get("long"),
+ "dirty": None,
+ "error": pieces["error"],
+ "date": None,
+ }
if not style or style == "default":
style = "pep440" # the default
@@ -1390,9 +1423,13 @@ def render(pieces, style):
else:
raise ValueError("unknown style '%s'" % style)
- return {"version": rendered, "full-revisionid": pieces["long"],
- "dirty": pieces["dirty"], "error": None,
- "date": pieces.get("date")}
+ return {
+ "version": rendered,
+ "full-revisionid": pieces["long"],
+ "dirty": pieces["dirty"],
+ "error": None,
+ "date": pieces.get("date"),
+ }
class VersioneerBadRootError(Exception):
@@ -1415,8 +1452,9 @@ def get_versions(verbose=False):
handlers = HANDLERS.get(cfg.VCS)
assert handlers, "unrecognized VCS '%s'" % cfg.VCS
verbose = verbose or cfg.verbose
- assert cfg.versionfile_source is not None, \
- "please set versioneer.versionfile_source"
+ assert (
+ cfg.versionfile_source is not None
+ ), "please set versioneer.versionfile_source"
assert cfg.tag_prefix is not None, "please set versioneer.tag_prefix"
versionfile_abs = os.path.join(root, cfg.versionfile_source)
@@ -1470,9 +1508,13 @@ def get_versions(verbose=False):
if verbose:
print("unable to compute version")
- return {"version": "0+unknown", "full-revisionid": None,
- "dirty": None, "error": "unable to compute version",
- "date": None}
+ return {
+ "version": "0+unknown",
+ "full-revisionid": None,
+ "dirty": None,
+ "error": "unable to compute version",
+ "date": None,
+ }
def get_version():
@@ -1521,6 +1563,7 @@ def run(self):
print(" date: %s" % vers.get("date"))
if vers["error"]:
print(" error: %s" % vers["error"])
+
cmds["version"] = cmd_version
# we override "build_py" in both distutils and setuptools
@@ -1553,14 +1596,15 @@ def run(self):
# now locate _version.py in the new build/ directory and replace
# it with an updated value
if cfg.versionfile_build:
- target_versionfile = os.path.join(self.build_lib,
- cfg.versionfile_build)
+ target_versionfile = os.path.join(self.build_lib, cfg.versionfile_build)
print("UPDATING %s" % target_versionfile)
write_to_version_file(target_versionfile, versions)
+
cmds["build_py"] = cmd_build_py
if "cx_Freeze" in sys.modules: # cx_freeze enabled?
from cx_Freeze.dist import build_exe as _build_exe
+
# nczeczulin reports that py2exe won't like the pep440-style string
# as FILEVERSION, but it can be used for PRODUCTVERSION, e.g.
# setup(console=[{
@@ -1581,17 +1625,21 @@ def run(self):
os.unlink(target_versionfile)
with open(cfg.versionfile_source, "w") as f:
LONG = LONG_VERSION_PY[cfg.VCS]
- f.write(LONG %
- {"DOLLAR": "$",
- "STYLE": cfg.style,
- "TAG_PREFIX": cfg.tag_prefix,
- "PARENTDIR_PREFIX": cfg.parentdir_prefix,
- "VERSIONFILE_SOURCE": cfg.versionfile_source,
- })
+ f.write(
+ LONG
+ % {
+ "DOLLAR": "$",
+ "STYLE": cfg.style,
+ "TAG_PREFIX": cfg.tag_prefix,
+ "PARENTDIR_PREFIX": cfg.parentdir_prefix,
+ "VERSIONFILE_SOURCE": cfg.versionfile_source,
+ }
+ )
+
cmds["build_exe"] = cmd_build_exe
del cmds["build_py"]
- if 'py2exe' in sys.modules: # py2exe enabled?
+ if "py2exe" in sys.modules: # py2exe enabled?
try:
from py2exe.distutils_buildexe import py2exe as _py2exe # py3
except ImportError:
@@ -1610,13 +1658,17 @@ def run(self):
os.unlink(target_versionfile)
with open(cfg.versionfile_source, "w") as f:
LONG = LONG_VERSION_PY[cfg.VCS]
- f.write(LONG %
- {"DOLLAR": "$",
- "STYLE": cfg.style,
- "TAG_PREFIX": cfg.tag_prefix,
- "PARENTDIR_PREFIX": cfg.parentdir_prefix,
- "VERSIONFILE_SOURCE": cfg.versionfile_source,
- })
+ f.write(
+ LONG
+ % {
+ "DOLLAR": "$",
+ "STYLE": cfg.style,
+ "TAG_PREFIX": cfg.tag_prefix,
+ "PARENTDIR_PREFIX": cfg.parentdir_prefix,
+ "VERSIONFILE_SOURCE": cfg.versionfile_source,
+ }
+ )
+
cmds["py2exe"] = cmd_py2exe
# we override different "sdist" commands for both environments
@@ -1643,8 +1695,10 @@ def make_release_tree(self, base_dir, files):
# updated value
target_versionfile = os.path.join(base_dir, cfg.versionfile_source)
print("UPDATING %s" % target_versionfile)
- write_to_version_file(target_versionfile,
- self._versioneer_generated_versions)
+ write_to_version_file(
+ target_versionfile, self._versioneer_generated_versions
+ )
+
cmds["sdist"] = cmd_sdist
return cmds
@@ -1699,11 +1753,13 @@ def do_setup():
root = get_root()
try:
cfg = get_config_from_root(root)
- except (EnvironmentError, configparser.NoSectionError,
- configparser.NoOptionError) as e:
+ except (
+ EnvironmentError,
+ configparser.NoSectionError,
+ configparser.NoOptionError,
+ ) as e:
if isinstance(e, (EnvironmentError, configparser.NoSectionError)):
- print("Adding sample versioneer config to setup.cfg",
- file=sys.stderr)
+ print("Adding sample versioneer config to setup.cfg", file=sys.stderr)
with open(os.path.join(root, "setup.cfg"), "a") as f:
f.write(SAMPLE_CONFIG)
print(CONFIG_ERROR, file=sys.stderr)
@@ -1712,15 +1768,18 @@ def do_setup():
print(" creating %s" % cfg.versionfile_source)
with open(cfg.versionfile_source, "w") as f:
LONG = LONG_VERSION_PY[cfg.VCS]
- f.write(LONG % {"DOLLAR": "$",
- "STYLE": cfg.style,
- "TAG_PREFIX": cfg.tag_prefix,
- "PARENTDIR_PREFIX": cfg.parentdir_prefix,
- "VERSIONFILE_SOURCE": cfg.versionfile_source,
- })
-
- ipy = os.path.join(os.path.dirname(cfg.versionfile_source),
- "__init__.py")
+ f.write(
+ LONG
+ % {
+ "DOLLAR": "$",
+ "STYLE": cfg.style,
+ "TAG_PREFIX": cfg.tag_prefix,
+ "PARENTDIR_PREFIX": cfg.parentdir_prefix,
+ "VERSIONFILE_SOURCE": cfg.versionfile_source,
+ }
+ )
+
+ ipy = os.path.join(os.path.dirname(cfg.versionfile_source), "__init__.py")
if os.path.exists(ipy):
try:
with open(ipy, "r") as f:
@@ -1762,8 +1821,10 @@ def do_setup():
else:
print(" 'versioneer.py' already in MANIFEST.in")
if cfg.versionfile_source not in simple_includes:
- print(" appending versionfile_source ('%s') to MANIFEST.in" %
- cfg.versionfile_source)
+ print(
+ " appending versionfile_source ('%s') to MANIFEST.in"
+ % cfg.versionfile_source
+ )
with open(manifest_in, "a") as f:
f.write("include %s\n" % cfg.versionfile_source)
else: