diff java-genomics-toolkit/src/edu/unc/genomics/wigmath/GaussianSmooth.java @ 0:1daf3026d231

Upload alpha version
author timpalpant
date Mon, 13 Feb 2012 21:55:55 -0500
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/java-genomics-toolkit/src/edu/unc/genomics/wigmath/GaussianSmooth.java	Mon Feb 13 21:55:55 2012 -0500
@@ -0,0 +1,74 @@
+package edu.unc.genomics.wigmath;
+
+import java.io.IOException;
+import java.util.Iterator;
+
+import org.apache.log4j.Logger;
+import org.broad.igv.bbfile.WigItem;
+
+import com.beust.jcommander.Parameter;
+
+import edu.unc.genomics.io.WigFile;
+import edu.unc.genomics.io.WigFileException;
+
+public class GaussianSmooth extends WigMathTool {
+
+	private static final Logger log = Logger.getLogger(GaussianSmooth.class);
+
+	@Parameter(names = {"-i", "--input"}, description = "Input file", required = true)
+	public WigFile inputFile;
+	@Parameter(names = {"-s", "--stdev"}, description = "Standard deviation of Gaussian (bp)")
+	public int stdev = 20;
+	
+	float[] filter;
+
+	@Override
+	public void setup() {
+		inputs.add(inputFile);
+		
+		// Use a window size equal to +/- 3 SD's
+		filter = new float[6*stdev+1];
+		float sum = 0;
+		for (int i = 0; i < filter.length; i++) {
+			float x = i - 3*stdev;
+			float value = (float) Math.exp(-(x*x) / (2*stdev*stdev));
+			filter[i] = value;
+			sum += value;
+		}
+		for (int i = 0; i < filter.length; i++) {
+			filter[i] /= sum;
+		}
+	}
+	
+	@Override
+	public float[] compute(String chr, int start, int stop) throws IOException, WigFileException {
+		log.debug("Smoothing chunk "+chr+":"+start+"-"+stop);
+		
+		// Pad the query for smoothing
+		int paddedStart = Math.max(start-3*stdev, inputFile.getChrStart(chr));
+		int paddedStop = Math.min(stop+3*stdev, inputFile.getChrStop(chr));
+		Iterator<WigItem> result = inputFile.query(chr, paddedStart, paddedStop);
+		float[] data = WigFile.flattenData(result, start-3*stdev, stop+3*stdev);
+		
+		// Convolve the data with the filter
+		float[] smoothed = new float[stop-start+1];
+		for (int i = 0; i < smoothed.length; i++) {
+			for (int j = 0; j < filter.length; j++) {
+				smoothed[i] += data[i+j] * filter[j];
+			}
+		}
+		
+		return smoothed;
+	}
+	
+	
+	/**
+	 * @param args
+	 * @throws WigFileException 
+	 * @throws IOException 
+	 */
+	public static void main(String[] args) throws IOException, WigFileException {
+		new GaussianSmooth().instanceMain(args);
+	}
+
+}