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.net.URI;
  30 import java.net.URISyntaxException;
  31 import java.util.Arrays;
  32 import java.util.Enumeration;
  33 import java.util.Hashtable;
  34 import java.util.Vector;
  35
  36 import javax.swing.JApplet;
  37
  38 import org.forester.archaeopteryx.AptxUtil;
  39 import org.forester.archaeopteryx.TreePanel;
  40 import org.forester.phylogeny.PhylogenyNode;
  41 import org.forester.phylogeny.data.Accession;
  42 import org.forester.phylogeny.data.Identifier;
  43 import org.forester.util.ForesterUtil;
  44 import org.forester.util.SequenceIdParser;
  45 import org.forester.ws.wabi.RestUtil;
  46
  47 public final class Blast {
  48
  49     final public static void openNcbiBlastWeb( final String query,
  50                                                final boolean is_nucleic_acids,
  51                                                final JApplet applet,
  52                                                final TreePanel p ) {
  53         //http://www.ncbi.nlm.nih.gov/blast/Blast.cgi?CMD=Web&PAGE=Proteins&DATABASE=swissprot&QUERY=gi|163848401
  54         final StringBuilder uri_str = new StringBuilder();
  55         uri_str.append( "http://www.ncbi.nlm.nih.gov/blast/Blast.cgi?CMD=Web&DATABASE=nr&PAGE=" );
  56         if ( is_nucleic_acids ) {
  57             uri_str.append( "Nucleotide" );
  58         }
  59         else {
  60             uri_str.append( "Proteins" );
  61         }
  62         uri_str.append( "&QUERY=" );
  63         uri_str.append( query );
  64         try {
  65             AptxUtil.launchWebBrowser( new URI( uri_str.toString() ), applet != null, applet, "_aptx_blast" );
  66         }
  67         catch ( final IOException e ) {
  68             AptxUtil.showErrorMessage( p, e.toString() );
  69             e.printStackTrace();
  70         }
  71         catch ( final URISyntaxException e ) {
  72             AptxUtil.showErrorMessage( p, e.toString() );
  73             e.printStackTrace();
  74         }
  75     }
  76
  77     final public static String obtainQueryForBlast( final PhylogenyNode node ) {
  78         String query = "";
  79         if ( node.getNodeData().isHasSequence() ) {
  80             if ( !ForesterUtil.isEmpty( node.getNodeData().getSequence().getMolecularSequence() ) ) {
  81                 query = node.getNodeData().getSequence().getMolecularSequence();
  82             }
  83             else if ( ( node.getNodeData().getSequence().getAccession() != null )
  84                     && !ForesterUtil.isEmpty( node.getNodeData().getSequence().getAccession().getValue() ) ) {
  85                 if ( !ForesterUtil.isEmpty( node.getNodeData().getSequence().getAccession().getSource() ) ) {
  86                     query = node.getNodeData().getSequence().getAccession().getSource() + "%7C";
  87                 }
  88                 query += node.getNodeData().getSequence().getAccession().getValue();
  89             }
  90             else if ( !ForesterUtil.isEmpty( node.getNodeData().getSequence().getName() ) ) {
  91                 final Accession acc = AptxUtil.obtainSequenceAccessionFromName( node.getNodeData().getSequence()
  92                         .getName() );
  93                 if ( acc != null ) {
  94                     query = acc.getSource() + "%7C" + acc.getValue();
  95                 }
  96             }
  97         }
  98         else if ( !ForesterUtil.isEmpty( node.getName() ) ) {
  99             final Accession acc = AptxUtil.obtainSequenceAccessionFromName( node.getName() );
 100             if ( acc != null ) {
 101                 query = acc.getSource() + "%7C" + acc.getValue();
 102             }
 103             else {
 104                 final Identifier id = SequenceIdParser.parse( node.getName() );
 105                 if ( id != null ) {
 106                     query = id.getValue();
 107                 }
 108             }
 109         }
 110         return query;
 111     }
 112
 113     final public static boolean isContainsQueryForBlast( final PhylogenyNode node ) {
 114         return !ForesterUtil.isEmpty( obtainQueryForBlast( node ) );
 115     }
 116
 117     final public void ddbjBlast( final String geneName ) {
 118         // Retrieve accession number list which has specified gene name from searchByXMLPath of ARSA. Please click here for details of ARSA.
 119         /*target: Sequence length is between 300bp and 1000bp.
 120         Feature key is CDS.
 121         Gene qualifire is same as specified gene name.*/
 122         String queryPath = "/ENTRY/DDBJ/division=='HUM' AND (/ENTRY/DDBJ/length>=300 AND "
 123                 + "/ENTRY/DDBJ/length<=1000) ";
 124         queryPath += "AND (/ENTRY/DDBJ/feature-table/feature{/f_key = 'CDS' AND ";
 125         queryPath += "/f_quals/qualifier{/q_name = 'gene' AND /q_value=='" + geneName + "'}})";
 126         String query = "service=ARSA&method=searchByXMLPath&queryPath=" + queryPath
 127                 + "&returnPath=/ENTRY/DDBJ/primary-accession&offset=1&count=100";
 128         //Execute ARSA
 129         String arsaResult = null;
 130         try {
 131             arsaResult = RestUtil.getResult( query );
 132         }
 133         catch ( final IOException e ) {
 134             // TODO Auto-generated catch block
 135             e.printStackTrace();
 136         }
 137         final String[] arsaResultLines = arsaResult.split( "\n" );
 138         //Get hit count
 139         final int arsaResultNum = Integer.parseInt( arsaResultLines[ 0 ].replaceAll( "hitscount       =", "" ).trim() );
 140         //If there is no hit, print a message and exit
 141         if ( arsaResultNum == 0 ) {
 142             System.out.println( "There is no entry for gene:" + geneName );
 143             return;
 144         }
 145         //Retrieve DNA sequence of top hit entry by using getFASTA_DDBJEntry of GetEntry.
 146         //Retrieve DNA sequence of first fit.
 147         final String repAccession = arsaResultLines[ 2 ];
 148         query = "service=GetEntry&method=getFASTA_DDBJEntry&accession=" + repAccession;
 149         String dnaSeq = null;
 150         try {
 151             dnaSeq = RestUtil.getResult( query );
 152         }
 153         catch ( final IOException e ) {
 154             // TODO Auto-generated catch block
 155             e.printStackTrace();
 156         }
 157         System.out.println( "Retrieved DNA sequence is: " + dnaSeq );
 158         //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).
 159         //Execute blastn with step3's sequence
 160         query = "service=Blast&method=searchParam&program=blastn&database=ddbjpri ddbjrod ddbjmam ddbjvrt "
 161                 + "ddbjinv&query=" + dnaSeq + "&param=-m 8 -b 50 -v 50 -e 0.0001";
 162         String blastResult = null;
 163         try {
 164             blastResult = RestUtil.getResult( query );
 165         }
 166         catch ( final IOException e ) {
 167             // TODO Auto-generated catch block
 168             e.printStackTrace();
 169         }
 170         // Extract both accession number and similarity score from BLAST result.
 171         // 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.
 172         final String blastResultLines[] = blastResult.split( "\n" );
 173         final Vector<String[]> parsedBlastResult = new Vector<String[]>();
 174         for( final String blastResultLine : blastResultLines ) {
 175             final String cols[] = blastResultLine.split( "\t" );
 176             final String accession = cols[ 1 ].substring( 0, cols[ 1 ].indexOf( "|" ) );
 177             final String[] result = { accession, cols[ 2 ] };
 178             parsedBlastResult.add( result );
 179         }
 180         // Retrieve species name by using searchByXMLPath of ARSA. If the plural subjects whose species
 181         // name are same are in the result, the highest similarity score is used.
 182         //Retrieve species from accession number.
 183         final Hashtable<String, String> organismAccession = new Hashtable<String, String>();
 184         for( int i = 0; i < parsedBlastResult.size(); i++ ) {
 185             final String[] parsed = parsedBlastResult.elementAt( i );
 186             query = "service=ARSA&method=searchByXMLPath&queryPath=/ENTRY/DDBJ/primary-accession=='" + parsed[ 0 ]
 187                     + "'&returnPath=/ENTRY/DDBJ/organism&offset=1&count=100";
 188             String organism = null;
 189             try {
 190                 organism = RestUtil.getResult( query );
 191             }
 192             catch ( final IOException e ) {
 193                 // TODO Auto-generated catch block
 194                 e.printStackTrace();
 195             }
 196             final String[] organismLines = organism.split( "\n" );
 197             organism = organismLines[ 2 ];
 198             //If same organism name hits, use first hit.
 199             if ( !organismAccession.containsKey( organism ) ) {
 200                 organismAccession.put( organism, parsed[ 0 ] + "\t" + parsed[ 1 ] );
 201             }
 202         }
 203         // Print result.
 204         // Print Result
 205         System.out.println( "DDBJ entries: " + arsaResultNum );
 206         System.out.println( "Representative accession: " + repAccession );
 207         System.out.println( "Organism name\tDDBJ accession number\tSequence similarity" );
 208         final String[] keys = new String[ organismAccession.size() ];
 209         final Enumeration<String> enu = organismAccession.keys();
 210         int count = 0;
 211         while ( enu.hasMoreElements() ) {
 212             keys[ count ] = enu.nextElement();
 213             ++count;
 214         }
 215         Arrays.sort( keys );
 216         for( final String key : keys ) {
 217             System.out.println( key + "\t" + organismAccession.get( key ) );
 218         }
 219     }
 220 }