Mercurial > repos > timpalpant > java_genomics_toolkit

comparison java-genomics-toolkit/src/edu/unc/genomics/ngs/RollingReadLength.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.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.log4j.Logger;
11
12 import com.beust.jcommander.Parameter;
13
14 import edu.ucsc.genome.TrackHeader;
15 import edu.unc.genomics.Assembly;
16 import edu.unc.genomics.CommandLineTool;
17 import edu.unc.genomics.Interval;
18 import edu.unc.genomics.io.IntervalFile;
19
20 public class RollingReadLength extends CommandLineTool {
21
22 private static final Logger log = Logger.getLogger(RollingReadLength.class);
23
24 @Parameter(names = {"-i", "--input"}, description = "Input file (reads)", required = true)
25 public IntervalFile<? extends Interval> intervalFile;
26 @Parameter(names = {"-a", "--assembly"}, description = "Genome assembly", required = true)
27 public Assembly assembly;
28 @Parameter(names = {"-o", "--output"}, description = "Output file (Wig)", required = true)
29 public Path outputFile;
30
31 @Override
32 public void run() throws IOException {
33 log.debug("Initializing output file");
34 int mapped = 0;
35 try (BufferedWriter writer = Files.newBufferedWriter(outputFile, Charset.defaultCharset())) {
36 // Write the Wiggle track header to the output file
37 TrackHeader header = new TrackHeader("wiggle_0");
38 header.setName("Converted " + intervalFile.getPath().getFileName());
39 header.setDescription("Converted " + intervalFile.getPath().getFileName());
40 writer.write(header.toString());
41 writer.newLine();
42
43 // Process each chromosome in the assembly
44 for (String chr : assembly) {
45 log.debug("Processing chromosome " + chr);
46 // Write the contig header to the output file
47 writer.write("fixedStep chrom="+chr+" start=1 step=1 span=1");
48 writer.newLine();
49
50 int start = 1;
51 while (start < assembly.getChrLength(chr)) {
52 int stop = start + DEFAULT_CHUNK_SIZE - 1;
53 int length = stop - start + 1;
54 int[] sum = new int[length];
55 int[] count = new int[length];
56
57 Iterator<? extends Interval> it = intervalFile.query(chr, start, stop);
58 while (it.hasNext()) {
59 Interval entry = it.next();
60 for (int i = entry.getStart(); i <= entry.getStop(); i++) {
61 sum[i-start] += entry.length();
62 count[i-start]++;
63 }
64 mapped++;
65 }
66
67 // Write the average at each base pair to the output file
68 for (int i = 0; i < sum.length; i++) {
69 if (count[i] == 0) {
70 writer.write(String.valueOf(Float.NaN));
71 } else {
72 writer.write(String.valueOf(sum[i]/count[i]));
73 }
74 writer.newLine();
75 }
76
77 // Process the next chunk
78 start = stop + 1;
79 }
80 }
81 }
82
83 log.info("Mapped "+mapped+" reads");
84 }
85
86 public static void main(String[] args) {
87 new RollingReadLength().instanceMain(args);
88 }
89 }