Changeset 1912 for molgenis4xgap_importwizard

Timestamp:

04/07/09 15:40:03 (3 years ago)

Author:

jvelde

Message:

 

Location:

molgenis4xgap_importwizard/handwritten/java/cistransanalysis3

Files:

• Run.java (modified) (2 diffs)
• Step3456.java (modified) (6 diffs)
• Step78.java (modified) (8 diffs)
• makeUniqueSortAndPrint.java (modified) (1 diff)

molgenis4xgap_importwizard/handwritten/java/cistransanalysis3/Run.java

@@ -52 +52 @@
                 for(int index=0; index<intervalGenes.size();index++){
                         Marker m = Step3456.getMarkerInfo(db, intervalGenes.get(index).getEr().getMarkerName());
-                        pw.println("Interval of " + intervalGenes.get(index).getEr().getTraitname() + " on chr " + m.getChr() + " has " + intervalGenes.get(index).getGenesFromInterval().size() + " genes (cM = " + intervalGenes.get(index).getEr().getCiLower() + " to "
-                                        + intervalGenes.get(index).getEr().getCiUpper() + ")");
+                        pw.println("Interval of " + intervalGenes.get(index).getEr().getTraitname() + " on chr " + m.getChr() + " has " + intervalGenes.get(index).getGenesFromInterval().size() + " genes.");
 
                 }
@@ -119 +118 @@
                 String[] pathways = new String[]{"carotenoids","flavonols", "glucosinolate", "phenylpropanoid", "tocopherol", "folate"};
                 
-                //String[] pathways = new String[]{"phenylpropanoid"};
+                //String[] pathways = new String[]{"glucosinolate"};
                 
                 PrintWriter pw = new PrintWriter(new BufferedWriter(new FileWriter("D:/data/xgapdata/Frank_for_dbGG/cistransanalysis/results/output.txt")));

molgenis4xgap_importwizard/handwritten/java/cistransanalysis3/Step3456.java

@@ -35 +35 @@
         {
                 List<EQTLResultMap> intervalGenes = new ArrayList<EQTLResultMap>();
-                
-                System.out.println("getting all markers");
+
+                // System.out.println("getting all markers");
                 List<Marker> allMarkers = db.find(Marker.class);
 
                 for (IntervalPeakResult er : reducedList2)
                 {
-                        System.out.println("getting marker/gene info");
+                        // System.out.println("getting marker/gene info");
                         Marker m = getMarkerInfo(db, er.getMarkerName());
                         Gene g = getGeneInfo(db, er.getTraitname());
-                        System.out.println("\tdone");
+                        // System.out.println("\tdone");
 
                         // Step 3 & 4 - get confidence intervals, and select closest markers
-                        System.out.println("getting closest markers");
+                        // System.out.println("getting closest markers");
                         Marker ciLowerClosestMarker = findMarkerClosestToCM(allMarkers, er.getCiLower(), Integer.parseInt(m
                                         .getChr()));
                         Marker ciUpperClosestMarker = findMarkerClosestToCM(allMarkers, er.getCiUpper(), Integer.parseInt(m
                                         .getChr()));
-                        System.out.println("\tdone");
+                        // System.out.println("\tdone");
+                        
+                        if(ciLowerClosestMarker.getName().equals(ciUpperClosestMarker.getName())){
+                                String warning = "WARNING: the same markers ("+ciLowerClosestMarker.getName()+") mapped for interval positions";
+                                System.out.println(warning);
+                                pw.println(warning);
+                        }
 
                         // Step 5 - get the bp positions of thse markers
@@ -61 +67 @@
                         List<Gene> genesFromInterval = findGenesFromInterval(db, Integer.parseInt(m.getChr()), lowerBp, upperBp);
 
-                        boolean isCis = checkIfCis(er, m, g);
+                        boolean isCis = checkIfCis(m, g, ciLowerClosestMarker.getBpStart(), ciUpperClosestMarker.getBpStart());
 
                         if (!isCis)
@@ -67 +73 @@
                                 intervalGenes.add(new EQTLResultMap(er, genesFromInterval, isCis, (((double) lowerBp) / 1000000.0),
                                                 (((double) upperBp) / 1000000.0)));
-                                pw.println((isCis ? "Cis" : "Trans") + " QTL of " + er.getTraitname() + " at marker "
-                                                + er.getMarkerName() + " on chr " + g.getChr() + " has interval: " + lowerBp + "bp to "
-                                                + upperBp + "bp on chr " + Integer.parseInt(m.getChr()) + ". Closest markers: "
-                                                + ciLowerClosestMarker.getName() + " to " + ciUpperClosestMarker.getName());
-
-                                System.out.println((isCis ? "Cis" : "Trans") + " QTL of " + er.getTraitname() + " at marker "
-                                                + er.getMarkerName() + " on chr " + g.getChr() + " has interval: " + lowerBp + "bp to "
-                                                + upperBp + "bp on chr " + Integer.parseInt(m.getChr()) + ". Closest markers: "
-                                                + ciLowerClosestMarker.getName() + " to " + ciUpperClosestMarker.getName());
-                                
+
+                                String info = (isCis ? "Cis" : "Trans") + " QTL of " + er.getTraitname() + " (~" + meanGeneBpPos(g)
+                                                + "bp, chr " + g.getChr() + ") at marker " + er.getMarkerName() + " (" + m.getCM() + "cM, chr "
+                                                + m.getChr() + ")" + " has interval mapped to true markers: " + ciLowerClosestMarker.getName()
+                                                + " (" + ciLowerClosestMarker.getCM() + "cM, " + ciLowerClosestMarker.getBpStart() + "bp) and "
+                                                + ciUpperClosestMarker.getName() + " (" + ciUpperClosestMarker.getCM() + "cM, "
+                                                + ciUpperClosestMarker.getBpStart() + "bp).";
+
+                                pw.println(info);
+                                System.out.println(info);
+
                         }
 
@@ -94 +101 @@
         }
 
-        private boolean checkIfCis(IntervalPeakResult er, Marker m, Gene g)
-        {
-                if (g.getChr() == m.getChr())
-                {
-                        if (m.getCM() > er.getCiLower() && m.getCM() < er.getCiUpper())
+        private boolean checkIfCis(Marker m, Gene g, long ciLowerBp, long ciUpperBp)
+        {
+                if (g.getChr().equals(m.getChr()))
+                {
+                        long geneBpPos = meanGeneBpPos(g);
+
+                        if (geneBpPos > ciLowerBp && geneBpPos < ciUpperBp)
                         {
                                 return true;
@@ -104 +113 @@
                 }
                 return false;
+        }
+
+        static long meanGeneBpPos(Gene g)
+        {
+                long geneBpPos;
+                if (g.getBpStart() < g.getBpEnd())
+                {
+                        geneBpPos = g.getBpStart() + ((long) g.getBpEnd() - (long) g.getBpStart()) / (long) 2.0;
+                }
+                else if (g.getBpStart() > g.getBpEnd())
+                {
+                        geneBpPos = g.getBpEnd() + ((long) g.getBpStart() - (long) g.getBpEnd()) / (long) 2.0;
+                }
+                else
+                {
+                        geneBpPos = g.getBpStart();
+                }
+                return geneBpPos;
+        }
+        
+        static long meanGeneBpPos(BiosynthGene g)
+        {
+                long geneBpPos;
+                if (g.getBpStart() < g.getBpEnd())
+                {
+                        geneBpPos = g.getBpStart() + ((long) g.getBpEnd() - (long) g.getBpStart()) / (long) 2.0;
+                }
+                else if (g.getBpStart() > g.getBpEnd())
+                {
+                        geneBpPos = g.getBpEnd() + ((long) g.getBpStart() - (long) g.getBpEnd()) / (long) 2.0;
+                }
+                else
+                {
+                        geneBpPos = g.getBpStart();
+                }
+                return geneBpPos;
         }
 
@@ -128 +173 @@
         private Marker findMarkerClosestToCM(List<Marker> allMarkers, double cM, int chr)
         {
-                System.out.println("finding marker...");
                 double smallestDiff = -1;
                 Marker closestMarker = null;

molgenis4xgap_importwizard/handwritten/java/cistransanalysis3/Step78.java

@@ -17 +17 @@
 
 import org.apache.commons.math.MathException;
+import org.apache.commons.math.stat.descriptive.moment.Mean;
 import org.apache.commons.math.stat.regression.SimpleRegression;
 import org.molgenis.assets.data.DatabaseException;
@@ -44 +45 @@
                                 + " for Pearson's product moment correlation coefficient\n");
 
+                List<Marker> allMarkers = db.find(Marker.class);
+                
                 Map<IntervalPeakResult, List<Correlation>> correlatedGenes = new HashMap<IntervalPeakResult, List<Correlation>>();
 
@@ -57 +60 @@
                         if (intervalGenes.get(index).getGenesFromInterval().size() > 0)
                         {
-                                
+
                                 Marker mt = Step3456.getMarkerInfo(db, intervalGenes.get(index).getEr().getMarkerName());
                                 Gene gt = Step3456.getGeneInfo(db, intervalGenes.get(index).getEr().getTraitname());
 
-                                pw.println("Investigating all " + intervalGenes.get(index).getGenesFromInterval().size()
-                                                + " genes within the QTL confidence interval at chr "
-                                                + mt.getChr() + " of gene "
-                                                + intervalGenes.get(index).getEr().getTraitname() + ", located at chr "
-                                                + gt.getChr() + "");
-                                
-                                System.out.println("Investigating all " + intervalGenes.get(index).getGenesFromInterval().size()
-                                                + " genes within the QTL confidence interval at chr "
-                                                + mt.getChr() + " of gene "
-                                                + intervalGenes.get(index).getEr().getTraitname() + ", located at chr "
-                                                + gt.getChr() + "");
+                                String info = "Investigating all " + intervalGenes.get(index).getGenesFromInterval().size()
+                                                + " genes within the QTL confidence interval at chr " + mt.getChr() + " of gene "
+                                                + intervalGenes.get(index).getEr().getTraitname() + ", located at chr " + gt.getChr();
+
+                                pw.println(info);
+                                System.out.println(info);
 
                                 for (Gene g : intervalGenes.get(index).getGenesFromInterval())
                                 {
-
-                                        int GeneID = db.find(Gene.class, new QueryRule("name", Operator.EQUALS, g.getName())).get(0)
-                                                        .getId();
-                                        QueryRule row = new QueryRule("row", Operator.EQUALS, GeneID);
-                                        List<DecimalDataElement> genePVals = db.find(DecimalDataElement.class, data, row);
-
                                         for (BiosynthGene bg : biosynthGenes)
                                         {
+                                                // onderzochte trait moet een pathway gene zijn
                                                 if (bg.getAgiId().equals(intervalGenes.get(index).getEr().getTraitname()))
                                                 {
-
+                                                        // interval gene mag niet dezelfde zijn als dit gene
+                                                        // (correlation 1.0)
                                                         if (!g.getName().equals(bg.getAgiId()))
                                                         {
+                                                                int GeneID = db.find(Gene.class, new QueryRule("name", Operator.EQUALS, g.getName()))
+                                                                                .get(0).getId();
+                                                                QueryRule row = new QueryRule("row", Operator.EQUALS, GeneID);
+                                                                List<DecimalDataElement> genePVals = db.find(DecimalDataElement.class, data, row);
+
                                                                 double corr = correlate(genePVals, biosynthExpressionMap.get(bg.getAgiId()));
+
+                                                                double mean1 = mean(genePVals);
+                                                                double mean2 = mean(biosynthExpressionMap.get(bg.getAgiId()));
                                                                 
-                                                                
-                                                                
-//                                                              if (corr > thresholdCorrStep8 || -corr > thresholdCorrStep8)
-//                                                              {
-
-                                                                        boolean isCisCorrelation = isCisCorrelation(bg, intervalGenes.get(index), mt);
-                                                                        Correlation corrObj = new Correlation(g, bg, corr, isCisCorrelation);
-                                                                        correlatedGenesList.add(corrObj);
-//                                                              }
-                                                                
-                                                                
+                                                                // if (corr > thresholdCorrStep8 || -corr >
+                                                                // thresholdCorrStep8)
+                                                                // {
+
+                                                                boolean isCisCorrelation = isCisCorrelation(bg, intervalGenes.get(index), mt);
+                                                                Correlation corrObj = new Correlation(g, bg, corr, isCisCorrelation, mean1, mean2);
+                                                                correlatedGenesList.add(corrObj);
+                                                                // }
+
                                                         }
                                                 }
@@ -121 +121 @@
                         else
                         {
-//                              pw.println("\tA total of " + correlatedGenesList.size() + " genes correlated with known " + pathway
-//                                              + " genes located elsewhere.");
-                        }
-                        
+                                // pw.println("\tA total of " + correlatedGenesList.size() +
+                                // " genes correlated with known " + pathway
+                                // + " genes located elsewhere.");
+                        }
+
                         Collections.sort(correlatedGenesList, new SortCorrelationList());
-                        correlatedGenesList = correlatedGenesList.subList(0, correlatedGenesList.size() > 10 ? 10 : correlatedGenesList.size());
-                        
-                        for(Correlation c : correlatedGenesList){       
-                                if (c.isCis())
+                        correlatedGenesList = correlatedGenesList.subList(0, correlatedGenesList.size() > 20 ? 20
+                                        : correlatedGenesList.size());
+
+                        for (Correlation c : correlatedGenesList)
+                        {
+//                              if (c.isCis())
+//                              {
+//                                      pw.print("\tCorr. (near): ");
+//                              }
+//                              else
+//                              {
+//                                      pw.print("\tCorr. (far): ");
+//                              }
+                                
+                                Marker m1 = findMarkerClosestToBp(allMarkers, (int)Step3456.meanGeneBpPos(c.getGene()), Integer.parseInt(c.getGene().getChr()));
+                                Marker m2 = findMarkerClosestToBp(allMarkers, (int)Step3456.meanGeneBpPos(c.getBiosynthGene()), c.getBiosynthGene().getChromosome());
+                                
+                                double lod1 = getLOD(c.getGene().getName(), m1, db);
+                                double lod2 = getLOD(c.getBiosynthGene().getAgiId(), m2, db);
+                                
+                                double avgExpr1 = c.getAvgOfGene();
+                                double avgExpr2 = c.getAvgOfBiosynthGene();
+                                
+                                String info = c.getGene().getName() + " (chr " + c.getGene().getChr() + ", "
+                                                + Step3456.meanGeneBpPos(c.getGene()) + "bp, LOD at closest marker = "+lod1+", avg. expr. = "+avgExpr1+") correlates with " + c.getBiosynthGene().getAgiId()
+                                                + " (chr " + c.getBiosynthGene().getChromosome() + ", ~"
+                                                + Step3456.meanGeneBpPos(c.getBiosynthGene()) + "bp, LOD at closest marker = "+lod2+", avg. expr. = "+avgExpr2+"), strength = " + c.getCorrelation();
+
+                                pw.println(info);
+                                System.out.println(info);
+                        }
+
+                        correlatedGenes.put(intervalGenes.get(index).getEr(), correlatedGenesList);
+                }
+                this.setCorrelatedGenes(correlatedGenes);
+        }
+        
+        private double mean(List<DecimalDataElement> genePVals)
+        {
+                double total = 0;
+                for(DecimalDataElement d : genePVals){
+                        total += d.getValue();
+                }
+                return total / ((double)genePVals.size());
+        }
+
+        private double getLOD(String geneName, Marker m, JDBCDatabase db) throws DatabaseException{
+                int dataID = db.find(Data.class, new QueryRule("name", Operator.EQUALS, "Gene_lodscores")).get(0).getId();
+                int geneID = db.find(Gene.class, new QueryRule("name", Operator.EQUALS, geneName)).get(0).getId();
+                
+                QueryRule data = new QueryRule("data", Operator.EQUALS, dataID);
+                QueryRule row = new QueryRule("row", Operator.EQUALS, geneID);
+                QueryRule col = new QueryRule("col", Operator.EQUALS, m.getId());
+                
+                return db.find(DecimalDataElement.class, data, row, col).get(0).getValue();
+        }
+        
+        private Marker findMarkerClosestToBp(List<Marker> allMarkers, int bp, int chr)
+        {
+                double smallestDiff = -1;
+                Marker closestMarker = null;
+                for (Marker tryMarker : allMarkers)
+                {
+                        if (tryMarker.getBpStart() != (long) 0)
+                        {
+                                if (Integer.parseInt(tryMarker.getChr()) == chr)
                                 {
-                                        pw.print("\tCorr. (cis): ");
+                                        double diff = Math.abs(tryMarker.getBpStart() - bp);
+
+                                        if (diff < smallestDiff && smallestDiff != -1)
+                                        {
+                                                smallestDiff = diff;
+                                                closestMarker = tryMarker;
+                                        }
+                                        else if (smallestDiff == -1)
+                                        {
+                                                smallestDiff = diff;
+                                                closestMarker = tryMarker;
+                                        }
                                 }
-                                else
-                                {
-                                        pw.print("\tCorr. (trans): ");
-                                }
-                                pw.println(c.getGene().getName() + " (chr " + c.getGene().getChr()
-                                                + ") correlates with " + c.getBiosynthGene().getAgiId() + ", strength = " + c.getCorrelation());
-                                
-                                System.out.println(c.getGene().getName() + " (chr " + c.getGene().getChr()
-                                                + ") correlates with " + c.getBiosynthGene().getAgiId() + ", strength = " + c.getCorrelation());
-                        }
-                        
-                        
-                        
-                        correlatedGenes.put(intervalGenes.get(index).getEr(), correlatedGenesList);
-                }
-                this.setCorrelatedGenes(correlatedGenes);
+                        }
+                }
+                return closestMarker;
         }
 
         private boolean isCisCorrelation(BiosynthGene bg, EQTLResultMap resultMap, Marker m)
         {
-                if (bg.getChromosome() == (int)Integer.parseInt(m.getChr()))
-                {
+                // System.out.println("*** isCisCorrelation?");
+                // System.out.println("\tbg.getChromosome() = " + bg.getChromosome() +
+                // ", m.getChr() = " + m.getChr());
+                if (bg.getChromosome() == (int) Integer.parseInt(m.getChr()))
+                {
+                        // System.out.println("\tYES");
+                        // System.out.println("\tbg.getMbPos() = " + bg.getMbPos()+
+                        // ", IntervalMbStart() = " + resultMap.getIntervalMbStart() +
+                        // ", IntervalMbEnd() = " + resultMap.getIntervalMbEnd());
+                        // System.out.println("\tbg > start && bg < end?");
                         if (bg.getMbPos() > resultMap.getIntervalMbStart() && bg.getMbPos() < resultMap.getIntervalMbEnd())
                         {
+                                // System.out.println("\tYES");
                                 return true;
                         }
@@ -181 +251 @@
         }
 
-        private String printGene(Gene g)
-        {
-                String print = "";
-                print += "chr = " + g.getChr() + ", ";
-                print += "bpstr = " + g.getBpStart() + ", ";
-                print += "bpend = " + g.getBpEnd();
-                return print;
-        }
-
-        private String printBiosynthGene(BiosynthGene bg)
-        {
-                String print = "";
-                print += "abbr = " + bg.getGeneAbbr() + ", ";
-                print += "prot = " + bg.getProtein() + ", ";
-                print += "pathw = " + bg.getPathway() + ", ";
-                print += "orien = " + bg.getOrientation() + ", ";
-                print += "chr = " + bg.getChromosome() + ", ";
-                print += "bpstr = " + bg.getBpStart() + ", ";
-                print += "bpend = " + bg.getBpEnd() + ", ";
-                print += "mbpos = " + bg.getMbPos();
-                return print;
-        }
-
-        private BiosynthGene getBiosynthGene(Gene gene, List<BiosynthGene> biosynthGenes)
-        {
-                for (BiosynthGene bg : biosynthGenes)
-                {
-                        if (bg.getAgiId().equals(gene.getName()))
-                        {
-                                return bg;
-                        }
-                }
-                return null;
-        }
-
         private double correlate(List<DecimalDataElement> genePVals, List<DecimalDataElement> corrGenePVals)
                         throws MathException
@@ -244 +279 @@
 }
 
-
-class SortCorrelationList implements Comparator<Correlation>{
-    public int compare(Correlation c1, Correlation c2) {
-        return ((Double)(Math.abs(c2.getCorrelation()))).compareTo((Double)(Math.abs(c1.getCorrelation())));
-    }
+class SortCorrelationList implements Comparator<Correlation>
+{
+        public int compare(Correlation c1, Correlation c2)
+        {
+                return ((Double) (Math.abs(c2.getCorrelation()))).compareTo((Double) (Math.abs(c1.getCorrelation())));
+        }
 }
 
 class Correlation
 {
-        public Correlation(Gene gene, BiosynthGene biosynthGene, double correlation, boolean isCis)
+        public Correlation(Gene gene, BiosynthGene biosynthGene, double correlation, boolean isCis, double avgOfGene, double avgOfBiosynthGene)
         {
                 this.gene = gene;
@@ -259 +295 @@
                 this.correlation = correlation;
                 this.isCis = isCis;
+                this.avgOfGene = avgOfGene;
+                this.avgOfBiosynthGene = avgOfBiosynthGene;
         }
 
@@ -266 +304 @@
         boolean isCis;
         
-        public boolean isCis(){
+        double avgOfGene;
+        double avgOfBiosynthGene;
+        
+        
+
+        public double getAvgOfGene()
+        {
+                return avgOfGene;
+        }
+
+        public double getAvgOfBiosynthGene()
+        {
+                return avgOfBiosynthGene;
+        }
+
+        public boolean isCis()
+        {
                 return isCis;
         }

molgenis4xgap_importwizard/handwritten/java/cistransanalysis3/makeUniqueSortAndPrint.java

@@ -24 +24 @@
                 
                 Collections.sort(uniqCorr, new SortSimpleCorrelation());
-                pw.println("UnknownGene_name\tUnknownGene_chr\tUnknownGene_isPathwayGene\tPathwayGene_name\tPathwayGene_chr\tCorrelationStrength");
+                pw.println("UnknownGene_name\tUnknownGene_chr\tPathwayGene_name\tPathwayGene_chr\tCorrelationStrength");
                 for (SimpleCorrelation sc : uniqCorr) {
-                        pw.println(sc.getGene1() + "\t" + sc.getChrOfGene1() + "\t" + sc.isPathwayGene1() + "\t" + sc.getGene2()  + "\t" + sc.getChrOfGene2() + "\t" + sc.getCorr());
+                        pw.println(sc.getGene1() + "\t" + sc.getChrOfGene1() + "\t" + sc.getGene2()  + "\t" + sc.getChrOfGene2() + "\t" + sc.getCorr());
                 }
                 
                 pw.println("___________\n");
-        }
-        
-        private String cisOrTrans(int chrOfGene1, int chrOfGene2)
-        {
-                if(chrOfGene1==chrOfGene2){
-                        return "cis";
-                }else{
-                        return "trans";
-                }
         }