Mercurial > repos > timpalpant > java_genomics_toolkit

comparison java-genomics-toolkit/src/edu/unc/genomics/visualization/IntervalAverager.java @ 0:1daf3026d231

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Upload alpha version
author timpalpant
date Mon, 13 Feb 2012 21:55:55 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:1daf3026d231
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 }