jpayne@68: package icecream;
jpayne@68: 
jpayne@68: import java.io.PrintStream;
jpayne@68: import java.util.ArrayList;
jpayne@68: import java.util.Arrays;
jpayne@68: import java.util.Random;
jpayne@68: import java.util.concurrent.atomic.AtomicLong;
jpayne@68: 
jpayne@68: import dna.AminoAcid;
jpayne@68: import fileIO.ByteFile;
jpayne@68: import fileIO.ByteStreamWriter;
jpayne@68: import fileIO.FileFormat;
jpayne@68: import fileIO.ReadWrite;
jpayne@68: import shared.Parse;
jpayne@68: import shared.Parser;
jpayne@68: import shared.PreParser;
jpayne@68: import shared.ReadStats;
jpayne@68: import shared.Shared;
jpayne@68: import shared.Timer;
jpayne@68: import shared.Tools;
jpayne@68: import stream.ConcurrentReadInputStream;
jpayne@68: import stream.ConcurrentReadOutputStream;
jpayne@68: import stream.FASTQ;
jpayne@68: import stream.FastaReadInputStream;
jpayne@68: import stream.Read;
jpayne@68: import stream.SamLine;
jpayne@68: import structures.ByteBuilder;
jpayne@68: import structures.IntList;
jpayne@68: import structures.ListNum;
jpayne@68: 
jpayne@68: /**
jpayne@68:  * Generates chimeric PacBio reads containing inverted repeats
jpayne@68:  * due to missing adapters.
jpayne@68:  * 
jpayne@68:  * @author Brian Bushnell
jpayne@68:  * @date June 8, 2019
jpayne@68:  *
jpayne@68:  */
jpayne@68: public class IceCreamMaker {
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Initialization        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Code entrance from the command line.
jpayne@68: 	 * @param args Command line arguments
jpayne@68: 	 */
jpayne@68: 	public static void main(String[] args){
jpayne@68: 		//Start a timer immediately upon code entrance.
jpayne@68: 		Timer t=new Timer();
jpayne@68: 		
jpayne@68: 		//Create an instance of this class
jpayne@68: 		IceCreamMaker x=new IceCreamMaker(args);
jpayne@68: 		
jpayne@68: 		//Run the object
jpayne@68: 		x.process(t);
jpayne@68: 		
jpayne@68: 		//Close the print stream if it was redirected
jpayne@68: 		Shared.closeStream(x.outstream);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/**
jpayne@68: 	 * Constructor.
jpayne@68: 	 * @param args Command line arguments
jpayne@68: 	 */
jpayne@68: 	public IceCreamMaker(String[] args){
jpayne@68: 		
jpayne@68: 		{//Preparse block for help, config files, and outstream
jpayne@68: 			PreParser pp=new PreParser(args, getClass(), false);
jpayne@68: 			args=pp.args;
jpayne@68: 			outstream=pp.outstream;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Set shared static variables prior to parsing
jpayne@68: 		ReadWrite.USE_PIGZ=ReadWrite.USE_UNPIGZ=true;
jpayne@68: 		ReadWrite.MAX_ZIP_THREADS=Shared.threads();
jpayne@68: 		FASTQ.TEST_INTERLEAVED=FASTQ.FORCE_INTERLEAVED=false;
jpayne@68: 		Shared.FAKE_QUAL=8;
jpayne@68: //		Shared.FASTA_WRAP=511;
jpayne@68: 		SamLine.SET_FROM_OK=true;
jpayne@68: 		
jpayne@68: 		{//Parse the arguments
jpayne@68: 			final Parser parser=parse(args);
jpayne@68: 			Parser.processQuality();
jpayne@68: 			
jpayne@68: 			maxReads=parser.maxReads;
jpayne@68: 			overwrite=ReadStats.overwrite=parser.overwrite;
jpayne@68: 			append=ReadStats.append=parser.append;
jpayne@68: 			
jpayne@68: 			in1=parser.in1;
jpayne@68: 			extin=parser.extin;
jpayne@68: 
jpayne@68: 			out1=parser.out1;
jpayne@68: 			qfout1=parser.qfout1;
jpayne@68: 			extout=parser.extout;
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		validateParams();
jpayne@68: 		fixExtensions(); //Add or remove .gz or .bz2 as needed
jpayne@68: 		checkFileExistence(); //Ensure files can be read and written
jpayne@68: 		checkStatics(); //Adjust file-related static fields as needed for this program 
jpayne@68: 		
jpayne@68: 		//Create output FileFormat objects
jpayne@68: 		ffout1=FileFormat.testOutput(out1, FileFormat.FASTQ, extout, true, overwrite, append, false);
jpayne@68: 
jpayne@68: 		//Create input FileFormat objects
jpayne@68: 		ffin1=FileFormat.testInput(in1, FileFormat.FASTQ, extin, true, true);
jpayne@68: 		
jpayne@68: 		ffIdHist=FileFormat.testOutput(outIdHist, FileFormat.TXT, null, false, overwrite, append, false);
jpayne@68: 		
jpayne@68: 		insThresh=insFraction;
jpayne@68: 		delThresh=delFraction+insThresh;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------    Initialization Helpers    ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** Parse arguments from the command line */
jpayne@68: 	private Parser parse(String[] args){
jpayne@68: 		
jpayne@68: 		//Create a parser object
jpayne@68: 		Parser parser=new Parser();
jpayne@68: 		
jpayne@68: 		//Set any necessary Parser defaults here
jpayne@68: 		//parser.foo=bar;
jpayne@68: 		
jpayne@68: 		//Parse each argument
jpayne@68: 		for(int i=0; i<args.length; i++){
jpayne@68: 			String arg=args[i];
jpayne@68: 			
jpayne@68: 			//Break arguments into their constituent parts, in the form of "a=b"
jpayne@68: 			String[] split=arg.split("=");
jpayne@68: 			String a=split[0].toLowerCase();
jpayne@68: 			String b=split.length>1 ? split[1] : null;
jpayne@68: 			if(b!=null && b.equalsIgnoreCase("null")){b=null;}
jpayne@68: 			
jpayne@68: 			if(a.equals("verbose")){
jpayne@68: 				verbose=Parse.parseBoolean(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("idhist")){
jpayne@68: 				outIdHist=b;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("minlength") || a.equals("minlen")){
jpayne@68: 				minMoleculeLength=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("maxlength") || a.equals("maxlen")){
jpayne@68: 				maxMoleculeLength=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("length") || a.equals("len")){
jpayne@68: 				minMoleculeLength=maxMoleculeLength=Parse.parseIntKMG(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("minmovie") || a.equals("minmov")){
jpayne@68: 				minMovie=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("maxmovie") || a.equals("maxmov")){
jpayne@68: 				maxMovie=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("movie") || a.equals("mov")){
jpayne@68: 				minMovie=maxMovie=Parse.parseIntKMG(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("missingrate") || a.equals("missing")){
jpayne@68: 				missingRate=Double.parseDouble(b);
jpayne@68: 				assert(missingRate<=1);
jpayne@68: 			}else if(a.equals("hiddenrate") || a.equals("hidden")){
jpayne@68: 				hiddenRate=Double.parseDouble(b);
jpayne@68: 				assert(hiddenRate<=1);
jpayne@68: 			}else if(a.equals("bothends")){
jpayne@68: 				allowBothEndsBad=Parse.parseBoolean(b);
jpayne@68: 				assert(false) : "TODO";
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("gc")){
jpayne@68: 				genomeGC=(float)Double.parseDouble(b);
jpayne@68: 				assert(genomeGC>=0 && genomeGC<=1);
jpayne@68: 			}else if(a.equals("genomesize")){
jpayne@68: 				genomeSize=Parse.parseKMG(b);
jpayne@68: 			}else if(a.equals("addns") || a.equals("ns")){
jpayne@68: 				addNs=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("seed")){
jpayne@68: 				seed=Long.parseLong(b);
jpayne@68: 			}else if(a.equals("zmws") || a.equals("maxzmws") || a.equals("reads")){
jpayne@68: 				maxZMWs=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("ccs")){
jpayne@68: 				makeCCS=Parse.parseBoolean(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("invertedrepeatrate") || a.equals("invertrepeatrate") || a.equals("irrate")){
jpayne@68: 				invertedRepeatRate=Double.parseDouble(b);
jpayne@68: 				assert(invertedRepeatRate>=0);
jpayne@68: 			}else if(a.equals("invertedrepeatminlen") || a.equals("invertrepeatminlen") || a.equals("irminlen")){
jpayne@68: 				invertedRepeatMinLength=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("invertedrepeatmaxlen") || a.equals("invertrepeatmaxlen") || a.equals("irmaxlen")){
jpayne@68: 				invertedRepeatMaxLength=Parse.parseIntKMG(b);
jpayne@68: 			}else if(a.equals("invertedrepeatlen") || a.equals("invertrepeatlen") || a.equals("irlen")){
jpayne@68: 				invertedRepeatMinLength=invertedRepeatMaxLength=Parse.parseIntKMG(b);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("miner") || a.equals("minerrorrate")){
jpayne@68: 				minErrorRate=(float)Double.parseDouble(b);
jpayne@68: 				assert(minErrorRate>=0 && minErrorRate<=1);
jpayne@68: 			}else if(a.equals("maxer") || a.equals("maxerrorrate")){
jpayne@68: 				maxErrorRate=(float)Double.parseDouble(b);
jpayne@68: 				assert(maxErrorRate>=0 && maxErrorRate<=1);
jpayne@68: 			}else if(a.equals("er") || a.equals("errorrate")){
jpayne@68: 				minErrorRate=maxErrorRate=(float)Double.parseDouble(b);
jpayne@68: 				assert(minErrorRate>=0 && minErrorRate<=1);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("minid") || a.equals("minidentity")){
jpayne@68: 				maxErrorRate=1-(float)Double.parseDouble(b);
jpayne@68: 				assert(maxErrorRate>=0 && maxErrorRate<=1);
jpayne@68: 			}else if(a.equals("maxid") || a.equals("maxidentity")){
jpayne@68: 				minErrorRate=1-(float)Double.parseDouble(b);
jpayne@68: 				assert(minErrorRate>=0 && minErrorRate<=1);
jpayne@68: 			}else if(a.equals("id") || a.equals("identity")){
jpayne@68: 				minErrorRate=maxErrorRate=1-(float)Double.parseDouble(b);
jpayne@68: 				assert(minErrorRate>=0 && minErrorRate<=1);
jpayne@68: 			}else if(a.equals("adderrors")){
jpayne@68: 				addErrors=Parse.parseBoolean(b);
jpayne@68: 			}else if(a.equals("ref")){
jpayne@68: 				parser.in1=b;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			else if(a.equals("parse_flag_goes_here")){
jpayne@68: 				long fake_variable=Parse.parseKMG(b);
jpayne@68: 				//Set a variable here
jpayne@68: 			}else if(parser.parse(arg, a, b)){//Parse standard flags in the parser
jpayne@68: 				//do nothing
jpayne@68: 			}else{
jpayne@68: 				outstream.println("Unknown parameter "+args[i]);
jpayne@68: 				assert(false) : "Unknown parameter "+args[i];
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		return parser;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Add or remove .gz or .bz2 as needed */
jpayne@68: 	private void fixExtensions(){
jpayne@68: 		in1=Tools.fixExtension(in1);
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Ensure files can be read and written */
jpayne@68: 	private void checkFileExistence(){
jpayne@68: 		//Ensure output files can be written
jpayne@68: 		if(!Tools.testOutputFiles(overwrite, append, false, out1, outIdHist)){
jpayne@68: 			throw new RuntimeException("\n\noverwrite="+overwrite+"; Can't write to output file "+out1+"\n");
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Ensure input files can be read
jpayne@68: 		if(!Tools.testInputFiles(false, true, in1)){
jpayne@68: 			throw new RuntimeException("\nCan't read some input files.\n");  
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Ensure that no file was specified multiple times
jpayne@68: 		if(!Tools.testForDuplicateFiles(true, in1, out1, outIdHist)){
jpayne@68: 			throw new RuntimeException("\nSome file names were specified multiple times.\n");
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Adjust file-related static fields as needed for this program */
jpayne@68: 	private static void checkStatics(){
jpayne@68: 		//Adjust the number of threads for input file reading
jpayne@68: 		if(!ByteFile.FORCE_MODE_BF1 && !ByteFile.FORCE_MODE_BF2 && Shared.threads()>2){
jpayne@68: 			ByteFile.FORCE_MODE_BF2=true;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		assert(FastaReadInputStream.settingsOK());
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Ensure parameter ranges are within bounds and required parameters are set */
jpayne@68: 	private boolean validateParams(){
jpayne@68: 		assert(minMoleculeLength>0 && minMoleculeLength<=maxMoleculeLength) : minMoleculeLength+", "+maxMoleculeLength;
jpayne@68: 		assert(minMovie>0 && minMovie<=maxMovie) : minMovie+", "+maxMovie;
jpayne@68: 
jpayne@68: 		assert(missingRate>=0 && missingRate<=1) : missingRate;
jpayne@68: 		assert(hiddenRate>=0 && hiddenRate<=1) : hiddenRate;
jpayne@68: 		assert(genomeGC>=0 && genomeGC<=1) : genomeGC;
jpayne@68: 		assert(in1!=null || genomeSize>maxMoleculeLength) : genomeSize;
jpayne@68: 		assert(in1!=null || invertedRepeatMaxLength*2<genomeSize) : genomeSize;
jpayne@68: 		assert(maxZMWs>0) : "zmsw="+maxZMWs+"; please set to a positive number.";
jpayne@68: 
jpayne@68: 		assert(invertedRepeatRate>=0) : invertedRepeatRate;
jpayne@68: 		assert(invertedRepeatMinLength>0 && invertedRepeatMinLength<=invertedRepeatMaxLength) : invertedRepeatMinLength+", "+invertedRepeatMaxLength;
jpayne@68: 
jpayne@68: 		assert(minErrorRate>=0 && minErrorRate<=maxErrorRate) : minErrorRate+", "+maxErrorRate;
jpayne@68: 		return true;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Outer Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/** Create read streams and process all data */
jpayne@68: 	void process(Timer t){
jpayne@68: 		
jpayne@68: 		//Turn off read validation in the input threads to increase speed
jpayne@68: 		final boolean vic=Read.VALIDATE_IN_CONSTRUCTOR;
jpayne@68: 		Read.VALIDATE_IN_CONSTRUCTOR=true;
jpayne@68: 		
jpayne@68: 		//Create a read input stream
jpayne@68: 		final ConcurrentReadInputStream cris=makeCris();
jpayne@68: 		
jpayne@68: 		//Optionally create a read output stream
jpayne@68: 		final ConcurrentReadOutputStream ros=makeCros(false);
jpayne@68: 		
jpayne@68: 		//Reset counters
jpayne@68: 		readsProcessed=readsOut=0;
jpayne@68: 		basesProcessed=basesOut=0;
jpayne@68: 		
jpayne@68: 		Random randy=Shared.threadLocalRandom(seed);
jpayne@68: 		if(cris==null){
jpayne@68: 			ref=genSynthGenome(randy);
jpayne@68: 		}else{
jpayne@68: 			ref=loadData(cris, randy);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		if(invertedRepeatRate>0){
jpayne@68: 			addInvertedRepeats(ref, randy);
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Process the reads in separate threads
jpayne@68: 		spawnThreads(cris, ros);
jpayne@68: 		
jpayne@68: 		if(verbose){outstream.println("Finished; closing streams.");}
jpayne@68: 		
jpayne@68: 		//Write anything that was accumulated by ReadStats
jpayne@68: 		errorState|=ReadStats.writeAll();
jpayne@68: 		//Close the read streams
jpayne@68: 		errorState|=ReadWrite.closeStreams(cris, ros);
jpayne@68: 		
jpayne@68: 		writeIdHist();
jpayne@68: 		
jpayne@68: 		//Reset read validation
jpayne@68: 		Read.VALIDATE_IN_CONSTRUCTOR=vic;
jpayne@68: 		
jpayne@68: 		//Report timing and results
jpayne@68: 		t.stop();
jpayne@68: 		outstream.println(Tools.timeReadsBasesProcessed(t, readsOut, basesOut, 8));
jpayne@68: //		outstream.println(Tools.readsBasesOut(readsProcessed, basesProcessed, readsOut, basesOut, 8, false));
jpayne@68: 		
jpayne@68: 		//Throw an exception of there was an error in a thread
jpayne@68: 		if(errorState){
jpayne@68: 			throw new RuntimeException(getClass().getName()+" terminated in an error state; the output may be corrupt.");
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private void writeIdHist(){
jpayne@68: 		if(ffIdHist==null){return;}
jpayne@68: 		ByteStreamWriter bsw=new ByteStreamWriter(ffIdHist);
jpayne@68: 		bsw.start();
jpayne@68: 		bsw.print("#Identity\tCount\n".getBytes());
jpayne@68: 		for(int i=0; i<idHist.length; i++){
jpayne@68: 			bsw.print(i*100.0/(ID_BINS-1), 3).print('\t').println(idHist[i]);
jpayne@68: 		}
jpayne@68: 		errorState|=bsw.poisonAndWait();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Create a Read Input Stream */
jpayne@68: 	private ConcurrentReadInputStream makeCris(){
jpayne@68: 		if(ffin1==null){return null;}
jpayne@68: 		ConcurrentReadInputStream cris=ConcurrentReadInputStream.getReadInputStream(maxReads, true, ffin1, null);
jpayne@68: 		cris.start(); //Start the stream
jpayne@68: 		if(verbose){outstream.println("Started cris");}
jpayne@68: 		return cris;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Create a Read Output Stream */
jpayne@68: 	private ConcurrentReadOutputStream makeCros(boolean pairedInput){
jpayne@68: 		if(ffout1==null){return null;}
jpayne@68: 
jpayne@68: 		//Set output buffer size
jpayne@68: 		final int buff=4;
jpayne@68: 
jpayne@68: 		final ConcurrentReadOutputStream ros=ConcurrentReadOutputStream.getStream(
jpayne@68: 				ffout1, null, qfout1, null, buff, null, false);
jpayne@68: 		ros.start(); //Start the stream
jpayne@68: 		return ros;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Spawn process threads */
jpayne@68: 	private void spawnThreads(final ConcurrentReadInputStream cris, final ConcurrentReadOutputStream ros){
jpayne@68: 		
jpayne@68: 		//Do anything necessary prior to processing
jpayne@68: 		
jpayne@68: 		//Determine how many threads may be used
jpayne@68: 		final int threads=Shared.threads();
jpayne@68: 		
jpayne@68: 		//Fill a list with ProcessThreads
jpayne@68: 		ArrayList<ProcessThread> alpt=new ArrayList<ProcessThread>(threads);
jpayne@68: 		for(int i=0; i<threads; i++){
jpayne@68: 			alpt.add(new ProcessThread(ros, i, nextZmwID, seed));
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Start the threads
jpayne@68: 		for(ProcessThread pt : alpt){
jpayne@68: 			pt.start();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Wait for threads to finish
jpayne@68: 		waitForThreads(alpt);
jpayne@68: 		
jpayne@68: 		//Do anything necessary after processing
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/** Wait until all worker threads are finished, then return */
jpayne@68: 	private void waitForThreads(ArrayList<ProcessThread> alpt){
jpayne@68: 		
jpayne@68: 		//Wait for completion of all threads
jpayne@68: 		boolean success=true;
jpayne@68: 		for(ProcessThread pt : alpt){
jpayne@68: 			
jpayne@68: 			//Wait until this thread has terminated
jpayne@68: 			while(pt.getState()!=Thread.State.TERMINATED){
jpayne@68: 				try {
jpayne@68: 					//Attempt a join operation
jpayne@68: 					pt.join();
jpayne@68: 				} catch (InterruptedException e) {
jpayne@68: 					//Potentially handle this, if it is expected to occur
jpayne@68: 					e.printStackTrace();
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			//Accumulate per-thread statistics
jpayne@68: 			readsOut+=pt.readsOutT;
jpayne@68: 			basesOut+=pt.basesOutT;
jpayne@68: 			Tools.add(idHist, pt.idHistT);
jpayne@68: 			
jpayne@68: 			success&=pt.success;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Track whether any threads failed
jpayne@68: 		if(!success){errorState=true;}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Inner Methods        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	private byte randomBase(Random randy) {
jpayne@68: 		float rGC=randy.nextFloat();
jpayne@68: 		if(rGC>=genomeGC){//AT
jpayne@68: 			return (byte)(randy.nextBoolean() ? 'A' : 'T');
jpayne@68: 		}else{//GC
jpayne@68: 			return (byte)(randy.nextBoolean() ? 'G' : 'C');
jpayne@68: 		}
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private static int randomLength(int min, int max, Random randy) {
jpayne@68: 		if(min==max){return min;}
jpayne@68: 		int range=max-min+1;
jpayne@68: 		int x=min+Tools.min(randy.nextInt(range), randy.nextInt(range));
jpayne@68: //		System.err.println(x+", "+min+", "+max);
jpayne@68: //		new Exception().printStackTrace();
jpayne@68: //		System.err.println(randy.getClass());
jpayne@68: //		System.err.println(randy.nextLong());
jpayne@68: 		return x;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private static float randomRate(float min, float max, Random randy) {
jpayne@68: 		if(min==max){return min;}
jpayne@68: 		float range=max-min;
jpayne@68: //		float a=(randy.nextFloat()+randy.nextFloat());
jpayne@68: 		float b=(randy.nextFloat()+randy.nextFloat());
jpayne@68: 		float c=(1.6f*randy.nextFloat()+0.4f*randy.nextFloat());
jpayne@68: 		float x=min+range*0.5f*Tools.min(b, c);
jpayne@68: //		assert(false) : "x="+x+", a="+a+", b="+b+", c="+c+", min="+min+", max="+max;
jpayne@68: //		System.err.println(x+", "+min+", "+max);
jpayne@68: 		return x;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private byte[] genSynthGenome(Random randy){
jpayne@68: 		assert(genomeSize<=MAX_GENOME_LENGTH) : genomeSize;
jpayne@68: 		byte[] array=new byte[(int)genomeSize];
jpayne@68: 		for(int i=0; i<genomeSize; i++){
jpayne@68: 			array[i]=randomBase(randy);
jpayne@68: 		}
jpayne@68: 		return array;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private byte[] loadData(ConcurrentReadInputStream cris, Random randy){
jpayne@68: 
jpayne@68: 		ByteBuilder bb=new ByteBuilder();
jpayne@68: 		
jpayne@68: 		//Grab the first ListNum of reads
jpayne@68: 		ListNum<Read> ln=cris.nextList();
jpayne@68: 
jpayne@68: 		//As long as there is a nonempty read list...
jpayne@68: 		while(ln!=null && ln.size()>0){
jpayne@68: //			if(verbose){outstream.println("Fetched "+reads.size()+" reads.");} //Disabled due to non-static access
jpayne@68: 			
jpayne@68: 			for(Read r : ln){
jpayne@68: 				
jpayne@68: //				System.err.println("Fetched read len="+r.length());//123
jpayne@68: 				
jpayne@68: 				//Increment counters
jpayne@68: 				readsProcessed+=r.pairCount();
jpayne@68: 				basesProcessed+=r.pairLength();
jpayne@68: 				
jpayne@68: 				
jpayne@68: 
jpayne@68: //				System.err.println("Filter: addNs="+addNs+", (r.length()>maxMoleculeLength && (invertedRepeatRate<=0 || r.length()>2*invertedRepeatMaxLength);//123
jpayne@68: 				
jpayne@68: 				//Filter disabled; it causes short sequences to be ignored.
jpayne@68: 				//Optional filter criteria
jpayne@68: //				if(!addNs || (r.length()>maxMoleculeLength && (invertedRepeatRate<=0 || r.length()>2*invertedRepeatMaxLength))){
jpayne@68: 					final byte[] bases=r.bases;
jpayne@68: 					
jpayne@68: 					if(addNs){
jpayne@68: 						if(bb.length()>0){bb.append('N');}
jpayne@68: 					}else{
jpayne@68: 						for(int i=0; i<bases.length; i++){
jpayne@68: 							if(!AminoAcid.isFullyDefined(bases[i])){
jpayne@68: 								bases[i]=randomBase(randy);
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: 					
jpayne@68: 					if(bb.length()+bases.length<=MAX_GENOME_LENGTH){
jpayne@68: 						bb.append(bases);
jpayne@68: //						System.err.println("Appended "+r.length());//123
jpayne@68: 					}else{
jpayne@68: //						System.err.println("Appending partial.");//123
jpayne@68: 						for(byte b : bases){
jpayne@68: 							bb.append(b);
jpayne@68: 							if(bb.length>=MAX_GENOME_LENGTH){
jpayne@68: 								cris.returnList(ln.id, true);
jpayne@68: //								System.err.println("Returning partial "+bb.length);//123
jpayne@68: 								return bb.toBytes();
jpayne@68: 							}
jpayne@68: 						}
jpayne@68: 					}
jpayne@68: //				}
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			//Notify the input stream that the list was used
jpayne@68: 			cris.returnList(ln);
jpayne@68: //			if(verbose){outstream.println("Returned a list.");} //Disabled due to non-static access
jpayne@68: 			
jpayne@68: 			//Fetch a new list
jpayne@68: 			ln=cris.nextList();
jpayne@68: 		}
jpayne@68: 
jpayne@68: 		//Notify the input stream that the final list was used
jpayne@68: 		if(ln!=null){
jpayne@68: 			cris.returnList(ln.id, ln.list==null || ln.list.isEmpty());
jpayne@68: 		}
jpayne@68: 
jpayne@68: //		System.err.println("Returning full "+bb.length);//123
jpayne@68: 		return bb.toBytes();
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	private void addInvertedRepeats(byte[] bases, Random randy){
jpayne@68: 		
jpayne@68: 		long added=0;
jpayne@68: 		long toAdd=(long) (bases.length*invertedRepeatRate);
jpayne@68: 		while(added<toAdd){
jpayne@68: 			int len=randomLength(invertedRepeatMinLength, invertedRepeatMaxLength, randy);
jpayne@68: 			int start=randy.nextInt(bases.length-2*len);
jpayne@68: 			int stop=start+len;
jpayne@68: 			boolean OK=true;
jpayne@68: 			for(int i=0; i<len && OK; i++){
jpayne@68: 				OK&=(bases[start+i]!='N' && bases[stop+i]!='N');
jpayne@68: 			}
jpayne@68: 			if(OK){
jpayne@68: 				for(int i=0; i<len; i++){
jpayne@68: 					byte b=bases[stop-i-1];
jpayne@68: 					bases[stop+i]=AminoAcid.baseToComplementExtended[b];
jpayne@68: 				}
jpayne@68: 				added+=2*len;
jpayne@68: //				System.err.println("added="+added+"/"+toAdd);
jpayne@68: 			}else{
jpayne@68: //				
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Inner Classes        ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** */
jpayne@68: 	private class ProcessThread extends Thread {
jpayne@68: 		
jpayne@68: 		//Constructor
jpayne@68: 		ProcessThread(final ConcurrentReadOutputStream ros_, final int tid_, 
jpayne@68: 				final AtomicLong nextZmwID_, final long seed){
jpayne@68: 			ros=ros_;
jpayne@68: 			tid=tid_;
jpayne@68: 			atomicZmwID=nextZmwID_;
jpayne@68: //			assert(false) : randy.getClass()+", "+randy.nextLong();
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		//Called by start()
jpayne@68: 		@Override
jpayne@68: 		public void run(){
jpayne@68: 			//Do anything necessary prior to processing
jpayne@68: 			randy=Shared.threadLocalRandom(seed<0 ? seed : seed+(tid+1)*tid*1000L);
jpayne@68: 			
jpayne@68: 			//Process the reads
jpayne@68: 			processInner();
jpayne@68: 			
jpayne@68: 			//Do anything necessary after processing
jpayne@68: 			
jpayne@68: 			//Indicate successful exit status
jpayne@68: 			success=true;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/** Iterate through the reads */
jpayne@68: 		void processInner(){
jpayne@68: 
jpayne@68: 			//As long as there is a nonempty read list...
jpayne@68: 			for(long generated=atomicZmwID.getAndAdd(readsPerList); generated<maxZMWs; 
jpayne@68: 					generated=atomicZmwID.getAndAdd(readsPerList)){
jpayne@68: //				if(verbose){outstream.println("Fetched "+reads.size()+" reads.");} //Disabled due to non-static access
jpayne@68: 
jpayne@68: 				long toGenerate=Tools.min(readsPerList, maxZMWs-generated);
jpayne@68: 				
jpayne@68: 				ArrayList<Read> reads=generateList((int)toGenerate, generated);
jpayne@68: 
jpayne@68: 				//Output reads to the output stream
jpayne@68: 				if(ros!=null){ros.add(reads, 0);}
jpayne@68: 			}
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/** Generate the next list of reads */
jpayne@68: 		private ArrayList<Read> generateList(int toGenerate, long nextID){
jpayne@68: 
jpayne@68: 			//Grab the actual read list from the ListNum
jpayne@68: 			final ArrayList<Read> reads=new ArrayList<Read>(toGenerate);
jpayne@68: 			
jpayne@68: 			//Loop through each read in the list
jpayne@68: 			for(int i=0; i<toGenerate; i++, nextID++){
jpayne@68: 				ArrayList<Read> zmw=generateZMW(nextID);
jpayne@68: 				if(zmw==null){
jpayne@68: 					i--;
jpayne@68: 					nextID--;
jpayne@68: 				}else{reads.addAll(zmw);}
jpayne@68: 			}
jpayne@68: 
jpayne@68: 			return reads;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		private ReadBuilder median(ArrayList<ReadBuilder> list){
jpayne@68: 			if(list.size()<3){return null;}
jpayne@68: 			IntList lengths=new IntList(list.size()-2);
jpayne@68: 			
jpayne@68: 			for(int i=1; i<list.size()-1; i++){
jpayne@68: 				lengths.add(list.get(i).length());
jpayne@68: 			}
jpayne@68: 			lengths.sort();
jpayne@68: 			int median=lengths.get((lengths.size-1)/2);
jpayne@68: 			
jpayne@68: 			for(int i=1; i<list.size()-1; i++){
jpayne@68: 				ReadBuilder rb=list.get(i);
jpayne@68: 				if(rb.length()==median){
jpayne@68: 					return rb;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			assert(false);
jpayne@68: 			return null;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/**
jpayne@68: 		 * Generate a single read.
jpayne@68: 		 */
jpayne@68: 		private ArrayList<Read> generateZMW(final long nextID){
jpayne@68: 			
jpayne@68: 			final int movieLength=randomLength(minMovie, maxMovie, randy);
jpayne@68: 			final float errorRate=randomRate(minErrorRate, maxErrorRate, randy);
jpayne@68: 			
jpayne@68: 			ArrayList<ReadBuilder> baseCalls=baseCallAllPasses(movieLength, errorRate, nextID);
jpayne@68: 			if(baseCalls==null){
jpayne@68: //				System.err.println(movieLength+", "+errorRate+", "+nextID);//123
jpayne@68: 				return null;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			final boolean missing=randy.nextFloat()<missingRate;
jpayne@68: 			if(missing){
jpayne@68: 				final int missingMod=randy.nextInt(2);
jpayne@68: 				ArrayList<ReadBuilder> temp=new ArrayList<ReadBuilder>();
jpayne@68: 				
jpayne@68: 				ReadBuilder current=null;
jpayne@68: 				for(int i=0; i<baseCalls.size(); i++){
jpayne@68: 					ReadBuilder rb=baseCalls.get(i);
jpayne@68: 					assert(rb.subreads==1);
jpayne@68: 					assert(rb.missing==0);
jpayne@68: 					assert(rb.adapters==0);
jpayne@68: 					if(current!=null && (i&1)==missingMod){
jpayne@68: 						current.add(rb);
jpayne@68: 						current.missing++;
jpayne@68: 						assert(current.subreads>1);
jpayne@68: 						assert(current.missing>0);
jpayne@68: 						temp.add(current);
jpayne@68: 						current=null;
jpayne@68: 					}else{
jpayne@68: 						if(current!=null){temp.add(current);}
jpayne@68: 						current=rb;
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 				if(current!=null){temp.add(current);}
jpayne@68: 				baseCalls=temp;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(makeCCS){
jpayne@68: 				ReadBuilder median=median(baseCalls);
jpayne@68: 				if(median==null){return null;}
jpayne@68: 				baseCalls.clear();
jpayne@68: 				baseCalls.add(median);
jpayne@68: 			}else if(hiddenRate>0){
jpayne@68: 				ArrayList<ReadBuilder> temp=new ArrayList<ReadBuilder>();
jpayne@68: 				
jpayne@68: 				ReadBuilder current=null;
jpayne@68: 				for(int i=0; i<baseCalls.size(); i++){
jpayne@68: 					ReadBuilder rb=baseCalls.get(i);
jpayne@68: 					assert(rb.adapters==0);
jpayne@68: 					if(current!=null && randy.nextFloat()<hiddenRate){
jpayne@68: 						current.add(baseCallAdapter(errorRate));
jpayne@68: 						current.add(rb);
jpayne@68: 						assert(current.adapters>0);
jpayne@68: 						assert(current.subreads>1);
jpayne@68: 					}else{
jpayne@68: 						if(current!=null){temp.add(current);}
jpayne@68: 						current=rb;
jpayne@68: 						assert(current.adapters==0);
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 				if(current!=null){temp.add(current);}
jpayne@68: 				baseCalls=temp;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			
jpayne@68: 			ArrayList<Read> reads=new ArrayList<Read>();
jpayne@68: 			for(ReadBuilder rb : baseCalls){
jpayne@68: 				Read r=rb.toRead();
jpayne@68: 				readsOutT+=r.pairCount();
jpayne@68: 				basesOutT+=r.length();
jpayne@68: 				reads.add(r);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			idHistT[(int)((1-errorRate)*(ID_BINS-1)+0.5f)]++;
jpayne@68: 			return reads;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		private ArrayList<ReadBuilder> baseCallAllPasses(final int movieLength, final float errorRate, long zmw){
jpayne@68: 			byte[] frag=null;
jpayne@68: 			
jpayne@68: 			for(int i=0; i<10 && frag==null; i++){//retry several times
jpayne@68: 				frag=fetchBases(ref);
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			if(frag==null){
jpayne@68: //				System.err.println("Failed baseCallAllPasses("+movieLength+", "+errorRate+", "+zmw);//123
jpayne@68: 				return null;
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			ArrayList<ReadBuilder> list=new ArrayList<ReadBuilder>();
jpayne@68: 			int movieRemaining=movieLength;
jpayne@68: 			int moviePos=0;
jpayne@68: 			assert(frag.length>0) : frag.length;
jpayne@68: 			int start=randy.nextInt(frag.length);
jpayne@68: 			while(movieRemaining>0){
jpayne@68: 				ReadBuilder rb=baseCallOnePass(frag, errorRate, start, moviePos, movieRemaining, zmw);
jpayne@68: 				list.add(rb);
jpayne@68: 				start=0;
jpayne@68: 				int elapsed=rb.length()+adapterLen;
jpayne@68: 				moviePos+=elapsed;
jpayne@68: 				movieRemaining-=elapsed;
jpayne@68: 				AminoAcid.reverseComplementBasesInPlace(frag);
jpayne@68: 			}
jpayne@68: 			return list;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/** Call bases for one pass */
jpayne@68: 		private ReadBuilder baseCallOnePass(final byte[] frag, final float errorRate, final int start, final int movieStart, int movieRemaining, long zmw){
jpayne@68: 			final float mult=1/(1-errorRate);
jpayne@68: 			ByteBuilder bb=new ByteBuilder();
jpayne@68: 			int fpos=start;
jpayne@68: 			for(; fpos<frag.length && movieRemaining>0; fpos++){
jpayne@68: 				float f=randy.nextFloat();
jpayne@68: 				byte b=frag[fpos];
jpayne@68: 				if(f>=errorRate){
jpayne@68: 					bb.append(b);
jpayne@68: 					movieRemaining--;
jpayne@68: 				}else{
jpayne@68: 					f=mult*(1-f);
jpayne@68: 					if(f<insThresh){//ins
jpayne@68: 						b=AminoAcid.numberToBase[randy.nextInt(4)];
jpayne@68: 						bb.append(b);
jpayne@68: 						fpos--;
jpayne@68: 						movieRemaining--;
jpayne@68: 					}else if(f<delThresh){//del
jpayne@68: 						
jpayne@68: 					}else{//sub
jpayne@68: 						int x=AminoAcid.baseToNumber[b]+randy.nextInt(3)+1;
jpayne@68: 						b=AminoAcid.numberToBase[x&3];
jpayne@68: 						bb.append(b);
jpayne@68: 						movieRemaining--;
jpayne@68: 					}
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			
jpayne@68: 			float passes=(fpos-start)*1.0f/frag.length;
jpayne@68: 			ReadBuilder rb=new ReadBuilder(bb, passes, movieStart, zmw);
jpayne@68: 			rb.errorRate=errorRate;
jpayne@68: 			return rb;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/** Call bases for an adapter sequence pass */
jpayne@68: 		private ReadBuilder baseCallAdapter(final float errorRate){
jpayne@68: 			ReadBuilder rb=baseCallOnePass(pacbioAdapter, errorRate, 0, 0, 999999, 0);
jpayne@68: 			rb.passes=0;
jpayne@68: 			rb.fullPasses=0;
jpayne@68: 			rb.subreads=0;
jpayne@68: 			rb.adapters=1;
jpayne@68: 			return rb;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		private byte[] fetchBases(byte[] source){
jpayne@68: 			
jpayne@68: 			final int len=randomLength(minMoleculeLength, maxMoleculeLength, randy);
jpayne@68: 			int start=len>=source.length ? 0 : randy.nextInt(source.length-len);
jpayne@68: 			int stop=Tools.min(start+len, source.length);
jpayne@68: 			
jpayne@68: //			System.err.println("fetchBases(len="+len+", slen="+source.length+", start="+start+", stop="+stop);//123
jpayne@68: 			
jpayne@68: 			for(int i=start; i<stop; i++){
jpayne@68: 				if(!AminoAcid.isFullyDefined(source[i])){
jpayne@68: 					return null;
jpayne@68: 				}
jpayne@68: 			}
jpayne@68: 			if(stop-start<1){return null;}
jpayne@68: 			byte[] frag=Arrays.copyOfRange(source, start, stop);
jpayne@68: 			if(randy.nextBoolean()){AminoAcid.reverseComplementBasesInPlace(frag);}
jpayne@68: 			return frag;
jpayne@68: 		}
jpayne@68: 		
jpayne@68: 		/** Number of reads retained by this thread */
jpayne@68: 		protected long readsOutT=0;
jpayne@68: 		/** Number of bases retained by this thread */
jpayne@68: 		protected long basesOutT=0;
jpayne@68: 		
jpayne@68: 		/** True only if this thread has completed successfully */
jpayne@68: 		boolean success=false;
jpayne@68: 		
jpayne@68: 		private final AtomicLong atomicZmwID;
jpayne@68: 		private final int readsPerList=Shared.bufferLen();
jpayne@68: 		
jpayne@68: 		/** Random number source */
jpayne@68: 		private Random randy;
jpayne@68: 		
jpayne@68: 		/** Random number source */
jpayne@68: 		private final long[] idHistT=new long[ID_BINS];
jpayne@68: 		
jpayne@68: 		/** Shared output stream */
jpayne@68: 		private final ConcurrentReadOutputStream ros;
jpayne@68: 		/** Thread ID */
jpayne@68: 		final int tid;
jpayne@68: 	}
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------            Fields            ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** Primary input file path */
jpayne@68: 	private String in1=null;
jpayne@68: 
jpayne@68: 	/** Primary output file path */
jpayne@68: 	private String out1=null;
jpayne@68: 
jpayne@68: 	/** Primary output file path */
jpayne@68: 	private String outIdHist=null;
jpayne@68: 
jpayne@68: 	private String qfout1=null;
jpayne@68: 	
jpayne@68: 	/** Override input file extension */
jpayne@68: 	private String extin=null;
jpayne@68: 	/** Override output file extension */
jpayne@68: 	private String extout=null;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/**  */
jpayne@68: 	private int minMoleculeLength=500;
jpayne@68: 	/**  */
jpayne@68: 	private int maxMoleculeLength=10000;
jpayne@68: 	/**  */
jpayne@68: 	private int minMovie=500;
jpayne@68: 	/**  */
jpayne@68: 	private int maxMovie=40000;
jpayne@68: 	
jpayne@68: 	/**  */
jpayne@68: 	private double missingRate=0.0;
jpayne@68: 	/**  */
jpayne@68: 	private double hiddenRate=0.0;
jpayne@68: 	/**  */
jpayne@68: 	private boolean allowBothEndsBad=false;
jpayne@68: 	
jpayne@68: 	/**  */
jpayne@68: 	private float genomeGC=0.6f;
jpayne@68: 	/**  */
jpayne@68: 	private long genomeSize=10000000;
jpayne@68: 	/**  */
jpayne@68: 	private boolean addNs=true;
jpayne@68: 	
jpayne@68: 	/**  */
jpayne@68: 	private double invertedRepeatRate=0.0;
jpayne@68: 	/**  */
jpayne@68: 	private int invertedRepeatMinLength=100;
jpayne@68: 	/**  */
jpayne@68: 	private int invertedRepeatMaxLength=10000;
jpayne@68: 	
jpayne@68: 	/**  */
jpayne@68: 	private float minErrorRate=0.05f;
jpayne@68: 	/**  */
jpayne@68: 	private float maxErrorRate=0.25f;
jpayne@68: 	/**  */
jpayne@68: 	private boolean addErrors=true;
jpayne@68: 	/** One read per ZMW */
jpayne@68: 	private boolean makeCCS=false;
jpayne@68: 	
jpayne@68: 	/** */
jpayne@68: 	private long seed=-1;
jpayne@68: 	
jpayne@68: 	private long[] idHist=new long[ID_BINS]; 
jpayne@68: 	
jpayne@68: 	//These should add to 1
jpayne@68: 	private float insFraction=0.40f;
jpayne@68: 	private float delFraction=0.35f;
jpayne@68: 	private float subFraction=0.25f;
jpayne@68: 
jpayne@68: 	private final float insThresh;
jpayne@68: 	private final float delThresh;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/** Number of reads processed */
jpayne@68: 	protected long readsProcessed=0;
jpayne@68: 	/** Number of bases processed */
jpayne@68: 	protected long basesProcessed=0;
jpayne@68: 
jpayne@68: 	/** Number of reads retained */
jpayne@68: 	protected long readsOut=0;
jpayne@68: 	/** Number of bases retained */
jpayne@68: 	protected long basesOut=0;
jpayne@68: 
jpayne@68: 	/** Quit after processing this many INPUT reads */
jpayne@68: 	private long maxReads=-1;
jpayne@68: 	
jpayne@68: 	/** Quit after generating this many OUTPUT zmws */
jpayne@68: 	private long maxZMWs=-1;
jpayne@68: 	
jpayne@68: 	/** Reference genome, max 2Gbp */
jpayne@68: 	private byte[] ref;
jpayne@68: 	
jpayne@68: 	private AtomicLong nextZmwID=new AtomicLong(0);
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------         Final Fields         ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 
jpayne@68: 	/** Primary input file */
jpayne@68: 	private final FileFormat ffin1;
jpayne@68: 	
jpayne@68: 	/** Primary output file */
jpayne@68: 	private final FileFormat ffout1;
jpayne@68: 	
jpayne@68: 	/** Primary output file */
jpayne@68: 	private final FileFormat ffIdHist;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------          Constants           ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	private static final int ID_BINS=201;
jpayne@68: 	
jpayne@68: 	private static final long MAX_GENOME_LENGTH=2000000000;
jpayne@68: 
jpayne@68: 	public static final byte[] pacbioAdapter="ATCTCTCTCAACAACAACAACGGAGGAGGAGGAAAAGAGAGAGAT".getBytes();
jpayne@68: 	public static final int adapterLen=pacbioAdapter.length;
jpayne@68: 	
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	/*----------------        Common Fields         ----------------*/
jpayne@68: 	/*--------------------------------------------------------------*/
jpayne@68: 	
jpayne@68: 	/** Print status messages to this output stream */
jpayne@68: 	private PrintStream outstream=System.err;
jpayne@68: 	/** Print verbose messages */
jpayne@68: 	public static boolean verbose=false;
jpayne@68: 	/** True if an error was encountered */
jpayne@68: 	public boolean errorState=false;
jpayne@68: 	/** Overwrite existing output files */
jpayne@68: 	private boolean overwrite=false;
jpayne@68: 	/** Append to existing output files */
jpayne@68: 	private boolean append=false;
jpayne@68: 	
jpayne@68: }