← Back to branch summary

~ubuntu-branches/ubuntu/natty/python-cogent/natty

« back to all changes in this revision

Viewing changes to doc/cookbook/blast.rst

  • Committer: Bazaar Package Importer
  • Author(s): Steffen Moeller
  • Date: 2010-12-04 22:30:35 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20101204223035-j11kinhcrrdgg2p2
Tags: 1.5-1
* Bumped standard to 3.9.1, no changes required.
* New upstream version.
  - major additions to Cookbook
  - added AlleleFreqs attribute to ensembl Variation objects.
  - added getGeneByStableId method to genome objects.
  - added Introns attribute to Transcript objects and an Intron class.
  - added Mann-Whitney test and a Monte-Carlo version
  - exploratory and confirmatory period estimation techniques (suitable for
    symbolic and continuous data)
  - Information theoretic measures (AIC and BIC) added
  - drawing of trees with collapsed nodes
  - progress display indicator support for terminal and GUI apps
  - added parser for illumina HiSeq2000 and GAiix sequence files as 
    cogent.parse.illumina_sequence.MinimalIlluminaSequenceParser.
  - added parser to FASTQ files, one of the output options for illumina's
    workflow, also added cookbook demo.
  - added functionality for parsing of SFF files without the Roche tools in
    cogent.parse.binary_sff
  - thousand fold performance improvement to nmds
  - >10-fold performance improvements to some Table operations

Show diffs side-by-side

added added

removed removed

Lines of Context:
doc/cookbook/blast.rst
4
4
Controlling BLAST
5
5
*****************
6
6
 
7
 
.. authors, Gavin Huttley, Tom Elliott
 
7
.. authors, Gavin Huttley, Tom Elliott, Jeremy Widmann
8
8
 
9
9
Preliminaries
10
10
-------------
11
11
 
12
12
In order to run BLAST locally (from a program running on your computer) you will need to do three things:
13
13
 
14
 
- download the BLAST "executables" from NCBI
15
 
- make sure these programs are available on your ``PATH``
16
 
- construct and format a database to search against
 
14
- Download the BLAST "executables" from NCBI
 
15
- Make sure these programs are available on your ``PATH``
 
16
- Construct and format a database to search against
17
17
 
18
18
NCBI has recently changed the BLAST programs, and as yet PyCogent does not support the new versions. The "legacy" versions are available from `here <http://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download>`_ (login as guest).
19
19
 
20
 
Detailed instructions are beyond the scope of this example, but after downloading the programs and setting up your ``PATH``, you should test BLAST by doing this from the command line:
 
20
 
 
21
Detailed installation instructions are beyond the scope of this example, but are available at `NCBI's website <http://www.ncbi.nlm.nih.gov/staff/tao/URLAPI/unix_setup.html>`_ .
 
22
After downloading the programs and setting up your ``PATH``, you should test BLAST by doing this from the command line:
21
23
 
22
24
::
23
25
    
33
35
      -d  Database [String]
34
36
        default = nr...
35
37
 
36
 
The file ``refseqs.fasta`` contains some short sequences for use in the following examples. It is available from :download:`here <../data/refseqs.fasta>`.
37
 
 
38
 
.. TODO add to data_file_links.rst
 
38
The file ``refseqs.fasta`` contains some short sequences for use in the following examples.
39
39
 
40
40
.. doctest::
41
41
 
140
140
    e-value              6e-05
141
141
    bit score            26.3
142
142
 
 
143
BLAST with XML output
 
144
---------------------
 
145
 
 
146
In this example, we load a DNA sequence from a file in the data directory and BLAST against our formatted database as above.
 
147
 
143
148
NCBI recommends that you use XML as the output for BLAST. (They reserve the right to change the format for other output types). XML is the output when we pass '-m':'7'.
144
149
 
145
150
.. doctest::
209
214
    >>> remove_files(['data/blast_test.txt', 'data/blast_test.xml'],
210
215
    ...              error_on_missing=False)
211
216
 
 
217
BLAST with protein sequences
 
218
----------------------------
 
219
 
 
220
In this example, we load a protein sequence from a file in the data directory and BLAST against a new protein database we will create.  Since we want to BLAST protein sequences instead of DNA, we will have to construct a new BLAST database.
 
221
 
 
222
The file ``refseqs_protein.fasta`` contains some short sequences for use in the following examples.
 
223
 
 
224
.. doctest::
 
225
    
 
226
    >>> from cogent.app.formatdb import build_blast_db_from_fasta_path
 
227
    >>> result = build_blast_db_from_fasta_path('data/refseqs_protein.fasta', is_protein=True)
 
228
    >>> result[0]
 
229
    'data/refseqs_protein.fasta'...
 
230
 
 
231
Notice that we set the parameter ``is_protein`` to ``True`` since our database consists of protein sequences this time.  This was not necessary in the previous example, because ``is_protein`` is set to ``False`` by default.
 
232
 
 
233
Now that we have built our protein BLAST database, we can load our sequence and BLAST against this database.
 
234
 
 
235
.. doctest::
 
236
 
 
237
    >>> from cogent import LoadSeqs, PROTEIN
 
238
    >>> from cogent.app.blast import blast_seqs, Blastall
 
239
    >>> seqs = LoadSeqs('data/inseqs_protein.fasta', moltype=PROTEIN, aligned=False)
 
240
    >>> seq = seqs.getSeq('1091044_fragment')
 
241
        >>> seq
 
242
        ProteinSequence(IPLDFDK... 26)
 
243
        
 
244
Notice we need to use '-p':'blastp' in the parameters dictionary, since ``blastp`` is used for protein.
 
245
 
 
246
.. doctest::
 
247
        
 
248
        >>> params={'-p':'blastp','-m':'9'}
 
249
        >>> result = blast_seqs([seq], 
 
250
        ...    Blastall, 
 
251
        ...    blast_db = 'data/refseqs_protein.fasta',
 
252
        ...    params = params)
 
253
        >>> data = result['StdOut'].read()
 
254
        >>> print data.split('\n')[:1]
 
255
        ['# BLASTP 2.2...
 
256
 
 
257
We save the results for further processing 
 
258
 
 
259
.. doctest::
 
260
    
 
261
    >>> outfile = open('data/blast_protein_test.txt','w')
 
262
    >>> outfile.write(data)
 
263
    >>> outfile.close()
 
264
 
 
265
Now we will explore some of the convenience methods of the ``BlastResult`` object.
 
266
 
 
267
.. doctest::
 
268
 
 
269
        >>> from cogent.parse.blast import BlastResult
 
270
        >>> blast_results = BlastResult(open('data/blast_protein_test.txt','r'))
 
271
 
 
272
Suppose we want to filter our results based on various criteria.  In many cases you may want to only keep the top '3' matches with the longest 'ALIGNMENT LENGTH' for the query sequence to the target.
 
273
 
 
274
.. doctest::
 
275
 
 
276
        >>> best_hits = dict(blast_results.bestHitsByQuery(field='ALIGNMENT LENGTH', n=3))
 
277
        >>> query_1_best_hits = best_hits['1']
 
278
        >>> for hit in query_1_best_hits:
 
279
        ...     for key, value in hit.items():
 
280
        ...             print key.ljust(20), value
 
281
        ...     print
 
282
        ... 
 
283
        MISMATCHES           0
 
284
        ALIGNMENT LENGTH     26
 
285
        Q. END               26
 
286
        BIT SCORE            56.2
 
287
        % IDENTITY           100.00
 
288
        Q. START             1
 
289
        S. START             30
 
290
        S. END               55
 
291
        GAP OPENINGS         0
 
292
        QUERY ID             1
 
293
        E-VALUE              5e-12
 
294
        SUBJECT ID           1091044
 
295
        <BLANKLINE>
 
296
        MISMATCHES           10
 
297
        ALIGNMENT LENGTH     27
 
298
        Q. END               25
 
299
        BIT SCORE            33.5
 
300
        % IDENTITY           55.56
 
301
        Q. START             1
 
302
        S. START             32
 
303
        S. END               58
 
304
        GAP OPENINGS         1
 
305
        QUERY ID             1
 
306
        E-VALUE              3e-05
 
307
        SUBJECT ID           5326864
 
308
        <BLANKLINE>
 
309
        MISMATCHES           16
 
310
        ALIGNMENT LENGTH     24
 
311
        Q. END               25
 
312
        BIT SCORE            22.3
 
313
        % IDENTITY           33.33
 
314
        Q. START             2
 
315
        S. START             19
 
316
        S. END               42
 
317
        GAP OPENINGS         0
 
318
        QUERY ID             1
 
319
        E-VALUE              0.077
 
320
        SUBJECT ID           14286173
 
321
        <BLANKLINE>
 
322
 
 
323
The fist of the top 3 hits for alignment length has 0 MISMATCHES and a % IDENTITY of 100.00.  The next 2 hits have many MISMATCHES and a much lower % IDENTITY.  Lets filter the results again, but by E-VALUE this time:
 
324
 
 
325
.. doctest::
 
326
 
 
327
        >>> best_hits = dict(blast_results.bestHitsByQuery(field='E-VALUE', n=3))
 
328
        >>> query_1_best_hits = best_hits['1']
 
329
        >>> for hit in query_1_best_hits:
 
330
        ...     for key, value in hit.items():
 
331
        ...             print key.ljust(20), value
 
332
        ...     print
 
333
        ... 
 
334
        MISMATCHES           0
 
335
        ALIGNMENT LENGTH     26
 
336
        Q. END               26
 
337
        BIT SCORE            56.2
 
338
        % IDENTITY           100.00
 
339
        Q. START             1
 
340
        S. START             30
 
341
        S. END               55
 
342
        GAP OPENINGS         0
 
343
        QUERY ID             1
 
344
        E-VALUE              5e-12
 
345
        SUBJECT ID           1091044
 
346
        <BLANKLINE>
 
347
        MISMATCHES           10
 
348
        ALIGNMENT LENGTH     27
 
349
        Q. END               25
 
350
        BIT SCORE            33.5
 
351
        % IDENTITY           55.56
 
352
        Q. START             1
 
353
        S. START             32
 
354
        S. END               58
 
355
        GAP OPENINGS         1
 
356
        QUERY ID             1
 
357
        E-VALUE              3e-05
 
358
        SUBJECT ID           5326864
 
359
        <BLANKLINE>
 
360
        MISMATCHES           6
 
361
        ALIGNMENT LENGTH     18
 
362
        Q. END               26
 
363
        BIT SCORE            30.4
 
364
        % IDENTITY           66.67
 
365
        Q. START             9
 
366
        S. START             31
 
367
        S. END               48
 
368
        GAP OPENINGS         0
 
369
        QUERY ID             1
 
370
        E-VALUE              3e-04
 
371
        SUBJECT ID           15964668
 
372
        <BLANKLINE>
 
373
 
 
374
You can filter the BLAST results by any of the fields you like.  You can also use the ``BlastResult`` object to do a quick assessment of your BLAST results looking only at the fields you like:
 
375
 
 
376
.. doctest::
 
377
 
 
378
        >>> fields = ['SUBJECT ID', 'BIT SCORE', 'E-VALUE']
 
379
        >>> for query, results in blast_results.items():
 
380
        ...     print ''.join([f.ljust(20) for f in fields])
 
381
        ...     for result in results[-1]:
 
382
        ...             print ''.join(map(str,[result[field].ljust(20) for field in fields]))
 
383
        SUBJECT ID          BIT SCORE           E-VALUE             
 
384
        1091044             56.2                5e-12               
 
385
        5326864             33.5                3e-05               
 
386
        15964668            30.4                3e-04               
 
387
        17229033            29.6                5e-04               
 
388
        21112072            28.1                0.001               
 
389
        4704732             25.8                0.007               
 
390
        13541117            24.6                0.016               
 
391
        15826629            24.3                0.020               
 
392
        14286173            22.3                0.077               
 
393
        6323138             21.9                0.10                
 
394
        18313548            20.8                0.22                
 
395
        21674812            20.0                0.38                
 
396
        14600438            20.0                0.38                
 
397
        4996210             18.5                1.1                 
 
398
        15605963            17.3                2.5                 
 
399
        15615431            16.5                4.2                 
 
400
 
 
401
.. doctest::
 
402
            :hide:
 
403
 
 
404
            >>> from cogent.util.misc import remove_files
 
405
            >>> remove_files(['data/blast_protein_test.txt'],
 
406
            ...              error_on_missing=False)