Mercurial > repos > timpalpant > java_genomics_toolkit
comparison src/edu/unc/genomics/wigmath/GaussianSmooth.java @ 2:e16016635b2a
Uploaded
author | timpalpant |
---|---|
date | Mon, 13 Feb 2012 22:12:06 -0500 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
1:a54db233ee3d | 2:e16016635b2a |
---|---|
1 package edu.unc.genomics.wigmath; | |
2 | |
3 import java.io.IOException; | |
4 import java.util.Iterator; | |
5 | |
6 import org.apache.log4j.Logger; | |
7 import org.broad.igv.bbfile.WigItem; | |
8 | |
9 import com.beust.jcommander.Parameter; | |
10 | |
11 import edu.unc.genomics.io.WigFile; | |
12 import edu.unc.genomics.io.WigFileException; | |
13 | |
14 public class GaussianSmooth extends WigMathTool { | |
15 | |
16 private static final Logger log = Logger.getLogger(GaussianSmooth.class); | |
17 | |
18 @Parameter(names = {"-i", "--input"}, description = "Input file", required = true) | |
19 public WigFile inputFile; | |
20 @Parameter(names = {"-s", "--stdev"}, description = "Standard deviation of Gaussian (bp)") | |
21 public int stdev = 20; | |
22 | |
23 float[] filter; | |
24 | |
25 @Override | |
26 public void setup() { | |
27 inputs.add(inputFile); | |
28 | |
29 // Use a window size equal to +/- 3 SD's | |
30 filter = new float[6*stdev+1]; | |
31 float sum = 0; | |
32 for (int i = 0; i < filter.length; i++) { | |
33 float x = i - 3*stdev; | |
34 float value = (float) Math.exp(-(x*x) / (2*stdev*stdev)); | |
35 filter[i] = value; | |
36 sum += value; | |
37 } | |
38 for (int i = 0; i < filter.length; i++) { | |
39 filter[i] /= sum; | |
40 } | |
41 } | |
42 | |
43 @Override | |
44 public float[] compute(String chr, int start, int stop) throws IOException, WigFileException { | |
45 log.debug("Smoothing chunk "+chr+":"+start+"-"+stop); | |
46 | |
47 // Pad the query for smoothing | |
48 int paddedStart = Math.max(start-3*stdev, inputFile.getChrStart(chr)); | |
49 int paddedStop = Math.min(stop+3*stdev, inputFile.getChrStop(chr)); | |
50 Iterator<WigItem> result = inputFile.query(chr, paddedStart, paddedStop); | |
51 float[] data = WigFile.flattenData(result, start-3*stdev, stop+3*stdev); | |
52 | |
53 // Convolve the data with the filter | |
54 float[] smoothed = new float[stop-start+1]; | |
55 for (int i = 0; i < smoothed.length; i++) { | |
56 for (int j = 0; j < filter.length; j++) { | |
57 smoothed[i] += data[i+j] * filter[j]; | |
58 } | |
59 } | |
60 | |
61 return smoothed; | |
62 } | |
63 | |
64 | |
65 /** | |
66 * @param args | |
67 * @throws WigFileException | |
68 * @throws IOException | |
69 */ | |
70 public static void main(String[] args) throws IOException, WigFileException { | |
71 new GaussianSmooth().instanceMain(args); | |
72 } | |
73 | |
74 } |