FastMap.java - cephb - CEPH library - Google Project Hosting

Source path: svn/ trunk/ src/ dbsnptools/ src/ fr/ cephb/ dbsnptools/ FastMap.java

‹r197

r204

package fr.cephb.dbsnptools;


import java.io.BufferedReader;
import java.io.ByteArrayOutputStream;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStream;
import java.io.PrintWriter;
import java.io.StringWriter;
import java.util.ArrayList;
import java.util.HashSet;
import java.util.Set;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.logging.Logger;
import java.util.regex.Pattern;

import org.xml.sax.InputSource;
import org.xml.sax.SAXException;

import fr.cephb.bio.AcidNucleicUtils;
import fr.cephb.bio.FastaMapping;
import fr.cephb.io.IOUtils;

/**
 * FastMapping of sequences (short reads) using a prefix array and perfect match
 * @author pierre
 *
 */
public class FastMap
	{
	 /** logger */
	private static final int START_FROM_INDEX_0=0;
	/** logger */
    protected static final Logger LOG= Logger.getLogger("fr.cephb");
    /** splitter for SNP context */
    private final Pattern slash=Pattern.compile("[/]");
    /** mapping chromosome /file */
	private FastaMapping mapping=new FastaMapping();
	/** input sources (file,url) containing the reads */
	private Set<String> shortReadFiles=new HashSet<String>();
	/** size of short read */
    private int flanking=60;
    /** max length of the SNP variation */
    private int max_snp_length=100;
    /** maximum number of hits of one SNP on a chromosome */
    private int upper_limit=1000;
    /** sequence of the genomic chromosome */
    private byte chrom_sequence[];
    /** prefix array */
    private int pointers[];
    /** number of threads */
    private int num_threads=2;
    /** number of short read per thread */
    private int num_reads_per_thread=10000;
    
    /**
     * A short read and the information about it (original mapping )
     *
     */
    private static class ShortRead
    	{	
    	String name=null;
    	String context=null;
    	@SuppressWarnings("unused")
		String chromosome=null;
    	Character strand=null;
    	@SuppressWarnings("unused")
    	Integer start=null;
    	@SuppressWarnings("unused")
    	Integer end=null;
    	byte[] seq5=null;
    	String seqMid=null;
    	byte[] seq3=null;
    	}
    
    /** */
    private static abstract class ShortReadsCallable
    	implements Callable<Boolean>
    	{
    	protected String chromName;
    	private List<ShortRead> array;
    	ShortReadsCallable(String chromName,List<ShortRead> array)
    		{
    		this.chromName=chromName;
    		this.array=new ArrayList<ShortRead>(array);
    		}
    	@Override
    	public Boolean call() throws Exception
    		{
    		while(!this.array.isEmpty())
    			{
    			int last=this.array.size()-1;
    			this.map(this.array.get(last));
    			this.array.remove(last);
    			}
    		return true;
    		}
    	public abstract void map(ShortRead shortRead) throws Exception;
    	}
    
    private class MapPairsCallable
	extends ShortReadsCallable
		{
		MapPairsCallable(String chromName,List<ShortRead> array)
			{
			super(chromName,array);
			}
		@Override
		public void map(ShortRead shortRead) throws Exception
			{
			mapPair(this.chromName,shortRead);
			}
		}
    
    
    /**
     * 
     * do something with the current named chromosome and
     * this input source 'f'
     *
     */
    private interface Handler
		{
		public void handle(String chromName,String f) throws IOException;
		}
    
   /** generic run */
    private void run() throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				runFastaInput(chromName, input);
				}
			});
	    }
    
    /** align pairs of short seq */
    private void runPairs() throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				runPairsOfSeq(chromName, input);
				}
			});
	    }
    
    /** align pairs of seq, more info in the input */
    private void runPairsWithContext() throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				runPairsOfSeqWithContext(chromName, input);
				}
			});
	    }
    
    /** align pairs of seq , input from XML-dbsnp */
    private void runDBSNP()  throws IOException
	    {
	    run(new Handler()
	    	{
			@Override
			public void handle(String chromName, String input) throws IOException
				{
				mapDBSNP(chromName, input);
				}
			});
	    }
    
    private void run(Handler handler) throws IOException
    	{
    	//loop over all the chromosomes file
    	for(File chromFile: this.mapping.getFiles())
    		{
    		LOG.info("reading "+chromFile);
    		//get the name for this chromosome
    		String chromName=this.mapping.getName(chromFile);
    		
    		//write genomic sequence in this dynamic byte array
    		ByteArrayOutputStream baos=new ByteArrayOutputStream(150000000);
    		//open the genomiic sequence
    		FileInputStream fin=new FileInputStream(chromFile);
    		
    		//ignore the fasta header
    		int c=fin.read();
    		if(c!='>') throw new IOException("not starting with '>' "+chromFile);
    		///skip fasta header
    		while((c=fin.read())!=-1)
    			{
    			if(c=='\n') break;
    			}
    		//read the fasta sequence
    		byte array[]=new byte[1000000];
    		int nRead;
    		while((nRead=fin.read(array))!=-1)
    			{
    			int j=0;
    			//remove blanks
    			for(int i=0;i< nRead;++i)
    				{
    				if( Character.isWhitespace(array[i])||
    					Character.isDigit(array[i])) continue;
    				array[j]=(byte)Character.toUpperCase(array[i]);
    				if(array[j]=='>') throw new IOException("expected only one sequence in  "+chromFile);
    				j++;
    				}
    			baos.write(array,0,j);
    			}
    		fin.close();
    		//get the genomic sequence as an array of bytes
    		this.chrom_sequence=baos.toByteArray();
    		baos=null;
    		LOG.info("OK in byte array. length= "+chrom_sequence.length);
    		
    		//buildt the array of pointer
    		this.pointers=new int[chrom_sequence.length];
    		for(int i=0;i+this.flanking<= this.pointers.length;++i)
    			{
    			this.pointers[i]=i;
    			}
    		
    		//remove short read containing degenerate bases
    		LOG.info("Removing degenerates");
    		int j=0;
    		for(int i=0;i< this.pointers.length;++i)
    			{
    			boolean ok=true;
    			int k=0;
    			while(k< flanking && this.pointers[i]+k< this.chrom_sequence.length)
    				{
    				byte base=this.chrom_sequence[this.pointers[i]+k];
    				if(!(AcidNucleicUtils.isATGC((char)base)))
    					{
    					ok=false;
    					break;
    					}
    				++k;
    				}
    			if(k==0 || !ok) continue;
    			this.pointers[j]=this.pointers[i];
    			j++;
    			}
    		
    		//compact the array of pointers
    		if(j!=this.pointers.length)
    			{
    			int copy[]=new int[j];
    			System.arraycopy(this.pointers, 0, copy, 0, j);
    			this.pointers=copy;
    			LOG.info("Done. Length is "+this.pointers.length);
    			}
    		
    		LOG.info("OK in int array. Now sorting "+this.pointers.length+" items");
    		//sort the suffix array
    		quicksort();
    		
    		//check sort ok
    		for(int i=1;i< this.pointers.length;++i)
    			{
    			if(compare(i-1, i, FastMap.START_FROM_INDEX_0)>0)
    				{
    				throw new RuntimeException("Not sorted at "+i);
    				}
    			}
    			
    		
    		LOG.info("Done. Sorted");
    		
    		/* echo the first bases
    		for(int i=0;i< 100 && i< this.pointers.length;i++)
    			{
    			int start=pointers[i];
    			StringBuilder x=new StringBuilder(flanking);
    			for(int k= 0;k< flanking && start+k< this.chrom_sequence.length;++k)
    				{
    				x.append((char)this.chrom_sequence[start+k]);
    				}
    			LOG.info("["+i+"] "+x);
    			} */
    		
    		for(String input:this.shortReadFiles)
    			{
    			handler.handle(chromName, input);
    			}
    		}
    	}
    
    private static class Hit
    	{
    	char strand;
    	int position;
		public String toString()
			{
			return "("+strand+") at "+position;
			}
    	}
    
    
    private void runPairsOfSeq(String chromName,String input)
    	throws IOException
	
    	{
    	LOG.info("reading short read file: "+input);
    	String line;
    	BufferedReader in=null;
    	Pattern pat=Pattern.compile("[\t]");
    	//future result, tells where are the available thread
    	Future<?> futures[]=null;
    	//thread manager
    	ExecutorService service=null;
    	//index in futures, where are we looking next
    	int index_in_futures=-1;
    	//the short-reads cuumulation. if size >= num_reads_per_thread, submit to a new thread
    	List<ShortRead> array=null;
    	//if multithread, initialize the variables
    	if(this.num_threads>1)
    		{
    		array=new ArrayList<ShortRead>(this.num_reads_per_thread);
    		futures=new Future[this.num_threads];
    		service=Executors.newFixedThreadPool(futures.length);
    		}
    	try{
	    	 in=IOUtils.openReader(input);
	    	while((line=in.readLine())!=null)
	    		{
	    		if(line.startsWith("#")) continue;
	    		String tokens[]=pat.split(line);
	    		if(tokens.length<3) continue;
	    		ShortRead shortRead=new ShortRead();
	    		shortRead.name=tokens[0];
	    		shortRead.seq5=tokens[1].toUpperCase().getBytes();
	    		shortRead.seq3=tokens[2].toUpperCase().getBytes();
	    		if(service!=null)
	    			{
	    			array.add(shortRead);
	    			
	    			if(array.size()==this.num_reads_per_thread)
	    				{
	    				//find a place for this new thread
	    				while(true)
	    					{
	    					index_in_futures++;
	    					if(index_in_futures>=futures.length)
	    						{
	    						index_in_futures=0;
	    						}
	    					if( futures[index_in_futures]!=null &&
	    						!futures[index_in_futures].isDone())
	    						{
	    						continue;
	    						}
	    					//place was found, submit new thread
	    					futures[index_in_futures]= service.submit(
	    						new MapPairsCallable(chromName, array)
	    						);
	    					break;
	    					}
	    				array.clear();
	    				}
	    			}
	    		/* monothread version */
	    		else
	    			{
	    			mapPair(chromName,shortRead);
	    			}
	    		}
	    	
	    	if(array!=null && !array.isEmpty())
	    		{
	    		new MapPairsCallable(chromName, array).call();
	    		}
    		}
    	catch(IOException err)
    		{
    		throw err;
    		}
    	catch(Exception err)
    		{
    		throw new IOException(err);
    		}
    	finally
    		{
    		if(in!=null) in.close();
    		if(service!=null) service.shutdown();
    		}
    	}
    
    private void runPairsOfSeqWithContext(String chromName,String input)
    throws IOException
	
		{
		LOG.info("reading short read file: "+input);
		String line;
		Pattern pat=Pattern.compile("[\t]");
		BufferedReader in=null;
		Future<?> futures[]=null;
    	ExecutorService service=null;
    	int index_in_futures=-1;
    	List<ShortRead> array=null;
    	if(this.num_threads>1)
    		{
    		array=new ArrayList<ShortRead>(this.num_reads_per_thread);
    		futures=new Future[this.num_threads];
    		service=Executors.newFixedThreadPool(futures.length);
    		}
		
		try
			{
			in=IOUtils.openReader(input);

			while((line=in.readLine())!=null)
				{
				if(line.startsWith("#")) continue;
				String tokens[]=pat.split(line);
				if(tokens.length<9)
					{
					System.err.println("Illegal line: "+line );
					continue;
					}
			
				ShortRead shortRead=new ShortRead();
				shortRead.name=tokens[0];
				shortRead.context=tokens[1];
				shortRead.chromosome=tokens[2];
				shortRead.strand=tokens[3].charAt(0);
				shortRead.start=Integer.parseInt(tokens[4]);
				shortRead.end=Integer.parseInt(tokens[5]);
				shortRead.seq5=tokens[6].toUpperCase().getBytes();
				shortRead.seqMid=tokens[7];
				shortRead.seq3=tokens[8].toUpperCase().getBytes();
				
				if(service!=null)
					{
	    			array.add(shortRead);
	    			
	    			if(array.size()==this.num_reads_per_thread)
	    				{
	    				//find a place for this new thread
	    				while(true)
	    					{
	    					index_in_futures++;
	    					if(index_in_futures>=futures.length)
	    						{
	    						index_in_futures=0;
	    						}
	    					if( futures[index_in_futures]!=null &&
	    						!futures[index_in_futures].isDone())
	    						{
	    						continue;
	    						}
	    					futures[index_in_futures]= service.submit(new MapPairsCallable(chromName, array));
	    					break;
	    					}
	    				array.clear();
	    				}
	    			}
				else //monothreaded
					{
					mapPair(chromName,shortRead);
					}	
				}
			
			if(array!=null && !array.isEmpty())
	    		{
	    		new MapPairsCallable(chromName, array).call();
	    		}
			}
    	catch(IOException err)
    		{
    		err.printStackTrace();
    		throw err;
    		}
    	catch(Exception err)
    		{
    		err.printStackTrace();
    		throw new IOException(err);
    		}
    	finally
    		{
    		if(in!=null) in.close();
    		if(service!=null) service.shutdown();
    		}
		}
    
    
    private void mapDBSNP(final String chromName,String input)
	throws IOException

	{
	InputStream in=IOUtils.openInputStream(input);
	SnpHandler handler=new SnpHandler()
		{
		@Override
		public void found(SnpRecord rec)
			{
			ShortRead shortRead=new ShortRead();
			shortRead.name="rs"+rec.getRsId();
			shortRead.context=rec.getObserved();
			shortRead.seq5=rec.getSeq5(flanking).toUpperCase().getBytes();
			shortRead.seq3=rec.getSeq3(flanking).toUpperCase().getBytes();
			
			mapPair(
					chromName,
					shortRead
					);
			}
		};
	try
		{
		handler.parse(new InputSource(in));
		}
	catch(SAXException err)
		{
		throw new IOException(err);
		}
	in.close();
	}
    
    private void mapPair(
    		String chromName,
    		ShortRead shortRead
    		)
    	{
    	final int max_gap=Math.max(shortRead.seqMid!=null?shortRead.seqMid.length():this.max_snp_length,this.max_snp_length);
    	byte leftSeq[]=shortRead.seq5;
		byte rightSeq[]=shortRead.seq3;
    	List<Hit> leftHits=new ArrayList<Hit>();
		List<Hit> rightHits=new ArrayList<Hit>();
		call(chromName,shortRead.name,leftSeq,leftHits);
		call(chromName,shortRead.name,rightSeq,rightHits);
		if(rightHits.isEmpty() && leftHits.isEmpty()) return;
		
		System.out.println("#c\t"+shortRead.name+"\t"+chromName+"\t"+leftHits.size()+"\t"+rightHits.size());
		if(rightHits.isEmpty() || rightHits.size()>=upper_limit) return;
		if(leftHits.isEmpty() || leftHits.size()>=upper_limit) return;
		
		StringWriter str=new StringWriter();
		PrintWriter w=new PrintWriter(str);
		
		for(Hit left:leftHits)
			{
			for(Hit right:rightHits)
				{

				if(right.strand!=left.strand) continue;
				if(right.strand=='-')
					{
					Hit tmp=right;
					right=left;
					left=tmp;
					}
				
				
				if(left.position>= right.position) continue;
				if(left.position+leftSeq.length+max_gap <  right.position) continue;
				
				
				
				w.print(
					shortRead.name+"\t"+
					chromName+"\t"+left.strand+"\t"+
					(left.position+leftSeq.length)+"\t"+
					right.position+"\t"
					);
				StringBuilder seqGenome=new StringBuilder(max_gap);
				for(int i=left.position+leftSeq.length;i<right.position;++i)
					{
					seqGenome.append((char)chrom_sequence[i]);
					}
				final String seqGenomic=seqGenome.toString();
				w.print(seqGenomic);
				w.println();
				
				
				
				
				if(shortRead.seqMid!=null &&
						shortRead.strand!=null &&
						shortRead.context!=null
						)
					{
					String orientedSeq=seqGenomic;
					
					if(right.strand=='-')
						{
						orientedSeq=AcidNucleicUtils.reverseComplement(orientedSeq);
						}
					if(!orientedSeq.equalsIgnoreCase(shortRead.seqMid))
						{
						w.println("#!d\t"+shortRead.name+"\t"+shortRead.seqMid+"\t"+seqGenomic);
						}
					
					String contexts[]=slash.split(shortRead.context);
					boolean ok=false;
					for(String ctx:contexts)
						{
						if( (ctx.equals("-") && orientedSeq.isEmpty()) ||
							ctx.equalsIgnoreCase(orientedSeq))
							{
							ok=true;
							break;
							}
						}
					if(!ok)
						{
						w.println("#!x\t"+shortRead.name+"\t"+
								shortRead.strand+"\t"+
								shortRead.context+"\t"+
								seqGenomic);
						}
						
					}
				
				}
			}
		w.flush();
		String content=str.toString();
		if(!content.isEmpty())
			{
			System.out.print(str.toString());
			}
    	}
    
    
    private void runFastaInput(String chromName,String input)
    	throws IOException
        {
        int c;
		LOG.info("reading short read file: "+input);
		BufferedReader in=IOUtils.openReader(input);
		ByteArrayOutputStream shortSeq=new ByteArrayOutputStream(flanking);
		StringBuilder shortName=new StringBuilder();
		byte dna[];
		while((c=in.read())!=-1)
			{
			if(c=='>')
				{
				dna=shortSeq.toByteArray();
				if(dna.length!=0)
					{
					call(chromName,shortName.toString(),dna,null);
					}
				
				shortName.setLength(0);
				shortSeq=new ByteArrayOutputStream(flanking);
				while((c=in.read())!=-1)
					{
					if(c=='\n') break;
					shortName.append((char)c);
					}
				}
			else if(Character.isLetter(c))
				{
				shortSeq.write((byte)Character.toUpperCase(c));
				}
			}
		dna=shortSeq.toByteArray();
		if(dna.length!=0)
			{
			shortSeq=new ByteArrayOutputStream(flanking);
			call(chromName,shortName.toString(),dna,null);
			}
		in.close();
		
        }
    
    private void call(String chromName,String shortName,byte shortSeq[],List<Hit> hits)
    	{
	//System.err.println("Calling "+(new String(shortSeq)));
    	if(shortSeq.length==0) return;
    	call(chromName, shortName, shortSeq,'+',hits);
    	String s=new String(shortSeq);
    	s=AcidNucleicUtils.reverseComplement(s);
    	shortSeq=s.getBytes();
    	call(chromName, shortName, shortSeq,'-',hits);
    	}
    
   
    

    private void call(
    	final String chromName,
    	final String shortName,
    	final byte shortSeq[],
    	final char strand,
    	final List<Hit> hits)
    	{
    	int i=lower_bound(0,this.pointers.length,shortSeq);
    	if(i>= this.pointers.length ) return;

	/*
    	for(int k=0;k< flanking && this.pointers[i]+k< chrom_sequence.length;++k)
		{
		System.err.print((char)(chrom_sequence[this.pointers[i]+k]));
		}
	System.err.println();
	*/

	while(i< this.pointers.length &&
			compare(i, shortSeq,START_FROM_INDEX_0)==0)
    		{
    		if(hits==null)
    			{
	    		System.out.println(
	    			shortName+"\t"+
	    			chromName+"\t"+
	    			strand+"\t"+
	    			this.pointers[i]+"\t"+
	    			(this.pointers[i]+shortSeq.length)
	    			);
    			}
    		else if(hits.size()>upper_limit)
    			{
    			return;
    			}
    		else
    			{
    			Hit h=new Hit();
    			h.position=this.pointers[i];
    			h.strand=strand;
    			hits.add(h);
    			}
    		i++;
    		}
    	}
    
    
    /** C+ lower_bound */
	private  int lower_bound(
				int first, int last,
	            byte shortSeq[]
	            )
	    {
		final /* rm later */ int index_in_short_read=START_FROM_INDEX_0;
	    int len = last - first;
	    while (len > 0)
            {
            int half = len / 2;
            int middle = first + half;
            
            
            if (compare(middle,shortSeq,index_in_short_read)<0)
                    {
                    first = middle + 1;
                    len -= half + 1;
                    }
            else
                    {
                    len = half;
                    }
            
            /* we don't have to compare the two short read from the beginning
             * after each while because we know that the strings start
             * with the same prefix
             */
            /*
             * check this LATER
             *
            while(
            	len>0 && 
            	first+len < this.pointers.length &&
            	(this.pointers[first+len]+index_in_short_read) < chrom_sequence.length &&
            	(chrom_sequence[this.pointers[first]+index_in_short_read] ==
            	 chrom_sequence[this.pointers[first+len]+index_in_short_read]) &&
            	index_in_short_read < this.flanking
            	)
            	{
            	index_in_short_read++;
            	}*/
            }
	    return first;
	    }
    
    
    private void quicksort()
    	{
    	quicksort(0,this.pointers.length-1);
    	}
    
    private void quicksort(final int low,final int high)
		{
		if(high<=low) return;
		int i= partition(low,high);
		quicksort(low,i-1);
		quicksort(i+1,high);
		}
	
	 private int partition(int left, int right)
		 {
	     int i = left - 1;
	     int j = right;
	     while (true)
	    	 {
	    	 int index_in_short_read= START_FROM_INDEX_0;
	         while(compare(++i, right,index_in_short_read)<0)      // find item on left to swap
	             {                               // a[right] acts as sentinel
	             //nothing
	             
	             /** TEST THIS LATER
	             while( this.pointers[i]+index_in_short_read < chrom_sequence.length &&
	                 	this.pointers[right]+index_in_short_read < chrom_sequence.length &&
	                 	(chrom_sequence[this.pointers[i]+index_in_short_read] ==
	                 	 chrom_sequence[this.pointers[right]+index_in_short_read]) &&
	                 	index_in_short_read < this.flanking
	                 	)
	                 	{
	                 	index_in_short_read++;
	                 	} */
	             
	             
	             }
	         
	         index_in_short_read= START_FROM_INDEX_0;
	         while (compare(right,--j,index_in_short_read)<0)      // find item on right to swap
	             {
	             if (j == left) break;           // don't go out-of-bounds
	             
	             /** TEST THIS LATER
	             while( this.pointers[right]+index_in_short_read < chrom_sequence.length &&
		                 	this.pointers[j]+index_in_short_read < chrom_sequence.length &&
		                 	(chrom_sequence[this.pointers[right]+index_in_short_read] ==
		                 	 chrom_sequence[this.pointers[j]+index_in_short_read]) &&
		                 	index_in_short_read < this.flanking
		                 	)
		                 	{
		                 	index_in_short_read++;
		                 	}*/
	             
	             
	             }
	         if (i >= j) break;                  // check if pointers cross
	        exchange (i, j);                      // swap two elements into place
	     	}
	     exchange(i, right);                      // swap with partition element
	     return i;
		 }
    
		private int compare(
			final int index1,
			final int index2,
			final int index_in_short_read
			)
			{
			final int base1= this.pointers[index1];
			final int base2= this.pointers[index2];
			
			int i=index_in_short_read,c1,c2;
			while(i<flanking)
				{
				c1=((base1+i)< chrom_sequence.length ?chrom_sequence[base1+i] :-1);
				c2=((base2+i)< chrom_sequence.length ?chrom_sequence[base2+i] :-1);
				if(c1==-1)
					{
					if(c2==-1) return 0;
					return -1;
					}
				else if(c2==-1)
					{
					return 1;
					}
				int diff =c1-c2;
				if(diff!=0) return diff;
				i++;
				}
			return 0;
			}
		
		private int compare(
				final int index1,
				final byte shortSeq[]    ,
				final int index_in_short_read
				)
				{
				final int genome_pos= this.pointers[index1];
				int i=index_in_short_read,c1,c2;
				while(i< shortSeq.length)
					{
					
					c1=((genome_pos+i)< chrom_sequence.length ?chrom_sequence[genome_pos+i] :-1);
					c2=(i< shortSeq.length ?shortSeq[i] :-1);
					if(c1==-1)
						{
						if(c2==-1) return 0;
						return -1;
						}
					else if(c2==-1)
						{
						return 1;
						}
					
					int diff =c1-c2;
					if(diff!=0) return diff;
					i++;
					}
				return 0;
				}
		
		
		
		private void exchange(int i, int j)
			{
			int tmp=this.pointers[i];
			this.pointers[i]=this.pointers[j];
			this.pointers[j]=tmp;
			}
		
		
	public static void main(String[] args) {
		try
			{
			FastMap app=new FastMap();
			int optind=0;
			int mode=0;
			
			while(optind< args.length)
				{
				if(args[optind].equals("-h") ||
				   args[optind].equals("-help") ||
				   args[optind].equals("--help"))
					{
					System.err.println("Options:");
					System.err.println("This program maps the flanking ends of a SNP on a chromosome.");
					System.err.println(" -h help; This screen.");
					System.err.println(" -t <int> number of thread default:"+app.num_threads);
					System.err.println(" -c <file-fasta> add this genomic chromosome fasta file");
					System.err.println(" -f <int> flanking size. default:"+app.flanking);
					System.err.println(" -L <int> max length of a SNP default:"+app.max_snp_length);
					System.err.println(" -H <int> max hit for a flanking end default:"+app.upper_limit);
					System.err.println(" -x <file> read xml mapping file");
					System.err.println(" -m");
					System.err.println("    0: read fasta");
					System.err.println("    1: read name<tab>seq5<tab>seq3");
					System.err.println("    2: read NCBI dbsnp.xml");
					System.err.println("    3: read name<tab>obs<tab>chrom<tab>strand<tab>start<tab>end<tab>seq5<tab>seqSNP<tab>seq3");
					System.err.println("<files for input>");
					return;
					}
				else if(args[optind].equals("-m"))
					{
					mode=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-t"))
					{
					app.num_threads=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-x"))
					{
					app.mapping.read(new File(args[++optind]));
					}
				else if(args[optind].equals("-f"))
					{
					app.flanking=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-L"))
					{
					app.max_snp_length=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-H"))
					{
					app.upper_limit=Integer.parseInt(args[++optind]);
					}
				else if(args[optind].equals("-c"))
					{
					File f=new File(args[++optind]);
					String name=f.getName();
					if(name.length()>3 && name.endsWith(".fa"))
						{
						name=name.substring(0, name.length()-3);
						}
					if(name.length()>3 && name.startsWith("chr"))
						{
						name=name.substring(3);
						}
					app.mapping.put(name,f);
					}
				else if(args[optind].equals("--"))
					{
					optind++;
					break;
					}
				else if(args[optind].startsWith("-"))
					{
					System.err.println("Unknown option "+args[optind]);
					return;
					}
				else 
					{
					break;
					}
				++optind;
				}
			
            if(optind==args.length)
                    {
                    System.err.println("Illegal number of arguments.");
                    return;
                    }
            else
                    {
                    while(optind< args.length)
                            {
                          	app.shortReadFiles.add(args[optind++]);
                            }
                    }
            if(app.shortReadFiles.isEmpty())
            	{
            	 System.err.println("no short read files");
            	 return;
            	}
            if(app.mapping.getFiles().isEmpty())
            	{
            	System.err.println("no genomic files");
            	return;
            	}
            switch(mode)
            	{
            	case 0:  app.run(); break;
            	case 1: app.runPairs(); break;
            	case 2: app.runDBSNP(); break;
            	case 3: app.runPairsWithContext(); break;
            	default: System.err.println("bad mode "+mode); break;
            	}
			} 
		catch(Throwable err)
			{
			err.printStackTrace();
			}
		}
	}

Show details Hide details

Change log

r198 by plindenbaum.ceph on Feb 1, 2010 Diff

cont

Go to:

Project members, sign in to write a code review

Older revisions

r197 by plindenbaum.ceph on Jan 30, 2010 Diff

cont

r196 by plindenbaum.ceph on Jan 29, 2010 Diff

cont

r195 by plindenbaum.ceph on Jan 29, 2010 Diff

cont

All revisions of this file

File info

Size: 27720 bytes, 1044 lines

View raw file