|
2
|
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 }
|
