New upstream version 2.18.17+dfsg (4b641b5f) · Commits · Debian Med / picard-tools

src/main/java/picard/analysis/TheoreticalSensitivity.java

+19 −5

Original line number	Diff line number	Diff line
		@@ -244,6 +244,11 @@ public class TheoreticalSensitivity {
		* @return Theoretical sensitivity for the given arguments at a constant depth.
		*/
		public static double sensitivityAtConstantDepth(final int depth, final Histogram<Integer> qualityHistogram, final double logOddsThreshold, final int sampleSize, final double alleleFraction, final long randomSeed) {
		// If the depth is 0 at a particular locus, the sensitivity is trivially 0.0.
		if (depth == 0) {
		return 0.0;
		}

		final RouletteWheel qualityRW = new RouletteWheel(trimDistribution(normalizeHistogram(qualityHistogram)));
		final Random randomNumberGenerator = new Random(randomSeed);
		final RandomGenerator rg = new Well19937c(randomSeed);
		@@ -350,13 +355,22 @@ public class TheoreticalSensitivity {
		/**
		* Removes trailing zeros in a distribution. The purpose of this function is to prevent other
		* functions from evaluating in regions where the distribution has zero probability.
		* If the array consists of only zeros or is empty, an empty array is returned.
		*
		* @param distribution Distribution of base qualities
		* @return Distribution of base qualities removing any trailing zeros
		*/
		static double[] trimDistribution(final double[] distribution) {
		int endOfDistribution = distribution.length - 1;
		while(distribution[endOfDistribution] == 0) {
		int endOfDistribution = distribution.length;
		if (distribution.length == 0) {
		return distribution;
		}

		// Locate last element in array that is non-zero if it exists.
		// If there are no non-zero elements, endOfDistribution should be 0.
		while (distribution[endOfDistribution - 1] == 0) {
		endOfDistribution--;
		if (endOfDistribution == 0) break;
		}

		// Remove trailing zeros and return.

src/main/java/picard/fingerprint/CheckFingerprint.java

+14 −7

Original line number	Diff line number	Diff line
		@@ -224,7 +224,9 @@ public class CheckFingerprint extends CommandLineProgram {
		outputDetailMetricsFile = DETAIL_OUTPUT;
		outputSummaryMetricsFile = SUMMARY_OUTPUT;
		} else {
		if (!OUTPUT.endsWith(".")) OUTPUT = OUTPUT + ".";
		if (!OUTPUT.endsWith(".")) {
		OUTPUT += ".";
		}
		outputDetailMetricsFile = new File(OUTPUT + FINGERPRINT_DETAIL_FILE_SUFFIX);
		outputSummaryMetricsFile = new File(OUTPUT + FINGERPRINT_SUMMARY_FILE_SUFFIX);
		}
		@@ -245,8 +247,7 @@ public class CheckFingerprint extends CommandLineProgram {
		List<FingerprintResults> results;

		String observedSampleAlias = null;
		final boolean isBamOrSamFile = isBamOrSam(inputPath);
		if (isBamOrSamFile) {
		if (isBamOrSam(inputPath)) {
		SequenceUtil.assertSequenceDictionariesEqual(SAMSequenceDictionaryExtractor.extractDictionary(inputPath), SAMSequenceDictionaryExtractor.extractDictionary(genotypesPath), true);
		SequenceUtil.assertSequenceDictionariesEqual(SAMSequenceDictionaryExtractor.extractDictionary(inputPath), checker.getHeader().getSequenceDictionary(), true);

		@@ -311,6 +312,7 @@ public class CheckFingerprint extends CommandLineProgram {
		final MetricsFile<FingerprintingSummaryMetrics, ?> summaryFile = getMetricsFile();
		final MetricsFile<FingerprintingDetailMetrics, ?> detailsFile = getMetricsFile();

		boolean allZeroLod = true;
		for (final FingerprintResults fpr : results) {
		final MatchResults mr = fpr.getMatchResults().first();

		@@ -365,11 +367,20 @@ public class CheckFingerprint extends CommandLineProgram {

		summaryFile.addMetric(metrics);
		log.info("Read Group: " + metrics.READ_GROUP + " / " + observedSampleAlias + " vs. " + metrics.SAMPLE + ": LOD = " + metrics.LOD_EXPECTED_SAMPLE);

		allZeroLod &= metrics.LOD_EXPECTED_SAMPLE == 0;
		}

		summaryFile.write(outputSummaryMetricsFile);
		detailsFile.write(outputDetailMetricsFile);

		if (allZeroLod) {
		log.error("No non-zero results found. This is likely an error. " +
		"Probable cause: EXPECTED_SAMPLE (if provided) or the sample name from INPUT (if EXPECTED_SAMPLE isn't provided)" +
		"isn't a sample in GENOTYPES file.");
		return 1;
		}

		return 0;
		}

		@@ -389,10 +400,6 @@ public class CheckFingerprint extends CommandLineProgram {
		return super.customCommandLineValidation();
		}

		static boolean isBamOrSam(final File f) {
		return isBamOrSam(f.toPath());
		}

		static boolean isBamOrSam(final Path p) {
		return (p.toUri().getRawPath().endsWith(BamFileIoUtils.BAM_FILE_EXTENSION) \|\| p.toUri().getRawPath().endsWith(IOUtil.SAM_FILE_EXTENSION));
		}

src/main/java/picard/illumina/IlluminaBasecallsToFastq.java

+1 −1

Original line number	Diff line number	Diff line
		@@ -281,7 +281,7 @@ public class IlluminaBasecallsToFastq extends CommandLineProgram {
		FIRST_TILE, TILE_LIMIT, queryNameComparator,
		new FastqRecordsForClusterCodec(readStructure.templates.length(),
		readStructure.sampleBarcodes.length(), readStructure.molecularBarcode.length()),
		FastqRecordsForCluster.class, bclQualityEvaluationStrategy, IGNORE_UNEXPECTED_BARCODES);
		FastqRecordsForCluster.class, bclQualityEvaluationStrategy, INCLUDE_NON_PF_READS, IGNORE_UNEXPECTED_BARCODES);
		} else {
		basecallsConverter = new IlluminaBasecallsConverter<>(BASECALLS_DIR, BARCODES_DIR, LANE, readStructure,
		sampleBarcodeFastqWriterMap, demultiplex, Math.max(1, MAX_READS_IN_RAM_PER_TILE / readsPerCluster), TMP_DIR, NUM_PROCESSORS,

src/main/java/picard/illumina/IlluminaBasecallsToSam.java

+1 −1

Original line number	Diff line number	Diff line
		@@ -312,7 +312,7 @@ public class IlluminaBasecallsToSam extends CommandLineProgram {
		TMP_DIR, NUM_PROCESSORS,
		FIRST_TILE, TILE_LIMIT, new QueryNameComparator(),
		new Codec(numOutputRecords),
		SAMRecordsForCluster.class, bclQualityEvaluationStrategy, IGNORE_UNEXPECTED_BARCODES);
		SAMRecordsForCluster.class, bclQualityEvaluationStrategy, INCLUDE_NON_PF_READS, IGNORE_UNEXPECTED_BARCODES);
		} else {
		basecallsConverter = new IlluminaBasecallsConverter<>(BASECALLS_DIR, BARCODES_DIR, LANE, readStructure,
		barcodeSamWriterMap, true, MAX_READS_IN_RAM_PER_TILE / numOutputRecords, TMP_DIR, NUM_PROCESSORS, FORCE_GC,

src/main/java/picard/illumina/NewIlluminaBasecallsConverter.java

+98 −46

Original line number	Diff line number	Diff line
		@@ -14,11 +14,9 @@ import picard.illumina.parser.ReadStructure;
		import picard.illumina.parser.readers.AbstractIlluminaPositionFileReader;
		import picard.illumina.parser.readers.BclQualityEvaluationStrategy;
		import picard.illumina.parser.readers.LocsFileReader;
		import picard.util.ThreadPoolExecutorUtil;
		import picard.util.ThreadPoolExecutorWithExceptions;

		import java.io.File;
		import java.time.Duration;
		import java.util.ArrayList;
		import java.util.Arrays;
		import java.util.Collections;
		@@ -26,7 +24,9 @@ import java.util.Comparator;
		import java.util.HashMap;
		import java.util.List;
		import java.util.Map;
		import java.util.concurrent.ConcurrentHashMap;
		import java.util.concurrent.ThreadPoolExecutor;
		import java.util.concurrent.TimeUnit;
		import java.util.regex.Matcher;
		import java.util.regex.Pattern;

		@@ -36,10 +36,12 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		private final List<AbstractIlluminaPositionFileReader.PositionInfo> locs = new ArrayList<>();
		private final File[] filterFiles;
		private final Map<String, ThreadPoolExecutorWithExceptions> barcodeWriterThreads = new HashMap<>();
		private final Map<Integer, List<RecordWriter>> completedWork = Collections.synchronizedMap(new HashMap<>());
		private final Map<Integer, File> barcodesFiles = new HashMap<>();
		private final boolean includeNonPfReads;

		/**
		* Constructor setting includeNonPfReads to true by default.
		*
		* @param basecallsDir Where to read basecalls from.
		* @param barcodesDir Where to read barcodes from (optional; use basecallsDir if not specified).
		* @param lane What lane to process.
		@@ -71,11 +73,52 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		final Class<CLUSTER_OUTPUT_RECORD> outputRecordClass,
		final BclQualityEvaluationStrategy bclQualityEvaluationStrategy,
		final boolean ignoreUnexpectedBarcodes) {
		this(basecallsDir, barcodesDir, lane, readStructure, barcodeRecordWriterMap,
		demultiplex, maxReadsInRamPerTile, tmpDirs, numProcessors, firstTile,
		tileLimit, outputRecordComparator, codecPrototype, outputRecordClass,
		bclQualityEvaluationStrategy, true, ignoreUnexpectedBarcodes);
		}

		/**
		* @param basecallsDir Where to read basecalls from.
		* @param barcodesDir Where to read barcodes from (optional; use basecallsDir if not specified).
		* @param lane What lane to process.
		* @param readStructure How to interpret each cluster.
		* @param barcodeRecordWriterMap Map from barcode to CLUSTER_OUTPUT_RECORD writer. If demultiplex is false, must contain
		* one writer stored with key=null.
		* @param demultiplex If true, output is split by barcode, otherwise all are written to the same output stream.
		* @param maxReadsInRamPerTile Configures number of reads each tile will store in RAM before spilling to disk.
		* @param tmpDirs For SortingCollection spilling.
		* @param numProcessors Controls number of threads. If <= 0, the number of threads allocated is
		* available cores - numProcessors.
		* @param firstTile (For debugging) If non-null, start processing at this tile.
		* @param tileLimit (For debugging) If non-null, process no more than this many tiles.
		* @param outputRecordComparator For sorting output records within a single tile.
		* @param codecPrototype For spilling output records to disk.
		* @param outputRecordClass Inconveniently needed to create SortingCollections.
		* @param includeNonPfReads If true, will include ALL reads (including those which do not have PF set)
		* @param ignoreUnexpectedBarcodes If true, will ignore reads whose called barcode is not found in barcodeRecordWriterMap,
		*/
		public NewIlluminaBasecallsConverter(final File basecallsDir, final File barcodesDir, final int lane,
		final ReadStructure readStructure,
		final Map<String, ? extends ConvertedClusterDataWriter<CLUSTER_OUTPUT_RECORD>> barcodeRecordWriterMap,
		final boolean demultiplex,
		final int maxReadsInRamPerTile,
		final List<File> tmpDirs, final int numProcessors,
		final Integer firstTile,
		final Integer tileLimit,
		final Comparator<CLUSTER_OUTPUT_RECORD> outputRecordComparator,
		final SortingCollection.Codec<CLUSTER_OUTPUT_RECORD> codecPrototype,
		final Class<CLUSTER_OUTPUT_RECORD> outputRecordClass,
		final BclQualityEvaluationStrategy bclQualityEvaluationStrategy,
		final boolean includeNonPfReads,
		final boolean ignoreUnexpectedBarcodes) {

		super(barcodeRecordWriterMap, maxReadsInRamPerTile, tmpDirs, codecPrototype, ignoreUnexpectedBarcodes,
		demultiplex, outputRecordComparator, bclQualityEvaluationStrategy,
		outputRecordClass, numProcessors, new IlluminaDataProviderFactory(basecallsDir,
		barcodesDir, lane, readStructure, bclQualityEvaluationStrategy));
		this.includeNonPfReads = includeNonPfReads;
		this.tiles = new ArrayList<>();

		barcodeRecordWriterMap.keySet().forEach(barcode -> barcodeWriterThreads.put(barcode, new ThreadPoolExecutorWithExceptions(1)));
		@@ -153,26 +196,31 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		completedWorkExecutor.shutdown();

		//thread by surface tile
		final ThreadPoolExecutorWithExceptions tileProcessingExecutor = new ThreadPoolExecutorWithExceptions(numThreads);
		final ThreadPoolExecutor tileProcessingExecutor = new ThreadPoolExecutorWithExceptions(numThreads);
		workChecker.setTileProcessingExecutor(tileProcessingExecutor);

		for (final Integer tile : tiles) {
		tileProcessingExecutor.submit(new TileProcessor(tile, barcodesFiles.get(tile)));
		tileProcessingExecutor.submit(new TileProcessor(tile, barcodesFiles.get(tile), workChecker));
		}

		tileProcessingExecutor.shutdown();

		ThreadPoolExecutorUtil.awaitThreadPoolTermination("Reading executor", tileProcessingExecutor, Duration.ofMinutes(5));
		awaitThreadPoolTermination("Reading executor", tileProcessingExecutor);
		awaitThreadPoolTermination("Tile completion executor", completedWorkExecutor);

		// if there was an exception reading then initiate an immediate shutdown.
		if (tileProcessingExecutor.exception != null) {
		int tasksStillRunning = completedWorkExecutor.shutdownNow().size();
		tasksStillRunning += barcodeWriterThreads.values().stream().mapToLong(executor -> executor.shutdownNow().size()).sum();
		throw new PicardException("Reading executor had exceptions. There were " + tasksStillRunning
		+ " tasks were still running or queued and have been cancelled.", tileProcessingExecutor.exception);
		} else {
		ThreadPoolExecutorUtil.awaitThreadPoolTermination("Tile completion executor", completedWorkExecutor, Duration.ofMinutes(5));
		barcodeWriterThreads.values().forEach(ThreadPoolExecutor::shutdown);
		barcodeWriterThreads.forEach((barcode, executor) -> ThreadPoolExecutorUtil.awaitThreadPoolTermination(barcode + " writer", executor, Duration.ofMinutes(5)));
		barcodeWriterThreads.forEach((barcode, executor) -> awaitThreadPoolTermination(barcode + " writer", executor));
		}

		private void awaitThreadPoolTermination(final String executorName, final ThreadPoolExecutor executorService) {
		try {
		while (!executorService.awaitTermination(1, TimeUnit.SECONDS)) {
		log.info(String.format("%s waiting for job completion. Finished jobs - %d : Running jobs - %d : Queued jobs - %d",
		executorName, executorService.getCompletedTaskCount(), executorService.getActiveCount(),
		executorService.getQueue().size()));
		}
		} catch (final InterruptedException e) {
		e.printStackTrace();
		}
		}

		@@ -201,30 +249,16 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		}
		}

		private class Closer implements Runnable {
		private final ConvertedClusterDataWriter<CLUSTER_OUTPUT_RECORD> writer;
		private final String barcode;

		private Closer(final ConvertedClusterDataWriter<CLUSTER_OUTPUT_RECORD> writer, final String barcode) {
		this.writer = writer;
		this.barcode = barcode;
		}

		@Override
		public void run() {
		log.debug("Closing writer for barcode " + barcode);
		this.writer.close();
		}
		}

		private class TileProcessor implements Runnable {
		private final int tileNum;
		private final Map<String, SortingCollection<CLUSTER_OUTPUT_RECORD>> barcodeToRecordCollection = new HashMap<>();
		private final File barcodeFile;
		private final CompletedWorkChecker workChecker;

		TileProcessor(final int tileNum, final File barcodeFile) {
		TileProcessor(final int tileNum, final File barcodeFile, CompletedWorkChecker workChecker) {
		this.tileNum = tileNum;
		this.barcodeFile = barcodeFile;
		this.workChecker = workChecker;
		}

		@Override
		@@ -234,9 +268,13 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		while (dataProvider.hasNext()) {
		final ClusterData cluster = dataProvider.next();
		readProgressLogger.record(null, 0);
		// If this cluster is passing, or we do NOT want to ONLY emit passing reads, then add it to
		// the next
		if (cluster.isPf() \|\| includeNonPfReads) {
		final String barcode = (demultiplex ? cluster.getMatchedBarcode() : null);
		addRecord(barcode, converter.convertClusterToOutputRecord(cluster));
		}
		}

		dataProvider.close();

		@@ -247,7 +285,8 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		writerList.add(new RecordWriter(writer, value, barcode));

		});
		completedWork.put(tileNum, writerList);

		workChecker.registerCompletedWork(tileNum, writerList);

		log.info("Finished processing tile " + tileNum);
		}
		@@ -275,19 +314,21 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		final int maxRecordsInRam =
		Math.max(1, maxReadsInRamPerTile /
		barcodeRecordWriterMap.size());
		return SortingCollection.newInstanceFromPaths(
		return SortingCollection.newInstance(
		outputRecordClass,
		codecPrototype.clone(),
		outputRecordComparator,
		maxRecordsInRam,
		IOUtil.filesToPaths(tmpDirs));
		tmpDirs);
		}
		}


		private class CompletedWorkChecker implements Runnable {

		private final Map<Integer, List<RecordWriter>> completedWork = new ConcurrentHashMap<>();
		private int currentTileIndex = 0;
		private ThreadPoolExecutor tileProcessingExecutor;

		@Override
		public void run() {
		@@ -296,19 +337,30 @@ public class NewIlluminaBasecallsConverter<CLUSTER_OUTPUT_RECORD> extends Baseca
		if (completedWork.containsKey(currentTile)) {
		log.info("Writing out tile " + currentTile);
		completedWork.get(currentTile).forEach(writer -> barcodeWriterThreads.get(writer.getBarcode()).submit(writer));
		completedWork.remove(currentTile);
		currentTileIndex++;
		} else {
		try {
		Thread.sleep(5000);
		} catch (final InterruptedException e) {
		throw new PicardException(e.getMessage(), e);
		}
		if(tileProcessingExecutor.isTerminated() && completedWork.size() == 0){
		break;
		}
		}

		//we are all done scheduling work.. now schedule the closes
		barcodeRecordWriterMap.forEach((barcode, writer) -> barcodeWriterThreads.get(barcode).submit(new Closer(writer, barcode)));
		barcodeRecordWriterMap.forEach((barcode, writer) -> {
		ThreadPoolExecutorWithExceptions executorService = barcodeWriterThreads.get(barcode);
		executorService.shutdown();
		awaitThreadPoolTermination( "Barcode writer (" + barcode + ")", executorService);
		log.debug("Closing writer for barcode " + barcode);
		writer.close();
		});
		}

		void registerCompletedWork(int tileNum, List<RecordWriter> writerList) {
		completedWork.put(tileNum, writerList);
		}

		void setTileProcessingExecutor(ThreadPoolExecutor tileProcessingExecutor) {
		this.tileProcessingExecutor = tileProcessingExecutor;
		}
		}
		}