Integrate Harmony

Use rJava instead of outdated packages from Maven Central.

Author: arteymix

environment.yml

@@ -0,0 +1,15 @@
+name: gemma
+dependencies:
+  - r=4.4
+  - r-rserve
+  - r-rjava
+  - repeatmasker
+  - blast
+  - sra-tools
+  - hdf5=1.12
+  - python
+  - pip
+  - pip:
+    - anndata
+    - scipy
+    - numpy

gemma-core/pom.xml

@@ -251,6 +251,38 @@
             <version>${hdf5.version}</version>
         </dependency>
 
+        <!-- rJava -->
+        <dependency>
+            <groupId>org.rosuda.REngine</groupId>
+            <artifactId>REngine</artifactId>
+            <version>${rJava.version}</version>
+            <optional>true</optional>
+        </dependency>
+        <dependency>
+            <groupId>org.rosuda.REngine</groupId>
+            <artifactId>JRI</artifactId>
+            <version>${rJava.version}</version>
+            <optional>true</optional>
+        </dependency>
+        <dependency>
+            <groupId>org.rosuda.REngine</groupId>
+            <artifactId>JRIEngine</artifactId>
+            <version>${rJava.version}</version>
+            <optional>true</optional>
+        </dependency>
+
+        <dependency>
+            <groupId>org.rosuda.REngine</groupId>
+            <artifactId>Rserve</artifactId>
+            <version>${rserve.version}</version>
+            <optional>true</optional>
+        </dependency>
+        <dependency>
+            <groupId>com.kohlschutter.junixsocket</groupId>
+            <artifactId>junixsocket-core</artifactId>
+            <version>2.10.1</version>
+        </dependency>
+
         <!-- Gemma Slack Bot -->
         <dependency>
             <groupId>com.slack.api</groupId>

gemma-core/src/main/java/ubic/gemma/core/analysis/service/ExpressionDataFileUtils.java

@@ -25,7 +25,9 @@ public class ExpressionDataFileUtils {
     private static final String SC_DATA_SUFFIX = ".scdata";
     public static final String MEX_SC_DATA_SUFFIX = SC_DATA_SUFFIX + ".mex";
     public static final String TABULAR_SC_DATA_SUFFIX = SC_DATA_SUFFIX + ".tsv.gz";
-    public static final String CELL_BROWSER_SC_DATA_SUFFIX = SC_DATA_SUFFIX + ".cellbrowser.tsv.gz";
+    private static final String SC_METADATA_SUFFIX = ".scmetadata";
+    public static final String TABULAR_SC_METADATA_SUFFIX = SC_METADATA_SUFFIX + ".tsv.gz";
+    public static final String CELL_BROWSER_SC_DATA_SUFFIX = SC_METADATA_SUFFIX + ".cellbrowser.tsv.gz";
 
     // for single-cell metadata
     private static final String SC_METADATA_SUFFIX = ".scmetadata";

gemma-core/src/main/java/ubic/gemma/core/analysis/singleCell/batcheffect/Harmony.java

@@ -1,16 +1,51 @@
 package ubic.gemma.core.analysis.singleCell.batcheffect;
 
+import org.rosuda.REngine.REXP;
+import org.rosuda.REngine.REXPNull;
 import org.springframework.util.Assert;
 import ubic.gemma.core.datastructure.matrix.SingleCellDesignMatrix;
 import ubic.gemma.core.datastructure.matrix.SingleCellExpressionDataMatrix;
+import ubic.gemma.core.util.r.RClient;
+import ubic.gemma.core.util.r.REngineFactory;
 
+/**
+ * Perform batch correction using the <a href="">Harmony algorithm</a>.
+ * <p>
+ * Requirements: an R engine with the Harmony R package installed.
+ * @author poirigui
+ */
 class Harmony implements BatchCorrection {
 
+    private final REngineFactory rEngineFactory;
+
+    Harmony( REngineFactory rEngineFactory ) {
+        this.rEngineFactory = rEngineFactory;
+    }
+
     @Override
     public SingleCellExpressionDataMatrix<?> perform( SingleCellExpressionDataMatrix<?> dataMatrix, SingleCellDesignMatrix singleCellDesignMatrix ) {
         Assert.isTrue( dataMatrix.getBioAssays().equals( singleCellDesignMatrix.getBioAssays() ),
                 "Assays in the data matrix must match exactly those of the design matrix." );
-        // TODO: serialize both matrices to disk and call Harmony R package
-        return dataMatrix;
+        try ( RClient rEngine = new RClient( rEngineFactory ) ) {
+            // TODO: serialize both matrices to disk and call Harmony R package
+            rEngine.parseAndEval( "library(harmony);" );
+            // rEngine.assignDataFrame( "dataMatrix", toDataFrame( dataMatrix ) );
+            // rEngine.assignDataFrame( "designMatrix", toDataFrame( singleCellDesignMatrix ) );
+            //language=R
+            return fromDataFrame( rEngine.parseAndEval( "harmony::HarmonyMatrix(dataMatrix, designMatrix);" ) );
+        }
+    }
+
+    private REXP toDataFrame( SingleCellExpressionDataMatrix<?> dataMatrix ) {
+        // Convert the SingleCellExpressionDataMatrix to an REXP object
+        return new REXPNull();
+    }
+
+    private REXP toDataFrame( SingleCellDesignMatrix singleCellDesignMatrix ) {
+        return new REXPNull();
+    }
+
+    private SingleCellExpressionDataMatrix<?> fromDataFrame( REXP rexp ) {
+        return null;
     }
 }

gemma-core/src/main/java/ubic/gemma/core/analysis/singleCell/batcheffect/SingleCellBatchCorrectionServiceImpl.java

@@ -0,0 +1,65 @@
+package ubic.gemma.core.analysis.singleCell.batcheffect;
+
+import org.springframework.beans.factory.annotation.Autowired;
+import org.springframework.stereotype.Service;
+import org.springframework.transaction.annotation.Transactional;
+import org.springframework.util.Assert;
+import ubic.gemma.core.datastructure.matrix.SingleCellDesignMatrix;
+import ubic.gemma.core.datastructure.matrix.SingleCellExpressionDataMatrix;
+import ubic.gemma.core.datastructure.matrix.SingleCellExpressionDataMatrixUtils;
+import ubic.gemma.core.util.r.REngineFactory;
+import ubic.gemma.model.common.quantitationtype.QuantitationType;
+import ubic.gemma.model.expression.bioAssayData.CellLevelCharacteristics;
+import ubic.gemma.model.expression.bioAssayData.SingleCellDimension;
+import ubic.gemma.model.expression.bioAssayData.SingleCellExpressionDataVector;
+import ubic.gemma.model.expression.experiment.ExpressionExperiment;
+import ubic.gemma.persistence.service.expression.experiment.SingleCellExpressionExperimentService;
+
+import java.util.ArrayList;
+import java.util.Collection;
+import java.util.List;
+
+@Service
+public class SingleCellBatchCorrectionServiceImpl implements SingleCellBatchCorrectionService {
+
+    @Autowired
+    private SingleCellExpressionExperimentService singleCellExpressionExperimentService;
+
+    @Autowired
+    private REngineFactory rEngineFactory;
+
+    @Override
+    @Transactional
+    public QuantitationType batchCorrect( ExpressionExperiment ee, QuantitationType qt, SingleCellBatchCorrectionMethod method ) {
+        Assert.notNull( ee.getExperimentalDesign(), ee + " does not have experimental design. It is required to perform batch correction." );
+        BatchCorrection m = createBatchCorrection( method );
+        SingleCellDimension dimension = singleCellExpressionExperimentService.getSingleCellDimension( ee, qt );
+        if ( dimension == null ) {
+            throw new IllegalArgumentException( qt + " does not have single cell dimension." );
+        }
+        List<SingleCellExpressionDataVector> vectors = new ArrayList<>( singleCellExpressionExperimentService.getSingleCellDataVectors( ee, qt ) );
+        SingleCellExpressionDataMatrix<?> dataMatrix = SingleCellExpressionDataMatrix.getMatrix( vectors );
+        Collection<CellLevelCharacteristics> clcs = new ArrayList<>();
+        // TODO: select relevant CTAs and CLCs
+        clcs.addAll( dimension.getCellTypeAssignments() );
+        clcs.addAll( dimension.getCellLevelCharacteristics() );
+        SingleCellDesignMatrix designMatrix = SingleCellDesignMatrix.from( dimension, ee.getExperimentalDesign(), clcs );
+        SingleCellExpressionDataMatrix<?> correctedMatrix = m.perform( dataMatrix, designMatrix );
+        QuantitationType correctedQt = correctedMatrix.getQuantitationType();
+        List<SingleCellExpressionDataVector> correctedVectors = SingleCellExpressionDataMatrixUtils.toVectors( correctedMatrix );
+        String details = "Batch correction using " + method + " for " + ee.getShortName() + " on quantitation type " + qt.getName();
+        singleCellExpressionExperimentService.addSingleCellDataVectors( ee, correctedQt, correctedVectors, details );
+        return correctedQt;
+    }
+
+    private BatchCorrection createBatchCorrection( SingleCellBatchCorrectionMethod method ) {
+        switch ( method ) {
+            case HARMONY:
+                return new Harmony( rEngineFactory );
+            case COMBAT:
+                return new ComBat();
+            default:
+                throw new IllegalArgumentException( "Unknown batch correction method: " + method );
+        }
+    }
+}

gemma-core/src/main/java/ubic/gemma/core/datastructure/matrix/SingleCellDesignMatrixImpl.java

@@ -1,5 +1,6 @@
 package ubic.gemma.core.datastructure.matrix;
 
+import ubic.gemma.core.util.ListUtils;
 import ubic.gemma.model.common.description.Characteristic;
 import ubic.gemma.model.expression.bioAssay.BioAssay;
 import ubic.gemma.model.expression.bioAssayData.CellLevelCharacteristics;
@@ -11,20 +12,33 @@
 import ubic.gemma.model.util.SparseRangeArrayList;
 
 import javax.annotation.Nullable;
-import java.util.ArrayList;
-import java.util.Collections;
-import java.util.List;
+import java.util.*;
 
 public class SingleCellDesignMatrixImpl implements SingleCellDesignMatrix {
 
+    // rows
     private final SparseRangeArrayList<BioAssay> assays;
     private final List<String> cellIds;
+    private final Map<BioAssay, Map<String, Integer>> index;
+
+    // columns
     private final List<ExperimentalFactor> factors;
+    private final Map<ExperimentalFactor, Integer> factorsIndex;
+
+    /**
+     * This is technically a matrix, but using {@link List} allows for sparse range array to be used for sample-level
+     * factors.
+     * <p>
+     * Also, this is transposed w.r.t. to rows/columns that the interface requires. This is due to the fact that
+     * sparsity is better handled along factors
+     */
+    private final List<List<FactorValue>> factorValues;
 
     public SingleCellDesignMatrixImpl( SingleCellDimension dimension, List<BioAssay> assays, List<ExperimentalFactor> factors, List<CellLevelCharacteristics> cellLevelCharacteristics ) {
         int[] bioAssayOffsets = new int[assays.size()];
         int k = 0;
         List<String> cellIdsL = new ArrayList<>( dimension.getNumberOfCells() );
+        Map<BioAssay, Map<String, Integer>> index = new HashMap<>( assays.size() );
         for ( int i = 0; i < assays.size(); i++ ) {
             BioAssay assay = assays.get( i );
             int sampleIndex = dimension.getBioAssays().indexOf( assay );
@@ -34,17 +48,26 @@ public SingleCellDesignMatrixImpl( SingleCellDimension dimension, List<BioAssay>
             List<String> sampleCellIds = dimension.getCellIdsBySample( sampleIndex );
             bioAssayOffsets[i] = k;
             cellIdsL.addAll( sampleCellIds );
+            Map<String, Integer> cellid2pos = new HashMap<>();
+            for ( int j = 0; j < sampleCellIds.size(); j++ ) {
+                cellid2pos.put( sampleCellIds.get( j ), k + j );
+            }
+            index.put( assay, cellid2pos );
             k += sampleCellIds.size();
         }
         this.assays = new SparseRangeArrayList<>( assays, bioAssayOffsets, k );
         this.cellIds = cellIdsL;
+        this.index = index;
         ArrayList<ExperimentalFactor> factorsL = new ArrayList<>( factors.size() + cellLevelCharacteristics.size() );
         factorsL.addAll( factors );
         for ( CellLevelCharacteristics clc : cellLevelCharacteristics ) {
             ExperimentalFactor factor = createFactorFromCellLevelCharacteristics( clc );
             factorsL.add( factor );
         }
         this.factors = Collections.unmodifiableList( factorsL );
+        this.factorsIndex = Collections.unmodifiableMap( ListUtils.indexOfElements( factorsL ) );
+        // TODO: fill the matrix
+        this.factorValues = new ArrayList<>( factors.size() );
     }
 
     @Override
@@ -69,13 +92,17 @@ public int columns() {
 
     @Override
     public List<FactorValue> getColumn( int column ) {
-        return new FactorValue[0];
+        return factorValues.get( column );
     }
 
     @Nullable
     @Override
     public List<FactorValue> getColumn( ExperimentalFactor factor ) {
-        return new FactorValue[0];
+        int index = factors.indexOf( factor );
+        if ( index == -1 ) {
+            return null;
+        }
+        return getColumn( index );
     }
 
     @Override
@@ -91,22 +118,30 @@ public ExperimentalFactor getFactorForColumn( int column ) {
     @Nullable
     @Override
     public List<FactorValue> getRow( BioAssay bioAssay, String cellId ) {
-        return new FactorValue[0];
+        int row = getRowIndex( bioAssay, cellId );
+        if ( row == -1 ) {
+            return null;
+        }
+        return getRow( row );
     }
 
     @Override
     public List<FactorValue> getRow( int row ) {
-        return new FactorValue[0];
+        List<FactorValue> fvs = new ArrayList<>( factors.size() );
+        for ( int i = 0; i < factors.size(); i++ ) {
+            fvs.add( factorValues.get( i ).get( row ) );
+        }
+        return fvs;
     }
 
     @Override
     public int rows() {
-        return 0;
+        return cellIds.size();
     }
 
     @Override
     public BioAssay getBioAssayForRow( int row ) {
-        return null;
+        return assays.get( row );
     }
 
     @Override
@@ -116,12 +151,16 @@ public BioMaterial getBioMaterialForRow( int row ) {
 
     @Override
     public String getCellIdForRow( int row ) {
-        return "";
+        return cellIds.get( row );
     }
 
     @Override
     public int getRowIndex( BioAssay bioAssay, String cellId ) {
-        return 0;
+        Map<String, Integer> cell2pos = index.get( bioAssay );
+        if ( cell2pos == null ) {
+            return -1;
+        }
+        return cell2pos.getOrDefault( cellId, -1 );
     }
 
     private ExperimentalFactor createFactorFromCellLevelCharacteristics( CellLevelCharacteristics characteristics ) {

gemma-core/src/main/java/ubic/gemma/core/datastructure/matrix/io/SingleCellMetadataWriter.java

@@ -0,0 +1,12 @@
+package ubic.gemma.core.datastructure.matrix.io;
+
+import ubic.gemma.model.expression.bioAssayData.SingleCellDimension;
+import ubic.gemma.model.expression.experiment.ExpressionExperiment;
+
+import java.io.IOException;
+import java.io.Writer;
+
+public interface SingleCellMetadataWriter {
+
+    void write( ExpressionExperiment ee, SingleCellDimension singleCellDimension, Writer writer ) throws IOException;
+}