Mercurial > repos > timpalpant > java_genomics_toolkit

comparison src/edu/unc/genomics/ngs/Autocorrelation.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.commons.lang3.StringUtils;
11 import org.apache.log4j.Logger;
12 import org.broad.igv.bbfile.WigItem;
13
14 import com.beust.jcommander.Parameter;
15
16 import edu.emory.mathcs.jtransforms.fft.FloatFFT_1D;
17 import edu.unc.genomics.CommandLineTool;
18 import edu.unc.genomics.Interval;
19 import edu.unc.genomics.PositiveIntegerValidator;
20 import edu.unc.genomics.io.IntervalFile;
21 import edu.unc.genomics.io.WigFile;
22 import edu.unc.genomics.io.WigFileException;
23
24 public class Autocorrelation extends CommandLineTool {
25
26 private static final Logger log = Logger.getLogger(Autocorrelation.class);
27
28 @Parameter(names = {"-i", "--input"}, description = "Input file", required = true)
29 public WigFile wig;
30 @Parameter(names = {"-l", "--loci"}, description = "Genomic loci (Bed format)", required = true)
31 public IntervalFile<? extends Interval> loci;
32 @Parameter(names = {"-o", "--output"}, description = "Output file", required = true)
33 public Path outputFile;
34 @Parameter(names = {"-m", "--max"}, description = "Autocorrelation limit (bp)", validateWith = PositiveIntegerValidator.class)
35 public int limit = 200;
36
37 private void abs2(float[] data) {
38 for (int i = 0; i < data.length; i+=2) {
39 data[i] = data[i]*data[i] + data[i+1]*data[i+1];
40 data[i+1] = 0;
41 }
42 }
43
44 @Override
45 public void run() throws IOException {
46 try (BufferedWriter writer = Files.newBufferedWriter(outputFile, Charset.defaultCharset())) {
47 log.debug("Computing autocorrelation for each window");
48 int skipped = 0;
49 for (Interval interval : loci) {
50 if (interval.length() < limit) {
51 log.debug("Skipping interval: " + interval.toString());
52 skipped++;
53 continue;
54 }
55
56 Iterator<WigItem> wigIter;
57 try {
58 wigIter = wig.query(interval);
59 } catch (IOException | WigFileException e) {
60 log.debug("Skipping interval: " + interval.toString());
61 skipped++;
62 continue;
63 }
64
65 float[] data = WigFile.flattenData(wigIter, interval.getStart(), interval.getStop());
66
67 // Compute the autocorrelation with the Wiener-Khinchin theorem
68 FloatFFT_1D fft = new FloatFFT_1D(data.length);
69 fft.realForward(data);
70 abs2(data);
71 fft.realInverse(data, true);
72
73 writer.write(StringUtils.join(data, "\t"));
74 writer.newLine();
75 }
76
77 log.info("Skipped " + skipped + " intervals");
78 }
79 }
80
81 public static void main(String[] args) {
82 new Autocorrelation().instanceMain(args);
83 }
84 }