|
2
|
1 package edu.unc.genomics.visualization;
|
|
|
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 import java.util.List;
|
|
|
10
|
|
|
11 import org.apache.commons.lang3.ArrayUtils;
|
|
|
12 import org.apache.log4j.Logger;
|
|
|
13 import org.broad.igv.bbfile.WigItem;
|
|
|
14
|
|
|
15 import com.beust.jcommander.Parameter;
|
|
|
16
|
|
|
17 import edu.unc.genomics.BedEntry;
|
|
|
18 import edu.unc.genomics.CommandLineTool;
|
|
|
19 import edu.unc.genomics.ReadablePathValidator;
|
|
|
20 import edu.unc.genomics.io.BedFile;
|
|
|
21 import edu.unc.genomics.io.WigFile;
|
|
|
22 import edu.unc.genomics.io.WigFileException;
|
|
|
23
|
|
|
24 public class IntervalAverager extends CommandLineTool {
|
|
|
25
|
|
|
26 private static final Logger log = Logger.getLogger(IntervalAverager.class);
|
|
|
27
|
|
|
28 @Parameter(names = {"-i", "--input"}, description = "Input file (Wig)", required = true)
|
|
|
29 public WigFile inputFile;
|
|
|
30 @Parameter(names = {"-l", "--loci"}, description = "Loci file (Bed)", required = true, validateWith = ReadablePathValidator.class)
|
|
|
31 public Path lociFile;
|
|
|
32 @Parameter(names = {"-o", "--output"}, description = "Output file (matrix2png format)", required = true)
|
|
|
33 public Path outputFile;
|
|
|
34
|
|
|
35 private List<BedEntry> loci;
|
|
|
36
|
|
|
37 @Override
|
|
|
38 public void run() throws IOException {
|
|
|
39 log.debug("Loading alignment intervals");
|
|
|
40 try (BedFile bed = new BedFile(lociFile)) {
|
|
|
41 loci = bed.loadAll();
|
|
|
42 }
|
|
|
43
|
|
|
44 // Compute the matrix dimensions
|
|
|
45 int leftMax = Integer.MIN_VALUE;
|
|
|
46 int rightMax = Integer.MIN_VALUE;
|
|
|
47 for (BedEntry entry : loci) {
|
|
|
48 int left = Math.abs(entry.getValue().intValue()-entry.getStart());
|
|
|
49 int right = Math.abs(entry.getValue().intValue()-entry.getStop());
|
|
|
50 if (left > leftMax) {
|
|
|
51 leftMax = left;
|
|
|
52 }
|
|
|
53 if (right > rightMax) {
|
|
|
54 rightMax = right;
|
|
|
55 }
|
|
|
56 }
|
|
|
57
|
|
|
58 int m = loci.size();
|
|
|
59 int n = leftMax + rightMax + 1;
|
|
|
60 int alignmentPoint = leftMax;
|
|
|
61 log.info("Intervals aligned into: " + m+"x"+n + " matrix");
|
|
|
62 log.info("Alignment point: " + alignmentPoint);
|
|
|
63
|
|
|
64 float[] sum = new float[n];
|
|
|
65 int[] counts = new int[n];
|
|
|
66 int count = 0, skipped = 0;
|
|
|
67 log.debug("Iterating over all intervals");
|
|
|
68 for (BedEntry entry : loci) {
|
|
|
69 Iterator<WigItem> result = null;
|
|
|
70 try {
|
|
|
71 result = inputFile.query(entry);
|
|
|
72 } catch (WigFileException e) {
|
|
|
73 skipped++;
|
|
|
74 continue;
|
|
|
75 }
|
|
|
76
|
|
|
77 float[] data = WigFile.flattenData(result, entry.getStart(), entry.getStop());
|
|
|
78 // Reverse if on the Crick strand
|
|
|
79 if (entry.isCrick()) {
|
|
|
80 ArrayUtils.reverse(data);
|
|
|
81 }
|
|
|
82
|
|
|
83 // Locus alignment point (entry value) should be positioned over the global alignment point
|
|
|
84 int n1 = alignmentPoint - Math.abs(entry.getValue().intValue()-entry.getStart());
|
|
|
85 int n2 = alignmentPoint + Math.abs(entry.getValue().intValue()-entry.getStop());
|
|
|
86 for (int bp = n1; bp <= n2; bp++) {
|
|
|
87 sum[bp] += data[bp-n1];
|
|
|
88 counts[bp]++;
|
|
|
89 }
|
|
|
90 }
|
|
|
91
|
|
|
92 inputFile.close();
|
|
|
93 log.info(count + " intervals processed");
|
|
|
94 log.info(skipped + " intervals skipped");
|
|
|
95
|
|
|
96 log.debug("Computing average");
|
|
|
97 float[] avg = new float[n];
|
|
|
98 for (int i = 0; i < n; i++) {
|
|
|
99 if (counts[i] == 0) {
|
|
|
100 avg[i] = Float.NaN;
|
|
|
101 } else {
|
|
|
102 avg[i] = sum[i] / counts[i];
|
|
|
103 }
|
|
|
104 }
|
|
|
105
|
|
|
106 log.debug("Writing average to output");
|
|
|
107 try (BufferedWriter writer = Files.newBufferedWriter(outputFile, Charset.defaultCharset())) {
|
|
|
108 for (int i = 0; i < n; i++) {
|
|
|
109 writer.write(i-alignmentPoint + "\t" + avg[i]);
|
|
|
110 writer.newLine();
|
|
|
111 }
|
|
|
112 }
|
|
|
113 }
|
|
|
114
|
|
|
115 public static void main(String[] args) {
|
|
|
116 new IntervalAverager().instanceMain(args);
|
|
|
117 }
|
|
|
118 } |
