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