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