forester/java/src/org/forester/archaeopteryx/tools/Blast.java

   1 // $Id:
   2 // FORESTER -- software libraries and applications
   3 // for evolutionary biology research and applications.
   4 //
   5 // Copyright (C) 2008-2009 Christian M. Zmasek
   6 // Copyright (C) 2008-2009 Burnham Institute for Medical Research
   7 // All rights reserved
   8 //
   9 // This library is free software; you can redistribute it and/or
  10 // modify it under the terms of the GNU Lesser General Public
  11 // License as published by the Free Software Foundation; either
  12 // version 2.1 of the License, or (at your option) any later version.
  13 //
  14 // This library is distributed in the hope that it will be useful,
  15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  17 // Lesser General Public License for more details.
  18 //
  19 // You should have received a copy of the GNU Lesser General Public
  20 // License along with this library; if not, write to the Free Software
  21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA
  22 //
  23 // Contact: phylosoft @ gmail . com
  24 // WWW: www.phylosoft.org/forester
  25
  26 package org.forester.archaeopteryx.tools;
  27
  28 import java.io.IOException;
  29 import java.util.Arrays;
  30 import java.util.Enumeration;
  31 import java.util.Hashtable;
  32 import java.util.Vector;
  33
  34 import org.forester.ws.wabi.RestUtil;
  35
  36 public class Blast {
  37
  38     public Blast() {
  39     }
  40
  41     public void go( final String geneName ) {
  42         // Retrieve accession number list which has specified gene name from searchByXMLPath of ARSA. Please click here for details of ARSA.
  43         /*target: Sequence length is between 300bp and 1000bp.
  44         Feature key is CDS.
  45         Gene qualifire is same as specified gene name.*/
  46         String queryPath = "/ENTRY/DDBJ/division=='HUM' AND (/ENTRY/DDBJ/length>=300 AND "
  47                 + "/ENTRY/DDBJ/length<=1000) ";
  48         queryPath += "AND (/ENTRY/DDBJ/feature-table/feature{/f_key = 'CDS' AND ";
  49         queryPath += "/f_quals/qualifier{/q_name = 'gene' AND /q_value=='" + geneName + "'}})";
  50         String query = "service=ARSA&method=searchByXMLPath&queryPath=" + queryPath
  51                 + "&returnPath=/ENTRY/DDBJ/primary-accession&offset=1&count=100";
  52         //Execute ARSA
  53         String arsaResult = null;
  54         try {
  55             arsaResult = RestUtil.getResult( query );
  56         }
  57         catch ( final IOException e ) {
  58             // TODO Auto-generated catch block
  59             e.printStackTrace();
  60         }
  61         final String[] arsaResultLines = arsaResult.split( "\n" );
  62         //Get hit count
  63         final int arsaResultNum = Integer.parseInt( arsaResultLines[ 0 ].replaceAll( "hitscount       =", "" ).trim() );
  64         //If there is no hit, print a message and exit
  65         if ( arsaResultNum == 0 ) {
  66             System.out.println( "There is no entry for gene:" + geneName );
  67             return;
  68         }
  69         //Retrieve DNA sequence of top hit entry by using getFASTA_DDBJEntry of GetEntry.
  70         //Retrieve DNA sequence of first fit.
  71         final String repAccession = arsaResultLines[ 2 ];
  72         query = "service=GetEntry&method=getFASTA_DDBJEntry&accession=" + repAccession;
  73         String dnaSeq = null;
  74         try {
  75             dnaSeq = RestUtil.getResult( query );
  76         }
  77         catch ( final IOException e ) {
  78             // TODO Auto-generated catch block
  79             e.printStackTrace();
  80         }
  81         System.out.println( "Retrieved DNA sequence is: " + dnaSeq );
  82         //Execute blastn by using searchParam of Blast with step2's sequence. Specified option is -e 0.0001 -m 8 -b 50 -v 50. It means "Extract top 50 hit which E-value is more than 0.0001.". The reference databases are specified as follows. ddbjpri(primates) ddbjrod(rodents) ddbjmam(mammals) ddbjvrt(vertebrates ) ddbjinv(invertebrates).
  83         //Execute blastn with step3's sequence
  84         query = "service=Blast&method=searchParam&program=blastn&database=ddbjpri ddbjrod ddbjmam ddbjvrt "
  85                 + "ddbjinv&query=" + dnaSeq + "&param=-m 8 -b 50 -v 50 -e 0.0001";
  86         String blastResult = null;
  87         try {
  88             blastResult = RestUtil.getResult( query );
  89         }
  90         catch ( final IOException e ) {
  91             // TODO Auto-generated catch block
  92             e.printStackTrace();
  93         }
  94         // Extract both accession number and similarity score from BLAST result.
  95         // This step does not use Web API and extract the part of result or edit the result. Please click here to see the details of each column in the BLAST tab delimited format which is generated by -m 8 option.
  96         final String blastResultLines[] = blastResult.split( "\n" );
  97         final Vector<String[]> parsedBlastResult = new Vector<String[]>();
  98         for( final String blastResultLine : blastResultLines ) {
  99             final String cols[] = blastResultLine.split( "\t" );
 100             final String accession = cols[ 1 ].substring( 0, cols[ 1 ].indexOf( "|" ) );
 101             final String[] result = { accession, cols[ 2 ] };
 102             parsedBlastResult.add( result );
 103         }
 104         // Retrieve species name by using searchByXMLPath of ARSA. If the plural subjects whose species
 105         // name are same are in the result, the highest similarity score is used.
 106         //Retrieve species from accession number.
 107         final Hashtable<String, String> organismAccession = new Hashtable<String, String>();
 108         for( int i = 0; i < parsedBlastResult.size(); i++ ) {
 109             final String[] parsed = parsedBlastResult.elementAt( i );
 110             query = "service=ARSA&method=searchByXMLPath&queryPath=/ENTRY/DDBJ/primary-accession=='" + parsed[ 0 ]
 111                     + "'&returnPath=/ENTRY/DDBJ/organism&offset=1&count=100";
 112             String organism = null;
 113             try {
 114                 organism = RestUtil.getResult( query );
 115             }
 116             catch ( final IOException e ) {
 117                 // TODO Auto-generated catch block
 118                 e.printStackTrace();
 119             }
 120             final String[] organismLines = organism.split( "\n" );
 121             organism = organismLines[ 2 ];
 122             //If same organism name hits, use first hit.
 123             if ( !organismAccession.containsKey( organism ) ) {
 124                 organismAccession.put( organism, parsed[ 0 ] + "\t" + parsed[ 1 ] );
 125             }
 126         }
 127         // Print result.
 128         // Print Result
 129         System.out.println( "DDBJ entries: " + arsaResultNum );
 130         System.out.println( "Representative accession: " + repAccession );
 131         System.out.println( "Organism name\tDDBJ accession number\tSequence similarity" );
 132         final String[] keys = new String[ organismAccession.size() ];
 133         final Enumeration<String> enu = organismAccession.keys();
 134         int count = 0;
 135         while ( enu.hasMoreElements() ) {
 136             keys[ count ] = enu.nextElement();
 137             ++count;
 138         }
 139         Arrays.sort( keys );
 140         for( final String key : keys ) {
 141             System.out.println( key + "\t" + organismAccession.get( key ) );
 142         }
 143     }
 144 }