From ea7e737441bea2a391ea9a44e2428d6a8411e03c Mon Sep 17 00:00:00 2001
From: jmaguire <jmaguire@348d0f76-0448-11de-a6fe-93d51630548a>
Date: Fri, 29 Jan 2010 23:23:00 +0000
Subject: [PATCH] Two new annotations:

	1. LowMQ: fraction of reads at MQ=0 or MQ<=10.
	2. Alignability: annotate SNPs with Heng's (or anyone else's) alignability mask.



git-svn-id: file:///humgen/gsa-scr1/gsa-engineering/svn_contents/trunk@2746 348d0f76-0448-11de-a6fe-93d51630548a
---
 .../gatk/walkers/annotator/Alignability.java  | 39 +++++++++++++++++++
 .../sting/gatk/walkers/annotator/LowMQ.java   | 36 +++++++++++++++++
 2 files changed, 75 insertions(+)
 create mode 100755 java/src/org/broadinstitute/sting/gatk/walkers/annotator/Alignability.java
 create mode 100755 java/src/org/broadinstitute/sting/gatk/walkers/annotator/LowMQ.java

diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/annotator/Alignability.java b/java/src/org/broadinstitute/sting/gatk/walkers/annotator/Alignability.java
new file mode 100755
index 000000000..6c01fc9c6
--- /dev/null
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/annotator/Alignability.java
@@ -0,0 +1,39 @@
+package org.broadinstitute.sting.gatk.walkers.annotator;
+
+import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
+import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
+import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
+import org.broadinstitute.sting.gatk.refdata.TabularROD;
+import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
+import org.broadinstitute.sting.utils.pileup.PileupElement;
+import org.broadinstitute.sting.utils.genotype.Variation;
+import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
+
+import java.util.Map;
+
+
+public class Alignability implements VariantAnnotation {
+
+    public String annotate(RefMetaDataTracker tracker, 
+						   ReferenceContext ref, 
+						   Map<String, StratifiedAlignmentContext> stratifiedContexts, 
+						   Variation variation) 
+	{
+		TabularROD record = (TabularROD)(tracker.lookup("alignability", null));
+		int value;
+		if (record == null) { value = 3; }
+		else 
+		{ 
+			if (record.get("alignability") == null) 
+			{ 
+				throw new RuntimeException("ERROR: alignability column not defined in alignability input.\n");
+			}
+			value = Integer.parseInt(record.get("alignability")); 
+		}
+		return String.format("%d", value);
+    }
+
+    public String getKeyName() { return "Alignability"; }
+
+    public VCFInfoHeaderLine getDescription() { return new VCFInfoHeaderLine(getKeyName(), 1, VCFInfoHeaderLine.INFO_TYPE.Integer, "Alignability according to a mask file (3 is best, 0 is worst)"); }
+}
diff --git a/java/src/org/broadinstitute/sting/gatk/walkers/annotator/LowMQ.java b/java/src/org/broadinstitute/sting/gatk/walkers/annotator/LowMQ.java
new file mode 100755
index 000000000..92bc87bda
--- /dev/null
+++ b/java/src/org/broadinstitute/sting/gatk/walkers/annotator/LowMQ.java
@@ -0,0 +1,36 @@
+package org.broadinstitute.sting.gatk.walkers.annotator;
+
+import org.broadinstitute.sting.gatk.contexts.ReferenceContext;
+import org.broadinstitute.sting.gatk.contexts.StratifiedAlignmentContext;
+import org.broadinstitute.sting.gatk.refdata.RefMetaDataTracker;
+import org.broadinstitute.sting.utils.pileup.ReadBackedPileup;
+import org.broadinstitute.sting.utils.pileup.PileupElement;
+import org.broadinstitute.sting.utils.genotype.Variation;
+import org.broadinstitute.sting.utils.genotype.vcf.VCFInfoHeaderLine;
+
+import java.util.Map;
+
+
+public class LowMQ implements VariantAnnotation {
+
+    public String annotate(RefMetaDataTracker tracker, ReferenceContext ref, Map<String, StratifiedAlignmentContext> stratifiedContexts, Variation variation) {
+        double mq0 = 0;
+		double mq10 = 0;
+		double total = 0;
+        for ( String sample : stratifiedContexts.keySet() ) 
+		{
+            ReadBackedPileup pileup = stratifiedContexts.get(sample).getContext(StratifiedAlignmentContext.StratifiedContextType.COMPLETE).getBasePileup();
+            for (PileupElement p : pileup ) 
+			{
+                if ( p.getMappingQual() == 0 )  { mq0 += 1; }
+                if ( p.getMappingQual() <= 10 ) { mq10 += 1; }
+				total += 1; 
+            }
+        }
+        return String.format("%.04f,%.04f,%.00f", mq0/total, mq10/total, total);
+    }
+
+    public String getKeyName() { return "LowMQ"; }
+
+    public VCFInfoHeaderLine getDescription() { return new VCFInfoHeaderLine(getKeyName(), 1, VCFInfoHeaderLine.INFO_TYPE.Integer, "3-tuple: <fraction of reads with MQ=0>,<fraction of reads with MQ<=10>,<total nubmer of reads>"); }
+}