Mercurial > repos > timpalpant > java_genomics_toolkit
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); + } + +}