gt: 					final int strand=r.strand();
context: package jgi;
  
  import java.io.File;
  import java.io.PrintStream;
  import java.util.ArrayList;
  import java.util.Arrays;
  import java.util.BitSet;
  import java.util.HashMap;
  import java.util.Locale;
  
  import dna.ChromosomeArray;
  import dna.Data;
  import dna.Scaffold;
  import fileIO.ByteFile;
  import fileIO.ByteStreamWriter;
  import fileIO.FileFormat;
  import fileIO.ReadWrite;
  import fileIO.TextFile;
  import fileIO.TextStreamWriter;
  import shared.KillSwitch;
  import shared.Parser;
  import shared.PreParser;
  import shared.ReadStats;
  import shared.Shared;
  import shared.Timer;
  import shared.Tools;
  import stream.Read;
  import stream.SamLine;
  import stream.SamLineStreamer;
  import stream.ScaffoldCoordinates;
  import stream.SiteScore;
  import structures.ByteBuilder;
  import structures.CoverageArray;
  import structures.CoverageArray2;
  import structures.CoverageArray3;
  import structures.ListNum;
  import structures.LongList;
  
  /**
   * @author Brian Bushnell
   * @date Jan 4, 2013
   *
   */
  public class CoveragePileup {
  	
  	/*--------------------------------------------------------------*/
  	/*----------------             Main             ----------------*/
  	/*--------------------------------------------------------------*/
  	
  	public static void main(String[] args){
  		Timer t=new Timer();
  		
  		CoveragePileup x=new CoveragePileup(args);
  		
  		x.process();
  		
  		t.stop();
  		if(!KEY_VALUE){
  			x.outstream.println();
  			x.outstream.println("Time: \t"+t);
  		}
  		
  		//Close the print stream if it was redirected
  		Shared.closeStream(x.outstream);
  	}
  	
  	/*--------------------------------------------------------------*/
  	/*----------------        Initialization        ----------------*/
  	/*--------------------------------------------------------------*/
  	
  	public CoveragePileup(String[] args){
  
  		{//Preparse block for help, config files, and outstream
  			PreParser pp=new PreParser(args, printCommand ? getClass() : null, false);
  			args=pp.args;
  			outstream=pp.outstream;
  		}
  		
  		int vectorMode=-1;
  		boolean outset=false;
  		ReadWrite.USE_UNPIGZ=true;
  //		SamLine.RNAME_AS_BYTES=false;
  		
  		for(int i=0; i<args.length; i++){
  			final String arg=args[i];
  			final String[] split=arg.split("=");
  			String a=split[0].toLowerCase();
  			String b=split.length>1 ? split[1] : null;
  			
  			if(a.equals("ref") || a.equals("reference") || a.equals("fasta")){
  				reference=b;
  			}else if(a.equals("addfromref")){
  				ADD_FROM_REF=Tools.parseBoolean(b);
  			}else if(a.equals("addfromreads")){
  				ADD_FROM_READS=Tools.parseBoolean(b);
  			}else if(a.equals("in") || a.equals("in1")){
  				in1=b;
  			}else if(a.equals("in2")){
  				in2=b;
  			}else if(a.equals("out") || a.equals("coveragestats") || a.equals("covstats") || a.equals("stats")){
  				covstats=b;
  				outset=true;
  			}else if(a.equals("minscaf") || a.equals("covminscaf")){
  				minscaf=Integer.parseInt(b);
  			}else if(a.equals("ss") || a.equals("samstreamer")){
  				if(b!=null && Tools.isDigit(b.charAt(0))){
  					useStreamer=true;
  					streamerThreads=Tools.max(1, Integer.parseInt(b));
  				}else{
  					useStreamer=Tools.parseBoolean(b);
  				}
  			}else if(a.equals("minq") || a.equals("minmapq")){
  				minMapq=Integer.parseInt(b);
  			}else if(a.equals("outsam")){
  				outsam=b;
  			}else if(a.equals("rpkm") || a.equals("fpkm") || a.equals("outrpkm")){
  				outrpkm=b;
  			}else if(a.equals("outorf")){
  				outorf=b;
  			}else if(a.equals("orffasta") || a.equals("fastaorf")){
  				orffasta=b;
  			}else if(a.equals("basecov") || a.equals("outcov")){
  				basecov=b;
  			}else if(a.equals("bincov") || a.equals("outbinned")){
  				bincov=b;
  			}else if(a.equals("normcov") || a.equals("outnormalized")){
  				normcov=b;
  			}else if(a.equals("normcovo") || a.equals("outnormalizedoverall")){
  				normcovOverall=b;
  			}else if(a.equals("delta")){
  				DELTA_ONLY=Tools.parseBoolean(b);
  			}else if(a.equals("physical") || a.equals("physicalcoverage") || a.equals("physcov")){
  				PHYSICAL_COVERAGE=Tools.parseBoolean(b);
  			}else if(a.equals("tlen")){
  				USE_TLEN=Tools.parseBoolean(b);
  			}else if(a.equals("hist") || a.equals("histogram") || a.equals("covhist")){
  				histogram=b;
  			}else if(a.equals("reads")){
  				maxReads=Tools.parseKMG(b);
  			}else if(a.equals("scafs") || a.equals("scaffolds")){
  				initialScaffolds=Tools.max(128, (int)(Tools.min(Long.parseLong(b),2000000000)));
  			}else if(a.equals("binsize")){
  				binsize=Integer.parseInt(b);
  			}else if(a.equals("32bit")){
  				bits32=Tools.parseBoolean(b);
  			}else if(a.equals("bitset") || a.equals("usebitset") || a.equals("bitsets") || a.equals("usebitsets")){
  //				if(Tools.parseBoolean(b)){arrayMode=BITSET_MODE;}
  				vectorMode=Tools.parseBoolean(b) ? BITSET_MODE : NOTHING_MODE;
  			}else if(a.equals("array") || a.equals("arrays") || a.equals("usearrays")){
  				vectorMode=Tools.parseBoolean(b) ? ARRAY_MODE : NOTHING_MODE;
  			}else if(a.equals("median") || a.equals("calcmedian")){
  				if(Tools.parseBoolean(b)){
  					vectorMode=ARRAY_MODE;
  				}
  			}else if(a.startsWith("nonzero") || a.equals("nzo")){
  				NONZERO_ONLY=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set NONZERO_ONLY to "+NONZERO_ONLY);}
  			}else if(a.equals("append") || a.equals("app")){
  				append=ReadStats.append=Tools.parseBoolean(b);
  			}else if(a.equals("overwrite") || a.equals("ow")){
  				overwrite=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set overwrite to "+overwrite);}
  			}else if(a.equalsIgnoreCase("twocolumn")){
  				TWOCOLUMN=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set TWOCOLUMN to "+TWOCOLUMN);}
  			}else if(a.equalsIgnoreCase("keyvalue") || a.equalsIgnoreCase("machineout")){
  				KEY_VALUE=Tools.parseBoolean(b);
  			}else if(a.equalsIgnoreCase("countgc")){
  				COUNT_GC=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set COUNT_GC to "+COUNT_GC);}
  			}else if(a.equals("secondary") || a.equals("usesecondary")){
  				USE_SECONDARY=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set USE_SECONDARY_ALIGNMENTS to "+USE_SECONDARY);}
  			}else if(a.equals("softclip") || a.equals("includesoftclip")){
  				INCLUDE_SOFT_CLIP=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set INCLUDE_SOFT_CLIP to "+INCLUDE_SOFT_CLIP);}
  			}else if(a.equals("keepshortbins") || a.equals("ksb")){
  				KEEP_SHORT_BINS=Tools.parseBoolean(b);
  				if(verbose){outstream.println("Set KEEP_SHORT_BINS to "+KEEP_SHORT_BINS);}
  			}else if(a.equals("strandedcoverage") || a.equals("strandedcov") || a.equals("covstranded") || a.equals("stranded")){
  				STRANDED=Tools.parseBoolean(b);
  			}else if(a.equals("startcov") || a.equals("covstart") || a.equals("startonly")){
  				START_ONLY=Tools.parseBoolean(b);
  			}else if(a.equals("stopcov") || a.equals("covstop") || a.equals("stoponly")){
  				STOP_ONLY=Tools.parseBoolean(b);
  			}else if(a.equals("concise")){
  				CONCISE=Tools.parseBoolean(b);
  			}else if(a.equals("verbose")){
  				verbose=Tools.parseBoolean(b);
  			}else if(a.equals("normc") || a.equals("normalizecoverage")){
  				NORMALIZE_COVERAGE=Tools.parseBoolean(b);
  			}else if(a.equals("header") || a.equals("hdr")){
  				printHeader=Tools.parseBoolean(b);
  			}else if(a.equals("headerpound") || a.equals("#")){
  				headerPound=Tools.parseBoolean(b);
  			}else if(a.equals("stdev")){
  				calcCovStdev=Tools.parseBoolean(b);
  			}else if(a.equals("delcov") || a.equals("dels") || a.equals("includedels") || a.equals("includedeletions") || a.equals("delcoverage")){
  				INCLUDE_DELETIONS=Tools.parseBoolean(b);
  			}else if(a.equals("dupecoverage") || a.equals("dupecov") || a.equals("dupes") || a.equals("duplicates") || a.equals("includeduplicates")){
  				INCLUDE_DUPLICATES=Tools.parseBoolean(b);
  			}else if(a.equals("ignoredupes") || a.equals("ignoreduplicates")){
  				INCLUDE_DUPLICATES=!Tools.parseBoolean(b);
  			}else if(a.equals("normb") || a.equals("normalizebins")){
  				try {
  					NORMALIZE_LENGTH_BINS=Integer.parseInt(b);
  				} catch (NumberFormatException e) {
  					boolean x=Tools.parseBoolean(b);
  					NORMALIZE_LENGTH_BINS=x ? 100 : -1;
  				}
  			}else if(a.equals("covwindow")){
  				if(b==null || b.length()<1 || Character.isLetter(b.charAt(0))){
  					USE_WINDOW=Tools.parseBoolean(b);
  				}else{
  					LOW_COV_WINDOW=Integer.parseInt(b);
  					USE_WINDOW=(LOW_COV_WINDOW>0);
  				}
  			}else if(a.equals("covwindowavg")){
  				LOW_COV_DEPTH=Double.parseDouble(b);
  			}else if(a.equals("k")){
  				k=Integer.parseInt(b);
  			}else if(Parser.parseCommonStatic(arg, a, b)){
  				//do nothing
  			}else if(Parser.parseZip(arg, a, b)){
  				//do nothing
  			}else if(in1==null && arg.indexOf('=')<0 && new File(arg).exists()){
  				in1=arg;
  			}else{
  				throw new RuntimeException("Unknown parameter: "+args[i]);
  			}
  			
  		}
  		
  		ReadWrite.SAMTOOLS_IGNORE_FLAG|=ReadWrite.SAM_UNMAPPED;
  		if(!USE_SECONDARY){ReadWrite.SAMTOOLS_IGNORE_FLAG|=ReadWrite.SAM_SECONDARY;}
  		if(!INCLUDE_DUPLICATES){ReadWrite.SAMTOOLS_IGNORE_FLAG|=ReadWrite.SAM_DUPLICATE;}
  		
  		if(outsam==null){
  			SamLine.PARSE_0=false;
  			SamLine.PARSE_6=false;
  			SamLine.PARSE_7=false;
  			SamLine.PARSE_8=false;
  			if(k<1){SamLine.PARSE_10=false;}
  			SamLine.PARSE_OPTIONAL=false;
  		}
  		
  		if(vectorMode>-1){
  			USE_BITSETS=(vectorMode==BITSET_MODE);
  			USE_COVERAGE_ARRAYS=(vectorMode==ARRAY_MODE);
  		}else{
  			if(histogram==null && basecov==null && bincov==null && normcov==null &&
  					normcovOverall==null && outorf==null && !calcCovStdev){//No need for coverage array!
  				USE_COVERAGE_ARRAYS=false;
  				if(TWOCOLUMN){//No need for bitset, either!
  					USE_BITSETS=false;
  				}else{
  					USE_BITSETS=true;
  				}
  			}
  		}
  		
  		if(verbose){
  			outstream.println("Set USE_COVERAGE_ARRAYS to "+USE_COVERAGE_ARRAYS);
  			outstream.println("Set USE_BITSETS to "+USE_BITSETS);
  		}
  		
  		if(maxReads<0){maxReads=Long.MAX_VALUE;}
  		{
  			final String a=(args.length>0 ? args[0] : null);
  //			final String b=(args.length>1 ? args[1] : null);
  			if(in1==null && a!=null && a.indexOf('=')<0 && (a.startsWith("stdin") || new File(a).exists())){in1=a;}
  //			if(covstats==null && b!=null && b.indexOf('=')<0){covstats=b;}
  			if(in1==null){in1="stdin";}
  //			if(covstats==null && !outset){
  ////				out="stdout";
  ////				outstream.println("Warning: output destination not set; producing no output.  To print to standard out, set 'out=stdout'");
  //				outstream=System.err;
  //			}
  		}
  		assert(in1!=null);
  //		assert(out!=null || outset) : "Output file was not set.";
  		
  		if(STRANDED){
  			assert(basecov==null || basecov.indexOf('#')>=0) : "Output filenames must contain '#' symbol for strand-specific output.";
  			assert(bincov==null || bincov.indexOf('#')>=0) : "Output filenames must contain '#' symbol for strand-specific output.";
  			assert(normcov==null || normcov.indexOf('#')>=0) : "Output filenames must contain '#' symbol for strand-specific output.";
  			assert(normcovOverall==null || normcovOverall.indexOf('#')>=0) : "Output filenames must contain '#' symbol for strand-specific output.";
  			assert(histogram==null || histogram.indexOf('#')>=0) : "Output filenames must contain '#' symbol for strand-specific output.";
  			assert(covstats==null || covstats.indexOf('#')>=0) : "Output filenames must contain '#' symbol for strand-specific output.";
  		}
  		
  		if(!Tools.testOutputFiles(overwrite, append, false, basecov, bincov, normcov, normcovOverall, histogram, covstats, outrpkm)){
  			throw new RuntimeException("\n\noverwrite="+overwrite+"; Can't write to output files "+
  					basecov+", "+bincov+", "+normcov+", "+normcovOverall+", "+histogram+", "+covstats+", "+outrpkm+"\n");
  		}
  	}
  
  	
  	/*--------------------------------------------------------------*/
  	/*----------------       Data Structures        ----------------*/
  	/*--------------------------------------------------------------*/
  	
  	
  	/** The goal if this is to garbage-collect unnecessary objects, not really for reusing the object */
  	public void clear(){
  		list=null;
  		table=null;
  		pairTable=null;
  		
  		program=null;
  		version=null;
  		
  		in1=null;
  		covstats=null;
  		outsam=null;
  		outorf=null;
  		outrpkm=null;
  		reference=null;
  		histogram=null;
  		basecov=null;
  		bincov=null;
  		normcov=null;
  		normcovOverall=null;
  		orffasta=null;
  		
  		error=false;
  
  		refBases=0;
  		mappedBases=0;
  		mappedNonClippedBases=0;
  		mappedReads=0;
  		properPairs=0;
  		readsProcessed=0;
  		basesProcessed=0;
  		kmersProcessed=0;
  		mappedKmers=0;
  		totalCoveredBases1=0;
  		totalCoveredBases2=0;
  		scaffoldsWithCoverage1=0;
  		scaffoldsWithCoverage2=0;
  		totalScaffolds=0;
  	}
  	
  	public void createDataStructures(){
  		refBases=0;
  		mappedBases=0;
  		mappedNonClippedBases=0;
  		mappedReads=0;
  		properPairs=0;
  		readsProcessed=0;
  		basesProcessed=0;
  		kmersProcessed=0;
  		mappedKmers=0;
  		totalCoveredBases1=0;
  		totalCoveredBases2=0;
  		scaffoldsWithCoverage1=0;
  		scaffoldsWithCoverage2=0;
  		totalScaffolds=0;
  		error=false;
  		list=new ArrayList<Scaffold>(initialScaffolds);
  		table=new HashMap<String, Scaffold>(initialScaffolds);
  		
  		if(PHYSICAL_COVERAGE){
  			pairTable=new HashMap<String, Object>();
  			if(COUNT_GC){
  				COUNT_GC=false;
  				outstream.println("COUNT_GC disabled for physical coverage mode.");
  			}
  			if(USE_SECONDARY){
  				USE_SECONDARY=false;
  				outstream.println("USE_SECONDARY disabled for physical coverage mode.");
  			}
  			
  			SamLine.PARSE_0=true;
  			SamLine.PARSE_6=true;
  			SamLine.PARSE_7=true;
  			SamLine.PARSE_8=true;
  			SamLine.PARSE_10=false;
  			SamLine.PARSE_OPTIONAL=false;
  		}
  	}
  	
  	
  	/*--------------------------------------------------------------*/
  	/*----------------         Outer Methods        ----------------*/
  	/*--------------------------------------------------------------*/
  	
  	/** Read and process all input data. */
  	public void process(){
  		createDataStructures();
  
  //		System.err.println("A");
  		if(Data.SAMTOOLS() && useStreamer){ReadWrite.USE_SAMBAMBA=false;} //Disable because it takes forever to read the header
  		ByteFile tf=ByteFile.makeByteFile(in1, false);
  
  //		System.err.println("B");
  		
  		final ByteStreamWriter tsw=(outsam==null ? null : new ByteStreamWriter(outsam, overwrite, false, true));
  		if(tsw!=null){tsw.start();}
  		ReadWrite.USE_SAMBAMBA=true;
  
  //		System.err.println("C");
  		processHeader(tf, tsw);
  
  //		System.err.println("D");
  		processReference();
  //		System.err.println("E");
  		if(maxReads<0){maxReads=Long.MAX_VALUE;}
  		if(useStreamer && !FileFormat.isStdio(in1) && maxReads==Long.MAX_VALUE){
  			tf.close();
  //			System.err.println("F");
  			processViaStreamer(tsw);
  //			System.err.println("G");
  		}else{
  //			processViaByteFile(tf, line, tsw);
  			processViaByteFile(tf, tsw);
  //			System.err.println("H");
  		}
  //		System.err.println("I");
  		
  		printOutput();
  		
  		if(orffasta!=null){
  			processOrfsFasta(orffasta, outorf, table);
  		}
  		
  		if(tsw!=null){tsw.waitForFinish();}
  	}
  	
  	private void processViaStreamer(ByteStreamWriter tsw){
  		SamLineStreamer ss=new SamLineStreamer(in1, streamerThreads, false);
  		ss.start();
  		ByteBuilder bb=new ByteBuilder(33000);
  		for(ListNum<SamLine> ln=ss.nextLines(); ln!=null && ln.size()>0; ln=ss.nextLines()){
  			ArrayList<SamLine> list=(ln==null ? null : ln.list);
  			for(SamLine sl : list){
  				if(tsw!=null){
  					sl.toBytes(bb);
  					bb.nl();
  					if(bb.length>=16384){
  						tsw.print(bb);
  						bb.clear();
  					}
  				}
  				processSamLine(sl);
  			}
  		}
  		if(tsw!=null){
  			if(bb.length>0){tsw.print(bb);}
  			tsw.poison();
  		}
  	}
  	
  	private void processViaByteFile(ByteFile tf, ByteStreamWriter tsw){
  		for(byte[] line=tf.nextLine(); line!=null && readsProcessed<maxReads; line=tf.nextLine()){
  			if(tsw!=null){tsw.println(line);}
  			processSamLine(line);
  		}
  		
  		tf.close();
  		if(tsw!=null){tsw.poison();}
  	}
  	
  	
  	/*--------------------------------------------------------------*/
  	/*----------------             Setup            ----------------*/
  	/*--------------------------------------------------------------*/
  	
  	
  	/** Process all sam header lines from the tf.
  	 * Once a non-header line is encountered, return it.
  	 * If non-null, print all lines to the tsw. */
  	public void processHeader(ByteFile tf, ByteStreamWriter tsw){
  		byte[] line=null;
  		for(line=tf.nextLine(); line!=null && (line.length==0 || line[0]=='@'); line=tf.nextLine()){
  //			System.err.println(new String(line));
  			if(tsw!=null){tsw.println(line);}
  
  			if(line.length>2){
  				final byte a=line[1], b=line[2];
  
  				if(a=='S' && b=='Q'){
  					Scaffold scaf=new Scaffold(line);
  					if(COUNT_GC){scaf.basecount=new long[8];}
  					assert(!table.containsKey(scaf.name)) : "\nDuplicate scaffold name!\n"+scaf+"\n\n"+table.get(scaf.name);
  					table.put(scaf.name, scaf);
  					list.add(scaf);
  					refBases+=scaf.length;
  //					sc.obj=new CoverageArray2(table.size(), sc.length+1);
  //					outstream.println("Made scaffold "+sc.name+" of length "+sc.length);
  				}else if(a=='P' && b=='G'){
  					String[] split=new String(line).split("\t");
  					for(String s : split){
  						if(s.startsWith("PN:")){
  							if(program==null){program=Data.forceIntern(s.substring(3));}
  						}else if(s.startsWith("VN:")){
  							if(version==null){version=Data.forceIntern(s.substring(3));}
  						}
  					}
  				}else if(a=='R' && b=='G'){
  					//Do nothing
  				}else if(a=='H' && b=='D'){
  					//Do nothing
  				}else if(a=='C' && b=='O'){
  					//Do nothing
  				}else{
  					//				assert(false) : line;
  				}
  			}
  		}
  		if(line!=null){tf.pushBack(line);}
  //		return line;
  	}
  	
  	
  	public void loadScaffoldsFromIndex(int minChrom, int maxChrom){
  		
  		final int[][] lengths=Data.scaffoldLengths;
  		final int[][] locs=Data.scaffoldLocs;
  		final byte[][][] names=Data.scaffoldNames;
  		final int[] counts=new int[8];
  		for(int chrom=minChrom; chrom<=maxChrom; chrom++){
  			final ChromosomeArray ca=Data.getChromosome(chrom);
  //			assert(false) : lengths[chrom]+", "+lengths.length+", "+names[chrom]+", "+locs[chrom]+", "+(ca==null);
  			if(lengths[chrom]!=null){
  				final int[] clengths=lengths[chrom];
  				final int[] clocs=locs[chrom];
  				final byte[][] cnames=names[chrom];
  				for(int idx=0; idx<clengths.length; idx++){
  					final int length=clengths[idx];
  					final int loc=clocs[idx];
  					final String name=new String(cnames[idx]);
  					final Scaffold scaf=new Scaffold(name, length);
  					if(ca!=null){
  						scaf.gc=ca.calcGC(loc, length, counts);
  					}
  					if(COUNT_GC){scaf.basecount=new long[8];}
  					assert(!table.containsKey(scaf.name)) : "\nDuplicate scaffold name!\n"+scaf+"\n\n"+table.get(scaf.name);
  					table.put(scaf.name, scaf);
  					list.add(scaf);
  					refBases+=scaf.length;
  				}
  			}
  		}
  	}
  	
  	
  	public void processReference(){
  		if(reference==null){return;}
  
  		ByteFile bf=ByteFile.makeByteFile(reference, false);
  		Scaffold scaf=null;
  		int len=0;
  		final long[] acgtn=new long[8];
  		boolean addLen=false;
  		for(byte[] s=bf.nextLine(); s!=null; s=bf.nextLine()){
  			if(s.length>0 && s[0]=='>'){
  				if(scaf!=null){
  					scaf.length=len;
  					if(addLen){
  						refBases+=scaf.length;
  						addLen=false;
  					}
  					scaf.gc=(float)((acgtn[1]+acgtn[2])*1d/Data.max(1, acgtn[0]+acgtn[1]+acgtn[2]+acgtn[3]));
  					scaf=null;
  					len=0;
  					Arrays.fill(acgtn, 0);
  				}
  				
  				String name=new String(s, 1, s.length-1);
  				scaf=table.get(name);
  				if(ADD_FROM_REF && scaf==null){
  					scaf=new Scaffold(name, 0);
  					if(!warned){
  						outstream.println("Warning - SAM header did not include "+name+"\nAbsent scaffolds will be added; further warnings will be suppressed.");
  						warned=true;
  					}
  					if(COUNT_GC){scaf.basecount=new long[8];}
  					table.put(name, scaf);
  					list.add(scaf);
  					addLen=true;
  				}
  			}else{
  				len+=s.length;
  				for(int i=0; i<s.length; i++){
  					acgtn[charToNum[s[i]]]++;
  				}
  			}
  		}
  		if(scaf!=null){
  			scaf.length=len;
  			if(addLen){
  				refBases+=scaf.length;
  				addLen=false;
  			}
  			scaf.gc=(float)((acgtn[1]+acgtn[2])*1d/Data.max(1, acgtn[0]+acgtn[1]+acgtn[2]+acgtn[3]));
  			scaf=null;
  			len=0;
  			Arrays.fill(acgtn, 0);
  		}
  	}
  	
  	
  	public void processOrfsFasta(String fname_in, String fname_out, HashMap<String, Scaffold> map){
  		TextFile tf=new TextFile(fname_in, false);
  		assert(!fname_in.equalsIgnoreCase(fname_out));
  		TextStreamWriter tsw=new TextStreamWriter(fname_out, overwrite, false, true);
  		tsw.start();
  		
  		if(printHeader){
  			String pound=(headerPound ? "#" : "");
  			tsw.print(pound+"mappedBases="+mappedBases+"\n");
  			tsw.print(pound+"mappedNonClippedBases="+mappedNonClippedBases+"\n");
  			tsw.print(pound+"mappedReads="+mappedReads+"\n");
  			tsw.print(pound+"name\tlength\tdepthSum\tavgDepth\tavgDepth/mappedBases\tminDepth\tmaxDepth\tmedianDepth\tstdDevDepth\tfractionCovered\n");
  		}
  		
  		String line;
  		final StringBuilder sb=new StringBuilder(500);
  //		Formatter formatter=new Formatter(sb);
  		
  		while((line=tf.nextLine())!=null){
  			if(line.length()>1 && line.charAt(0)=='>'){
  				
  				String[] split=line.split(" # "); //' # ' used as delimiters
  				
  				String orfname=split[0].substring(1).trim(); //In case there are spaces around the ' # ' delimiters
  				String scafname=orfname;
  				if(scafname.contains("_")){//PRODIGAL pads _1 to the name of the first orf of a scaffold, and etc
  					int last=scafname.lastIndexOf('_');
  					boolean numeric=false;
  					for(int i=last+1; i<scafname.length(); i++){
  						if(Tools.isDigit(scafname.charAt(i))){numeric=true;}
  						else{numeric=false; break;}
  					}
  					if(numeric){scafname=scafname.substring(0, last);}
  				}
  				
  				int start=Integer.parseInt(split[1].trim());
  				int stop=Integer.parseInt(split[2].trim());
  				int strand=Integer.parseInt(split[3].trim());
  				if(strand==1){strand=Shared.PLUS;}else{strand=Shared.MINUS;}
  				Orf orf=new Orf(orfname, start, stop, (byte)strand);
  				
  				Scaffold scaf=map.get(scafname);
  //				if(scaf==null){scaf=map.get(orfname);}
  				
  //				assert(scaf!=null) : "\nCan't find scaffold for ("+orf+")\nfrom line\n"+line+"\n";
  //				assert(orf.start>=0 && orf.stop<scaf.length) : "\norf goes out of scaffold bounds.\n"+orf+"\n"+scaf+"\n";
  				
  				if(scaf==null){
  					outstream.println("Can't find scaffold for ("+orf+")\nfrom line\n"+line+"\nscafname='"+scafname+"'\norfname='"+orfname+"'");
  					if(ABORT_ON_ERROR){
  						tsw.poison();
  						throw new RuntimeException("Aborting.");
  					}
  				}else{
  					if(orf.start<0 && orf.stop>=scaf.length){
  						outstream.println("orf goes out of scaffold bounds.\n"+orf+"\n"+scaf);
  						if(ABORT_ON_ERROR){
  							tsw.poison();
  							throw new RuntimeException("Aborting.");
  						}
  					}
  					
  					CoverageArray ca=(CoverageArray)(STRANDED && strand==1 ? scaf.obj2 : scaf.obj1); //TODO:  Strand logic here depends on stranding protocol.
  					orf.readCoverageArray(ca);
  				}
  				
  //				{
  //					tsw.print(String.format(Locale.ROOT, "%s\t%d\t", args));
  //				}
  				
  				sb.append(orf.name).append('\t');
  				sb.append(orf.length()).append('\t');
  				sb.append(orf.baseDepth).append('\t');
  				sb.append(String.format(Locale.ROOT, "%.4f", orf.avgCoverage())).append('\t');
  				sb.append(orf.avgCoverage()/mappedNonClippedBases);
  
  				sb.append('\t');
  				sb.append(orf.minDepth).append('\t');
  				sb.append(orf.maxDepth).append('\t');
  				sb.append(orf.medianDepth).append('\t');
  				sb.append(String.format(Locale.ROOT, "%.4f",orf.stdevDepth)).append('\t');
  				sb.append(String.format(Locale.ROOT, "%.4f",orf.fractionCovered()));
  
  				sb.append('\n');
  				tsw.print(sb.toString());
  				sb.setLength(0);
  			}
  		}
  		
  		tsw.poison();
  		tsw.waitForFinish();
  	}
  	
  	
  	/*--------------------------------------------------------------*/
  	/*----------------         Inner Methods        ----------------*/
  	/*--------------------------------------------------------------*/
  	
  	public boolean addCoverage(final String scafName, final byte[] seq, byte[] match, final int start0, final int stop0, final int readlen, final int strand, int incrementFrags,
  			boolean properPair, SamLine sl){//sl is optional
  		Scaffold scaf=table.get(scafName);
  		if(scaf==null){
  			if(ADD_FROM_READS){
  				if(!warned){
  					outstream.println("Warning: A read was mapped to unknown reference sequence "+scafName+
  							"\nFurther warnings will be suppressed.  These scaffolds will be included, but \n"
  							+ "some coverage statistics will be inaccurate as scaffold lengths are not known.");
  					warned=true;
  				}
  				scaf=new Scaffold(scafName, 0);
  				if(COUNT_GC){scaf.basecount=new long[8];}
  				table.put(scafName, scaf);
  				list.add(scaf);
  				return addCoverage(scaf, seq, match, start0, stop0, readlen, strand, incrementFrags, properPair, sl);
  			}else if(EA){
  				KillSwitch.kill("ERROR: A read was mapped to unknown reference sequence "+scafName);
  			}else if(!warned){
  				outstream.println("Warning: Can't find "+scafName+"\nThis scaffold will not be included.  Further warnings will be suppressed.");
  				warned=true;
  				error=true;
  			}
  			return false;
  		}
  		return addCoverage(scaf, seq, match, start0, stop0, readlen, strand, incrementFrags, properPair, sl);
  	}
  	
  	public boolean addCoverage(final Scaffold scaf, final byte[] seq, byte match[], final int start0, final int stop0, final int readlen, final int strand, int incrementFrags,
  			boolean properPair, SamLine sl){//sl is optional
  		if(scaf==null){
  			assert(false) : "Adding coverage to a null Scaffold.";
  			return false;
  		}
  		
  		if(ADD_FROM_READS){
  			int x=scaf.length;
  			scaf.length=Tools.max(scaf.length, stop0+1);
  			if(scaf.length!=x){
  				refBases=refBases+scaf.length-x;
  			}
  		}
  		
  		final int start=Tools.max(start0, 0);
  		final int stop=Tools.min(stop0, scaf.length-1);
  		
  		assert(start>=0 && stop>=0) : "\nAn error was encountered when processing a read. Output will not be valid.\n"+
  			"\nscafName="+scaf.name+"\nseq="+new String(seq)+"\nstart="+start+
  			"\nstop="+stop+"\nreadlen="+readlen+"\nstrand="+strand+"\nscaf.length="+scaf.length+"\nscaf="+scaf;
  		
  		mappedBases+=readlen;
  		mappedNonClippedBases+=(stop-start+1);
  		mappedReads++;
  		if(properPair){properPairs++;}
  
  		scaf.readhits++;
  		scaf.fraghits+=incrementFrags;
  		if(strand==1){scaf.readhitsMinus++;}
  		
  		if(seq!=null && scaf.basecount!=null){
  			final long[] counts=scaf.basecount;
  			for(int i=0; i<seq.length; i++){
  				counts[charToNum[seq[i]]]++;
  			}
  		}
  		
  		if(!INCLUDE_DELETIONS && !START_ONLY && !STOP_ONLY){
  			assert(match!=null) : "Coverage excluding deletions cannot be calculated without a match string.";
  			return addCoverageIgnoringDeletions(scaf, seq, match, start, stop, readlen, strand, incrementFrags);
  		}
  		
  		final int basehits=stop-start+1;
  		scaf.basehits+=basehits;
  		
  		if(USE_COVERAGE_ARRAYS){
  			if(scaf.obj1==null){
  				scaf.obj1=(bits32 ? new CoverageArray3(table.size(), scaf.length+1) : new CoverageArray2(table.size(), scaf.length+1));
  				if(STRANDED){
  					scaf.obj2=(bits32 ? new CoverageArray3(table.size(), scaf.length+1) : new CoverageArray2(table.size(), scaf.length+1));
  				}
  			}
  			CoverageArray ca=(CoverageArray)(STRANDED && strand==1 ? scaf.obj2 : scaf.obj1);
  			if(START_ONLY){
  				ca.increment(start);
  			}else if(STOP_ONLY){
  				ca.increment(stop);
  			}else{
  				ca.incrementRange(start, stop);
  //				mappedBases+=(stop-start+1);
  //				mappedBases+=(readlen);
  //				assert(readlen==(stop-start+1)) : "readlen="+readlen+", (stop-start+1)="+(stop-start+1)+", start="+start+", stop="+stop+", start0="+start0+", stop0="+stop0+
  //				"\n"+sl+"\n";
  				//123
  			}
  		}else if(USE_BITSETS){
  			if(scaf.obj1==null){
  				scaf.obj1=new BitSet(scaf.length+1);
  				if(STRANDED){
  					scaf.obj2=new BitSet(scaf.length+1);
  				}
  			}
  			BitSet bs=(BitSet)(STRANDED && strand==1 ? scaf.obj2 : scaf.obj1);
  			if(START_ONLY){
  				bs.set(start);
  			}else if(STOP_ONLY){
  				bs.set(stop);
  			}else{
  				bs.set(start, stop+1);
  			}
  		}
  		
  		return true;
  	}
  	
  	private boolean addCoverageIgnoringDeletions(final Scaffold scaf, final byte[] seq, byte match[], final int start, final int stop, final int readlen, final int strand, int incrementFrags){
  		assert(!INCLUDE_DELETIONS && !START_ONLY && !STOP_ONLY);
  		assert(match!=null) : "Coverage excluding deletions cannot be calculated without a match string.";
  		
  		if(Read.isShortMatchString(match)){
  			match=Read.toLongMatchString(match);
  		}
  		
  		int basehits=0;
  		if(USE_COVERAGE_ARRAYS){
  			if(scaf.obj1==null){
  				scaf.obj1=(bits32 ? new CoverageArray3(table.size(), scaf.length+1) : new CoverageArray2(table.size(), scaf.length+1));
  				if(STRANDED){
  					scaf.obj2=(bits32 ? new CoverageArray3(table.size(), scaf.length+1) : new CoverageArray2(table.size(), scaf.length+1));
  				}
  			}
  			CoverageArray ca=(CoverageArray)(STRANDED && strand==1 ? scaf.obj2 : scaf.obj1);
  			for(int rpos=start, mpos=0; mpos<match.length && rpos<=stop; mpos++){
  				byte m=match[mpos];
  				if(m=='m' || m=='S' || m=='N'){
  					ca.increment(rpos, 1);
  					basehits++;
  					rpos++;
  				}else if(m=='X' || m=='Y' || m=='C' || m=='I'){
  					//do nothing
  				}else if(m=='D'){
  					rpos++;
  				}else{
  					assert(false) : "Unhandled symbol "+m;
  				}
  			}
  		}else if(USE_BITSETS){
  			if(scaf.obj1==null){
  				scaf.obj1=new BitSet(scaf.length+1);
  				if(STRANDED){
  					scaf.obj2=new BitSet(scaf.length+1);
  				}
  			}
  			BitSet bs=(BitSet)(STRANDED && strand==1 ? scaf.obj2 : scaf.obj1);
  			for(int rpos=start, mpos=0; mpos<match.length && rpos<=stop; mpos++){
  				byte m=match[mpos];
  				if(m=='m' || m=='S' || m=='N'){
  					bs.set(rpos);
  					basehits++;
  					rpos++;
  				}else if(m=='X' || m=='Y' || m=='C' || m=='I'){
  					//do nothing
  				}else if(m=='D'){
  					rpos++;
  				}else{
  					assert(false) : "Unhandled symbol "+m;
  				}
  			}
  		}
  		scaf.basehits+=basehits;
  		
  		return true;
  	}
  	
  	
  	public boolean processSamLine(byte[] line){
  		if(line==null || line.length<3){
  			return false;
  		}else if(line[0]=='@'){
  			if(!error){
  				outstream.println("Unexpected header line: "+line);
  				outstream.println("This should not cause problems, and is probably due to concatenated sam files.\n" +
  						"Supressing future unexpected header warnings.");
  				error=true;
  			}
  			
  			if(line[1]=='S' && line[2]=='Q'){
  				Scaffold scaf=new Scaffold(line);
  				if(!table.containsKey(scaf.name)){
  					if(COUNT_GC){scaf.basecount=new long[8];}
  					table.put(scaf.name, scaf);
  					list.add(scaf);
  					refBases+=scaf.length;
  				}
  			}
  		}else{
  			SamLine sl=new SamLine(line);
  			return processSamLine(sl);
  		}
  		return false;
  	}
  	
  	public boolean processSamLine(SamLine sl){
  		readsProcessed++;
  		basesProcessed+=sl.length();
  		if(sl.duplicate() && !INCLUDE_DUPLICATES){return false;}
  		final int kmers=Tools.countKmers(sl.seq, k);
  		final double cKmers=sl.qual==null ? kmers : Tools.countCorrectKmers(sl.qual, k);
  		kmersProcessed+=kmers;
  		correctKmers+=cKmers;
  		final boolean properPair=(sl.hasMate() && sl.mapped() && sl.primary() && sl.properPair());
  		if(PHYSICAL_COVERAGE && properPair){
  			SamLine mate=(SamLine)pairTable.remove(sl.qname);
  			if(mate==null){pairTable.put(sl.qname, sl);}
  			else{
  				final int start1=sl.start(INCLUDE_SOFT_CLIP, false);
  				final int stop1=sl.stop(start1, INCLUDE_SOFT_CLIP, false);
  				final int start2=mate.start(INCLUDE_SOFT_CLIP, false);
  				final int stop2=mate.stop(start2, INCLUDE_SOFT_CLIP, false);
  				final int strand=(sl.pairnum()==0 ? sl.strand() : mate.strand());
  				final int length=USE_TLEN ? sl.tlen : Tools.max(stop1, stop2)-Tools.min(start1, start2)+1;
  				mappedKmers+=kmers;
  				addCoverage(sl.rnameS(), null, null, Tools.min(start1, start2), Tools.max(stop1, stop2), length, strand, 2, sl.properPair(), sl);
  			}
  		}else if(sl.mapped() && (USE_SECONDARY || sl.primary()) && sl.mapq>=minMapq){
  			assert(sl.seq!=null || sl.cigar!=null) : "This program requires bases or a cigar string for every sam line.  Problem line:\n"+sl+"\n";
  //			assert(sl.seq!=null) : sl.toString();
  			final int length=sl.length();
  			final int start=sl.start(INCLUDE_SOFT_CLIP, false);
  			final int stop=sl.stop(start, INCLUDE_SOFT_CLIP, false);
  //			assert(false && length==stop-start+1) : length+", "+start+", "+stop+", "+(stop-start+1);
  //			assert(false) : "'"+new String(sl.rname())+"', '"+sl.rnameS()+"'";
  //			assert(false) : "'"+sl.rnameS()+"'";
  			final byte[] match=(INCLUDE_DELETIONS ? null : sl.toShortMatch(true));
  			mappedKmers+=kmers;
  			return addCoverage(sl.rnameS(), sl.seq, match, start, stop, length, sl.strand(), sl.hasMate() ? 1 : 2, sl.properPair(), sl);
  		}
  		return false;
  	}
  	
  	
  	public boolean processRead(Read r){
  		readsProcessed++;
  		final int kmers=Tools.countKmers(r.bases, k);
  		final double cKmers=r.quality==null ? kmers : Tools.countCorrectKmers(r.quality, k);
  		kmersProcessed+=kmers;
  		correctKmers+=cKmers;
  		if(r.mapped() && r.bases!=null){
  			if(USE_SECONDARY && r.sites!=null && r.sites.size()>0){
  				boolean b=false;
  				for(SiteScore ss : r.sites){
  					b=processRead(r, ss) || b;
  				}
  				return b;
  			}else{
  				final Read mate=r.mate;
  				final boolean set1=coords.set(r);
  				final boolean set2=(PHYSICAL_COVERAGE && r.paired() && r.pairnum()==0 && coords2.set(mate));
  				if(set1 && set2 && Tools.equals(coords.name, coords2.name)){
  					final int start1=coords.start;
  					final int stop1=coords.stop;
  					final int start2=coords2.start;
  					final int stop2=coords2.stop;
  