#175209 - Pastie

Theme:
package berkeleydb;



import java.io.BufferedReader;
import java.io.File;

import java.io.IOException;
import java.io.InputStreamReader;
import java.net.URL;

import java.util.HashMap;
import java.util.regex.Pattern;
import java.util.zip.GZIPInputStream;



import com.sleepycat.bind.tuple.TupleBinding;
import com.sleepycat.bind.tuple.TupleInput;
import com.sleepycat.bind.tuple.TupleOutput;
import com.sleepycat.je.Cursor;
import com.sleepycat.je.Database;
import com.sleepycat.je.DatabaseConfig;
import com.sleepycat.je.DatabaseEntry;
import com.sleepycat.je.DatabaseException;
import com.sleepycat.je.Environment;
import com.sleepycat.je.EnvironmentConfig;
import com.sleepycat.je.EnvironmentMutableConfig;
import com.sleepycat.je.LockMode;
import com.sleepycat.je.OperationStatus;
import com.sleepycat.je.SecondaryConfig;
import com.sleepycat.je.SecondaryCursor;
import com.sleepycat.je.SecondaryDatabase;
import com.sleepycat.je.SecondaryKeyCreator;
import com.sleepycat.je.StatsConfig;
import com.sleepycat.je.Transaction;



public class Test01
	{
	private Environment env=null;
	private Database markerDB=null;
	private Database individualDB=null;
	private Database genotypeDB=null;
	private SecondaryDatabase individualIdDB=null;
	private SecondaryDatabase positionDB=null;
	private Database databases[]=new Database[0];
	
	
	//http://cardiogenomics.med.harvard.edu/groups/comp08/pages/CEPH_pedigrees.html

	private static class Position
		{
		String chromosome;
		int position;
		Position(String chromosome,int position)
			{
			this.chromosome=chromosome;
			this.position=position;
			}
		}
	
	
	private  static class Genotype
		{
		char a1,a2;
		Genotype(char a1,char a2)
			{
			this.a1=a1;
			this.a2=a2;
			}
		
		public boolean same(char a1,char a2)
			{
			return (this.a1==a1 && this.a2==a2) ||
				   (this.a1==a2 && this.a2==a1)
				   ;
			}
		public boolean same(Genotype g)
			{
			return same(g.a1,g.a2);
			}
		
		@Override
		public boolean equals(Object obj) {
			if(obj==this) return true;
			if(obj==null || !(obj instanceof Genotype )) return false;
			return same(Genotype.class.cast(obj));
			}
		
		@Override
		public String toString() {
			return ""+a1+a2;
			}
		}
	
	
	
	private  static class Marker
		{
		int rsId;
		String alleles;
		Position pos;
		char strand;
		Marker(int rsId)
			{
			this.rsId=rsId;
			}
		public int getRSId()
			{
			return rsId;
			}
		}
	
	private  static class Individual
		{
		/**Coriell Repository Number*/
		String name;
		/** id in corriel */
		String individualID=null;
		/** father id */
		String fatherID=null;
		/** mother id */
		String motherID=null;
		
		@Override
		public String toString() {
			return name+" "+individualID+" "+fatherID+" "+motherID;
			}
		}
	
	private static class MarkerIndividual
		{
		private int rsId;
		/**Coriell Repository Number*/
		private String individualId;
		public MarkerIndividual( int rsId,String individualId)
			{
			this.rsId=rsId;
			this.individualId=individualId.toUpperCase().trim();
			}
		int getRS() { return this.rsId;}
		String getIndividual() {return this.individualId;}
		}
	
	private TupleBinding markerIndividualTupleBinding=new TupleBinding()
		{
		@Override
		public Object entryToObject(TupleInput input)
			{
			int markerId= input.readInt();
			String individualId= input.readString();
			return new MarkerIndividual(markerId,individualId);
			}
		
		@Override
		public void objectToEntry(Object object, TupleOutput output) {
			MarkerIndividual mi= MarkerIndividual.class.cast(object);
			output.writeInt(mi.getRS());
			output.writeString(mi.getIndividual());
			}
		};
		
		
		private TupleBinding positionTupleBinding = new TupleBinding()
			{
			@Override
			public Object entryToObject(TupleInput input)
				{
				String chrom= input.readString();
				int position= input.readInt();
				return new Position(chrom,position);
				}
			
			@Override
			public void objectToEntry(Object object, TupleOutput output) {
				Position p= Position.class.cast(object);
				output.writeString(p.chromosome);
				output.writeInt(p.position);
				}
			};
		
		
	private TupleBinding individualTupleBinding=new TupleBinding()
		{
		@Override
		public Object entryToObject(TupleInput input)
			{
			Individual indi= new Individual();
			indi.name= input.readString();
			indi.individualID= input.readString();
			indi.fatherID=input.readString();
			indi.motherID= input.readString();
			return indi;
			}
		
		@Override
		public void objectToEntry(Object object, TupleOutput output) {
			Individual indi= Individual.class.cast(object);
			output.writeString(indi.name);
			output.writeString(indi.individualID);
			output.writeString(indi.fatherID);
			output.writeString(indi.motherID);
			}
		};	
	
		public class IndividualCorrielNameCreator implements SecondaryKeyCreator
			{
		    public IndividualCorrielNameCreator()
		    	{

		    	}
	
		    public boolean createSecondaryKey(SecondaryDatabase secDb,
	              DatabaseEntry keyEntry, 
	              DatabaseEntry dataEntry,
	              DatabaseEntry resultEntry) {
	
		        try {
		        	Individual indi= (Individual)Test01.this.individualTupleBinding.entryToObject(dataEntry);
		            resultEntry.setData(indi.individualID==null?null:indi.individualID.getBytes("UTF-8"));
		        } catch (IOException willNeverOccur) {}
		        return true;
		    }
		} 	
		
		public class MarkerPositionCreator implements SecondaryKeyCreator
			{
		    public MarkerPositionCreator()
		    	{

		    	}
	
		    public boolean createSecondaryKey(SecondaryDatabase secDb,
	              DatabaseEntry keyEntry, 
	              DatabaseEntry dataEntry,
	              DatabaseEntry resultEntry)
		    		{
		        	Marker marker= (Marker)Test01.this.markerTupleBinding.entryToObject(dataEntry);
		           Test01.this.positionTupleBinding.objectToEntry(marker.pos, resultEntry);
		           return true;
		    	   }
		} 	
		
		
		
	private TupleBinding markerTupleBinding=new TupleBinding()
		{
		@Override
		public Object entryToObject(TupleInput input)
			{
			int rsId= input.readInt();
			Marker m= new Marker(rsId);
			m.alleles= input.readString();
			String chrom=input.readString();
			int position= input.readInt();
			m.pos= new Position(chrom,position);
			m.strand= input.readChar();
			return m;
			}
		
		@Override
		public void objectToEntry(Object object, TupleOutput output) {
			Marker m= Marker.class.cast(object);
			output.writeInt(m.getRSId());
			output.writeString(m.alleles);
			output.writeString(m.pos.chromosome);
			output.writeInt(m.pos.position);
			output.writeChar(m.strand);
			}
		};	
	
	private TupleBinding genotypeTupleBinding =new TupleBinding()
		{
		@Override
		public Object entryToObject(TupleInput input)
			{
			char a1= input.readChar();
			char a2= input.readChar();
			return new Genotype(a1,a2);
			}
		
		@Override
		public void objectToEntry(Object object, TupleOutput output) {
			Genotype g= Genotype.class.cast(object);
			output.writeChar(g.a1);
			output.writeChar(g.a2);
			}
		};		
		
		
	private Test01()
		{
		
		}
	
	void open(File file, boolean readOnly) throws DatabaseException
		{
		if(this.env!=null) close();
		EnvironmentConfig envConf= new EnvironmentConfig();
		envConf.setAllowCreate(!readOnly);
		envConf.setReadOnly(readOnly);
		envConf.setTransactional(false);
		this.env= new Environment(
				file,
			envConf
			);
		DatabaseConfig dbConfig= new DatabaseConfig();
		/* allow create if database does not exists */
		dbConfig.setAllowCreate(!readOnly);
		dbConfig.setReadOnly(readOnly);
		this.markerDB= this.env.openDatabase(null, "marker", dbConfig);
		this.individualDB= this.env.openDatabase(null, "individual", dbConfig);
		this.genotypeDB= this.env.openDatabase(null, "genotype", dbConfig);
		this.databases= new Database[]{this.markerDB,this.individualDB,this.genotypeDB};
		
		SecondaryConfig secondConfig= new SecondaryConfig();
		secondConfig.setAllowCreate(!readOnly);
		/* on open, if the secondary database is empty then the primary database is read in its entirety and additions/modifications to the secondary's records occur automatically */
		secondConfig.setAllowPopulate(true);
		secondConfig.setSortedDuplicates(true);
		IndividualCorrielNameCreator indiCreator = new IndividualCorrielNameCreator();
		/* Identifies the key creator object to be used for secondary key creation.  */
		secondConfig.setKeyCreator(indiCreator);
		individualIdDB = this.env.openSecondaryDatabase(null, "corriel", this.individualDB,secondConfig);
		
		
	 	secondConfig= new SecondaryConfig();
		secondConfig.setAllowCreate(!readOnly);
		/* on open, if the secondary database is empty then the primary database is read in its entirety and additions/modifications to the secondary's records occur automatically */
		secondConfig.setAllowPopulate(true);
		secondConfig.setSortedDuplicates(true);
		MarkerPositionCreator posCreator = new MarkerPositionCreator();
		/* Identifies the key creator object to be used for secondary key creation.  */
		secondConfig.setKeyCreator(posCreator);
		positionDB = this.env.openSecondaryDatabase(null, "position", this.markerDB,secondConfig);
		}
	
	Environment getEnv()
		{
		return this.env;
		}
	
	Database getMarkerDB(){ return this.markerDB;}
	Database getIndividualDB(){ return this.individualDB;}
	Database getGenotypeDB(){ return this.genotypeDB;}
	
	Individual findIndividualByCorrielId(String individualId) throws DatabaseException
		{
		if(individualId==null) return null;
		DatabaseEntry searchKey =  new DatabaseEntry(individualId.getBytes());
		DatabaseEntry primaryKey = new DatabaseEntry();
		DatabaseEntry primaryData = new DatabaseEntry();
		if(this.individualIdDB.get(null,searchKey,primaryKey,primaryData,LockMode.DEFAULT)!=OperationStatus.SUCCESS) return null;
		return (Individual)individualTupleBinding.entryToObject(primaryData);
		}
	
	/** close the environement */
	void close() throws DatabaseException
		{
		if(this.individualIdDB!=null)
			{
			this.individualIdDB.close();
			}
		if(this.positionDB!=null)
			{
			this.positionDB.close();
			}
		for(Database db:this.databases)
			{
			if(db!=null) db.close();
			}
		this.databases= new Database[0];
		if(env!=null)
			{
			/*cleaner thread is responsible for cleaning, or compacting, your log files for you.*/
			env.cleanLog();
			/* close Environment */
			env.close();
			env=null;
			}
		}
	
	/** truncate all databases */
	void clear() throws DatabaseException
		{
		for(Object o:getEnv().getDatabaseNames())
			{
			this.getEnv().truncateDatabase(null, o.toString(), false);
			}
		}
	
	@Override
	protected void finalize() throws Throwable {
		try {close();}catch(Throwable err) {}
		super.finalize();
		}
	
	public void readHapMap(URL url) throws IOException,DatabaseException
		{
		Pattern space=Pattern.compile("[ ]");
		String HEADER[]=new String[]{"rs#","SNPalleles","chrom","pos","strand","genome_build","center","protLSID","assayLSID","panelLSID","QC_code"};
		BufferedReader in= new  BufferedReader(new InputStreamReader(new GZIPInputStream(url.openStream())));
		
		String line= in.readLine();
		if(line==null) throw new IOException("Empty file");
		String header[]=space.split(line);
		for(int i=0;i< HEADER.length;++i)
			{
			if(!header[i].equals(HEADER[i])) throw new IOException("Bad header "+header[i]+" expected "+HEADER[i]);
			}
		for(int i=HEADER.length;i< header.length;++i)
			{
			Individual individual= new Individual();
			individual.name=header[i];
			DatabaseEntry key= new DatabaseEntry(individual.name.getBytes());
			DatabaseEntry data= new DatabaseEntry();
			this.individualTupleBinding.objectToEntry(individual, data);
			getIndividualDB().put(null
					,key
					,data
					);
			}
		TupleBinding INT_BINDING=TupleBinding.getPrimitiveBinding(Integer.class);
		while((line=in.readLine())!=null)
			{
			if(!line.startsWith("rs")) continue;
			String tokens[]=space.split(line);
			//System.err.println(line);
			Marker marker= new Marker(Integer.parseInt( tokens[0].substring(2)));
			marker.alleles= tokens[1];
			marker.pos= new Position(tokens[2],Integer.parseInt(tokens[3]));
			marker.strand= tokens[4].charAt(0);
			
			DatabaseEntry key= new DatabaseEntry();
			INT_BINDING.objectToEntry(marker.getRSId(), key);
			DatabaseEntry data= new DatabaseEntry();
			this.markerTupleBinding.objectToEntry(marker, data);
			getMarkerDB().put(null
					,key
					,data
					);
			for(int i=HEADER.length;i<header.length;++i)
				{
				if(tokens[i].length()!=2 || tokens[i].equals("NN")) continue;
				MarkerIndividual mi= new MarkerIndividual(marker.rsId,header[i]);
				this.markerIndividualTupleBinding.objectToEntry(mi, key);
				this.genotypeTupleBinding.objectToEntry(
						new Genotype(tokens[i].charAt(0),tokens[i].charAt(1)),
						data
						);
				getGenotypeDB().put(null
						,key
						,data
						);
				}
			}
		in.close();
		System.err.println("Individuals count:="+getIndividualDB().count());
		System.err.println("Marker count:="+getMarkerDB().count());
		System.err.println("Genotypes count:="+getGenotypeDB().count());
		
		}
	
	public void makePedigree() throws DatabaseException,IOException
		{
		URL url=new URL("http://cardiogenomics.med.harvard.edu/groups/comp08/pages/CEPH_pedigrees.html");
		String line;
		Pattern removeTag = Pattern.compile("<[/]?[^>]*>");
		BufferedReader r= new BufferedReader(new InputStreamReader(url.openStream()));
		DatabaseEntry key=new DatabaseEntry();
		DatabaseEntry value=new DatabaseEntry();
		while((line=r.readLine())!=null)
			{
			if(!line.startsWith("<tr>")) continue;
			line=line.replaceAll("</td><td>", "\t");
			line=removeTag.matcher(line).replaceAll("");
			String tokens[]=line.split("[\t]");
			if(tokens.length!=6) continue;
			key.setData(tokens[5].getBytes());
			if(individualDB.get(null, key, value, LockMode.DEFAULT)!=OperationStatus.SUCCESS) continue;

			Individual indi= Individual.class.cast(individualTupleBinding.entryToObject(value));
			indi.individualID= tokens[1];
			indi.fatherID= (tokens[2].equals("0")?null:tokens[2]);
			indi.motherID= (tokens[3].equals("0")?null:tokens[3]);
			
			this.individualTupleBinding.objectToEntry(indi,value);
			this.individualDB.put(null
					,key
					,value
					);

			}
		r.close();
		}
	
	Genotype getGenotypeAt(Marker marker,Individual indi) throws DatabaseException
		{
		if(marker==null || indi==null) return null;
		DatabaseEntry key=new DatabaseEntry();
		DatabaseEntry value=new DatabaseEntry();
		this.markerIndividualTupleBinding.objectToEntry(new MarkerIndividual(marker.rsId,indi.name),key);
		if(this.genotypeDB.get(null, key, value, LockMode.DEFAULT)!=OperationStatus.SUCCESS) return null;
		Genotype g= Genotype.class.cast(this.genotypeTupleBinding.entryToObject(value));
		return g;
		}
	
	
	void frequences() throws DatabaseException
		{
		SecondaryCursor cM= this.positionDB.openSecondaryCursor(null, null);
		DatabaseEntry key=new DatabaseEntry();
		DatabaseEntry value=new DatabaseEntry();
		while(cM.getNext(key, value,LockMode.DEFAULT)==OperationStatus.SUCCESS)
			{
			Marker m= Marker.class.cast(this.markerTupleBinding.entryToObject(value));
			HashMap<Character, Integer> allele2count= new HashMap<Character, Integer>();
			int totalGenotyped=0;
			int totalFailures=0;
			System.out.print("rs"+m.rsId+" "+m.alleles+" "+m.pos.chromosome+" "+m.pos.position+" "+m.strand);
			Cursor cI= this.individualDB.openCursor(null, null);
			while(cI.getNext(key, value,LockMode.DEFAULT)==OperationStatus.SUCCESS)
				{
				Individual indi=Individual.class.cast(this.individualTupleBinding.entryToObject(value));
				Genotype g= getGenotypeAt(m, indi);
				if(g!=null)
					{
					totalGenotyped++;
					for(int i=0;i< 2;++i)
						{
						char c= (i==0?g.a1:g.a2);
						Integer count= allele2count.get(c);
						if(count==null) count=0;
						allele2count.put(c,count+1);
						}
					}
				else
					{
					totalFailures++;
					}
				}
			cI.close();
			System.out.print(" genotyped:"+(int)(100.0*(totalGenotyped-totalFailures)/(float)totalGenotyped)+"%");
			for(Character allele: allele2count.keySet())
				{
				System.out.print(" f("+allele+")="+allele2count.get(allele)/(2.0*totalGenotyped));
				}
			System.out.println();
			}
		cM.close();
		}
	
	
	void incompats() throws DatabaseException
		{
		DatabaseEntry key=new DatabaseEntry();
		DatabaseEntry value=new DatabaseEntry();
		
		Cursor cI= this.individualDB.openCursor(null, null);
		while(cI.getNext(key, value,LockMode.DEFAULT)==OperationStatus.SUCCESS)
			{
			Individual indi=Individual.class.cast(this.individualTupleBinding.entryToObject(value));

			if(indi.fatherID==null && indi.motherID==null)
				{
				continue;
				}
			
			Individual father= findIndividualByCorrielId(indi.fatherID);
			Individual mother= findIndividualByCorrielId(indi.motherID);
			
			Cursor cM= this.markerDB.openCursor(null, null);
			while(cM.getNext(key, value,LockMode.DEFAULT)==OperationStatus.SUCCESS)
				{
				Marker m= Marker.class.cast(this.markerTupleBinding.entryToObject(value));
				Genotype gChild= getGenotypeAt(m, indi);
				Genotype gFather= getGenotypeAt(m, father);
				if(isIncompat(gChild,gFather))
					{
					System.out.println("Incompat: for rs"+m.getRSId()+"("+m.alleles+") "+
						indi.individualID+" is "+gChild+" and his father "+indi.fatherID+" is "+
						gFather
						);
					continue;
					}
				Genotype gMother= getGenotypeAt(m, mother);
				if(isIncompat(gChild,gMother))
					{
					System.out.println("Incompat: for rs"+m.getRSId()+"("+m.alleles+") "+
						indi.individualID+" is "+gChild+" and his mother "+indi.motherID+" is "+
						gMother
						);
					continue;
					}
				if(isIncompat(gChild,gFather,gMother))
					{
					System.out.println("Incompat: for rs"+m.getRSId()+"("+m.alleles+") "+
						indi.individualID+" is "+gChild+
						" and his father "+indi.fatherID+" is "+ gFather+
						" and his mother "+indi.motherID+" is "+ gMother
						);
					}
				}
			cM.close();
			}
		cI.close();
		
			
		
		}
	
	static boolean isIncompat(Genotype child,Genotype parent)
		{
		if(child==null || parent==null) return false;
		return !(
			child.a1==parent.a1 ||
			child.a1==parent.a2 ||
			child.a2==parent.a1 ||
			child.a2==parent.a2
			);
		}
	
	static boolean isIncompat(Genotype child,Genotype father,Genotype mother)
		{
		if(child==null || father==null || mother==null) return false;
		return !(
			   child.same(father.a1,mother.a1) ||
			   child.same(father.a1,mother.a2) ||
			   child.same(father.a2,mother.a1) ||
			   child.same(father.a2,mother.a2)
			   )
			   ;
		}
	
	/**
	 * main
	 * @param args arguments
	 */
	public static void main(String[] args) {
		try {
			File f= new File("/home/pierre/tmp/berkeleydb");
			Test01 test= new Test01();
			test.open(f, false);
			//test.clear();
			
			//List databases in test
			for(Object name: test.getEnv().getDatabaseNames())
				{
				System.out.println("Database:"+name);
				}
			test.readHapMap( new URL("http://www.hapmap.org/genotypes/latest/rs_strand/non-redundant/genotypes_chrMT_CEU_r23a_nr.b36.txt.gz"));
			test.makePedigree();
			test.frequences();
			test.incompats();
			test.close();
			} 
		catch (Exception e)
			{
			e.printStackTrace();
			}
		}
}
Pastie Show Pastie some love
