OUTPUT_DIRECTORY.mkdirs();
}
IOUtil.assertDirectoryIsWritable(OUTPUT_DIRECTORY);
// Load up the targets and the reference
final IntervalList targets;
final IntervalList originalTargets = IntervalList.fromFile(TARGETS);
{
// Apply padding
final IntervalList padded = new IntervalList(originalTargets.getHeader());
final SAMSequenceDictionary dict = padded.getHeader().getSequenceDictionary();
for (final Interval i : originalTargets.getIntervals()) {
padded.add(new Interval(i.getSequence(),
Math.max(i.getStart() - PADDING, 1),
Math.min(i.getEnd() + PADDING, dict.getSequence(i.getSequence()).getSequenceLength()),
i.isNegativeStrand(),
i.getName()));
}
log.info("Starting with " + padded.size() + " targets.");
padded.unique();
log.info("After uniquing " + padded.size() + " targets remain.");
if (MERGE_NEARBY_TARGETS) {
final ListIterator<Interval> iterator = padded.getIntervals().listIterator();
Interval previous = iterator.next();
targets = new IntervalList(padded.getHeader());
while (iterator.hasNext()) {
final Interval next = iterator.next();
if (previous.getSequence().equals(next.getSequence()) &&
estimateBaits(previous.getStart(), previous.getEnd()) + estimateBaits(next.getStart(), next.getEnd()) >=
estimateBaits(previous.getStart(), next.getEnd())) {
previous = new Interval(previous.getSequence(),
previous.getStart(),
Math.max(previous.getEnd(), next.getEnd()),
previous.isNegativeStrand(),
previous.getName());
}
else {
targets.add(previous);
previous = next;
}
}
if (previous != null) targets.add(previous);
log.info("After collapsing nearby targets " + targets.size() + " targets remain.");
}
else {
targets = padded;
}
}
final ReferenceSequenceFileWalker referenceWalker = new ReferenceSequenceFileWalker(REFERENCE_SEQUENCE);
// Check that the reference and the target list have matching headers
SequenceUtil.assertSequenceDictionariesEqual(referenceWalker.getSequenceDictionary(),
targets.getHeader().getSequenceDictionary());
// Design the baits!
int discardedBaits = 0;
final IntervalList baits = new IntervalList(targets.getHeader());
for (final Interval target : targets) {
final int sequenceIndex = targets.getHeader().getSequenceIndex(target.getSequence());
final ReferenceSequence reference = referenceWalker.get(sequenceIndex);
for (final Bait bait : DESIGN_STRATEGY.design(this, target, reference)) {
if (bait.length() != BAIT_SIZE) {
throw new PicardException("Bait designed at wrong length: " + bait);
}
if (bait.getMaskedBaseCount() <= REPEAT_TOLERANCE) {
baits.add(bait);
for (final byte b : bait.getBases()) {
final byte upper = StringUtil.toUpperCase(b);
if (upper != 'A' && upper != 'C' && upper != 'G' && upper != 'T') {
log.warn("Bait contains non-synthesizable bases: " + bait);
}
}
}
else {
log.debug("Discarding bait: " + bait);
discardedBaits++;
}
}
}
calculateStatistics(targets, baits);
log.info("Designed and kept " + baits.size() + " baits, discarded " + discardedBaits);
// Write out some files!
originalTargets.write(new File(OUTPUT_DIRECTORY, DESIGN_NAME + ".targets.interval_list"));
baits.write(new File(OUTPUT_DIRECTORY, DESIGN_NAME + ".baits.interval_list"));
writeParametersFile(new File(OUTPUT_DIRECTORY, DESIGN_NAME + ".design_parameters.txt"));
writeDesignFastaFile(new File(OUTPUT_DIRECTORY, DESIGN_NAME + ".design.fasta"), baits);
if (POOL_SIZE > 0) writePoolFiles(OUTPUT_DIRECTORY, DESIGN_NAME, baits);
return 0;