Merge branch 'master' into csweaver/release1

MANIFEST.in

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

client/src/actions/index.js

@@ -90,16 +90,12 @@ async function _doRequestExpressionData(dispatch, getState, genes) {
   const state = getState();
   const { universe } = state.controls;
   /* preload data already in cache */
-  let expressionData = _.transform(
+  let expressionData = _.transform(genes, (expData, g) => {
-    genes,
+    const data = kvCache.get(universe.varDataCache, g);
-    (expData, g) => {
+    if (data) {
-      const data = kvCache.get(universe.varDataCache, g);
+      expData[g] = data;
-      if (data) {
+    }
-        expData[g] = data;
+  }); // --> { gene: data }
-      }
-    },
-    {}
-  ); // --> { gene: data }
   /* make a list of genes for which we do not have data */
   const genesToFetch = _.filter(genes, g => expressionData[g] === undefined);
 
@@ -123,6 +119,7 @@ async function _doRequestExpressionData(dispatch, getState, genes) {
           }),
           headers: new Headers({
             accept: "application/json",
+            "Accept-Encoding": "gzip, deflate, br",
             "Content-Type": "application/json"
           })
         }
@@ -242,6 +239,7 @@ 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({
@@ -299,9 +297,6 @@ 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,8 +11,7 @@ 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) {
@@ -55,7 +54,7 @@ class App extends React.Component {
   }
 
   render() {
-    const { loading, error, graphRenderCounter } = this.props;
+    const { loading } = this.props;
     return (
       <Container>
         <Helmet title="cellxgene" />
@@ -80,8 +79,10 @@ class App extends React.Component {
               marginLeft: 350 /* but responsive */
             }}
           >
-            {loading ? null : <Graph key={graphRenderCounter} />}
+            {loading ? null : <Graph />}
+
             <Legend />
+            {}
           </div>
         </div>
       </Container>

client/src/components/brushableHistogram/index.js

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

client/src/components/graph/drawPointsRegl.js

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

client/src/components/graph/graph.js

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

client/src/components/graph/setupSVGandBrush.js

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

client/src/components/scatterplot/drawPointsRegl.js

@@ -38,7 +38,14 @@ export default function(regl) {
     uniforms: {
       distance: regl.prop("distance"),
       view: regl.prop("view"),
-      projection: () => mat4.perspective([], Math.PI / 2, 1, 0.01, 1000)
+      projection: (context, props) =>
+        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,6 +82,12 @@ 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);
 
@@ -92,35 +98,43 @@ class Scatterplot extends React.Component {
     const colorBuffer = regl.buffer();
     const sizeBuffer = regl.buffer();
 
-    const reglRender = regl.frame(() => {
+    regl.frame(({ viewportWidth, viewportHeight }) => {
-      this.reglDraw(
+      regl.clear({
-        regl,
+        depth: 1,
-        drawPoints,
+        color: [1, 1, 1, 1]
-        sizeBuffer,
+      });
-        colorBuffer,
+
-        pointBuffer,
+      drawPoints({
-        camera
+        distance: camera.distance,
-      );
+        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,
@@ -130,18 +144,6 @@ class Scatterplot extends React.Component {
       expressionY,
       colorRGB
     } = this.props;
-    const {
-      reglRender,
-      xScale,
-      yScale,
-      regl,
-      pointBuffer,
-      colorBuffer,
-      sizeBuffer,
-      svg,
-      drawPoints,
-      camera
-    } = this.state;
 
     if (
       world &&
@@ -157,11 +159,6 @@ class Scatterplot extends React.Component {
       this.drawAxesSVG(xScale, yScale, svg);
     }
 
-    if (reglRender && this.reglRenderState === "rendering") {
-      reglRender.cancel();
-      this.reglRenderState = "paused";
-    }
-
     if (
       world &&
       regl &&
@@ -201,16 +198,6 @@ 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 (
@@ -240,22 +227,6 @@ 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,7 +54,6 @@ 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
   },
@@ -411,13 +410,7 @@ 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,7 +31,8 @@ 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 allowed within GET URL query parameter filters
+- Index filters are not be 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:
 
-- `topN`: return top N differentially expressed variables (across all variables)
+- Return top N differentially expressed variables (genes)
-- `varFilter`: return DE for caller-provided variable filter (_future_)
+- 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`). These filters must not include a variable filter.
+Both modes perform calculations using a subset of observations, where each subset is defined by an observation filter (`set1` and `set2`).
 
 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,22 +568,24 @@ 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 varindex, where each contains the following values:
+- For 200 Success, differential expression statistics returned as array of arrays sorted by obs index, where each contains the following values:
 
   - **varIndex**: variable index for the computed results
-  - **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,
+  - **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,
   - **pVal**: unadjusted p-value,
-  - **pValAdj**: adjusted p-value
+  - **pValAdj**: Adjusted p-value, based on bonferroni correction using all genes in the original dataset),
+  - **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_, _logfoldchange_, _pVal_, _pValAdj_
+  _varIndex_, _avgDiff_, _pVal_, _pValAdj_, _set1AvgExp_, _set2AvgExp_
 
   For example:
 
   ```
   [
-    [ 1720, 2.4679039, 2.3124478092035228e-175, 4.250279073316075e-172 ]
+    [ 328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9 ],
     // ...
   ]
   ```
@@ -614,8 +616,8 @@ POST /diffexp/obs
 200 - Success
 {
   "diffexp": [
-    [ 328, -2.569489, 2.655706e-63, 3.642036e-57 ],
+    [ 328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9 ],
-    // [ varIdx, logfoldchange, pVal, pValAdj ],
+    // [ varIdx, avgDiff, pVal, pValAdj, set1AvgExp, set2AvgExp ],
     // ...
   ]
 }
@@ -688,9 +690,10 @@ Routes:
 - `GET /schema`
 - `GET /annotations/obs`
 - `GET /annotations/var`
-- `GET /layout/obs` - get the default layout
+- `GET /layout/obs`
 - `PUT /data/obs` - request will contain a filter by var `name`
-- `POST /diffexp/obs` - mode `topN`, typically with a `count` of 10, and two sets defined by an obs index filter (`{ filter: { obs: { index: [...] } } }`)
+- `POST /diffexp/obs` - mode `topN`, typically with a couple 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/ --extra-index-url https://pypi.org/simple cellxgene`
+      `pip install --index-url https://test.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,18 +83,16 @@ class CXGDriver(metaclass=ABCMeta):
         pass
 
     @abstractmethod
-    def diffexp_topN(self, obsFilter1, obsFilter2, top_n=None, interactive_limit=None):
+    def diffexp(self, filter1, filter2, top_n=None, interactive_limit=None):
         """
-        Computes the top N differentially expressed variables between two observation sets. If mode
+        Computes the top differentially expressed variables between two observation sets. If dataframes
-        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 obsFilter1: filter: dictionary with filter params for first set of observations
+        :param filter1: filter: dictionary with filter params for first set of observations
-        :param obsFilter2: filter: dictionary with filter params for second 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 N genes and corresponding stats
+        :return: top genes, stats and expression values for variables
         """
         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, logfoldchange,  pVal, pValAdj",
+                               "varIndex, avgDiff,  pVal, pValAdj, set1AvgExp, set2AvgExp",
                 "examples": {
                     "application/json": [
-                        [328, -2.569489, 2.655706e-63, 3.642036e-57],
+                        [328, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9],
-                        [1250, -2.569489, 2.655706e-63, 3.642036e-57],
+                        [1250, -2.569489, 2.655706e-63, 3.642036e-57, 383.393, 583.9],
                     ]
                 }
             },
@@ -584,12 +584,11 @@ 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 or "varFilter" in args:
+        if mode == DiffExpMode.VAR_FILTER:
-            # not NOT_IMPLEMENTED
+            if "varFilter" not in args:
-            return make_response("mode=varfilter not implemented", HTTPStatus.NOT_IMPLEMENTED)
+                return make_response("varFilter is required when mode is set to varFilter ", HTTPStatus.BAD_REQUEST)
-        if mode == DiffExpMode.TOP_N and "count" not in args:
+            if Axis.OBS in args["varFilter"]["filter"]:
-            return make_response("mode=topN requires a count parameter", HTTPStatus.BAD_REQUEST)
+                return make_response("Obs filter not allowed in varFilter", HTTPStatus.BAD_REQUEST)
-
         if "set1" not in args:
             return make_response("set1 is required.", HTTPStatus.BAD_REQUEST)
         if Axis.VAR in args["set1"]["filter"]:
@@ -599,17 +598,16 @@ 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:
-        # TODO: implement varfilter mode
+            set1_filter[Axis.VAR] = args["varFilter"]["filter"][Axis.VAR]
-
+            set2_filter[Axis.VAR] = args["varFilter"]["filter"][Axis.VAR]
-        # mode=topN
+        # mode
         count = args.get("count", None)
         try:
-            diffexp = current_app.data.diffexp_topN(set1_filter, set2_filter, count,
+            diffexp = current_app.data.diffexp(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

@@ -1,83 +0,0 @@
-
-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,12 +4,11 @@ import numpy as np
 from pandas import DataFrame
 from pandas.core.dtypes.dtypes import CategoricalDtype
 import scanpy.api as sc
-from scipy import sparse
+from scipy import stats, sparse
 
 from server.app.driver.driver import CXGDriver
-from server.app.util.constants import Axis, DEFAULT_TOP_N
+from server.app.util.constants import Axis, DEFAULT_TOP_N, DiffExpMode
 from server.app.util.errors import FilterError, InteractiveError, PrepareError, ScanpyFileError
-from server.app.scanpy_engine.diffexp import diffexp_ttest
 
 """
 Sort order for methods
@@ -115,6 +114,35 @@ 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. "
@@ -152,7 +180,7 @@ class ScanpyEngine(CXGDriver):
                 f"`cellxgene prepare --layout {self.layout_method} <datafile>` "
                 f"to solve this problem. ")
 
-    def filter_dataframe(self, filter):
+    def filter_dataframe(self, filter, include_uns=False):
         """
          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.
@@ -161,68 +189,70 @@ 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
-        obs_selector, var_selector = self._filter_to_mask(filter, use_slices=False)
+        cells_idx = np.ones((self.cell_count,), dtype=bool)
-        data = self._slice(self.data, obs_selector, var_selector)
+        genes_idx = np.ones((self.gene_count,), dtype=bool)
+        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
 
-    @staticmethod
+    def _filter_index(self, filter, index, axis):
-    def _annotation_filter_to_mask(filter, d_axis, count):
+        """
-        mask = np.ones((count, ), dtype=bool)
+        Filter data based on index. ex. [1, 3, [111:200]]
+        :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"])
-                mask = np.logical_and(mask, key_idx)
+                index = np.logical_and(index, 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()
-                    mask = np.logical_and(mask, key_idx)
+                    index = np.logical_and(index, key_idx)
                 if max_ is not None:
                     key_idx = (getattr(d_axis, v["name"]) <= max_).ravel()
-                    mask = np.logical_and(mask, key_idx)
+                    index = np.logical_and(index, key_idx)
-        return mask
+        return index
-
-    @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):
@@ -263,25 +293,16 @@ class ScanpyEngine(CXGDriver):
         [observation ids, val1, val2...]
         """
         try:
-            obs_selector, var_selector = self._filter_to_mask(filter)
+            df = self.filter_dataframe(filter)
-        except (KeyError, IndexError) as e:
+        except KeyError as e:
             raise FilterError(f"Error parsing filter: {e}") from e
-        if axis == Axis.OBS:
+        df_axis = getattr(df, axis)
-            obs = self.data.obs[obs_selector]
+        if not fields:
-            if not fields:
+            fields = df_axis.columns.tolist()
-                fields = obs.columns.tolist()
+        result = {
-            result = {
+            "names": fields,
-                "names": fields,
+            "data": DataFrame(df_axis[fields]).to_records(index=True).tolist()
-                "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):
@@ -295,38 +316,85 @@ class ScanpyEngine(CXGDriver):
         }
         """
         try:
-            obs_selector, var_selector = self._filter_to_mask(filter)
+            slice = self.filter_dataframe(filter)
-        except (KeyError, IndexError) as e:
+        except KeyError as e:
             raise FilterError(f"Error parsing filter: {e}") from e
-        _X = self.data._X[obs_selector, var_selector]
+        # convert sparse slice to dense
-        if sparse.issparse(_X):
+        X = slice._X.toarray() if sparse.issparse(slice._X) else slice._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": var_index_sliced.tolist(),
+                "var": slice.var.index.tolist(),
-                "obs": DataFrame(_X, index=obs_index_sliced).to_records(index=True).tolist()
+                "obs": DataFrame(X, index=slice.obs.index).to_records(index=True).tolist()
             }
         else:
             result = {
-                "obs": obs_index_sliced.tolist(),
+                "obs": slice.obs.index.tolist(),
-                "var": DataFrame(_X.T, index=var_index_sliced).to_records(index=True).tolist()
+                "var": DataFrame(X.T, index=slice.var.index).to_records(index=True).tolist()
             }
         return result
 
-    def diffexp_topN(self, obsFilterA, obsFilterB, top_n=None, interactive_limit=None):
+    def diffexp(self, filter1, filter2, top_n=None, interactive_limit=None):
-        if Axis.VAR in obsFilterA or Axis.VAR in obsFilterB:
+        """
-            raise FilterError("Observation filters may not contain vaiable conditions")
+        Computes the top differentially expressed variables between two observation sets. If 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 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:
-            obs_mask_A = self._axis_filter_to_mask(obsFilterA["obs"], self.data.obs, self.data.n_obs)
+            df1 = self.filter_dataframe(filter1)
-            obs_mask_B = self._axis_filter_to_mask(obsFilterB["obs"], self.data.obs, self.data.n_obs)
+        except KeyError as e:
-        except (KeyError, IndexError) as e:
+            raise FilterError(f"Error parsing filter for set 1: {e}") from e
-            raise FilterError(f"Error parsing filter: {e}") from e
+        # TODO df2 should be inverse if not filter2 provided
-        if top_n is None:
+        try:
-            top_n = DEFAULT_TOP_N
+            df2 = self.filter_dataframe(filter2)
-        result = diffexp_ttest(self.data, obs_mask_A, obs_mask_B, top_n)
+        except KeyError as e:
-        return sorted(result, key=lambda r: r[0])
+            raise FilterError(f"Error parsing filter for set 2: {e}") from e
+        # 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):
         """
@@ -336,8 +404,8 @@ class ScanpyEngine(CXGDriver):
         :return:  [cellid, x, y, ...]
         """
         try:
-            df = self.filter_dataframe(filter)
+            df = self.filter_dataframe(filter, include_uns=True)
-        except (KeyError, IndexError) as e:
+        except KeyError 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.13
+anndata>=0.6.12
 click>=6.7
 Flask>=1.0.2
 Flask-Caching>=1.4.0

server/test/test_api.py

@@ -198,7 +198,6 @@ 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"])
+        data = self.data.filter_dataframe(filter_["filter"], include_uns=False)
         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_topN(self):
+    def test_diffexp(self):
         f1 = {
             "filter": {
                 "obs": {
@@ -199,11 +199,11 @@ class UtilTest(unittest.TestCase):
                 }
             }
         }
-        result = self.data.diffexp_topN(f1["filter"], f2["filter"])
+        result = self.data.diffexp(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_topN(f1["filter"], f2["filter"], 20)
+        result = self.data.diffexp(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", "rb") as fh:
+with open("README.md", "r") as fh:
-    long_description = fh.read().decode()
+    long_description = fh.read()
 
 with open("server/requirements.txt") as fh:
     requirements = fh.read().splitlines()
@@ -16,7 +16,6 @@ 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,