simulate PPKTs

Author: pnrobinson

pom.xml

@@ -4,7 +4,7 @@
     <modelVersion>4.0.0</modelVersion>
     <groupId>org.monarchinitiative.hpotools</groupId>
     <artifactId>hpotools</artifactId>
-    <version>0.1.3</version>
+    <version>0.1.4</version>
     <packaging>jar</packaging>
     <name>hpotools</name>
 
@@ -167,7 +167,11 @@
             <artifactId>commons-math3</artifactId>
             <version>3.6.1</version>
         </dependency>
-
+        <dependency>
+            <groupId>com.google.protobuf</groupId>
+            <artifactId>protobuf-java-util</artifactId>
+            <version>3.21.12</version> <!-- or compatible version -->
+        </dependency>
     </dependencies>
 
     <build>

src/main/java/org/monarchinitiative/hpotools/Main.java

@@ -24,9 +24,8 @@ public static void main(String[] args) {
                     .addSubcommand("maxo", new MaxoCommand())
                     .addSubcommand("mondo", new MondoCommand())
                     .addSubcommand("onset", new OnsetCommand())
-                    .addSubcommand("simhpo", new SimHpoCommand())
+                    .addSubcommand("sim", new SimHpoCommand())
                     .addSubcommand("stats", new StatsCommand())
-                    .addSubcommand("translate", new DiseaseTranslateCommand())
                     .addSubcommand("tsv", new Hpo2TsvCommand())
                     .addSubcommand("word", new WordCommand())
                    ;

src/main/java/org/monarchinitiative/hpotools/analysis/simhpo/SimulatedHpoDiseaseGenerator.java

@@ -68,13 +68,13 @@ public Optional<Phenopacket > generateSimulatedPhenopacket(TermId omimId, int nT
                 // choose a random onset from the range
                 int start = onsetRange.start().days();
                 int end = onsetRange.end().days();
-                onset = random.nextInt(start, end + 1);
+                onset = (int)(start+end)/2;
             } else {
                 LOGGER.debug("No onset information available for disease {}", omimId.getValue());
             }
 
             // Add some age to the phenopacket by adding a few years to the onset
-            if (onset != 0) {
+            if (onset > 0) {
                 age = onset + random.nextInt(0, 10 * 365); // add up to 10 years
             }
 

src/main/java/org/monarchinitiative/hpotools/cmd/DiseaseTranslateCommand.java

@@ -12,10 +12,9 @@
 
 
 /**
- * Implementation of download in HpoWorkbench. The command is intended to download
- * both the OBO file and the association file. For HPO, this is {@code hp.obo} and
- * {@code phenotype_annotation.tab}.
- * Code modified from Download command in Jannovar.
+ * Implementation of download in HpoTools. The command is intended to download
+ * both the JSON file and the association file. For HPO, this is {@code hp.json} and
+ * {@code phenotype.hpoa}.
  * @author <a href="mailto:manuel.holtgrewe@charite.de">Manuel Holtgrewe</a>
  * @author <a href="mailto:peter.robinson@jax.org">Peter Robinson</a>
  * @version 0.0.1 (May 10, 2017)

src/main/java/org/monarchinitiative/hpotools/cmd/DownloadCommand.java

@@ -13,6 +13,9 @@
 import java.util.concurrent.Callable;
 
 
+/**
+ * Create TSV file with a target term and all of its descendants.
+ */
 @CommandLine.Command(name = "word",
         mixinStandardHelpOptions = true,
         description = "Output subontology as word file (experimental)")
@@ -35,14 +38,10 @@ public Integer call() {
         TermId hpoId = TermId.of(startTermId);
         Optional<Term> opt = hpOntology.termForTermId(hpoId);
         if (opt.isEmpty()) {
-            System.err.printf("[ERROR] No HPO term found for %s.\n", startTermId);
-        }
-        if (opt.isEmpty()) {
-            LOGGER.error("[ERROR] No term found for {}.", startTermId);
+            LOGGER.error("[ERROR] No HPO term found for {}.", startTermId);
             return 1;
         }
         Term targetTerm = opt.get();
-
         String name = targetTerm.getName().replaceAll(" ", "_");
         String id = targetTerm.id().getValue().replaceAll(":", "_");
         String outfilename = String.format("%s_%s.tsv", name, id);

src/main/java/org/monarchinitiative/hpotools/cmd/Hpo2TsvCommand.java

@@ -13,7 +13,10 @@
 import java.util.*;
 import java.util.concurrent.Callable;
 
-
+/**
+ * Test the distribution of terms in the various top-level HPO categories. Are the categories with
+ * significantly more terms than we would expect by chance?
+ */
 @CommandLine.Command(name = "dist",
         mixinStandardHelpOptions = true,
         description = "Output subontology as word file (experimental)")
@@ -33,12 +36,6 @@ public Integer call() {
         double probability = ExactMultinomial.exactMultinomialTest(orderedObservedCounts, expectedProportions);
         System.out.println("Exact multinomial p value " + probability);
         outputProprotions(topLevelHpoTermCounts, orderedObservedCounts, expectedProportions, topLevelHpoLabels);
-        /*
-        ChiSquareTest chiSquareTest = new ChiSquareTest();
-        double pValue = chiSquareTest.chiSquareTest(expected, observedCounts);
-        double chiSquareStatistic = chiSquareTest.chiSquare(expected, observedCounts);
-    */
-
         return 0;
     }
 
@@ -74,7 +71,7 @@ private void outputProprotions(List<Map.Entry<TermId, Integer>> topLevelHpoTermC
      * Create an array with observed counts and arrange it in the same order
      * @param topLevelHpoTermCounts List with top-level HPO terms in order
      * @param observedHpoCounts map with counts observed in an experiment
-     * @return
+     * @return A list of observed counts
      */
     private int[] getOrderedObservedCounts(List<Map.Entry<TermId, Integer>> topLevelHpoTermCounts, Map<TermId, Integer> observedHpoCounts) {
         int[] observed_counts = new int[topLevelHpoTermCounts.size()];
@@ -108,7 +105,7 @@ private double[] getExpectedProportions(List<Map.Entry<TermId, Integer>> topLeve
 
     /**
      * Get a list of the top level terms, ordered by total number of children
-     * @param hpo
+     * @param hpo the Ontology
      * @return A list of Map Entries, with key the TermId and value the count
      */
     private List<Map.Entry<TermId, Integer>> getTopLevelTermList(Ontology hpo) {
@@ -154,21 +151,21 @@ private Map<TermId, Integer> observedHpoCounts() {
         // Abnormality of the musculoskeletal system
         hpoCounts.put(TermId.of("HP:0033127"), 38);
         // Abnormality of the nervous system
-        hpoCounts.put(TermId.of("HP:0000707"), 37);
+        hpoCounts.put(TermId.of("HP:0000707"), 40);
         //Abnormality of limbs
         hpoCounts.put(TermId.of("HP:0040064"), 19);
         //Abnormality of the cardiovascular system
         hpoCounts.put(TermId.of("HP:0001626"), 11);
         //Neoplasm
-        hpoCounts.put(TermId.of("HP:0002664"), 13);
+        hpoCounts.put(TermId.of("HP:0002664"), 16);
         // Abnormality of head or neck
         hpoCounts.put(TermId.of("HP:0000152"), 8);
         // Abnormality of the eye
         hpoCounts.put(TermId.of("HP:0000478"), 6);
         // Abnormality of the integument
-        hpoCounts.put(TermId.of("HP:0001574"), 8);
+        hpoCounts.put(TermId.of("HP:0001574"), 10);
         // Growth abnormality
-        hpoCounts.put(TermId.of("HP:0001507"), 6);
+        hpoCounts.put(TermId.of("HP:0001507"), 5);
         // Abnormality of blood and blood-forming tissues
         hpoCounts.put(TermId.of("HP:0001871"), 6);
         // Abnormality of the immune system

src/main/java/org/monarchinitiative/hpotools/cmd/HpoDistCommand.java

@@ -16,26 +16,35 @@
 import org.slf4j.LoggerFactory;
 import picocli.CommandLine;
 
+import com.google.protobuf.util.JsonFormat;
+import org.phenopackets.schema.v2.Phenopacket;
+
 import java.io.File;
+import java.io.IOException;
+import java.nio.file.Files;
 import java.nio.file.Path;
-import java.util.Optional;
-import java.util.Set;
+import java.nio.file.Paths;
+import java.util.*;
 import java.util.concurrent.Callable;
 
-@CommandLine.Command(name = "onset",
+@CommandLine.Command(name = "sim",
         mixinStandardHelpOptions = true,
-        description = "Calculate number of diseases with onset data")
+        description = "Simulate phenopackets")
 public class SimHpoCommand extends HPOCommand implements Callable<Integer> {
     private final static Logger LOGGER = LoggerFactory.getLogger(SimHpoCommand.class);
 
     /** default Noonan syndrome 1	163950 */
-    @CommandLine.Option(names={"--disease"}, description = "OMIM identifer of disease to be simulated", required = true)
-    private String omimIdentifier = "163950";
+    @CommandLine.Option(names={"-c","--cases"}, description = "Number of cases to be simulated", required = false)
+    private int nCases = 100;
 
     /** default simulate 5 HPO terms */
     @CommandLine.Option(names={"-n","--nterms"}, description = "number of HPO terms to be simulated", required = false)
     private int nterms = 5;
 
+    /** output directory for simulated phenopackets */
+    @CommandLine.Option(names={"--outdir"}, description = "Output directory", required = false)
+    private File outdir = new File("sim_data");
+
     @Override
     public Integer call() throws Exception {
         if (hpopath==null) {
@@ -58,18 +67,48 @@ public Integer call() throws Exception {
         HpoDiseaseLoader loader = HpoDiseaseLoaders.defaultLoader(ontology, options);
         Path annotpath = annotFile.toPath();
         HpoDiseases diseases = loader.load(annotpath);
-
-        SimulatedHpoDiseaseGenerator generator = new SimulatedHpoDiseaseGenerator(diseases, ontology);
-        TermId diseaseId = TermId.of("OMIM", omimIdentifier);
-        Optional<Phenopacket> opt = generator.generateSimulatedPhenopacket(diseaseId);
-        // When we get here, we want to output the simulated phenopackets to file
-        // now, the generator is just a skeleton and it always returns Optional.empty()!!!
-        if (opt.isPresent()) {
-            System.out.println(opt.get());
+        // make directory if needed
+        if (outdir.exists()) {
+            System.out.println("[WARN] Output directory already exists: " + outdir.getAbsolutePath());
         } else {
-            System.out.println("Could not retrieve phenopacket for \"" + omimIdentifier + "\"");
-            System.out.println("NEED TO IMPLEMENT generateSimulatedPhenopacket()");
+            boolean created = outdir.mkdirs();
+            if (!created) {
+                System.err.println("[ERROR] Could not create directory at: " + outdir.getAbsolutePath());
+                return 1;
+            }
         }
+        SimulatedHpoDiseaseGenerator generator = new SimulatedHpoDiseaseGenerator(diseases, ontology);
+        List<TermId> diseaseIds = new ArrayList<>(diseases.diseaseIds());
+        int count = 0;
+        while (count < nCases) {
+            Random rand = new Random();
+            int idx = rand.nextInt(diseaseIds.size());
+            TermId randomOmimId = diseaseIds.get(idx);
+            Optional<Phenopacket> opt = generator.generateSimulatedPhenopacket(randomOmimId);
+            // When we get here, we want to output the simulated phenopackets to file
+            // now, the generator is just a skeleton and it always returns Optional.empty()!!!
+            if (opt.isPresent()) {
+                Phenopacket ppkt = opt.get();
+                String jsonString = JsonFormat.printer().print(ppkt);
+                System.out.println(jsonString);
+                String ppkt_id = ppkt.getId();
+                String cleaned = ppkt_id.replaceAll("[^\\x00-\\x7F]", "_");
+                String outname = String.format("%s.json", cleaned);
+                try {
+                    Path filePath = Paths.get(String.valueOf(outdir), outname);
+                    Files.write(filePath, jsonString.getBytes());
+                    System.out.println("File written successfully");
+                } catch (IOException e) {
+                    System.err.println("Error writing file: " + e.getMessage());
+                }
+            } else {
+                System.out.println("Could not retrieve phenopacket for \"" + count + "\"");
+                System.out.println("NEED TO IMPLEMENT generateSimulatedPhenopacket()");
+            }
+            count++;
+        }
+
+
 
         return 0;
     }