|
0
|
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.Arrays;
|
|
|
9 import java.util.Iterator;
|
|
|
10 import java.util.List;
|
|
|
11
|
|
|
12 import org.apache.commons.lang3.ArrayUtils;
|
|
|
13 import org.apache.log4j.Logger;
|
|
|
14 import org.broad.igv.bbfile.WigItem;
|
|
|
15
|
|
|
16 import com.beust.jcommander.Parameter;
|
|
|
17
|
|
|
18 import edu.unc.genomics.BedEntry;
|
|
|
19 import edu.unc.genomics.CommandLineTool;
|
|
|
20 import edu.unc.genomics.ReadablePathValidator;
|
|
|
21 import edu.unc.genomics.io.BedFile;
|
|
|
22 import edu.unc.genomics.io.WigFile;
|
|
|
23 import edu.unc.genomics.io.WigFileException;
|
|
|
24
|
|
|
25 public class MatrixAligner extends CommandLineTool {
|
|
|
26
|
|
|
27 private static final Logger log = Logger.getLogger(MatrixAligner.class);
|
|
|
28
|
|
|
29 @Parameter(names = {"-i", "--input"}, description = "Input file (Wig)", required = true)
|
|
|
30 public WigFile inputFile;
|
|
|
31 @Parameter(names = {"-l", "--loci"}, description = "Loci file (Bed)", required = true, validateWith = ReadablePathValidator.class)
|
|
|
32 public Path lociFile;
|
|
|
33 @Parameter(names = {"-m", "--max"}, description = "Truncate width (base pairs)")
|
|
|
34 public Integer maxWidth;
|
|
|
35 @Parameter(names = {"-o", "--output"}, description = "Output file (matrix2png format)", required = true)
|
|
|
36 public Path outputFile;
|
|
|
37
|
|
|
38 private List<BedEntry> loci;
|
|
|
39
|
|
|
40 @Override
|
|
|
41 public void run() throws IOException {
|
|
|
42 log.debug("Loading alignment intervals");
|
|
|
43 try (BedFile bed = new BedFile(lociFile)) {
|
|
|
44 loci = bed.loadAll();
|
|
|
45 }
|
|
|
46
|
|
|
47 // Compute the matrix dimensions
|
|
|
48 int leftMax = Integer.MIN_VALUE;
|
|
|
49 int rightMax = Integer.MIN_VALUE;
|
|
|
50 for (BedEntry entry : loci) {
|
|
|
51 int left = Math.abs(entry.getValue().intValue()-entry.getStart());
|
|
|
52 int right = Math.abs(entry.getValue().intValue()-entry.getStop());
|
|
|
53 if (left > leftMax) {
|
|
|
54 leftMax = left;
|
|
|
55 }
|
|
|
56 if (right > rightMax) {
|
|
|
57 rightMax = right;
|
|
|
58 }
|
|
|
59 }
|
|
|
60
|
|
|
61 int m = loci.size();
|
|
|
62 int n = leftMax + rightMax + 1;
|
|
|
63 int alignmentPoint = leftMax;
|
|
|
64 log.info("Intervals aligned into: " + m+"x"+n + " matrix");
|
|
|
65 log.info("Alignment point: " + alignmentPoint);
|
|
|
66
|
|
|
67 int leftBound = 0;
|
|
|
68 int rightBound = n-1;
|
|
|
69 if (maxWidth != null && maxWidth < n) {
|
|
|
70 log.info("Truncated to: " + m+"x"+maxWidth);
|
|
|
71 int leftAlignDistance = alignmentPoint;
|
|
|
72 int rightAlignDistance = n - alignmentPoint - 1;
|
|
|
73 int halfMax = maxWidth / 2;
|
|
|
74
|
|
|
75 if (halfMax < leftAlignDistance && halfMax < rightAlignDistance) {
|
|
|
76 leftBound = alignmentPoint - halfMax;
|
|
|
77 rightBound = alignmentPoint + halfMax;
|
|
|
78 } else {
|
|
|
79 if (leftAlignDistance <= rightAlignDistance) {
|
|
|
80 rightBound = maxWidth;
|
|
|
81 } else {
|
|
|
82 leftBound = n - maxWidth;
|
|
|
83 }
|
|
|
84 }
|
|
|
85 }
|
|
|
86
|
|
|
87 log.debug("Initializing output file");
|
|
|
88 int count = 0, skipped = 0;
|
|
|
89 try (BufferedWriter writer = Files.newBufferedWriter(outputFile, Charset.defaultCharset())) {
|
|
|
90 writer.write("ID");
|
|
|
91 for (int i = leftBound-alignmentPoint; i <= rightBound-alignmentPoint; i++) {
|
|
|
92 writer.write("\t"+i);
|
|
|
93 }
|
|
|
94 writer.newLine();
|
|
|
95
|
|
|
96 log.debug("Iterating over all intervals");
|
|
|
97 String[] row = new String[n];
|
|
|
98 for (BedEntry entry : loci) {
|
|
|
99 Iterator<WigItem> result = null;
|
|
|
100 try {
|
|
|
101 result = inputFile.query(entry);
|
|
|
102 } catch (WigFileException e) {
|
|
|
103 skipped++;
|
|
|
104 continue;
|
|
|
105 }
|
|
|
106
|
|
|
107 float[] data = WigFile.flattenData(result, entry.getStart(), entry.getStop());
|
|
|
108 // Reverse if on the crick strand
|
|
|
109 if (entry.isCrick()) {
|
|
|
110 ArrayUtils.reverse(data);
|
|
|
111 }
|
|
|
112
|
|
|
113 // Position the data in the matrix
|
|
|
114 // Locus alignment point (entry value) should be positioned over the matrix alignment point
|
|
|
115 int n1 = alignmentPoint - Math.abs(entry.getValue().intValue()-entry.getStart());
|
|
|
116 int n2 = alignmentPoint + Math.abs(entry.getValue().intValue()-entry.getStop());
|
|
|
117 assert data.length == n2-n1+1;
|
|
|
118
|
|
|
119 Arrays.fill(row, "-");
|
|
|
120 for (int i = 0; i < data.length; i++) {
|
|
|
121 if (!Float.isNaN(data[i])) {
|
|
|
122 row[n1+i] = String.valueOf(data[i]);
|
|
|
123 }
|
|
|
124 }
|
|
|
125
|
|
|
126 // Write to output
|
|
|
127 String id = ((entry.getId() == null) ? entry.getId() : "Row "+(count++));
|
|
|
128 writer.write(id);
|
|
|
129 for (int i = leftBound; i <= rightBound; i++) {
|
|
|
130 writer.write("\t"+row[i]);
|
|
|
131 }
|
|
|
132 writer.newLine();
|
|
|
133 }
|
|
|
134 }
|
|
|
135
|
|
|
136 inputFile.close();
|
|
|
137 log.info(count + " intervals processed");
|
|
|
138 log.info(skipped + " intervals skipped");
|
|
|
139 }
|
|
|
140
|
|
|
141 public static void main(String[] args) {
|
|
|
142 new MatrixAligner().instanceMain(args);
|
|
|
143 }
|
|
|
144
|
|
|
145 }
|
