Mercurial > repos > timpalpant > java_genomics_toolkit

comparison src/edu/unc/genomics/ngs/IntervalStats.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.nio.file.Paths;
9 import java.util.ArrayList;
10 import java.util.Iterator;
11 import java.util.List;
12
13 import org.apache.commons.lang3.StringUtils;
14 import org.apache.commons.math.stat.descriptive.SummaryStatistics;
15 import org.apache.log4j.Logger;
16 import org.broad.igv.bbfile.WigItem;
17
18 import com.beust.jcommander.Parameter;
19
20 import edu.unc.genomics.CommandLineTool;
21 import edu.unc.genomics.Interval;
22 import edu.unc.genomics.io.IntervalFile;
23 import edu.unc.genomics.io.WigFile;
24 import edu.unc.genomics.io.WigFileException;
25
26 public class IntervalStats extends CommandLineTool {
27
28 private static final Logger log = Logger.getLogger(IntervalStats.class);
29
30 @Parameter(description = "Input files", required = true)
31 public List<String> inputFiles = new ArrayList<String>();
32 @Parameter(names = {"-l", "--loci"}, description = "Loci file (Bed)", required = true)
33 public IntervalFile<? extends Interval> lociFile;
34 @Parameter(names = {"-o", "--output"}, description = "Output file", required = true)
35 public Path outputFile;
36
37 private List<WigFile> wigs = new ArrayList<>();
38
39 @Override
40 public void run() throws IOException {
41 log.debug("Initializing input Wig file(s)");
42 for (String inputFile : inputFiles) {
43 try {
44 WigFile wig = WigFile.autodetect(Paths.get(inputFile));
45 wigs.add(wig);
46 } catch (WigFileException e) {
47 log.error("Error initializing Wig input file: " + inputFile);
48 e.printStackTrace();
49 throw new RuntimeException("Error initializing Wig input file: " + inputFile);
50 }
51 }
52
53 log.debug("Initializing output file");
54 int count = 0, skipped = 0;
55 try (BufferedWriter writer = Files.newBufferedWriter(outputFile, Charset.defaultCharset())) {
56 writer.write("#Chr\tStart\tStop\tID\tValue\tStrand");
57 for (String inputFile : inputFiles) {
58 Path p = Paths.get(inputFile);
59 writer.write("\t" + p.getFileName().toString());
60 }
61 writer.newLine();
62
63 log.debug("Iterating over all intervals and computing statistics");
64 SummaryStatistics stats = new SummaryStatistics();
65 for (Interval interval : lociFile) {
66 List<Double> means = new ArrayList<>(wigs.size());
67 for (WigFile wig : wigs) {
68 stats.clear();
69 try {
70 Iterator<WigItem> result = wig.query(interval);
71 while(result.hasNext()) {
72 WigItem item = result.next();
73 for (int i = item.getStartBase(); i <= item.getEndBase(); i++) {
74 stats.addValue(item.getWigValue());
75 }
76 }
77 means.add(stats.getMean());
78 } catch (WigFileException e) {
79 means.add(Double.NaN);
80 skipped++;
81 }
82 }
83
84 writer.write(interval.toBed() + "\t" + StringUtils.join(means, "\t"));
85 writer.newLine();
86 count++;
87 }
88 }
89
90 lociFile.close();
91 for (WigFile wig : wigs) {
92 wig.close();
93 }
94 log.info(count + " intervals processed");
95 log.info(skipped + " interval skipped");
96 }
97
98 public static void main(String[] args) {
99 new IntervalStats().instanceMain(args);
100 }
101
102 }