|
2
|
1 package edu.unc.genomics.ngs;
|
|
|
2
|
|
|
3 import java.io.BufferedWriter;
|
|
|
4 import java.io.IOException;
|
|
|
5 import java.nio.charset.Charset;
|
|
|
6 import java.nio.file.Files;
|
|
|
7 import java.nio.file.Path;
|
|
|
8 import java.util.Iterator;
|
|
|
9
|
|
|
10 import org.apache.log4j.Logger;
|
|
|
11 import org.broad.igv.bbfile.WigItem;
|
|
|
12
|
|
|
13 import com.beust.jcommander.Parameter;
|
|
|
14
|
|
|
15 import edu.emory.mathcs.jtransforms.fft.FloatFFT_1D;
|
|
|
16 import edu.unc.genomics.CommandLineTool;
|
|
|
17 import edu.unc.genomics.Interval;
|
|
|
18 import edu.unc.genomics.io.IntervalFile;
|
|
|
19 import edu.unc.genomics.io.WigFile;
|
|
|
20 import edu.unc.genomics.io.WigFileException;
|
|
|
21
|
|
|
22 public class PowerSpectrum extends CommandLineTool {
|
|
|
23
|
|
|
24 private static final Logger log = Logger.getLogger(PowerSpectrum.class);
|
|
|
25
|
|
|
26 @Parameter(names = {"-i", "--input"}, description = "Input file (Wig)", required = true)
|
|
|
27 public WigFile inputFile;
|
|
|
28 @Parameter(names = {"-l", "--loci"}, description = "Genomic loci (Bed format)", required = true)
|
|
|
29 public IntervalFile<? extends Interval> loci;
|
|
|
30 @Parameter(names = {"-o", "--output"}, description = "Output file (tabular)", required = true)
|
|
|
31 public Path outputFile;
|
|
|
32
|
|
|
33 /**
|
|
|
34 * Computes the power spectrum from FFT data
|
|
|
35 * taking into accound even/odd length arrays
|
|
|
36 * refer to JTransforms documentation for layout of the FFT data
|
|
|
37 * @param f
|
|
|
38 * @return
|
|
|
39 */
|
|
|
40 private float[] abs2(float[] f) {
|
|
|
41 int n = f.length;
|
|
|
42 float[] ps = new float[n/2+1];
|
|
|
43 // DC component
|
|
|
44 ps[0] = (f[0]*f[0]) / (n*n);
|
|
|
45
|
|
|
46 // Even
|
|
|
47 if (n % 2 == 0) {
|
|
|
48 for (int k = 1; k < n/2; k++) {
|
|
|
49 ps[k] = f[2*k]*f[2*k] + f[2*k+1]*f[2*k+1];
|
|
|
50 }
|
|
|
51 ps[n/2] = f[1]*f[1];
|
|
|
52 // Odd
|
|
|
53 } else {
|
|
|
54 for (int k = 1; k < (n-1)/2; k++) {
|
|
|
55 ps[k] = f[2*k]*f[2*k] + f[2*k+1]*f[2*k+1];
|
|
|
56 }
|
|
|
57
|
|
|
58 ps[(n-1)/2] = f[n-1]*f[n-1] + f[1]*f[1];
|
|
|
59 }
|
|
|
60
|
|
|
61 return ps;
|
|
|
62 }
|
|
|
63
|
|
|
64 public void run() throws IOException {
|
|
|
65 try (BufferedWriter writer = Files.newBufferedWriter(outputFile, Charset.defaultCharset())) {
|
|
|
66 log.debug("Computing power spectrum for each window");
|
|
|
67 int skipped = 0;
|
|
|
68 for (Interval interval : loci) {
|
|
|
69 Iterator<WigItem> wigIter;
|
|
|
70 try {
|
|
|
71 wigIter = inputFile.query(interval);
|
|
|
72 } catch (IOException | WigFileException e) {
|
|
|
73 log.debug("Skipping interval: " + interval.toString());
|
|
|
74 skipped++;
|
|
|
75 continue;
|
|
|
76 }
|
|
|
77
|
|
|
78 float[] data = WigFile.flattenData(wigIter, interval.getStart(), interval.getStop());
|
|
|
79 // Compute the power spectrum
|
|
|
80 FloatFFT_1D fft = new FloatFFT_1D(data.length);
|
|
|
81 fft.realForward(data);
|
|
|
82 float[] ps = abs2(data);
|
|
|
83 // and normalize the power spectrum
|
|
|
84 float sum = 0;
|
|
|
85 for (int i = 1; i < ps.length; i++) {
|
|
|
86 sum += ps[i];
|
|
|
87 }
|
|
|
88 for (int i = 1; i < ps.length; i++) {
|
|
|
89 ps[i] /= sum;
|
|
|
90 }
|
|
|
91
|
|
|
92 writer.write(interval.toBed());
|
|
|
93 for (int i = 1; i < Math.min(ps.length, 40); i++) {
|
|
|
94 writer.write("\t"+ps[i]);
|
|
|
95 }
|
|
|
96 writer.newLine();
|
|
|
97 }
|
|
|
98
|
|
|
99 log.info("Skipped " + skipped + " intervals");
|
|
|
100 }
|
|
|
101 }
|
|
|
102
|
|
|
103 public static void main(String[] args) {
|
|
|
104 new PowerSpectrum().instanceMain(args);
|
|
|
105 }
|
|
|
106 } |
