Merge pull request #373 from AdamaJava/qcoverage_blacklist

Qcoverage blacklist

Author: newellf

qcoverage/src/org/qcmg/coverage/Algorithm.java

@@ -10,6 +10,7 @@
 
 interface Algorithm {
 	public String getName();
+	public CoverageType getCoverageType();
 	public void applyTo(final SAMRecord read, Object coverageCounter);
 	public void applyTo(final SAMRecord read, Object coverageCounter, boolean fullyPopulated);
 //	public void applyTo(final SAMRecord read, final int[] perBaseCoverages);

qcoverage/src/org/qcmg/coverage/Configuration.java

@@ -13,7 +13,7 @@
 import org.qcmg.qbamfilter.query.QueryExecutor;
 
 public final class Configuration {
-	private final boolean perFeatureFlag;
+	private boolean perFeatureFlag;
 	private final int numberThreads;
 	private final String type;
 	private final String outputFileName;
@@ -47,6 +47,9 @@ public Configuration(final Options options) throws Exception {
 		} else if (type.equals("physical") || type.equals("phys")) {
 			coverageType = CoverageType.PHYSICAL;
 			algorithm = new PhysicalCoverageAlgorithm();
+		} else if (type.equals("low_readdepth")) {
+			coverageType = CoverageType.LOW_READDEPTH;
+			algorithm = new LowReadDepthAlgorithm(options.getLowReadDepthCutoff());
 		} else {
 			throw new Exception("Unknown coverage type: '" + type + "'");
 		}
@@ -63,6 +66,12 @@ public Configuration(final Options options) throws Exception {
 		inferMissingBaiFileNames();
 
 		perFeatureFlag = options.hasPerFeatureOption();
+
+		//Set to by feature for low read depth option
+		if (type.equals("low_readdepth")) {
+			perFeatureFlag = true;
+		}
+
 		if (options.hasNumberThreadsOption()) {
 			numberThreads = options.getNumberThreads()[0];
 		} else {

qcoverage/src/org/qcmg/coverage/Coverage.java

@@ -25,9 +25,7 @@
 import java.math.BigInteger;
 import java.text.DateFormat;
 import java.text.SimpleDateFormat;
-import java.util.ArrayList;
-import java.util.Calendar;
-import java.util.List;
+import java.util.*;
 
 public final class Coverage {
 	private final Options options;
@@ -38,7 +36,7 @@ public Coverage(final Options options) throws Exception {
 		this.options = options;
 		Configuration invariants = new Configuration(options);
 		jobQueue = new JobQueue(invariants);
-		saveCoverageReport();
+		saveCoverageReport(invariants.getCoverageType());
 	}
 
 	/**
@@ -200,24 +198,57 @@ private static VcfRecord convertCoverageToVCFRecord(org.qcmg.coverage.CoverageRe
 		return vcf;
 	}
 
-	private void saveCoverageReport() throws Exception {
-		final QCoverageStats stats = new QCoverageStats();
-		for (final CoverageReport report : jobQueue.getCoverageReport()) {
-			stats.getCoverageReport().add(report);
-		}
-		
-		if (options.hasVcfFlag() && options.hasPerFeatureOption()) {
-			writeVCFReport(stats);
-		}
-		
-		if (options.hasXmlFlag()) {
-			writeXMLCoverageReport(stats);
+	private void saveCoverageReport(CoverageType coverageType) throws Exception {
+
+		//save low read depth bed file
+		if (coverageType.equals(CoverageType.LOW_READDEPTH)) {
+
+			String outfile = options.getOutputFileNames()[0];
+			//rename .bed extension if present to add the mnin coverage value
+			if (outfile.endsWith(".bed")) {
+				outfile = outfile.substring(0, outfile.length() - 4);
+			}
+			outfile = String.format("%s.low_read_depth.%s.bed", outfile, options.getLowReadDepthCutoff());
+
+			final HashMap<String, List<LowReadDepthRegion>> lowReadDepthResultsFinalMap = jobQueue.getLowReadDepthResultsFinalMap();
+			LinkedHashSet<String> refNamesOrdered = jobQueue.getRefNamesOrdered();
+
+			try (final BufferedWriter out = new BufferedWriter(new FileWriter(outfile))) {
+				for (String refName : refNamesOrdered) {
+					List<LowReadDepthRegion> lowReadDepthValues = lowReadDepthResultsFinalMap.get(refName);
+					if (lowReadDepthValues == null) {
+						continue;
+					}
+					//start and end so that bed is in correct order
+					lowReadDepthValues.sort(new LowReadDepthRegionComparator(refNamesOrdered));
+						for (final LowReadDepthRegion region : lowReadDepthValues) {
+							out.write(region.toBedString() + StringUtils.RETURN);
+						}
+				}
+			}
+		} else {
+			//save coverage report
+			final QCoverageStats stats = new QCoverageStats();
+			for (final CoverageReport report : jobQueue.getCoverageReport()) {
+				stats.getCoverageReport().add(report);
+			}
+			if (options.hasVcfFlag() && options.hasPerFeatureOption()) {
+				writeVCFReport(stats);
+			}
+
+			if (options.hasXmlFlag()) {
+				writeXMLCoverageReport(stats);
+			}
+
+			if (options.hasTxtFlag()) {
+				if( options.hasPerFeatureOption()) writePerFeatureTabDelimitedCoverageReport(stats);
+				else writePerTypeTabDelimitedCoverageReport(stats);
+			}
 		}
 
-		if (options.hasTxtFlag()) {			
-			if( options.hasPerFeatureOption()) writePerFeatureTabDelimitedCoverageReport(stats);
-			else writePerTypeTabDelimitedCoverageReport(stats);
-		}
+
+
+
 
 	}
 

qcoverage/src/org/qcmg/coverage/CoverageJob.java

@@ -6,11 +6,10 @@
  */
 package org.qcmg.coverage;
 
+
 import java.io.File;
-import java.util.Arrays;
-import java.util.HashMap;
-import java.util.HashSet;
-import java.util.Iterator;
+import java.io.IOException;
+import java.util.*;
 import java.util.concurrent.atomic.AtomicLong;
 
 import htsjdk.samtools.SamReader;
@@ -29,15 +28,20 @@ class CoverageJob implements Job {
 	private final HashSet<Gff3Record> features;
 	private int[] perBaseCoverages; // Uses 0-based coordinate indexing
 	private final HashMap<String, HashMap<Integer, AtomicLong>> idToCoverageToBaseCountMap = new HashMap<String, HashMap<Integer, AtomicLong>>();
+	private final HashMap<String, List<LowReadDepthRegion>> lowReadDepthMap = new HashMap<>();
 	private final QLogger logger;
 	private final QueryExecutor filter;
 	private final boolean perFeatureFlag;
-	private final HashSet<SamReader> fileReaders = new HashSet<SamReader>();
+
+	private final HashSet<SamReader> fileReaders = new HashSet<>();
 	private final Algorithm alg;
 	private final ReadsNumberCounter counterIn;
 	private final ReadsNumberCounter counterOut;
 	private boolean fullyPopulated;
 
+
+
+
 	CoverageJob(final String refName, final int refLength, final HashMap<String, HashSet<Gff3Record>> refToFeaturesMap,
 			final HashSet<Pair<File, File>> filePairs, final QueryExecutor filter,
 			final boolean perFeatureFlag, final Algorithm algorithm, final ReadsNumberCounter counterIn,final ReadsNumberCounter counterOut) throws Exception {
@@ -72,6 +76,11 @@ synchronized public HashMap<String, HashMap<Integer, AtomicLong>> getResults() {
 		return idToCoverageToBaseCountMap;
 	}
 
+	@Override
+	synchronized public HashMap<String, List<LowReadDepthRegion>> getLowReadDepthResults() {
+		return lowReadDepthMap;
+	}
+
 	@Override
 	public String toString() {
 		return refName + " coverage";
@@ -86,7 +95,7 @@ synchronized public void run() throws Exception{
 			logger.info("performing coverage for: " + refName);
 			performCoverage();
 			logger.info("assembling results for: " + refName);
-			assembleResults();
+			assembleResultsByAlgorithm();
 			logger.debug("assembled results for: " + refName + " are: " + getResults());
 			logger.info("ending job for: " + refName);
 		} catch (Exception ex) {
@@ -97,6 +106,7 @@ synchronized public void run() throws Exception{
 
 	void constructCoverageMap() {
 		perBaseCoverages = new int[refLength]; // All elements default to zero
+
 		boolean isArrayFull = true;
 		// Initially set all values to -1 for no coverage at that coordinate
 		Arrays.fill(perBaseCoverages, -1);
@@ -122,6 +132,14 @@ void constructCoverageMap() {
 		logger.info("fully populated: " + isArrayFull);
 	}
 
+	private void assembleResultsByAlgorithm() throws IOException {
+		if (alg.getCoverageType().equals(CoverageType.LOW_READDEPTH)) {
+			assembleLowReadDepthResults();
+		} else {
+			assembleResults( );
+		}
+	}
+
 	private void performCoverage() throws Exception {
  		for (final SamReader fileReader : fileReaders) {
 
@@ -180,17 +198,56 @@ private void performCoverage() throws Exception {
 
 	private void assembleResults() {
 		for (Gff3Record feature : features) {
-			String id = null;
+			String id;
 			if (perFeatureFlag) {
 				id = feature.getRawData();
 			} else {
 				id = feature.getType();
 			}
-			HashMap<Integer, AtomicLong> covToBaseCountMap = idToCoverageToBaseCountMap.get(id);
-			if (null == covToBaseCountMap) {
-				covToBaseCountMap = new HashMap<Integer, AtomicLong>();
-				idToCoverageToBaseCountMap.put(id, covToBaseCountMap);
+            HashMap<Integer, AtomicLong> covToBaseCountMap = idToCoverageToBaseCountMap.computeIfAbsent(id, k -> new HashMap<>());
+            for (int pos = feature.getStart(); pos <= feature.getEnd(); pos++) {
+				// GFF3 format uses 1-based feature coordinates; avoid problem
+				// of GFF3 accidentally containing 0 coordinate
+				if (pos > 0 && (pos - 1) < perBaseCoverages.length) {
+					// Adjust from 1-based to 0-based indexing
+					int cov = perBaseCoverages[pos - 1];
+					if (-1 >= cov) {
+						throw new IllegalStateException(
+								"Malformed internal state. -1 coverage values are invalid. Report this bug.");
+					}
+					covToBaseCountMap.computeIfAbsent(cov, v -> new AtomicLong()).incrementAndGet();
+				}
+			}
+		}
+		// Attempt to release coverage memory by nullifying
+		perBaseCoverages = null;
+	}
+
+	private int addLowReadDepthRegionIfNeeded(int cov, int pos, int coverageLimit, int startPos, HashMap<String, List<LowReadDepthRegion>> lowRDepthMap) {
+		if (cov < coverageLimit) {
+			if (startPos == -1) {
+				startPos = pos;
+			}
+		} else {
+			//Already a low read depth position previously, but now is higher coverage, so time
+			//to create the low read depth region and reset the startPos
+			if (startPos != -1) {
+				int endPos = pos - 1;//the end is the pos - 1
+				lowRDepthMap.get(refName).add(new LowReadDepthRegion(refName, startPos, endPos, coverageLimit));
+				startPos = -1;
 			}
+		}
+		return(startPos);
+	}
+
+	private void assembleLowReadDepthResults() throws IOException {
+		for (Gff3Record feature : features) {
+            //If low read depth flag is being requested, then we need to find regions with <=8 and <=12 coverage
+			LowReadDepthAlgorithm lowRdepthAlg = (LowReadDepthAlgorithm) alg;
+            lowReadDepthMap.computeIfAbsent(refName, k -> new ArrayList<>());
+
+			int lowReadDepthStart = -1;
+
 			for (int pos = feature.getStart(); pos <= feature.getEnd(); pos++) {
 				// GFF3 format uses 1-based feature coordinates; avoid problem
 				// of GFF3 accidentally containing 0 coordinate
@@ -201,7 +258,15 @@ private void assembleResults() {
 						throw new IllegalStateException(
 								"Malformed internal state. -1 coverage values are invalid. Report this bug.");
 					}
-					covToBaseCountMap.computeIfAbsent(cov, v -> new AtomicLong()).incrementAndGet();
+
+					lowReadDepthStart = addLowReadDepthRegionIfNeeded(cov,pos, lowRdepthAlg.getReaddepthCutoff(), lowReadDepthStart, lowReadDepthMap);
+
+					//add final low read depth region if we are at the end of the feature
+					if (pos == feature.getEnd()) {
+						if (lowReadDepthStart != -1) {
+							lowReadDepthMap.get(refName).add(new LowReadDepthRegion(refName, lowReadDepthStart, pos, lowRdepthAlg.getReaddepthCutoff()));
+						}
+					}
 				}
 			}
 		}

qcoverage/src/org/qcmg/coverage/CoverageType.java

@@ -6,7 +6,10 @@
 public enum CoverageType {
 
     SEQUENCE("sequence"),
-    PHYSICAL("physical");
+    PHYSICAL("physical"),
+    LOW_READDEPTH("low_readdepth");
+
+
     private final String value;
 
     CoverageType(String v) {

qcoverage/src/org/qcmg/coverage/Job.java

@@ -4,10 +4,12 @@
 package org.qcmg.coverage;
 
 import java.util.HashMap;
+import java.util.List;
 import java.util.concurrent.atomic.AtomicLong;
 
 interface Job {
 	public HashMap<String, HashMap<Integer, AtomicLong>> getResults();
+	public HashMap<String, List<LowReadDepthRegion>> getLowReadDepthResults();
 	public void run() throws Exception;
 	public String toString();
 }