SamTranscriptomeParser & DefinedRegionDifferentialSeq

* Replaced the IH tag with NH tag for reporting the number of alignments present per read to be compatible with DEXSeq's HTSeq app.

SamReadDepthSubSampler
* New app for reducing extreme read depths in amplicon based sequencing datasets

More mods to the SpliceJunction annotator for zscore estimations

Source/edu/utah/seq/analysis/BisStat.java

@@ -319,6 +319,10 @@ else if (plusPD != null && minusPD != null) {
 
 	/**Only call when one of the two stranded datasets is missing.*/
 	private void processStrand(PointData[] pd, PointData nonCon, PointData con, String strand){	
+			if (pd == null || pd[0] == null){
+				System.out.println(chromosome+"+\tskipping, too few observations");
+				return;
+			}
 			//save fractions
 			baseFractionMethylation.add(pd[0].getScores());
 

Source/edu/utah/seq/analysis/DefinedRegionDifferentialSeq.java

@@ -45,7 +45,7 @@ public class DefinedRegionDifferentialSeq {
 	private int minimumSpliceCounts = 10;
 	private boolean performStrandedAnalysis = false;
 	private boolean performReverseStrandedAnalysis = false;
-	private int maxRepeats = 0;
+	private int maxNumAlignments = 1;
 	private boolean verbose = true;
 	private int maxAlignmentsDepth = 50000;
 	private boolean secondStrandFlipped = false;
@@ -627,11 +627,11 @@ private boolean alignmentFails(SAMRecord sam){
 		//aligned?
 		if (sam.getReadUnmappedFlag()) return true;
 		//limit to max matches?
-		if (maxRepeats !=0){
-			Object o = sam.getAttribute("IH");
+		if (maxNumAlignments !=0){
+			Object o = sam.getAttribute("NH");
 			if (o != null)  {
-				int numRepeats = (Integer)o;
-				if (numRepeats > maxRepeats) return true;
+				int num = (Integer)o;
+				if (num > maxNumAlignments) return true;
 			}
 		}
 
@@ -1854,7 +1854,7 @@ public void processArgs(String[] args){
 					case 'l': minLog2Ratio = Float.parseFloat(args[++i]); break;
 					case 'e': minimumCounts = Integer.parseInt(args[++i]); break;
 					case 'm': removeOverlappingRegions = true; break;
-					case 'n': maxRepeats = Integer.parseInt(args[++i]); break;
+					case 'n': maxNumAlignments = Integer.parseInt(args[++i]); break;
 					case 'x': maxAlignmentsDepth = Integer.parseInt(args[++i]); break;
 					case 'p': performStrandedAnalysis = true; break;
 					case 'i': scoreIntrons = true; break;
@@ -1975,8 +1975,8 @@ public static void printDocs(){
 				"-m Mask overlapping gene annotations, recommended for well annotated genomes.\n"+
 				"-x Max per base alignment depth, defaults to 50000. Genes containing such high\n"+
 				"       density coverage are ignored.\n"+
-				"-n Max number repeat alignments. Defaults to all.  Assumes 'IH' tags have been set by\n" +
-				"       processing raw alignments with the SamTranscriptomeProcessor.\n"+
+				"-n Max number alignments per read. Defaults to 1, unique.  Assumes 'NH' tags have\n"+
+				"      been set by processing raw alignments with the SamTranscriptomeProcessor.\n"+
 				"-e Minimum number alignments per gene-region per replica, defaults to 20.\n"+
 				"-i Score introns instead of exons.\n"+
 				"-p Perform a stranded analysis. Only collect reads from the same strand as the\n" +

Source/edu/utah/seq/base/NMer.java

@@ -192,7 +192,7 @@ public String toString(){
 			return sb.toString();
 		}
 
-		public long getNumberObservations(){
+		public double getNumberObservations(){
 			return sdA[0].getNumberObservations();
 		}
 

Source/edu/utah/seq/data/Sam2USeq.java

@@ -217,7 +217,7 @@ public void splitSamBamFiles(){
 					//add number of repeats to cigar string for subsequent repeat scaling?
 					double forScaling = 1;
 					if (scaleRepeats){
-						Object o = sam.getAttribute("IH");
+						Object o = sam.getAttribute("NH");
 						if (o != null)  {
 							int numRepeats = (Integer)o;
 							cigar = cigar + numRepeats;