Mercurial > repos > timpalpant > java_genomics_toolkit

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