Add requirements to manifest (#445)

This _should_ make it look more attractive on pypi.

MANIFEST.in

@@ -1,3 +1,4 @@
 recursive-include server/app/web/templates *
 recursive-include server/app/web/static *
 
+include server/requirements.txt

client/src/actions/index.js

@@ -90,12 +90,16 @@ async function _doRequestExpressionData(dispatch, getState, genes) {
   const state = getState();
   const { universe } = state.controls;
   /* preload data already in cache */
-  let expressionData = _.transform(genes, (expData, g) => {
+  let expressionData = _.transform(
-    const data = kvCache.get(universe.varDataCache, g);
+    genes,
-    if (data) {
+    (expData, g) => {
-      expData[g] = data;
+      const data = kvCache.get(universe.varDataCache, g);
-    }
+      if (data) {
-  }); // --> { gene: data }
+        expData[g] = data;
+      }
+    },
+    {}
+  ); // --> { gene: data }
   /* make a list of genes for which we do not have data */
   const genesToFetch = _.filter(genes, g => expressionData[g] === undefined);
 
@@ -119,7 +123,6 @@ async function _doRequestExpressionData(dispatch, getState, genes) {
           }),
           headers: new Headers({
             accept: "application/json",
-            "Accept-Encoding": "gzip, deflate, br",
             "Content-Type": "application/json"
           })
         }
@@ -239,7 +242,6 @@ const requestDifferentialExpression = (set1, set2, num_genes = 10) => async (
         method: "POST",
         headers: new Headers({
           Accept: "application/json",
-          "Accept-Encoding": "gzip, deflate, br",
           "Content-Type": "application/json"
         }),
         body: JSON.stringify({
@@ -297,6 +299,9 @@ const resetInterface = () => (dispatch, getState) => {
     type: "reset World to eq Universe",
     universe
   });
+  dispatch({
+    type: "increment graph render counter"
+  });
 };
 
 export default {

client/src/components/app.js

@@ -11,7 +11,8 @@ import actions from "../actions";
 
 @connect(state => ({
   loading: state.controls.loading,
-  error: state.controls.error
+  error: state.controls.error,
+  graphRenderCounter: state.controls.graphRenderCounter
 }))
 class App extends React.Component {
   constructor(props) {
@@ -54,7 +55,7 @@ class App extends React.Component {
   }
 
   render() {
-    const { loading } = this.props;
+    const { loading, error, graphRenderCounter } = this.props;
     return (
       <Container>
         <Helmet title="cellxgene" />
@@ -79,10 +80,8 @@ class App extends React.Component {
               marginLeft: 350 /* but responsive */
             }}
           >
-            {loading ? null : <Graph />}
+            {loading ? null : <Graph key={graphRenderCounter} />}
-
             <Legend />
-            {}
           </div>
         </div>
       </Container>

client/src/components/brushableHistogram/index.js

@@ -215,7 +215,12 @@ class HistogramBrush extends React.Component {
     d3.select(svgRef)
       .append("g")
       .attr("class", "brush")
-      .call(d3.brushX().on("end", this.onBrush(field, x.invert).bind(this)));
+      .call(
+        d3
+          .brushX()
+          .on("brush", this.onBrush(field, x.invert).bind(this))
+          .on("end", this.onBrush(field, x.invert).bind(this))
+      );
 
     /* AXIS */
     d3.select(svgRef)
@@ -243,9 +248,9 @@ class HistogramBrush extends React.Component {
       colorAccessor,
       isUserDefined,
       isDiffExp,
-      avgDiff,
+      logFoldChange,
-      set1AvgExp,
+      pval,
-      set2AvgExp,
+      pvalAdj,
       scatterplotXXaccessor,
       scatterplotYYaccessor,
       zebra
@@ -332,25 +337,17 @@ class HistogramBrush extends React.Component {
             }}
           >
             <span>
-              <strong>1:</strong>
+              <strong>log fold change:</strong>
-              {` ${set1AvgExp.toPrecision(2)}`}
+              {` ${logFoldChange.toPrecision(4)}`}
             </span>
             <span
               style={{
                 marginLeft: 7,
-                backgroundColor: globals.lighterGrey,
                 padding: 2
               }}
             >
-              <strong>2:</strong>
+              <strong>p-value (adj):</strong>
-              {` ${set2AvgExp.toPrecision(2)}`}
+              {pvalAdj < 0.0001 ? " < 0.0001" : ` ${pvalAdj.toFixed(4)}`}
-            </span>
-            <span
-              style={{
-                marginLeft: 7
-              }}
-            >
-              {`Av. Diff: ${avgDiff.toFixed(2)}`}
             </span>
           </div>
         ) : null}

client/src/components/geneExpression/index.js

@@ -152,9 +152,9 @@ class GeneExpression extends React.Component {
                     zebra={index % 2 === 0}
                     ranges={d3.extent(values)}
                     isDiffExp
-                    avgDiff={value[1]}
+                    logFoldChange={value[1]}
-                    set1AvgExp={value[4]}
+                    pval={value[2]}
-                    set2AvgExp={value[5]}
+                    pvalAdj={value[3]}
                   />
                 );
               })

client/src/components/graph/drawPointsRegl.js

@@ -38,7 +38,7 @@ export default function(regl) {
     uniforms: {
       distance: regl.prop("distance"),
       view: regl.prop("view"),
-      projection: () => mat4.perspective([], Math.PI / 2, 1, 0.01, 1000)
+      projection: ({viewportWidth, viewportHeight}) => mat4.perspective([], Math.PI / 2, viewportWidth / viewportHeight, 0.01, 1000)
     },
 
     count: regl.prop("count"),

client/src/components/graph/graph.js

@@ -6,8 +6,7 @@ import { connect } from "react-redux";
 import mat4 from "gl-mat4";
 import _regl from "regl";
 import { Button, AnchorButton, Tooltip } from "@blueprintjs/core";
-import { worldEqUniverse } from "../../util/stateManager/world";
+import * as globals from "../../globals";
-
 import setupSVGandBrushElements from "./setupSVGandBrush";
 import actions from "../../actions";
 import _camera from "../../util/camera";
@@ -30,9 +29,9 @@ class Graph extends React.Component {
     super(props);
     this.count = 0;
     this.inverse = mat4.identity([]);
-    this.graphPaddingTop = 100;
+    this.graphPaddingTop = 0;
     this.graphPaddingBottom = 45;
-    this.graphPaddingRight = 10;
+    this.graphPaddingRight = globals.leftSidebarWidth;
     this.renderCache = {
       positions: null,
       colors: null
@@ -126,18 +125,24 @@ class Graph extends React.Component {
         const glScaleX = scaleLinear([0, 1], [-1, 1]);
         const glScaleY = scaleLinear([0, 1], [1, -1]);
 
+        const offset = [d3.mean(obsLayout.X) - 0.5, d3.mean(obsLayout.Y) - 0.5];
+
         for (
           let i = 0, { positions } = this.renderCache;
           i < cellCount;
           i += 1
         ) {
-          positions[2 * i] = glScaleX(obsLayout.X[i]);
+          positions[2 * i] = glScaleX(obsLayout.X[i] - offset[0]);
-          positions[2 * i + 1] = glScaleY(obsLayout.Y[i]);
+          positions[2 * i + 1] = glScaleY(obsLayout.Y[i] - offset[1]);
         }
         pointBuffer({
           data: this.renderCache.positions,
           dimension: 2
         });
+
+        this.setState({
+          offset
+        });
       }
 
       // Colors for each point - a cached value that only changes when
@@ -196,7 +201,7 @@ class Graph extends React.Component {
         this.handleBrushSelectAction.bind(this),
         this.handleBrushDeselectAction.bind(this),
         responsive,
-        this.graphPaddingTop
+        this.graphPaddingRight
       );
       this.setState({ svg: newSvg, brush });
     }
@@ -251,7 +256,7 @@ class Graph extends React.Component {
     an event on procedural deselect because it is move: null
     */
 
-    const { camera } = this.state;
+    const { camera, offset } = this.state;
     const { dispatch, responsive } = this.props;
 
     if (d3.event.sourceEvent !== null) {
@@ -262,6 +267,7 @@ class Graph extends React.Component {
       https://bl.ocks.org/EfratVil/0e542f5fc426065dd1d4b6daaa345a9f
     */
       const s = d3.event.selection;
+      const gl = this.state.regl._gl;
       /*
       event describing brush position:
       @-------|
@@ -270,19 +276,23 @@ class Graph extends React.Component {
       |-------@
     */
 
+      // get aspect ratio
+      const aspect = gl.drawingBufferWidth / gl.drawingBufferHeight;
+
       // compute inverse view matrix
       const inverse = mat4.invert([], camera.view());
 
       // transform screen coordinates -> cell coordinates
       const invert = pin => {
-        const x = (2 * pin[0]) / (responsive.height - this.graphPaddingTop) - 1;
+        const x =
+          (2 * pin[0]) / (responsive.width - this.graphPaddingRight) - 1;
         const y =
           2 * (1 - pin[1] / (responsive.height - this.graphPaddingTop)) - 1;
         const pout = [
-          x * inverse[14] + inverse[12],
+          x * inverse[14] * aspect + inverse[12],
           y * inverse[14] + inverse[13]
         ];
-        return [(pout[0] + 1) / 2, (pout[1] + 1) / 2];
+        return [(pout[0] + 1) / 2 + offset[0], (pout[1] + 1) / 2 + offset[1]];
       };
 
       const brushCoords = {
@@ -366,6 +376,7 @@ class Graph extends React.Component {
                 style={{ marginRight: 10 }}
                 onClick={() => {
                   dispatch(actions.regraph());
+                  dispatch({ type: "increment graph render counter" });
                 }}
               >
                 subset to current selection
@@ -421,12 +432,10 @@ class Graph extends React.Component {
         </div>
         <div
           style={{
-            marginRight: 50,
-            marginTop: 50,
             zIndex: -9999,
             position: "fixed",
-            right: this.graphPaddingRight,
+            top: 0,
-            bottom: this.graphPaddingBottom
+            right: 0
           }}
         >
           <div
@@ -437,7 +446,7 @@ class Graph extends React.Component {
           />
           <div style={{ padding: 0, margin: 0 }}>
             <canvas
-              width={responsive.height - this.graphPaddingTop}
+              width={responsive.width - this.graphPaddingRight}
               height={responsive.height - this.graphPaddingTop}
               ref={canvas => {
                 this.reglCanvas = canvas;

client/src/components/graph/setupSVGandBrush.js

@@ -12,22 +12,18 @@ export default (
   handleBrushSelectAction,
   handleBrushDeselectAction,
   responsive,
-  graphPaddingTop
+  graphPaddingRight
 ) => {
-  const side = responsive.height - graphPaddingTop;
   const svg = d3
     .select("#graphAttachPoint")
     .append("svg")
-    .attr("width", side)
+    .attr("width", responsive.width - graphPaddingRight)
-    .attr("height", side)
+    .attr("height", responsive.height)
     .attr("class", `${styles.graphSVG}`);
 
   const brush = d3
     .brush()
-    .extent([
+    .extent([[0, 0], [responsive.width - graphPaddingRight, responsive.height]])
-      [0, 0],
-      [responsive.height - graphPaddingTop, responsive.height - graphPaddingTop]
-    ])
     .on("brush", handleBrushSelectAction)
     .on("end", handleBrushDeselectAction);
 

client/src/components/scatterplot/drawPointsRegl.js

@@ -38,14 +38,7 @@ export default function(regl) {
     uniforms: {
       distance: regl.prop("distance"),
       view: regl.prop("view"),
-      projection: (context, props) =>
+      projection: () => mat4.perspective([], Math.PI / 2, 1, 0.01, 1000)
-        mat4.perspective(
-          [],
-          Math.PI / 2,
-          (context.viewportWidth * props.scale) / context.viewportHeight,
-          0.01,
-          1000
-        )
     },
 
     count: regl.prop("count"),

client/src/components/scatterplot/scatterplot.js

@@ -82,12 +82,6 @@ class Scatterplot extends React.Component {
       this.drawAxesSVG(scales.xScale, scales.yScale, svg);
     }
 
-    this.setState({
-      svg,
-      xScale: scales ? scales.xScale : null,
-      yScale: scales ? scales.yScale : null
-    });
-
     const camera = _camera(this.reglCanvas, { scale: true, rotate: false });
     const regl = _regl(this.reglCanvas);
 
@@ -98,43 +92,35 @@ class Scatterplot extends React.Component {
     const colorBuffer = regl.buffer();
     const sizeBuffer = regl.buffer();
 
-    regl.frame(({ viewportWidth, viewportHeight }) => {
+    const reglRender = regl.frame(() => {
-      regl.clear({
+      this.reglDraw(
-        depth: 1,
+        regl,
-        color: [1, 1, 1, 1]
+        drawPoints,
-      });
+        sizeBuffer,
-
+        colorBuffer,
-      drawPoints({
+        pointBuffer,
-        distance: camera.distance,
+        camera
-        color: colorBuffer,
+      );
-        position: pointBuffer,
-        size: sizeBuffer,
-        count: this.count,
-        view: camera.view(),
-        scale: viewportHeight / viewportWidth
-      });
-
       camera.tick();
     });
 
+    this.reglRenderState = "rendering";
+
     this.setState({
       regl,
       sizeBuffer,
       pointBuffer,
-      colorBuffer
+      colorBuffer,
+      svg,
+      xScale: scales ? scales.xScale : null,
+      yScale: scales ? scales.yScale : null,
+      reglRender,
+      camera,
+      drawPoints
     });
   }
 
   componentDidUpdate(prevProps) {
-    const {
-      svg,
-      xScale,
-      yScale,
-      regl,
-      pointBuffer,
-      colorBuffer,
-      sizeBuffer
-    } = this.state;
     const {
       world,
       crossfilter,
@@ -144,6 +130,18 @@ class Scatterplot extends React.Component {
       expressionY,
       colorRGB
     } = this.props;
+    const {
+      reglRender,
+      xScale,
+      yScale,
+      regl,
+      pointBuffer,
+      colorBuffer,
+      sizeBuffer,
+      svg,
+      drawPoints,
+      camera
+    } = this.state;
 
     if (
       world &&
@@ -159,6 +157,11 @@ class Scatterplot extends React.Component {
       this.drawAxesSVG(xScale, yScale, svg);
     }
 
+    if (reglRender && this.reglRenderState === "rendering") {
+      reglRender.cancel();
+      this.reglRenderState = "paused";
+    }
+
     if (
       world &&
       regl &&
@@ -198,6 +201,16 @@ class Scatterplot extends React.Component {
       colorBuffer({ data: colorsBuf, dimension: 3 });
       sizeBuffer({ data: sizesBuf, dimension: 1 });
       this.count = cellCount;
+
+      regl._refresh();
+      this.reglDraw(
+        regl,
+        drawPoints,
+        sizeBuffer,
+        colorBuffer,
+        pointBuffer,
+        camera
+      );
     }
 
     if (
@@ -227,6 +240,22 @@ class Scatterplot extends React.Component {
     };
   }
 
+  reglDraw(regl, drawPoints, sizeBuffer, colorBuffer, pointBuffer, camera) {
+    regl.clear({
+      depth: 1,
+      color: [1, 1, 1, 1]
+    });
+
+    drawPoints({
+      size: sizeBuffer,
+      distance: camera.distance,
+      color: colorBuffer,
+      position: pointBuffer,
+      count: this.count,
+      view: camera.view()
+    });
+  }
+
   drawAxesSVG(xScale, yScale, svg) {
     const { scatterplotYYaccessor, scatterplotXXaccessor } = this.props;
     svg.selectAll("*").remove();

client/src/reducers/controls.js

@@ -54,6 +54,7 @@ const Controls = (
     scatterplotXXaccessor: null, // just easier to read
     scatterplotYYaccessor: null,
     axesHaveBeenDrawn: false,
+    graphRenderCounter: 0 /* integer as <Component key={graphRenderCounter} - a change in key forces a remount */,
     __storedStateForCelllist1__: null /* will need procedural control of brush ie., brush.extent https://bl.ocks.org/micahstubbs/3cda05ca68cba260cb81 */,
     __storedStateForCelllist2__: null
   },
@@ -410,7 +411,13 @@ const Controls = (
         ...state,
         opacityForDeselectedCells: action.data
       };
-
+    case "increment graph render counter": {
+      const c = state.graphRenderCounter + 1;
+      return {
+        ...state,
+        graphRenderCounter: c
+      };
+    }
     /*******************************
           Categorical metadata
     *******************************/

client/src/util/actionHelpers.js

@@ -31,8 +31,7 @@ export const doJsonRequest = async url => {
   const res = await fetch(url, {
     method: "get",
     headers: new Headers({
-      "Content-Type": "application/json",
+      "Content-Type": "application/json"
-      "Accept-Encoding": "gzip, deflate, br"
     })
   });
   if (res.ok && res.headers.get("Content-Type") === "application/json") {

docs/REST_API.md

@@ -75,7 +75,7 @@ For a GET URL query parameter:
 - Annotation name is encoded as `obs:name` or `var:name`<sup>[2](#endnote-2)</sup>.
 - Enumerated values (string, categorical, boolean) are encoded as option lists, ie, `var:tissue=lung, obs:tumor=true`
 - Scalar values (int32, float32) are encoded as ranges, ie, `obs:num_reads=1000,10000` where either min or max may be replaced with an asterisk to indicate a half-open range.
-- Index filters are not be allowed within GET URL query parameter filters
+- Index filters are not allowed within GET URL query parameter filters
 - Logically, filters are ANDed, except for repeated annotation names which are ORed. For example, `?X=A&X=B&Y=1` is evaluated as `((X==A or X==B) and Y==1)`
 
 Example selection for _lung_ and _heart_ tissue with more than 1000 reads:
@@ -505,10 +505,10 @@ Generate differential expression (DE) statistics for two specified subsets of da
 
 Two modes are provided:
 
-- Return top N differentially expressed variables (genes)
+- `topN`: return top N differentially expressed variables (across all variables)
-- Return DE for caller-provided variable filter (future)
+- `varFilter`: return DE for caller-provided variable filter (_future_)
 
-Both modes perform calculations using a subset of observations, where each subset is defined by an observation filter (`set1` and `set2`).
+Both modes perform calculations using a subset of observations, where each subset is defined by an observation filter (`set1` and `set2`). These filters must not include a variable filter.
 
 If differential expression is not supported by the server, must return an HTTP 501 response. If, in the view of the server, the request will exceed a reasonable interactive time period, must immediately return HTTP 403 error (error return _before_ attempting computation).
 
@@ -568,24 +568,22 @@ If differential expression is not supported by the server, must return an HTTP 5
 
 **Response body:**
 
-- For 200 Success, differential expression statistics returned as array of arrays sorted by obs index, where each contains the following values:
+- For 200 Success, differential expression statistics returned as array of arrays sorted by varindex, where each contains the following values:
 
   - **varIndex**: variable index for the computed results
-  - **avgDiff**: log fold-change of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group,
+  - **logfoldchange**: log fold-change of the average expression between the two groups. Positive values indicate that the gene is more highly expressed in the first group,
   - **pVal**: unadjusted p-value,
-  - **pValAdj**: Adjusted p-value, based on bonferroni correction using all genes in the original dataset),
+  - **pValAdj**: adjusted p-value
-  - **set1AvgExp:** average expression value for all observations in set 1,
-  - **set2AvgExp**: average expression value for all observations in set 2
 
   Statistics are encoded as an array of arrays, with fields ordered as:
 
-  _varIndex_, _avgDiff_, _pVal_, _pValAdj_, _set1AvgExp_, _set2AvgExp_
+  _varIndex_, _logfoldchange_, _pVal_, _pValAdj_
 
   For example:
 
   ```
   [
-    [ 328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9 ],
+    [ 1720, 2.4679039, 2.3124478092035228e-175, 4.250279073316075e-172 ]
     // ...
   ]
   ```
@@ -616,8 +614,8 @@ POST /diffexp/obs
 200 - Success
 {
   "diffexp": [
-    [ 328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9 ],
+    [ 328, -2.569489, 2.655706e-63, 3.642036e-57 ],
-    // [ varIdx, avgDiff, pVal, pValAdj, set1AvgExp, set2AvgExp ],
+    // [ varIdx, logfoldchange, pVal, pValAdj ],
     // ...
   ]
 }
@@ -690,10 +688,9 @@ Routes:
 - `GET /schema`
 - `GET /annotations/obs`
 - `GET /annotations/var`
-- `GET /layout/obs`
+- `GET /layout/obs` - get the default layout
 - `PUT /data/obs` - request will contain a filter by var `name`
-- `POST /diffexp/obs` - mode `topN`, typically with a couple of 10, and two sets defined by an obs index filter (`{ filter: { obs: { index: [...] } } }`)
+- `POST /diffexp/obs` - mode `topN`, typically with a `count` of 10, and two sets defined by an obs index filter (`{ filter: { obs: { index: [...] } } }`)
-- `PUT /layout/obs` - (_coming soon_) request will contain a filter by obs index
 
 Requests include the following content negotiation headers:
 

docs/release_process.md

@@ -38,7 +38,7 @@ Follow these steps to create a release.
     - [optional] upload the package to test pypi
       `twine upload --repository-url https://test.pypi.org/legacy/ dist/*`
     - [optional] test the test installation in a fresh virtual environment using
-      `pip install --index-url https://test.pypi.org/simple/ cellxgene`
+      `pip install --index-url https://test.pypi.org/simple/ --extra-index-url https://pypi.org/simple cellxgene`
     - upload the package to real pypi using `twine upload dist/*`
     - [optional] test the installation in a fresh virtual environment using
       `pip install cellxgene`

server/app/driver/driver.py

@@ -83,16 +83,18 @@ class CXGDriver(metaclass=ABCMeta):
         pass
 
     @abstractmethod
-    def diffexp(self, filter1, filter2, top_n=None, interactive_limit=None):
+    def diffexp_topN(self, obsFilter1, obsFilter2, top_n=None, interactive_limit=None):
         """
-        Computes the top differentially expressed variables between two observation sets. If dataframes
+        Computes the top N differentially expressed variables between two observation sets. If mode
+        is "TOP_N", then stats for the top N
+        dataframes
         contain a subset of variables, then statistics for all variables will be returned, otherwise
         only the top N vars will be returned.
-        :param filter1: filter: dictionary with filter params for first set of observations
+        :param obsFilter1: filter: dictionary with filter params for first set of observations
-        :param filter2: filter: dictionary with filter params for second set of observations
+        :param obsFilter2: filter: dictionary with filter params for second set of observations
         :param top_n: Limit results to top N (Top var mode only)
         :param interactive_limit: -- don't compute if total # genes in dataframes are larger than this
-        :return: top genes, stats and expression values for variables
+        :return: top N genes and corresponding stats
         """
         pass
 

server/app/rest_api/rest.py

@@ -555,11 +555,11 @@ class DiffExpObsAPI(Resource):
         "responses": {
             "200": {
                 "description": "Statistics are encoded as an array of arrays, with fields ordered as: "
-                               "varIndex, avgDiff,  pVal, pValAdj, set1AvgExp, set2AvgExp",
+                               "varIndex, logfoldchange,  pVal, pValAdj",
                 "examples": {
                     "application/json": [
-                        [328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9],
+                        [328, -2.569489, 2.655706e-63, 3.642036e-57],
-                        [1250, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9],
+                        [1250, -2.569489, 2.655706e-63, 3.642036e-57],
                     ]
                 }
             },
@@ -584,11 +584,12 @@ class DiffExpObsAPI(Resource):
         except ValueError:
             return make_response(f"Error: invalid mode option {args['mode']}", HTTPStatus.BAD_REQUEST)
         # Validate filters
-        if mode == DiffExpMode.VAR_FILTER:
+        if mode == DiffExpMode.VAR_FILTER or "varFilter" in args:
-            if "varFilter" not in args:
+            # not NOT_IMPLEMENTED
-                return make_response("varFilter is required when mode is set to varFilter ", HTTPStatus.BAD_REQUEST)
+            return make_response("mode=varfilter not implemented", HTTPStatus.NOT_IMPLEMENTED)
-            if Axis.OBS in args["varFilter"]["filter"]:
+        if mode == DiffExpMode.TOP_N and "count" not in args:
-                return make_response("Obs filter not allowed in varFilter", HTTPStatus.BAD_REQUEST)
+            return make_response("mode=topN requires a count parameter", HTTPStatus.BAD_REQUEST)
+
         if "set1" not in args:
             return make_response("set1 is required.", HTTPStatus.BAD_REQUEST)
         if Axis.VAR in args["set1"]["filter"]:
@@ -598,16 +599,17 @@ class DiffExpObsAPI(Resource):
             return make_response("Set2 as inverse of set1 is not implemented", HTTPStatus.NOT_IMPLEMENTED)
         if Axis.VAR in args["set2"]["filter"]:
             return make_response("Var filter not allowed for set2", HTTPStatus.BAD_REQUEST)
+
         set1_filter = args["set1"]["filter"]
         set2_filter = args.get("set2", {"filter": {}})["filter"]
-        if "varFilter" in args:
+
-            set1_filter[Axis.VAR] = args["varFilter"]["filter"][Axis.VAR]
+        # TODO: implement varfilter mode
-            set2_filter[Axis.VAR] = args["varFilter"]["filter"][Axis.VAR]
+
-        # mode
+        # mode=topN
         count = args.get("count", None)
         try:
-            diffexp = current_app.data.diffexp(set1_filter, set2_filter, count,
+            diffexp = current_app.data.diffexp_topN(set1_filter, set2_filter, count,
-                                               current_app.data.features["diffexp"]["interactiveLimit"])
+                                                    current_app.data.features["diffexp"]["interactiveLimit"])
         except (ValueError, FilterError) as e:
             return make_response(e.message, HTTPStatus.BAD_REQUEST)
         except InteractiveError:

server/app/scanpy_engine/diffexp.py

@@ -0,0 +1,83 @@
+
+import numpy as np
+from scipy import sparse, stats
+
+
+# Convenience function which handles sparse data
+def _mean_var_n(X):
+    """
+    Two-pass variance calculation.  Numerically (more) stable
+    than naive methods (and same method used by numpy.var())
+    https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Two-pass
+    """
+    n = X.shape[0]
+    if sparse.issparse(X):
+        mean = X.mean(axis=0).A1
+        dfm = X - mean
+        sumsq = np.sum(np.multiply(dfm, dfm), axis=0).A1
+        v = sumsq / (n - 1)
+    else:
+        mean = X.mean(axis=0)
+        dfm = X - mean
+        sumsq = np.sum(np.multiply(dfm, dfm), axis=0)
+        v = sumsq / (n - 1)
+
+    return mean, v, n
+
+
+def diffexp_ttest(adata, maskA, maskB, top_n=8):
+    """
+    Return differential expression statistics for top N variables, sorted by
+    t statistic.   Implemented as a unequal variance t-test.
+
+    :param adata: anndata dataframe
+    :param maskA: observation selection mask for set 1
+    :param maskB: observation selection mask for set 2
+    :param top_n: number of variables to return stats for
+    :return:  for top N genes, [ varindex, logfoldchange, pval, pval_adj ]
+    """
+
+    # mean, variance, N
+    meanA, vA, nA = _mean_var_n(adata._X[maskA])
+    meanB, vB, nB = _mean_var_n(adata._X[maskB])
+
+    # variance / N
+    vnA = vA / nA
+    vnB = vB / nB
+    sum_vn = vnA + vnB
+
+    # degrees of freedom for Welch's t-test
+    with np.errstate(divide='ignore', invalid='ignore'):
+        dof = sum_vn**2 / (vnA**2 / (nA - 1) + vnB**2 / (nB - 1))
+    dof[np.isnan(dof)] = 1
+
+    # Welch's t-test score calculation
+    with np.errstate(divide='ignore', invalid='ignore'):
+        tscores = (meanA - meanB) / np.sqrt(sum_vn)
+    tscores[np.isnan(tscores)] = 0
+
+    # p-value
+    pvals = stats.t.sf(np.abs(tscores), dof) * 2
+    pvals_adj = pvals * adata._X.shape[1]
+
+    # logfoldchanges: log2(meanA / meanB)
+    logfoldchanges = np.log2(np.abs((meanA + 1e-9) / (meanB + 1e-9)))
+
+    # top n sort
+    stats_to_sort = np.abs(tscores)
+    partition = np.argpartition(stats_to_sort, -top_n)[-top_n:]
+    rel_sort_order = np.argsort(stats_to_sort[partition])[::-1]
+    vars_indices = np.arange(adata.n_vars, dtype=int)
+    sort_order = vars_indices[partition][rel_sort_order]
+
+    # top n slice
+    logfoldchanges_top_n = logfoldchanges[sort_order]
+    pvals_top_n = pvals[sort_order]
+    pvals_adj_top_n = pvals_adj[sort_order]
+
+    # varIndex, logfoldchange, pval, pval_adj
+    result = [[sort_order[i],
+               logfoldchanges_top_n[i],
+               pvals_top_n[i],
+               pvals_adj_top_n[i]] for i in range(top_n)]
+    return result

server/app/scanpy_engine/scanpy_engine.py

@@ -4,11 +4,12 @@ import numpy as np
 from pandas import DataFrame
 from pandas.core.dtypes.dtypes import CategoricalDtype
 import scanpy.api as sc
-from scipy import stats, sparse
+from scipy import sparse
 
 from server.app.driver.driver import CXGDriver
-from server.app.util.constants import Axis, DEFAULT_TOP_N, DiffExpMode
+from server.app.util.constants import Axis, DEFAULT_TOP_N
 from server.app.util.errors import FilterError, InteractiveError, PrepareError, ScanpyFileError
+from server.app.scanpy_engine.diffexp import diffexp_ttest
 
 """
 Sort order for methods
@@ -114,35 +115,6 @@ class ScanpyEngine(CXGDriver):
                                   f"that your input and try again.")
         return result
 
-    @staticmethod
-    def _top_sort(values, sort_order, top_n=None):
-        """
-        Sorts an iterable in sort order limited by top_n
-        :param values: iterable of values to sort
-        :param sort_order: ndarray order to sort in
-        :param top_n: cutoff number to return
-        :return: values sorted by sort_order limited by top_n
-        """
-        return values[sort_order][:top_n]
-
-    @staticmethod
-    def _nan_to_one(values):
-        """
-        Replaces NaN values with 1
-        :param values: numpy ndarray
-        :return: ndarray
-        """
-        return np.where(np.isnan(values), 1, values)
-
-    @staticmethod
-    def _nan_to_zero(values):
-        """
-        Replaces NaN values with 0
-        :param values: numpy ndarray
-        :return: ndarray
-        """
-        return np.where(np.isnan(values), 0, values)
-
     def _validate_data_types(self):
         if self.data.X.dtype != "float32":
             warnings.warn(f"Scanpy data matrix is in {self.data.X.dtype} format not float32. "
@@ -180,7 +152,7 @@ class ScanpyEngine(CXGDriver):
                 f"`cellxgene prepare --layout {self.layout_method} <datafile>` "
                 f"to solve this problem. ")
 
-    def filter_dataframe(self, filter, include_uns=False):
+    def filter_dataframe(self, filter):
         """
          Filter cells from data and return a subset of the data. They can operate on both obs and var dimension with
          indexing and filtering by annotation value. Filters are combined with the and operator.
@@ -189,70 +161,68 @@ class ScanpyEngine(CXGDriver):
          https://docs.google.com/document/d/1Fxjp1SKtCk7l8QP9-7KAjGXL0eldi_qEnNT0NmlGzXI/edit#heading=h.8qc9q57amldx
 
         :param filter: dictionary with filter params
-        :param include_uns: bool, include unstructured annotations
         :return: View into scanpy object with cells/genes filtered
         """
         if not filter:
             return self.data
-        cells_idx = np.ones((self.cell_count,), dtype=bool)
+        obs_selector, var_selector = self._filter_to_mask(filter, use_slices=False)
-        genes_idx = np.ones((self.gene_count,), dtype=bool)
+        data = self._slice(self.data, obs_selector, var_selector)
-        if Axis.OBS in filter:
-            if "index" in filter["obs"]:
-                cells_idx = self._filter_index(filter["obs"]["index"], cells_idx, Axis.OBS)
-            if "annotation_value" in filter["obs"]:
-                cells_idx = self._filter_annotation(filter["obs"]["annotation_value"], cells_idx, Axis.OBS)
-        if Axis.VAR in filter:
-            if "index" in filter["var"]:
-                genes_idx = self._filter_index(filter["var"]["index"], genes_idx, Axis.VAR)
-            if "annotation_value" in filter["var"]:
-                genes_idx = self._filter_annotation(filter["var"]["annotation_value"], genes_idx, Axis.VAR)
-
-        data = self._slice(self.data, cells_idx, genes_idx)
         return data
 
-    def _filter_index(self, filter, index, axis):
+    @staticmethod
-        """
+    def _annotation_filter_to_mask(filter, d_axis, count):
-        Filter data based on index. ex. [1, 3, [111:200]]
+        mask = np.ones((count, ), dtype=bool)
-        :param filter: subset of filter dict for obs/var:index
-        :param index: np logical vector containing true for passing false for failing filter
-        :param axis: string obs or var
-        :return: np logical vector for whether the data passes the filter
-        """
-        if axis == Axis.OBS:
-            count_ = self.cell_count
-        elif axis == Axis.VAR:
-            count_ = self.gene_count
-        idx_filter = np.zeros((count_,), dtype=bool)
-        for i in filter:
-            if type(i) == list:
-                idx_filter[i[0]:i[1]] = True
-            else:
-                idx_filter[i] = True
-        return np.logical_and(index, idx_filter)
-
-    def _filter_annotation(self, filter, index, axis):
-        """
-        Filter data based on annotation value
-        :param filter: subset of filter dict for obs/var:annotation_value
-        :param index: np logical vector containing true for passing false for failing filter
-        :param axis: string obs or var
-        :return: np logical vector for whether the data passes the filter
-        """
-        d_axis = getattr(self.data, axis.value)
         for v in filter:
             if d_axis[v["name"]].dtype.name in ["boolean", "category", "object"]:
                 key_idx = np.in1d(getattr(d_axis, v["name"]), v["values"])
-                index = np.logical_and(index, key_idx)
+                mask = np.logical_and(mask, key_idx)
             else:
                 min_ = v.get("min", None)
                 max_ = v.get("max", None)
                 if min_ is not None:
                     key_idx = (getattr(d_axis, v["name"]) >= min_).ravel()
-                    index = np.logical_and(index, key_idx)
+                    mask = np.logical_and(mask, key_idx)
                 if max_ is not None:
                     key_idx = (getattr(d_axis, v["name"]) <= max_).ravel()
-                    index = np.logical_and(index, key_idx)
+                    mask = np.logical_and(mask, key_idx)
-        return index
+        return mask
+
+    @staticmethod
+    def _index_filter_to_mask(filter, count):
+        mask = np.zeros((count, ), dtype=bool)
+        for i in filter:
+            if type(i) == list:
+                mask[i[0]:i[1]] = True
+            else:
+                mask[i] = True
+        return mask
+
+    @staticmethod
+    def _axis_filter_to_mask(filter, d_axis, count):
+        mask = np.ones((count, ), dtype=bool)
+        if "index" in filter:
+            mask = np.logical_and(mask, ScanpyEngine._index_filter_to_mask(filter["index"], count))
+        if "annotation_value" in filter:
+            mask = np.logical_and(mask,
+                                  ScanpyEngine._annotation_filter_to_mask(filter["annotation_value"],
+                                                                          d_axis,
+                                                                          count))
+        return mask
+
+    def _filter_to_mask(self, filter, use_slices=True):
+        if use_slices:
+            obs_selector = slice(0, self.data.n_obs)
+            var_selector = slice(0, self.data.n_vars)
+        else:
+            obs_selector = None
+            var_selector = None
+
+        if filter is not None:
+            if Axis.OBS in filter:
+                obs_selector = self._axis_filter_to_mask(filter["obs"], self.data.obs, self.data.n_obs)
+            if Axis.VAR in filter:
+                var_selector = self._axis_filter_to_mask(filter["var"], self.data.var, self.data.n_vars)
+        return obs_selector, var_selector
 
     @staticmethod
     def _slice(data, obs_selector=None, vars_selector=None):
@@ -293,16 +263,25 @@ class ScanpyEngine(CXGDriver):
         [observation ids, val1, val2...]
         """
         try:
-            df = self.filter_dataframe(filter)
+            obs_selector, var_selector = self._filter_to_mask(filter)
-        except KeyError as e:
+        except (KeyError, IndexError) as e:
             raise FilterError(f"Error parsing filter: {e}") from e
-        df_axis = getattr(df, axis)
+        if axis == Axis.OBS:
-        if not fields:
+            obs = self.data.obs[obs_selector]
-            fields = df_axis.columns.tolist()
+            if not fields:
-        result = {
+                fields = obs.columns.tolist()
-            "names": fields,
+            result = {
-            "data": DataFrame(df_axis[fields]).to_records(index=True).tolist()
+                "names": fields,
-        }
+                "data": DataFrame(obs[fields]).to_records(index=True).tolist()
+            }
+        else:
+            var = self.data.var[var_selector]
+            if not fields:
+                fields = var.columns.tolist()
+            result = {
+                "names": fields,
+                "data": DataFrame(var[fields]).to_records(index=True).tolist()
+            }
         return result
 
     def data_frame(self, filter, axis):
@@ -316,85 +295,38 @@ class ScanpyEngine(CXGDriver):
         }
         """
         try:
-            slice = self.filter_dataframe(filter)
+            obs_selector, var_selector = self._filter_to_mask(filter)
-        except KeyError as e:
+        except (KeyError, IndexError) as e:
             raise FilterError(f"Error parsing filter: {e}") from e
-        # convert sparse slice to dense
+        _X = self.data._X[obs_selector, var_selector]
-        X = slice._X.toarray() if sparse.issparse(slice._X) else slice._X
+        if sparse.issparse(_X):
+            _X = _X.toarray()
+        var_index_sliced = self.data.var.index[var_selector]
+        obs_index_sliced = self.data.obs.index[obs_selector]
         if axis == Axis.OBS:
             result = {
-                "var": slice.var.index.tolist(),
+                "var": var_index_sliced.tolist(),
-                "obs": DataFrame(X, index=slice.obs.index).to_records(index=True).tolist()
+                "obs": DataFrame(_X, index=obs_index_sliced).to_records(index=True).tolist()
             }
         else:
             result = {
-                "obs": slice.obs.index.tolist(),
+                "obs": obs_index_sliced.tolist(),
-                "var": DataFrame(X.T, index=slice.var.index).to_records(index=True).tolist()
+                "var": DataFrame(_X.T, index=var_index_sliced).to_records(index=True).tolist()
             }
         return result
 
-    def diffexp(self, filter1, filter2, top_n=None, interactive_limit=None):
+    def diffexp_topN(self, obsFilterA, obsFilterB, top_n=None, interactive_limit=None):
-        """
+        if Axis.VAR in obsFilterA or Axis.VAR in obsFilterB:
-        Computes the top differentially expressed variables between two observation sets. If dataframes
+            raise FilterError("Observation filters may not contain vaiable conditions")
-        contain a subset of variables, then statistics for all variables will be returned, otherwise
-        only the top N vars will be returned.
-        :param filter1: filter: dictionary with filter params for first set of observations
-        :param filter2: filter: dictionary with filter params for second set of observations
-        :param top_n: Limit results to top N (Top var mode only)
-        :param interactive_limit: -- don't compute if total # genes in dataframes are larger than this
-        :return: top genes, stats and expression values for variables
-        """
         try:
-            df1 = self.filter_dataframe(filter1)
+            obs_mask_A = self._axis_filter_to_mask(obsFilterA["obs"], self.data.obs, self.data.n_obs)
-        except KeyError as e:
+            obs_mask_B = self._axis_filter_to_mask(obsFilterB["obs"], self.data.obs, self.data.n_obs)
-            raise FilterError(f"Error parsing filter for set 1: {e}") from e
+        except (KeyError, IndexError) as e:
-        # TODO df2 should be inverse if not filter2 provided
+            raise FilterError(f"Error parsing filter: {e}") from e
-        try:
+        if top_n is None:
-            df2 = self.filter_dataframe(filter2)
+            top_n = DEFAULT_TOP_N
-        except KeyError as e:
+        result = diffexp_ttest(self.data, obs_mask_A, obs_mask_B, top_n)
-            raise FilterError(f"Error parsing filter for set 2: {e}") from e
+        return sorted(result, key=lambda r: r[0])
-        # If not the same genes, test is wrong!
-        if np.any(df1.var.index != df2.var.index):
-            raise ValueError("Variables ares not the same in set1 and set2")
-        if interactive_limit and df1.shape[0] + df2.shape[0] > interactive_limit:
-            raise InteractiveError("Size of set 1 and 2 is too large for interactive computation")
-        # If not all genes, they used a var filter
-        if df1.var.shape[0] < self.gene_count:
-            mode = DiffExpMode.VAR_FILTER
-            if top_n:
-                raise Warning("Top N was specified but will not be used in 'Var Filter' mode")
-        else:
-            mode = DiffExpMode.TOP_N
-            if not top_n:
-                top_n = DEFAULT_TOP_N
-
-        genes_idx = df1.var.index
-        # ensure we are using a dense ndarray
-        X1 = df1._X.toarray() if sparse.issparse(df1._X) else df1._X
-        X2 = df2._X.toarray() if sparse.issparse(df2._X) else df2._X
-        diffexp_result = stats.ttest_ind(X1, X2)
-        tstats = self._nan_to_zero(diffexp_result.statistic)
-        pval = self._nan_to_one(diffexp_result.pvalue)
-        bonferroni_pval = 1 - (1 - pval) ** self.gene_count
-        ave_exp_set1 = np.mean(X1, axis=0)
-        ave_exp_set2 = np.mean(X2, axis=0)
-        ave_diff = ave_exp_set1 - ave_exp_set2
-        if mode == DiffExpMode.TOP_N:
-            sort_order = np.argsort(np.abs(tstats))[::-1]
-            # If top_n > length it will just return length
-            genes = self._top_sort(genes_idx, sort_order, top_n)
-            pval = self._top_sort(pval, sort_order, top_n)
-            bonferroni_pval = self._top_sort(bonferroni_pval, sort_order, top_n)
-            ave_exp_set1 = self._top_sort(ave_exp_set1, sort_order, top_n)
-            ave_exp_set2 = self._top_sort(ave_exp_set2, sort_order, top_n)
-            ave_diff = self._top_sort(ave_diff, sort_order, top_n)
-
-        # varIndex, avgDiff,  pVal, pValAdj, set1AvgExp, set2AvgExp
-        result = []
-        for i in range(len(genes)):
-            result.append([genes[i], ave_diff[i], pval[i], bonferroni_pval[i], ave_exp_set1[i], ave_exp_set2[i]])
-        # Results need to be returned in var index order
-        return sorted(result, key=lambda gene: gene[0])
 
     def layout(self, filter, interactive_limit=None):
         """
@@ -404,8 +336,8 @@ class ScanpyEngine(CXGDriver):
         :return:  [cellid, x, y, ...]
         """
         try:
-            df = self.filter_dataframe(filter, include_uns=True)
+            df = self.filter_dataframe(filter)
-        except KeyError as e:
+        except (KeyError, IndexError) as e:
             raise FilterError(f"Error parsing filter: {e}") from e
         if interactive_limit and len(df.obs.index) > interactive_limit:
             raise InteractiveError("Size data is too large for interactive computation")

server/requirements.txt

@@ -1,4 +1,4 @@
-anndata>=0.6.12
+anndata>=0.6.13
 click>=6.7
 Flask>=1.0.2
 Flask-Caching>=1.4.0

server/test/test_api.py

@@ -198,6 +198,7 @@ class EndPoints(unittest.TestCase):
         url = f"{URL_BASE}{endpoint}"
         params = {
             "mode": "topN",
+            "count": 10,
             "set1": {
                 "filter": {
                     "obs": {

server/test/test_scanpy_engine.py

@@ -85,7 +85,7 @@ class UtilTest(unittest.TestCase):
                 }
             }
         }
-        data = self.data.filter_dataframe(filter_["filter"], include_uns=False)
+        data = self.data.filter_dataframe(filter_["filter"])
         self.assertEqual(data.shape[1], 1)
 
     def test_filter_complex(self):
@@ -184,7 +184,7 @@ class UtilTest(unittest.TestCase):
         layout = self.data.layout(filter_["filter"])
         self.assertEqual(len(layout["coordinates"]), 497)
 
-    def test_diffexp(self):
+    def test_diffexp_topN(self):
         f1 = {
             "filter": {
                 "obs": {
@@ -199,11 +199,11 @@ class UtilTest(unittest.TestCase):
                 }
             }
         }
-        result = self.data.diffexp(f1["filter"], f2["filter"])
+        result = self.data.diffexp_topN(f1["filter"], f2["filter"])
         self.assertEqual(len(result), 10)
         var_idx = [i[0] for i in result]
         self.assertEqual(var_idx, sorted(var_idx))
-        result = self.data.diffexp(f1["filter"], f2["filter"], 20)
+        result = self.data.diffexp_topN(f1["filter"], f2["filter"], 20)
         self.assertEqual(len(result), 20)
 
     def test_data_frame(self):

setup.py

@@ -1,7 +1,7 @@
 from setuptools import setup, find_packages
 
-with open("README.md", "r") as fh:
+with open("README.md", "rb") as fh:
-    long_description = fh.read()
+    long_description = fh.read().decode()
 
 with open("server/requirements.txt") as fh:
     requirements = fh.read().splitlines()
@@ -16,6 +16,7 @@ setup(
     author_email="cweaver@chanzuckerberg.com",
     description="Web application for exploration of large scale scRNA-seq datasets",
     long_description=long_description,
+    long_description_content_type='text/markdown',
     install_requires=requirements,
     include_package_data=True,
     zip_safe=False,