File: chapter_searchio.tex

package info (click to toggle)
python-biopython 1.68%2Bdfsg-3~bpo8%2B1
  • links: PTS, VCS
  • area: main
  • in suites: jessie-backports
  • size: 46,856 kB
  • sloc: python: 160,306; xml: 93,216; ansic: 9,118; sql: 1,208; makefile: 155; sh: 63
file content (1184 lines) | stat: -rw-r--r-- 53,598 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
\chapter{BLAST and other sequence search tools (\textit{experimental code})}
\label{chapter:searchio}

\emph{WARNING}: This chapter of the Tutorial describes an \emph{experimental}
module in Biopython. It is being included in Biopython and documented
here in  the tutorial in a pre-final state to allow a period of feedback
and refinement before we declare it stable. Until then the details will
probably change, and any scripts using the current \verb|Bio.SearchIO|
would need to be updated. Please keep this in mind! For stable code
working with NCBI BLAST, please continue to use Bio.Blast described
in the preceding Chapter~\ref{chapter:blast}.

Biological sequence identification is an integral part of bioinformatics.
Several tools are available for this, each with their own algorithms and
approaches, such as BLAST (arguably the most popular), FASTA, HMMER, and many
more. In general, these tools usually use your sequence to search a database of
potential matches. With the growing number of known sequences (hence the
growing number of potential matches), interpreting the results becomes
increasingly hard as there could be hundreds or even thousands of potential
matches. Naturally, manual interpretation of these searches' results is out of
the question. Moreover, you often need to work with several sequence search
tools, each with its own statistics, conventions, and output format. Imagine how
daunting it would be when you need to work with multiple sequences using
multiple search tools.

We know this too well ourselves, which is why we created the \verb|Bio.SearchIO|
submodule in Biopython. \verb|Bio.SearchIO| allows you to extract information
from your search results in a convenient way, while also dealing with the
different standards and conventions used by different search tools.
The name \verb|SearchIO| is a homage to BioPerl's module of the same name.

In this chapter, we'll go through the main features of \verb|Bio.SearchIO| to
show what it can do for you. We'll use two popular search tools along the way:
BLAST and BLAT. They are used merely for illustrative purposes, and you should
be able to adapt the workflow to any other search tools supported by
\verb|Bio.SearchIO| in a breeze. You're very welcome to follow along with the
search output files we'll be using. The BLAST output file can be downloaded
\href{http://biopython.org/SRC/biopython/Doc/examples/my_blast.xml}{here},
and the BLAT output file
\href{http://biopython.org/SRC/biopython/Doc/examples/my_blat.psl}{here}
or are included with the Biopython source code under the \verb|Doc/examples/|
folder. Both output files were generated using this sequence:

\begin{verbatim}
>mystery_seq
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
\end{verbatim}

The BLAST result is an XML file generated using \verb|blastn| against the NCBI
\verb|refseq_rna| database. For BLAT, the sequence database was the February 2009
\verb|hg19| human genome draft and the output format is PSL.

We'll start from an introduction to the \verb|Bio.SearchIO| object model. The
model is the representation of your search results, thus it is core to
\verb|Bio.SearchIO| itself. After that, we'll check out the main functions in
\verb|Bio.SearchIO| that you may often use.

Now that we're all set, let's go to the first step: introducing the core
object model.

\section{The SearchIO object model}
\label{sec:searchio-model}

Despite the wildly differing output styles among many sequence search tools,
it turns out that their underlying concept is similar:

\begin{itemize}
\item The output file may contain results from one or more search queries.
\item In each search query, you will see one or more hits from the given
search database.
\item In each database hit, you will see one or more regions containing the
actual sequence alignment between your query sequence and the database
sequence.
\item Some programs like BLAT or Exonerate may further split these regions into
several alignment fragments (or blocks in BLAT and possibly exons in
exonerate). This is not something you always see, as programs like BLAST and
HMMER do not do this.
\end{itemize}

Realizing this generality, we decided use it as base for creating the
\verb|Bio.SearchIO| object model. The object model consists of a nested
hierarchy of Python objects, each one representing one concept outlined above.
These objects are:

\begin{itemize}
\item \verb|QueryResult|, to represent a single search query.
\item \verb|Hit|, to represent a single database hit. \verb|Hit| objects are
    contained within \verb|QueryResult| and in each \verb|QueryResult| there is
    zero or more \verb|Hit| objects.
\item \verb|HSP| (short for high-scoring pair), to represent region(s) of
    significant alignments between query and hit sequences. \verb|HSP| objects
    are contained within \verb|Hit| objects and each \verb|Hit| has one or more
    \verb|HSP| objects.
\item \verb|HSPFragment|, to represent a single contiguous alignment between
    query and hit sequences. \verb|HSPFragment| objects are contained within
    \verb|HSP| objects. Most sequence search tools like BLAST and HMMER unify
    \verb|HSP| and \verb|HSPFragment| objects as each \verb|HSP| will only have
    a single \verb|HSPFragment|. However there are tools like BLAT and Exonerate
    that produce \verb|HSP| containing multiple \verb|HSPFragment|. Don't worry
    if this seems a tad confusing now, we'll elaborate more on these two objects
    later on.
\end{itemize}

These four objects are the ones you will interact with when you use
\verb|Bio.SearchIO|. They are created using one of the main \verb|Bio.SearchIO|
methods: \verb|read|, \verb|parse|, \verb|index|, or \verb|index_db|. The
details of these methods are provided in later sections. For this section, we'll
only be using read and parse. These functions behave similarly to their
\verb|Bio.SeqIO| and \verb|Bio.AlignIO| counterparts:

\begin{itemize}
\item \verb|read| is used for search output files with a single query and
    returns a \verb|QueryResult| object
\item \verb|parse| is used for search output files with multiple queries and
    returns a generator that yields \verb|QueryResult| objects
\end{itemize}

With that settled, let's start probing each \verb|Bio.SearchIO| object,
beginning with \verb|QueryResult|.

\subsection{QueryResult}
\label{sec:searchio-qresult}

The QueryResult object represents a single search query and contains zero or
more Hit objects. Let's see what it looks like using the BLAST file we have:

%doctest ../Doc/examples
\begin{verbatim}
>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> print(blast_qresult)
Program: blastn (2.2.27+)
  Query: 42291 (61)
         mystery_seq
 Target: refseq_rna
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0      1  gi|262205317|ref|NR_030195.1|  Homo sapiens microRNA 52...
            1      1  gi|301171311|ref|NR_035856.1|  Pan troglodytes microRNA...
            2      1  gi|270133242|ref|NR_032573.1|  Macaca mulatta microRNA ...
            3      2  gi|301171322|ref|NR_035857.1|  Pan troglodytes microRNA...
            4      1  gi|301171267|ref|NR_035851.1|  Pan troglodytes microRNA...
            5      2  gi|262205330|ref|NR_030198.1|  Homo sapiens microRNA 52...
            6      1  gi|262205302|ref|NR_030191.1|  Homo sapiens microRNA 51...
            7      1  gi|301171259|ref|NR_035850.1|  Pan troglodytes microRNA...
            8      1  gi|262205451|ref|NR_030222.1|  Homo sapiens microRNA 51...
            9      2  gi|301171447|ref|NR_035871.1|  Pan troglodytes microRNA...
           10      1  gi|301171276|ref|NR_035852.1|  Pan troglodytes microRNA...
           11      1  gi|262205290|ref|NR_030188.1|  Homo sapiens microRNA 51...
           12      1  gi|301171354|ref|NR_035860.1|  Pan troglodytes microRNA...
           13      1  gi|262205281|ref|NR_030186.1|  Homo sapiens microRNA 52...
           14      2  gi|262205298|ref|NR_030190.1|  Homo sapiens microRNA 52...
           15      1  gi|301171394|ref|NR_035865.1|  Pan troglodytes microRNA...
           16      1  gi|262205429|ref|NR_030218.1|  Homo sapiens microRNA 51...
           17      1  gi|262205423|ref|NR_030217.1|  Homo sapiens microRNA 52...
           18      1  gi|301171401|ref|NR_035866.1|  Pan troglodytes microRNA...
           19      1  gi|270133247|ref|NR_032574.1|  Macaca mulatta microRNA ...
           20      1  gi|262205309|ref|NR_030193.1|  Homo sapiens microRNA 52...
           21      2  gi|270132717|ref|NR_032716.1|  Macaca mulatta microRNA ...
           22      2  gi|301171437|ref|NR_035870.1|  Pan troglodytes microRNA...
           23      2  gi|270133306|ref|NR_032587.1|  Macaca mulatta microRNA ...
           24      2  gi|301171428|ref|NR_035869.1|  Pan troglodytes microRNA...
           25      1  gi|301171211|ref|NR_035845.1|  Pan troglodytes microRNA...
           26      2  gi|301171153|ref|NR_035838.1|  Pan troglodytes microRNA...
           27      2  gi|301171146|ref|NR_035837.1|  Pan troglodytes microRNA...
           28      2  gi|270133254|ref|NR_032575.1|  Macaca mulatta microRNA ...
           29      2  gi|262205445|ref|NR_030221.1|  Homo sapiens microRNA 51...
           ~~~
           97      1  gi|356517317|ref|XM_003527287.1|  PREDICTED: Glycine ma...
           98      1  gi|297814701|ref|XM_002875188.1|  Arabidopsis lyrata su...
           99      1  gi|397513516|ref|XM_003827011.1|  PREDICTED: Pan panisc...
\end{verbatim}

We've just begun to scratch the surface of the object model, but you can see that
there's already some useful information. By invoking \verb|print| on the
\verb|QueryResult| object, you can see:

\begin{itemize}
\item The program name and version (blastn version 2.2.27+)
\item The query ID, description, and its sequence length (ID is 42291,
    description is `mystery\_seq', and it is 61 nucleotides long)
\item The target database to search against (refseq\_rna)
\item A quick overview of the resulting hits. For our query sequence, there are
    100 potential hits (numbered 0--99 in the table). For each hit, we can also see
    how many HSPs it contains, its ID, and a snippet of its description. Notice
    here that \verb|Bio.SearchIO| truncates the hit table overview, by showing
    only hits numbered 0--29, and then 97--99.
\end{itemize}

Now let's check our BLAT results using the same procedure as above:

%cont-doctest
\begin{verbatim}
>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> print(blat_qresult)
Program: blat (<unknown version>)
  Query: mystery_seq (61)
         <unknown description>
 Target: <unknown target>
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0     17  chr19  <unknown description>
\end{verbatim}

You'll immediately notice that there are some differences. Some of these are
caused by the way PSL format stores its details, as you'll see. The rest are
caused by the genuine program and target database differences between our BLAST
and BLAT searches:

\begin{itemize}
\item The program name and version. \verb|Bio.SearchIO| knows that the program
    is BLAT, but in the output file there is no information regarding the
    program version so it defaults to `<unknown version>'.
\item The query ID, description, and its sequence length. Notice here that these
    details are slightly different from the ones we saw in BLAST. The ID is
    `mystery\_seq' instead of 42991, there is no known description, but the query
    length is still 61. This is actually a difference introduced by the file
    formats themselves. BLAST sometimes creates its own query IDs and uses your
    original ID as the sequence description.
\item The target database is not known, as it is not stated in the BLAT output
    file.
\item And finally, the list of hits we have is completely different. Here, we
    see that our query sequence only hits the `chr19' database entry, but in it
    we see 17 HSP regions. This should not be surprising however, given that we
    are using a different program, each with its own target database.
\end{itemize}

All the details you saw when invoking the \verb|print| method can be accessed
individually using Python's object attribute access notation (a.k.a. the dot
notation). There are also other format-specific attributes that you can access
using the same method.

%cont-doctest
\begin{verbatim}
>>> print("%s %s" % (blast_qresult.program, blast_qresult.version))
blastn 2.2.27+
>>> print("%s %s" % (blat_qresult.program, blat_qresult.version))
blat <unknown version>
>>> blast_qresult.param_evalue_threshold    # blast-xml specific
10.0
\end{verbatim}

For a complete list of accessible attributes, you can check each format-specific
documentation. Here are the ones
\href{http://biopython.org/DIST/docs/api/Bio.SearchIO.BlastIO-module.html}{for BLAST}
and for
\href{http://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html}{BLAT}.

Having looked at using \verb|print| on \verb|QueryResult| objects, let's drill
down deeper. What exactly is a \verb|QueryResult|? In terms of Python objects,
\verb|QueryResult| is a hybrid between a list and a dictionary. In other words,
it is a container object with all the convenient features of lists and
dictionaries.

Like Python lists and dictionaries, \verb|QueryResult| objects are iterable.
Each iteration returns a \verb|Hit| object:

\begin{verbatim}
>>> for hit in blast_qresult:
...     hit
Hit(id='gi|262205317|ref|NR_030195.1|', query_id='42291', 1 hsps)
Hit(id='gi|301171311|ref|NR_035856.1|', query_id='42291', 1 hsps)
Hit(id='gi|270133242|ref|NR_032573.1|', query_id='42291', 1 hsps)
Hit(id='gi|301171322|ref|NR_035857.1|', query_id='42291', 2 hsps)
Hit(id='gi|301171267|ref|NR_035851.1|', query_id='42291', 1 hsps)
...
\end{verbatim}

To check how many items (hits) a \verb|QueryResult| has, you can simply invoke
Python's \verb|len| method:

%cont-doctest
\begin{verbatim}
>>> len(blast_qresult)
100
>>> len(blat_qresult)
1
\end{verbatim}

Like Python lists, you can retrieve items (hits) from a \verb|QueryResult| using
the slice notation:

%cont-doctest
\begin{verbatim}
>>> blast_qresult[0]        # retrieves the top hit
Hit(id='gi|262205317|ref|NR_030195.1|', query_id='42291', 1 hsps)
>>> blast_qresult[-1]       # retrieves the last hit
Hit(id='gi|397513516|ref|XM_003827011.1|', query_id='42291', 1 hsps)
\end{verbatim}

To retrieve multiple hits, you can slice \verb|QueryResult| objects using the
slice notation as well. In this case, the slice will return a new
\verb|QueryResult| object containing only the sliced hits:

%cont-doctest
\begin{verbatim}
>>> blast_slice = blast_qresult[:3]     # slices the first three hits
>>> print(blast_slice)
Program: blastn (2.2.27+)
  Query: 42291 (61)
         mystery_seq
 Target: refseq_rna
   Hits: ----  -----  ----------------------------------------------------------
            #  # HSP  ID + description
         ----  -----  ----------------------------------------------------------
            0      1  gi|262205317|ref|NR_030195.1|  Homo sapiens microRNA 52...
            1      1  gi|301171311|ref|NR_035856.1|  Pan troglodytes microRNA...
            2      1  gi|270133242|ref|NR_032573.1|  Macaca mulatta microRNA ...
\end{verbatim}

Like Python dictionaries, you can also retrieve hits using the hit's ID. This is
particularly useful if you know a given hit ID exists within a search query
results:

%cont-doctest
\begin{verbatim}
>>> blast_qresult['gi|262205317|ref|NR_030195.1|']
Hit(id='gi|262205317|ref|NR_030195.1|', query_id='42291', 1 hsps)
\end{verbatim}

You can also get a full list of \verb|Hit| objects using \verb|hits| and a full
list of \verb|Hit| IDs using \verb|hit_keys|:

\begin{verbatim}
>>> blast_qresult.hits
[...]       # list of all hits
>>> blast_qresult.hit_keys
[...]       # list of all hit IDs
\end{verbatim}

What if you just want to check whether a particular hit is present in the query
results? You can do a simple Python membership test using the \verb|in| keyword:

%cont-doctest
\begin{verbatim}
>>> 'gi|262205317|ref|NR_030195.1|' in blast_qresult
True
>>> 'gi|262205317|ref|NR_030194.1|' in blast_qresult
False
\end{verbatim}

Sometimes, knowing whether a hit is present is not enough; you also want to know
the rank of the hit. Here, the \verb|index| method comes to the rescue:

%cont-doctest
\begin{verbatim}
>>> blast_qresult.index('gi|301171437|ref|NR_035870.1|')
22
\end{verbatim}

Remember that we're using Python's indexing style here, which is zero-based.
This means our hit above is ranked at no. 23, not 22.

Also, note that the hit rank you see here is based on the native hit ordering
present in the original search output file. Different search tools may order
these hits based on different criteria.

If the native hit ordering doesn't suit your taste, you can use the \verb|sort|
method of the \verb|QueryResult| object. It is very similar to Python's
\verb|list.sort| method, with the addition of an option to create a new sorted
\verb|QueryResult| object or not.

Here is an example of using \verb|QueryResult.sort| to sort the hits based on
each hit's full sequence length. For this particular sort, we'll set the
\verb|in_place| flag to \verb|False| so that sorting will return a new
\verb|QueryResult| object and leave our initial object unsorted. We'll also set
the \verb|reverse| flag to \verb|True| so that we sort in descending order.

%cont-doctest
\begin{verbatim}
>>> for hit in blast_qresult[:5]:   # id and sequence length of the first five hits
...     print("%s %i" % (hit.id, hit.seq_len))
...
gi|262205317|ref|NR_030195.1| 61
gi|301171311|ref|NR_035856.1| 60
gi|270133242|ref|NR_032573.1| 85
gi|301171322|ref|NR_035857.1| 86
gi|301171267|ref|NR_035851.1| 80

>>> sort_key = lambda hit: hit.seq_len
>>> sorted_qresult = blast_qresult.sort(key=sort_key, reverse=True, in_place=False)
>>> for hit in sorted_qresult[:5]:
...     print("%s %i" % (hit.id, hit.seq_len))
...
gi|397513516|ref|XM_003827011.1| 6002
gi|390332045|ref|XM_776818.2| 4082
gi|390332043|ref|XM_003723358.1| 4079
gi|356517317|ref|XM_003527287.1| 3251
gi|356543101|ref|XM_003539954.1| 2936
\end{verbatim}

The advantage of having the \verb|in_place| flag here is that we're preserving
the native ordering, so we may use it again later. You should note that this is
not the default behavior of \verb|QueryResult.sort|, however, which is why we
needed to set the \verb|in_place| flag to \verb|True| explicitly.

At this point, you've known enough about \verb|QueryResult| objects to make it
work for you. But before we go on to the next object in the \verb|Bio.SearchIO|
model, let's take a look at two more sets of methods that could make it even
easier to work with \verb|QueryResult| objects: the \verb|filter| and \verb|map|
methods.

If you're familiar with Python's list comprehensions, generator expressions
or the built in \verb|filter| and \verb|map| functions,
you'll know how useful they are for working with list-like objects (if you're
not, check them out!). You can use these built in methods to manipulate
\verb|QueryResult| objects, but you'll end up with regular Python lists and lose
the ability to do more interesting manipulations.

That's why, \verb|QueryResult| objects provide its own flavor of
\verb|filter| and \verb|map| methods. Analogous to \verb|filter|, there are
\verb|hit_filter| and \verb|hsp_filter| methods. As their name implies, these
methods filter its \verb|QueryResult| object either on its \verb|Hit| objects
or \verb|HSP| objects. Similarly, analogous to \verb|map|, \verb|QueryResult|
objects also provide the \verb|hit_map| and \verb|hsp_map| methods. These
methods apply a given function to all hits or HSPs in a \verb|QueryResult|
object, respectively.

Let's see these methods in action, beginning with \verb|hit_filter|. This method
accepts a callback function that checks whether a given \verb|Hit| object passes
the condition you set or not. In other words, the function must accept as its
argument a single \verb|Hit| object and returns \verb|True| or \verb|False|.

Here is an example of using \verb|hit_filter| to filter out \verb|Hit| objects
that only have one HSP:

%cont-doctest
\begin{verbatim}
>>> filter_func = lambda hit: len(hit.hsps) > 1     # the callback function
>>> len(blast_qresult)      # no. of hits before filtering
100
>>> filtered_qresult = blast_qresult.hit_filter(filter_func)
>>> len(filtered_qresult)   # no. of hits after filtering
37
>>> for hit in filtered_qresult[:5]:    # quick check for the hit lengths
...     print("%s %i" % (hit.id, len(hit.hsps)))
gi|301171322|ref|NR_035857.1| 2
gi|262205330|ref|NR_030198.1| 2
gi|301171447|ref|NR_035871.1| 2
gi|262205298|ref|NR_030190.1| 2
gi|270132717|ref|NR_032716.1| 2
\end{verbatim}

\verb|hsp_filter| works the same as \verb|hit_filter|, only instead of looking
at the \verb|Hit| objects, it performs filtering on the \verb|HSP| objects in
each hits.

As for the \verb|map| methods, they too accept a callback function as their
arguments. However, instead of returning \verb|True| or \verb|False|, the
callback function must return the modified \verb|Hit| or \verb|HSP| object
(depending on whether you're using \verb|hit_map| or \verb|hsp_map|).

Let's see an example where we're using \verb|hit_map| to rename the hit IDs:

%cont-doctest
\begin{verbatim}
>>> def map_func(hit):
...     hit.id = hit.id.split('|')[3]   # renames 'gi|301171322|ref|NR_035857.1|' to 'NR_035857.1'
...     return hit
...
>>> mapped_qresult = blast_qresult.hit_map(map_func)
>>> for hit in mapped_qresult[:5]:
...     print(hit.id)
NR_030195.1
NR_035856.1
NR_032573.1
NR_035857.1
NR_035851.1
\end{verbatim}

Again, \verb|hsp_map| works the same as \verb|hit_map|, but on \verb|HSP|
objects instead of \verb|Hit| objects.

\subsection{Hit}
\label{sec:searchio-hit}

\verb|Hit| objects represent all query results from a single database entry.
They are the second-level container in the \verb|Bio.SearchIO| object hierarchy.
You've seen that they are contained by \verb|QueryResult| objects, but they
themselves contain \verb|HSP| objects.

Let's see what they look like, beginning with our BLAST search:

%doctest ../Doc/examples
\begin{verbatim}
>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> blast_hit = blast_qresult[3]    # fourth hit from the query result
>>> print(blast_hit)
Query: 42291
       mystery_seq
  Hit: gi|301171322|ref|NR_035857.1| (86)
       Pan troglodytes microRNA mir-520c (MIR520C), microRNA
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0   8.9e-20     100.47      60           [1:61]                [13:73]
          1   3.3e-06      55.39      60           [0:60]                [13:73]
\end{verbatim}

You see that we've got the essentials covered here:

\begin{itemize}
\item The query ID and description is present. A hit is always tied to a query,
    so we want to keep track of the originating query as well. These values can
    be accessed from a hit using the \verb|query_id| and
    \verb|query_description| attributes.
\item We also have the unique hit ID, description, and full sequence lengths.
    They can be accessed using \verb|id|, \verb|description|, and
    \verb|seq_len|, respectively.
\item Finally, there's a table containing quick information about the HSPs this
    hit contains. In each row, we've got the important HSP details listed: the
    HSP index, its e-value, its bit score, its span (the alignment length
    including gaps), its query coordinates, and its hit coordinates.
\end{itemize}

Now let's contrast this with the BLAT search. Remember that in the BLAT search we
had one hit with 17 HSPs.

%cont-doctest
\begin{verbatim}
>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> blat_hit = blat_qresult[0]      # the only hit
>>> print(blat_hit)
Query: mystery_seq
       <unknown description>
  Hit: chr19 (59128983)
       <unknown description>
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0         ?          ?       ?           [0:61]    [54204480:54204541]
          1         ?          ?       ?           [0:61]    [54233104:54264463]
          2         ?          ?       ?           [0:61]    [54254477:54260071]
          3         ?          ?       ?           [1:61]    [54210720:54210780]
          4         ?          ?       ?           [0:60]    [54198476:54198536]
          5         ?          ?       ?           [0:61]    [54265610:54265671]
          6         ?          ?       ?           [0:61]    [54238143:54240175]
          7         ?          ?       ?           [0:60]    [54189735:54189795]
          8         ?          ?       ?           [0:61]    [54185425:54185486]
          9         ?          ?       ?           [0:60]    [54197657:54197717]
         10         ?          ?       ?           [0:61]    [54255662:54255723]
         11         ?          ?       ?           [0:61]    [54201651:54201712]
         12         ?          ?       ?           [8:60]    [54206009:54206061]
         13         ?          ?       ?          [10:61]    [54178987:54179038]
         14         ?          ?       ?           [8:61]    [54212018:54212071]
         15         ?          ?       ?           [8:51]    [54234278:54234321]
         16         ?          ?       ?           [8:61]    [54238143:54238196]
\end{verbatim}

Here, we've got a similar level of detail as with the BLAST hit we saw earlier.
There are some differences worth explaining, though:

\begin{itemize}
\item The e-value and bit score column values. As BLAT HSPs do not have e-values
    and bit scores, the display defaults to `?'.
\item What about the span column? The span values is meant to display the
    complete alignment length, which consists of all residues and any gaps that
    may be present. The PSL format do not have this information readily available
    and \verb|Bio.SearchIO| does not attempt to try guess what it is, so we get a
    `?' similar to the e-value and bit score columns.
\end{itemize}

In terms of Python objects, \verb|Hit| behaves almost the same as Python lists,
but contain \verb|HSP| objects exclusively. If you're familiar with lists, you
should encounter no difficulties working with the \verb|Hit| object.

Just like Python lists, \verb|Hit| objects are iterable, and each iteration
returns one \verb|HSP| object it contains:

%cont-doctest
\begin{verbatim}
>>> for hsp in blast_hit:
...     hsp
HSP(hit_id='gi|301171322|ref|NR_035857.1|', query_id='42291', 1 fragments)
HSP(hit_id='gi|301171322|ref|NR_035857.1|', query_id='42291', 1 fragments)
\end{verbatim}

You can invoke \verb|len| on a \verb|Hit| to see how many \verb|HSP| objects it
has:

%cont-doctest
\begin{verbatim}
>>> len(blast_hit)
2
>>> len(blat_hit)
17
\end{verbatim}

You can use the slice notation on \verb|Hit| objects, whether to retrieve single
\verb|HSP| or multiple \verb|HSP| objects. Like \verb|QueryResult|, if you slice
for multiple \verb|HSP|, a new \verb|Hit| object will be returned containing
only the sliced \verb|HSP| objects:

%cont-doctest
\begin{verbatim}
>>> blat_hit[0]                 # retrieve single items
HSP(hit_id='chr19', query_id='mystery_seq', 1 fragments)
>>> sliced_hit = blat_hit[4:9]  # retrieve multiple items
>>> len(sliced_hit)
5
>>> print(sliced_hit)
Query: mystery_seq
       <unknown description>
  Hit: chr19 (59128983)
       <unknown description>
 HSPs: ----  --------  ---------  ------  ---------------  ---------------------
          #   E-value  Bit score    Span      Query range              Hit range
       ----  --------  ---------  ------  ---------------  ---------------------
          0         ?          ?       ?           [0:60]    [54198476:54198536]
          1         ?          ?       ?           [0:61]    [54265610:54265671]
          2         ?          ?       ?           [0:61]    [54238143:54240175]
          3         ?          ?       ?           [0:60]    [54189735:54189795]
          4         ?          ?       ?           [0:61]    [54185425:54185486]
\end{verbatim}

You can also sort the \verb|HSP| inside a \verb|Hit|, using the exact same
arguments like the sort method you saw in the \verb|QueryResult| object.

Finally, there are also the \verb|filter| and \verb|map| methods you can use
on \verb|Hit| objects. Unlike in the \verb|QueryResult| object, \verb|Hit|
objects only have one variant of \verb|filter| (\verb|Hit.filter|) and one
variant of \verb|map| (\verb|Hit.map|). Both of \verb|Hit.filter| and
\verb|Hit.map| work on the \verb|HSP| objects a \verb|Hit| has.

\subsection{HSP}
\label{sec:searchio-hsp}

\verb|HSP| (high-scoring pair) represents region(s) in the hit sequence that
contains significant alignment(s) to the query sequence. It contains the actual
match between your query sequence and a database entry. As this match is
determined by the sequence search tool's algorithms, the \verb|HSP| object
contains the bulk of the statistics computed by the search tool. This also makes
the distinction between \verb|HSP| objects from different search tools more
apparent compared to the differences you've seen in \verb|QueryResult| or
\verb|Hit| objects.

Let's see some examples from our BLAST and BLAT searches. We'll look at the
BLAST HSP first:

%doctest ../Doc/examples
\begin{verbatim}
>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> blast_hsp = blast_qresult[0][0]    # first hit, first hsp
>>> print(blast_hsp)
      Query: 42291 mystery_seq
        Hit: gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 520b (MIR520...
Query range: [0:61] (1)
  Hit range: [0:61] (1)
Quick stats: evalue 4.9e-23; bitscore 111.29
  Fragments: 1 (61 columns)
     Query - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
             |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
       Hit - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
\end{verbatim}

Just like \verb|QueryResult| and \verb|Hit|, invoking \verb|print| on an
\verb|HSP| shows its general details:
\begin{itemize}
\item There are the query and hit IDs and descriptions. We need these to
    identify our \verb|HSP|.
\item We've also got the matching range of the query and hit sequences. The
    slice notation we're using here is an indication that the range is displayed
    using Python's indexing style (zero-based, half open). The number inside the
    parenthesis denotes the strand. In this case, both sequences have the plus
    strand.
\item Some quick statistics are available: the e-value and bitscore.
\item There is information about the HSP fragments. Ignore this for now; it will
    be explained later on.
\item And finally, we have the query and hit sequence alignment itself.
\end{itemize}

These details can be accessed on their own using the dot notation, just like in
\verb|QueryResult| and \verb|Hit|:

%cont-doctest
\begin{verbatim}
>>> blast_hsp.query_range
(0, 61)
\end{verbatim}
%hack! since float display may be different across versions
\begin{verbatim}
>>> blast_hsp.evalue
4.91307e-23
\end{verbatim}

They're not the only attributes available, though. \verb|HSP| objects come with
a default set of properties that makes it easy to probe their various
details. Here are some examples:

%cont-doctest
\begin{verbatim}
>>> blast_hsp.hit_start         # start coordinate of the hit sequence
0
>>> blast_hsp.query_span        # how many residues in the query sequence
61
>>> blast_hsp.aln_span          # how long the alignment is
61
\end{verbatim}

Check out the \verb|HSP|
\href{http://biopython.org/DIST/docs/api/Bio.SearchIO._model.hsp-module.html}{documentation}
for a full list of these predefined properties.

Furthermore, each sequence search tool usually computes its own statistics /
details for its \verb|HSP| objects. For example, an XML BLAST search also
outputs the number of gaps and identical residues. These attributes can be
accessed like so:

%cont-doctest
\begin{verbatim}
>>> blast_hsp.gap_num       # number of gaps
0
>>> blast_hsp.ident_num     # number of identical residues
61
\end{verbatim}

These details are format-specific; they may not be present in other formats.
To see which details are available for a given sequence search tool, you
should check the format's documentation in \verb|Bio.SearchIO|. Alternatively,
you may also use \verb|.__dict__.keys()| for a quick list of what's available:

\begin{verbatim}
>>> blast_hsp.__dict__.keys()
['bitscore', 'evalue', 'ident_num', 'gap_num', 'bitscore_raw', 'pos_num', '_items']
\end{verbatim}

Finally, you may have noticed that the \verb|query| and \verb|hit| attributes
of our HSP are not just regular strings:

%cont-doctest
\begin{verbatim}
>>> blast_hsp.query
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG', DNAAlphabet()), id='42291', name='aligned query sequence', description='mystery_seq', dbxrefs=[])
>>> blast_hsp.hit
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG', DNAAlphabet()), id='gi|262205317|ref|NR_030195.1|', name='aligned hit sequence', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
\end{verbatim}

They are \verb|SeqRecord| objects you saw earlier in
Section~\ref{chapter:SeqRecord}! This means that you can do all sorts of
interesting things you can do with \verb|SeqRecord| objects on \verb|HSP.query|
and/or \verb|HSP.hit|.

It should not surprise you now that the \verb|HSP| object has an
\verb|alignment| property which is a \verb|MultipleSeqAlignment| object:

%cont-doctest
\begin{verbatim}
>>> print(blast_hsp.aln)
DNAAlphabet() alignment with 2 rows and 61 columns
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAG...GGG 42291
CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAG...GGG gi|262205317|ref|NR_030195.1|
\end{verbatim}

Having probed the BLAST HSP, let's now take a look at HSPs from our BLAT
results for a different kind of HSP. As usual, we'll begin by invoking
\verb|print| on it:

%cont-doctest
\begin{verbatim}
>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> blat_hsp = blat_qresult[0][0]       # first hit, first hsp
>>> print(blat_hsp)
      Query: mystery_seq <unknown description>
        Hit: chr19 <unknown description>
Query range: [0:61] (1)
  Hit range: [54204480:54204541] (1)
Quick stats: evalue ?; bitscore ?
  Fragments: 1 (? columns)
\end{verbatim}

Some of the outputs you may have already guessed. We have the query and hit IDs
and descriptions and the sequence coordinates. Values for evalue and bitscore is
`?' as BLAT HSPs do not have these attributes. But The biggest difference here
is that you don't see any sequence alignments displayed. If you look closer, PSL
formats themselves do not have any hit or query sequences, so
\verb|Bio.SearchIO| won't create any sequence or alignment objects. What happens
if you try to access \verb|HSP.query|, \verb|HSP.hit|, or \verb|HSP.aln|?
You'll get the default values for these attributes, which is \verb|None|:

%cont-doctest
\begin{verbatim}
>>> blat_hsp.hit is None
True
>>> blat_hsp.query is None
True
>>> blat_hsp.aln is None
True
\end{verbatim}

This does not affect other attributes, though. For example, you can still
access the length of the query or hit alignment. Despite not displaying any
attributes, the PSL format still have this information so \verb|Bio.SearchIO|
can extract them:

%cont-doctest
\begin{verbatim}
>>> blat_hsp.query_span     # length of query match
61
>>> blat_hsp.hit_span       # length of hit match
61
\end{verbatim}

Other format-specific attributes are still present as well:

%cont-doctest
\begin{verbatim}
>>> blat_hsp.score          # PSL score
61
>>> blat_hsp.mismatch_num   # the mismatch column
0
\end{verbatim}

So far so good? Things get more interesting when you look at another `variant'
of HSP present in our BLAT results. You might recall that in BLAT searches,
sometimes we get our results separated into `blocks'. These blocks are
essentially alignment fragments that may have some intervening sequence between
them.

Let's take a look at a BLAT HSP that contains multiple blocks to see how
\verb|Bio.SearchIO| deals with this:

%cont-doctest
\begin{verbatim}
>>> blat_hsp2 = blat_qresult[0][1]      # first hit, second hsp
>>> print(blat_hsp2)
      Query: mystery_seq <unknown description>
        Hit: chr19 <unknown description>
Query range: [0:61] (1)
  Hit range: [54233104:54264463] (1)
Quick stats: evalue ?; bitscore ?
  Fragments: ---  --------------  ----------------------  ----------------------
               #            Span             Query range               Hit range
             ---  --------------  ----------------------  ----------------------
               0               ?                  [0:18]     [54233104:54233122]
               1               ?                 [18:61]     [54264420:54264463]
\end{verbatim}

What's happening here? We still some essential details covered: the IDs and
descriptions, the coordinates, and the quick statistics are similar to what
you've seen before. But the fragments detail is all different. Instead of
showing `Fragments: 1', we now have a table with two data rows.

This is how \verb|Bio.SearchIO| deals with HSPs having multiple fragments. As
mentioned before, an HSP alignment may be separated by intervening sequences
into fragments. The intervening sequences are not part of the query-hit match,
so they should not be considered part of query nor hit sequence. However, they
do affect how we deal with sequence coordinates, so we can't ignore them.

Take a look at the hit coordinate of the HSP above. In the \verb|Hit range:| field,
we see that the coordinate is \verb|[54233104:54264463]|. But looking at the
table rows, we see that not the entire region spanned by this coordinate matches
our query. Specifically, the intervening region spans from \verb|54233122| to
\verb|54264420|.

Why then, is the query coordinates seem to be contiguous, you ask? This is
perfectly fine. In this case it means that the query match is contiguous (no
intervening regions), while the hit match is not.

All these attributes are accessible from the HSP directly, by the way:

%cont-doctest
\begin{verbatim}
>>> blat_hsp2.hit_range         # hit start and end coordinates of the entire HSP
(54233104, 54264463)
>>> blat_hsp2.hit_range_all     # hit start and end coordinates of each fragment
[(54233104, 54233122), (54264420, 54264463)]
>>> blat_hsp2.hit_span          # hit span of the entire HSP
31359
>>> blat_hsp2.hit_span_all      # hit span of each fragment
[18, 43]
>>> blat_hsp2.hit_inter_ranges  # start and end coordinates of intervening regions in the hit sequence
[(54233122, 54264420)]
>>> blat_hsp2.hit_inter_spans   # span of intervening regions in the hit sequence
[31298]
\end{verbatim}

Most of these attributes are not readily available from the PSL file we have,
but \verb|Bio.SearchIO| calculates them for you on the fly when you parse the
PSL file. All it needs are the start and end coordinates of each fragment.

What about the \verb|query|, \verb|hit|, and \verb|aln| attributes? If the
HSP has multiple fragments, you won't be able to use these attributes as they
only fetch single \verb|SeqRecord| or \verb|MultipleSeqAlignment| objects.
However, you can use their \verb|*_all| counterparts: \verb|query_all|,
\verb|hit_all|, and \verb|aln_all|. These properties will return a list containing
\verb|SeqRecord| or \verb|MultipleSeqAlignment| objects from each of the HSP
fragment. There are other attributes that behave similarly, i.e. they only work
for HSPs with one fragment. Check out the \verb|HSP| \href{http://biopython.org/DIST/docs/api/Bio.SearchIO._model.hsp-module.html}{documentation}
for a full list.

Finally, to check whether you have multiple fragments or not, you can use the
\verb|is_fragmented| property like so:

%cont-doctest
\begin{verbatim}
>>> blat_hsp2.is_fragmented     # BLAT HSP with 2 fragments
True
>>> blat_hsp.is_fragmented      # BLAT HSP from earlier, with one fragment
False
\end{verbatim}

Before we move on, you should also know that we can use the slice notation on
\verb|HSP| objects, just like \verb|QueryResult| or \verb|Hit| objects. When
you use this notation, you'll get an \verb|HSPFragment| object in return, the
last component of the object model.

\subsection{HSPFragment}
\label{sec:searchio-hspfragment}

\verb|HSPFragment| represents a single, contiguous match between the query and
hit sequences. You could consider it the core of the object model and search
result, since it is the presence of these fragments that determine whether your
search have results or not.

In most cases, you don't have to deal with \verb|HSPFragment| objects directly
since not that many sequence search tools fragment their HSPs. When you do have
to deal with them, what you should remember is that \verb|HSPFragment| objects
were written with to be as compact as possible. In most cases, they only contain
attributes directly related to sequences: strands, reading frames, alphabets,
coordinates, the sequences themselves, and their IDs and descriptions.

These attributes are readily shown when you invoke \verb|print| on an
\verb|HSPFragment|. Here's an example, taken from our BLAST search:

%doctest ../Doc/examples
\begin{verbatim}
>>> from Bio import SearchIO
>>> blast_qresult = SearchIO.read('my_blast.xml', 'blast-xml')
>>> blast_frag = blast_qresult[0][0][0]    # first hit, first hsp, first fragment
>>> print(blast_frag)
      Query: 42291 mystery_seq
        Hit: gi|262205317|ref|NR_030195.1| Homo sapiens microRNA 520b (MIR520...
Query range: [0:61] (1)
  Hit range: [0:61] (1)
  Fragments: 1 (61 columns)
     Query - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
             |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
       Hit - CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTTTAGAGGG
\end{verbatim}

At this level, the BLAT fragment looks quite similar to the BLAST fragment, save
for the query and hit sequences which are not present:

%cont-doctest
\begin{verbatim}
>>> blat_qresult = SearchIO.read('my_blat.psl', 'blat-psl')
>>> blat_frag = blat_qresult[0][0][0]    # first hit, first hsp, first fragment
>>> print(blat_frag)
      Query: mystery_seq <unknown description>
        Hit: chr19 <unknown description>
Query range: [0:61] (1)
  Hit range: [54204480:54204541] (1)
  Fragments: 1 (? columns)
\end{verbatim}

In all cases, these attributes are accessible using our favorite dot notation.
Some examples:

%cont-doctest
\begin{verbatim}
>>> blast_frag.query_start      # query start coordinate
0
>>> blast_frag.hit_strand       # hit sequence strand
1
>>> blast_frag.hit              # hit sequence, as a SeqRecord object
SeqRecord(seq=Seq('CCCTCTACAGGGAAGCGCTTTCTGTTGTCTGAAAGAAAAGAAAGTGCTTCCTTT...GGG', DNAAlphabet()), id='gi|262205317|ref|NR_030195.1|', name='aligned hit sequence', description='Homo sapiens microRNA 520b (MIR520B), microRNA', dbxrefs=[])
\end{verbatim}

\section{A note about standards and conventions}
\label{sec:searchio-standards}

Before we move on to the main functions, there is something you ought to know
about the standards \verb|Bio.SearchIO| uses. If you've worked with multiple
sequence search tools, you might have had to deal with the many different ways
each program deals with things like sequence coordinates. It might not have been
a pleasant experience as these search tools usually have their own standards.
For example, one tools might use one-based coordinates, while the other uses
zero-based coordinates. Or, one program might reverse the start and end
coordinates if the strand is minus, while others don't. In short, these often
creates unnecessary mess must be dealt with.

We realize this problem ourselves and we intend to address it in
\verb|Bio.SearchIO|. After all, one of the goals of \verb|Bio.SearchIO| is to
create a common, easy to use interface to deal with various search output files.
This means creating standards that extend beyond the object model you just saw.

Now, you might complain, "Not another standard!". Well, eventually we have to
choose one convention or the other, so this is necessary. Plus, we're not
creating something entirely new here; just adopting a standard we think is best
for a Python programmer (it is Biopython, after all).

There are three implicit standards that you can expect when working with
\verb|Bio.SearchIO|:

\begin{itemize}
\item The first one pertains to sequence coordinates. In \verb|Bio.SearchIO|,
    all sequence coordinates follows Python's coordinate style: zero-based and
    half open. For example, if in a BLAST XML output file the start and end
    coordinates of an HSP are 10 and 28, they would become 9 and 28 in
    \verb|Bio.SearchIO|. The start coordinate becomes 9 because Python indices
    start from zero, while the end coordinate remains 28 as Python slices omit
    the last item in an interval.
\item The second is on sequence coordinate orders. In \verb|Bio.SearchIO|, start
    coordinates are always less than or equal to end coordinates. This isn't
    always the case with all sequence search tools, as some of them have larger
    start coordinates when the sequence strand is minus.
\item The last one is on strand and reading frame values. For strands, there are
    only four valid choices: \verb|1| (plus strand), \verb|-1| (minus strand),
    \verb|0| (protein sequences), and \verb|None| (no strand). For reading
    frames, the valid choices are integers from \verb|-3| to \verb|3| and
    \verb|None|.
\end{itemize}

Note that these standards only exist in \verb|Bio.SearchIO| objects. If you
write \verb|Bio.SearchIO| objects into an output format, \verb|Bio.SearchIO|
will use the format's standard for the output. It does not force its standard
over to your output file.

\section{Reading search output files}
\label{sec:searchio-input}

There are two functions you can use for reading search output files into
\verb|Bio.SearchIO| objects: \verb|read| and \verb|parse|. They're essentially
similar to \verb|read| and \verb|parse| functions in other submodules like
\verb|Bio.SeqIO| or \verb|Bio.AlignIO|. In both cases, you need to supply the
search output file name and the file format name, both as Python strings. You
can check the documentation for a list of format names \verb|Bio.SearchIO|
recognizes.

\verb|Bio.SearchIO.read| is used for reading search output files with only one
query and returns a \verb|QueryResult| object. You've seen \verb|read| used in
our previous examples. What you haven't seen is that \verb|read| may also accept
additional keyword arguments, depending on the file format.

Here are some examples. In the first one, we use \verb|read| just like
previously to read a BLAST tabular output file. In the second one, we use a
keyword argument to modify so it parses the BLAST tabular variant with comments
in it:

%doctest ../Tests/Blast
\begin{verbatim}
>>> from Bio import SearchIO
>>> qresult = SearchIO.read('tab_2226_tblastn_003.txt', 'blast-tab')
>>> qresult
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)
>>> qresult2 = SearchIO.read('tab_2226_tblastn_007.txt', 'blast-tab', comments=True)
>>> qresult2
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)
\end{verbatim}

These keyword arguments differs among file formats. Check the format
documentation to see if it has keyword arguments that modifies its parser's
behavior.

As for the \verb|Bio.SearchIO.parse|, it is used for reading search output
files with any number of queries. The function returns a generator object that
yields a \verb|QueryResult| object in each iteration. Like
\verb|Bio.SearchIO.read|, it also accepts format-specific keyword arguments:

%doctest ../Tests/Blast
\begin{verbatim}
>>> from Bio import SearchIO
>>> qresults = SearchIO.parse('tab_2226_tblastn_001.txt', 'blast-tab')
>>> for qresult in qresults:
...     print(qresult.id)
gi|16080617|ref|NP_391444.1|
gi|11464971:4-101
>>> qresults2 = SearchIO.parse('tab_2226_tblastn_005.txt', 'blast-tab', comments=True)
>>> for qresult in qresults2:
...     print(qresult.id)
random_s00
gi|16080617|ref|NP_391444.1|
gi|11464971:4-101
\end{verbatim}

\section{Dealing with large search output files with indexing}
\label{sec:searchio-index}

Sometimes, you're handed a search output file containing hundreds or thousands
of queries that you need to parse. You can of course use
\verb|Bio.SearchIO.parse| for this file, but that would be grossly inefficient
if you need to access only a few of the queries. This is because \verb|parse|
will parse all queries it sees before it fetches your query of interest.

In this case, the ideal choice would be to index the file using
\verb|Bio.SearchIO.index| or \verb|Bio.SearchIO.index_db|. If the names sound
familiar, it's because you've seen them before in  Section~\ref{sec:SeqIO-index}.
These functions also behave similarly to their \verb|Bio.SeqIO| counterparts,
with the addition of format-specific keyword arguments.

Here are some examples. You can use \verb|index| with just the filename and
format name:

%doctest ../Tests/Blast
\begin{verbatim}
>>> from Bio import SearchIO
>>> idx = SearchIO.index('tab_2226_tblastn_001.txt', 'blast-tab')
>>> sorted(idx.keys())
['gi|11464971:4-101', 'gi|16080617|ref|NP_391444.1|']
>>> idx['gi|16080617|ref|NP_391444.1|']
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)
>>> idx.close()
\end{verbatim}

Or also with the format-specific keyword argument:

%cont-doctest
\begin{verbatim}
>>> idx = SearchIO.index('tab_2226_tblastn_005.txt', 'blast-tab', comments=True)
>>> sorted(idx.keys())
['gi|11464971:4-101', 'gi|16080617|ref|NP_391444.1|', 'random_s00']
>>> idx['gi|16080617|ref|NP_391444.1|']
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)
>>> idx.close()
\end{verbatim}

Or with the \verb|key_function| argument, as in \verb|Bio.SeqIO|:

%cont-doctest
\begin{verbatim}
>>> key_function = lambda id: id.upper()    # capitalizes the keys
>>> idx = SearchIO.index('tab_2226_tblastn_001.txt', 'blast-tab', key_function=key_function)
>>> sorted(idx.keys())
['GI|11464971:4-101', 'GI|16080617|REF|NP_391444.1|']
>>> idx['GI|16080617|REF|NP_391444.1|']
QueryResult(id='gi|16080617|ref|NP_391444.1|', 3 hits)
>>> idx.close()
\end{verbatim}

\verb|Bio.SearchIO.index_db| works like as \verb|index|, only it writes the
query offsets into an SQLite database file.

\section{Writing and converting search output files}
\label{sec:searchio-write}

It is occasionally useful to be able to manipulate search results from an output
file and write it again to a new file. \verb|Bio.SearchIO| provides a
\verb|write| function that lets you do exactly this. It takes as its arguments
an iterable returning \verb|QueryResult| objects, the output filename to write
to, the format name to write to, and optionally some format-specific keyword
arguments. It returns a four-item tuple, which denotes the number or
\verb|QueryResult|, \verb|Hit|, \verb|HSP|, and \verb|HSPFragment| objects that
were written.

\begin{verbatim}
>>> from Bio import SearchIO
>>> qresults = SearchIO.parse('mirna.xml', 'blast-xml')     # read XML file
>>> SearchIO.write(qresults, 'results.tab', 'blast-tab')    # write to tabular file
(3, 239, 277, 277)
\end{verbatim}

You should note different file formats require different attributes of the
\verb|QueryResult|, \verb|Hit|, \verb|HSP| and \verb|HSPFragment| objects. If
these attributes are not present, writing won't work. In other words, you can't
always write to the output format that you want. For example, if you read a
BLAST XML file, you wouldn't be able to write the results to a PSL file as PSL
files require attributes not calculated by BLAST (e.g. the number of repeat
matches). You can always set these attributes manually, if you really want to
write to PSL, though.

Like \verb|read|, \verb|parse|, \verb|index|, and \verb|index_db|, \verb|write|
also accepts format-specific keyword arguments. Check out the documentation for
a complete list of formats \verb|Bio.SearchIO| can write to and their arguments.

Finally, \verb|Bio.SearchIO| also provides a \verb|convert| function, which is
simply a shortcut for \verb|Bio.SearchIO.parse| and \verb|Bio.SearchIO.write|.
Using the convert function, our example above would be:

\begin{verbatim}
>>> from Bio import SearchIO
>>> SearchIO.convert('mirna.xml', 'blast-xml', 'results.tab', 'blast-tab')
(3, 239, 277, 277)
\end{verbatim}

As \verb|convert| uses \verb|write|, it is only limited to format conversions
that have all the required attributes. Here, the BLAST XML file provides all the
default values a BLAST tabular file requires, so it works just fine. However,
other format conversions are less likely to work since you need to manually
assign the required attributes first.