Mercurial > repos > timpalpant > java_genomics_toolkit

comparison src/edu/unc/genomics/ngs/PowerSpectrum.java @ 2:e16016635b2a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.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 }