Added a GCCalculatorWalker for a oneoff analysis for Mark Daly (GC content of agilent 1.1 targets)

git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@3395 348d0f76-0448-11de-a6fe-93d51630548a
2010-05-19 18:49:51 +00:00 · 2010-05-19 18:49:51 +00:00 · ffb1b46166
parent 0036df7b03
commit ffb1b46166
1 changed files with 69 additions and 0 deletions
--- a/java/src/org/broadinstitute/sting/oneoffprojects/walkers/GCCalculatorWalker.java
+++ b/java/src/org/broadinstitute/sting/oneoffprojects/walkers/GCCalculatorWalker.java
@ -0,0 +1,69 @@
+package org.broadinstitute.sting.oneoffprojects.walkers;
+
+import org.broadinstitute.sting.gatk.contexts.AlignmentContext;
+import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
+import org.broadinstitute.sting.gatk.refdata.IntervalRod;
+import org.broadinstitute.sting.gatk.refdata.ReadMetaDataTracker;
+import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
+import org.broadinstitute.sting.gatk.refdata.utils.GATKFeature;
+import org.broadinstitute.sting.gatk.walkers.RefWalker;
+import org.broadinstitute.sting.utils.BaseUtils;
+import org.broadinstitute.sting.utils.GenomeLoc;
+import org.broadinstitute.sting.utils.collections.Pair;
+
+import java.util.HashMap;
+import java.util.HashSet;
+import java.util.Map;
+import java.util.Set;
+
+/**
+ * IF THERE IS NO JAVADOC RIGHT HERE, YELL AT chartl
+ *
+ * @Author chartl
+ * @Date May 19, 2010
+ */
+public class GCCalculatorWalker extends RefWalker<Pair<Set<GenomeLoc>,Boolean>, Map<GenomeLoc,Pair<Long,Long>>> {
+
+    public Map<GenomeLoc,Pair<Long,Long>> reduceInit() {
+        return new HashMap<GenomeLoc,Pair<Long,Long>>();
+    }
+
+    public Pair<Set<GenomeLoc>,Boolean> map(RefMetaDataTracker tracker, ReferenceContext ref, AlignmentContext context) {
+        if ( tracker == null || tracker.getReferenceMetaData("interval_list") == null ) {
+            return null;
+        } else {
+            Set<GenomeLoc> overlappingIntervals = new HashSet<GenomeLoc>();
+            for ( GATKFeature f : tracker.getGATKFeatureMetaData("interval_list",true) ) {
+                overlappingIntervals.add( f.getLocation() );
+            }
+
+            return new Pair<Set<GenomeLoc>,Boolean>(overlappingIntervals, ref.getBaseIndex() == BaseUtils.cIndex || ref.getBaseIndex() == BaseUtils.gIndex );
+        }
+    }
+
+    public Map<GenomeLoc,Pair<Long,Long>> reduce(Pair<Set<GenomeLoc>,Boolean> map, Map<GenomeLoc,Pair<Long,Long>> prevReduce) {
+        if ( map == null ) {
+            return prevReduce;
+        }
+
+        for ( GenomeLoc loc : map.first ) {
+            if ( ! prevReduce.keySet().contains(loc) ) {
+                prevReduce.put(loc,new Pair<Long,Long>(0l,0l));
+            }
+
+            prevReduce.get(loc).first ++;
+            if ( map.second ) {
+                prevReduce.get(loc).second ++;
+            }
+        }
+
+        return prevReduce;
+    }
+
+    public void onTraversalDone(Map<GenomeLoc,Pair<Long,Long>> reduced ) {
+        for ( Map.Entry<GenomeLoc,Pair<Long,Long>> gcCounts : reduced.entrySet() ) {
+            double gc_content = ( (double) gcCounts.getValue().second )/( (double) gcCounts.getValue().first );
+            out.printf("%s\t%.2f%n",gcCounts.getKey().toString(),100*gc_content);
+        }
+    }
+}