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 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288
|
%% LyX 1.3 created this file. For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{article}
\usepackage{bookman}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{geometry}
\geometry{verbose,a4paper,tmargin=20mm,bmargin=20mm}
\usepackage{graphicx}
\IfFileExists{url.sty}{\usepackage{url}}
{\newcommand{\url}{\texttt}}
\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
%% Bold symbol macro for standard LaTeX users
\newcommand{\boldsymbol}[1]{\mbox{\boldmath $#1$}}
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
\newenvironment{lyxcode}
{\begin{list}{}{
\setlength{\rightmargin}{\leftmargin}
\setlength{\listparindent}{0pt}% needed for AMS classes
\raggedright
\setlength{\itemsep}{0pt}
\setlength{\parsep}{0pt}
\normalfont\ttfamily}%
\item[]}
{\end{list}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
% header
\usepackage{fancyhdr}
\pagestyle{fancy}
\lhead{Structural Biopython FAQ}
\rhead{}
% remove date
\date{}
% make everything have section numbers
% Make links between references
\usepackage{hyperref}
\newif\ifpdf
\ifx\pdfoutput\undefined
\pdffalse
\else
\pdfoutput=1
\pdftrue
\fi
\ifpdf
\hypersetup{colorlinks=true, hyperindex=true, citecolor=red, urlcolor=blue}
\fi
\usepackage{babel}
\makeatother
\begin{document}
\title{{\huge The Biopython}\\
{\huge Structural Bioinformatics FAQ}}
\author{Thomas Hamelryck}
\author{{\normalsize Bioinformatics center}\\
{\normalsize Institute of Molecular Biology}\\
{\normalsize University of Copenhagen }\\
{\normalsize Universitetsparken 15, Bygning 10}\\
{\normalsize DK-2100 Kbenhavn }\\
{\normalsize Denmark }\\
{\normalsize thamelry@binf.ku.dk}\\
\url{http://www.binf.ku.dk/users/thamelry/}}
\maketitle
\section{Introduction}
The Biopython Project is an international association of developers
of freely available Python (\url{http://www.python.org}) tools for
computational molecular biology. Python is an object oriented, interpreted,
flexible language that is becoming increasingly popular for scientific
computing. Python is easy to learn, has a very clear syntax and can
easily be extended with modules written in C, C++ or FORTRAN.
The Biopython web site (\url{http://www.biopython.org}) provides
an online resource for modules, scripts, and web links for developers
of Python-based software for bioinformatics use and research. Basically,
the goal of biopython is to make it as easy as possible to use python
for bioinformatics by creating high-quality, reusable modules and
classes. Biopython features include parsers for various Bioinformatics
file formats (BLAST, Clustalw, FASTA, Genbank,...), access to online
services (NCBI, Expasy,...), interfaces to common and not-so-common
programs (Clustalw, DSSP, MSMS...), a standard sequence class, various
clustering modules, a KD tree data structure etc. and even documentation.
Bio.PDB is a biopython module that focuses on working with crystal
structures of biological macromolecules. This document gives a fairly
complete overview of Bio.PDB.
\section{Bio.PDB's installation}
Bio.PDB is automatically installed as part of Biopython. Biopython
can be obtained from \url{http://www.biopython.org}. It runs on many
platforms (Linux/Unix, windows, Mac,...).
However, the \verb|Bio.PDB.mmCIF.MMCIFlex| module (used internally by
\verb|Bio.PDB.MMCIFParser| to parse mmCIF files) is \emph{not} currently
installed by default. This module relies on a third party tool called
flex (fast lexical analyzer generator).
At the time of writing, in order to parse mmCIF files you'll have to
install flex, then tweak your \verb|setup.py| file and (re)install
Biopython from source.
\section{Who's using Bio.PDB?}
Bio.PDB was used in the construction of DISEMBL, a web server that
predicts disordered regions in proteins (\url{http://dis.embl.de/}),
and COLUMBA, a website that provides annotated protein structures
(\url{http://www.columba-db.de/}). Bio.PDB has also been used to
perform a large scale search for active sites similarities between
protein structures in the PDB (see \textit{Proteins Struct. Func.
Gen.}, \textbf{2003}, 51, 96-108), and to develop a new algorithm
that identifies linear secondary structure elements (\emph{BMC Bioinformatics},
\textbf{2005}, 6, 202, \url{http://www.biomedcentral.com/1471-2105/6/202}).
Judging from requests for features and information, Bio.PDB is also
used by several LPCs (Large Pharmaceutical Companies :-).
\section{Is there a Bio.PDB reference?}
Yes, and I'd appreciate it if you would refer to Bio.PDB in publications
if you make use of it. The reference is:
\begin{quote}
Hamelryck, T., Manderick, B. (2003) PDB parser and structure class
implemented in Python. \textit{Bioinformatics}, \textbf{19}, 2308-2310.
\end{quote}
The article can be freely downloaded via the Bioinformatics journal
website (\url{http://www.binf.ku.dk/users/thamelry/references.html}).
I welcome e-mails telling me what you are using Bio.PDB for. Feature
requests are welcome too.
\section{How well tested is Bio.PDB?}
Pretty well, actually. Bio.PDB has been extensively tested on nearly
5500 structures from the PDB - all structures seemed to be parsed
correctly. More details can be found in the Bio.PDB Bioinformatics
article. Bio.PDB has been used/is being used in many research projects
as a reliable tool. In fact, I'm using Bio.PDB almost daily for research
purposes and continue working on improving it and adding new features.
\section{How fast is it?}
The \texttt{PDBParser} performance was tested on about 800 structures
(each belonging to a unique SCOP superfamily). This takes about 20
minutes, or on average 1.5 seconds per structure. Parsing the structure
of the large ribosomal subunit (1FKK), which contains about 64000
atoms, takes 10 seconds on a 1000 MHz PC. In short: it's more than
fast enough for many applications.
\section{Why should I use Bio.PDB?}
Bio.PDB might be exactly what you want, and then again it might not.
If you are interested in data mining the PDB header, you might want
to look elsewhere because there is only limited support for this.
If you look for a powerful, complete data structure to access the
atomic data Bio.PDB is probably for you.
\section{Usage}
\subsection{General questions}
\subsubsection*{Importing Bio.PDB}
That's simple:
\begin{lyxcode}
from~Bio.PDB~import~{*}
\end{lyxcode}
\subsubsection*{Is there support for molecular graphics?}
Not directly, mostly since there are quite a few Python based/Python
aware solutions already, that can potentially be used with Bio.PDB.
My choice is Pymol, BTW (I've used this successfully with Bio.PDB,
and there will probably be specific PyMol modules in Bio.PDB soon/some
day). Python based/aware molecular graphics solutions include:
\begin{itemize}
\item PyMol: \url{http://pymol.sourceforge.net/}
\item Chimera: \url{http://www.cgl.ucsf.edu/chimera/}
\item PMV: \url{http://www.scripps.edu/~sanner/python/}
\item Coot: \url{http://www.ysbl.york.ac.uk/~emsley/coot/}
\item CCP4mg: \url{http://www.ysbl.york.ac.uk/~lizp/molgraphics.html}
\item mmLib: \url{http://pymmlib.sourceforge.net/}
\item VMD: \url{http://www.ks.uiuc.edu/Research/vmd/}
\item MMTK: \url{http://starship.python.net/crew/hinsen/MMTK/}
\end{itemize}
I'd be crazy to write another molecular graphics application (been
there - done that, actually :-).
\subsection{Input/output}
\subsubsection*{How do I create a structure object from a PDB file?}
First, create a \texttt{PDBParser} object:
\begin{lyxcode}
parser=PDBParser()
\end{lyxcode}
Then, create a structure object from a PDB file in the following way
(the PDB file in this case is called '1FAT.pdb', 'PHA-L' is a user
defined name for the structure):
\begin{lyxcode}
structure=parser.get\_structure('PHA-L',~'1FAT.pdb')
\end{lyxcode}
\subsubsection*{How do I create a structure object from an mmCIF file?}
Similarly to the case the case of PDB files, first create an \texttt{MMCIFParser}
object:
\begin{lyxcode}
parser=MMCIFParser()
\end{lyxcode}
Then use this parser to create a structure object from the mmCIF file:
\begin{lyxcode}
structure=parser.get\_structure('PHA-L',~'1FAT.cif')
\end{lyxcode}
\subsubsection*{...and what about the new PDB XML format?}
That's not yet supported, but I'm definitely planning to support that
in the future (it's not a lot of work). Contact me if you need this,
it might encourage me :-).
\subsubsection*{I'd like to have some more low level access to an mmCIF file...}
You got it. You can create a python dictionary that maps all mmCIF
tags in an mmCIF file to their values. If there are multiple values
(like in the case of tag \texttt{\_atom\_site.Cartn\_y}, which holds
the y coordinates of all atoms), the tag is mapped to a list of values.
The dictionary is created from the mmCIF file as follows:
\begin{lyxcode}
mmcif\_dict=MMCIF2Dict('1FAT.cif')
\end{lyxcode}
Example: get the solvent content from an mmCIF file:
\begin{lyxcode}
sc=mmcif\_dict{[}'\_exptl\_crystal.density\_percent\_sol'{]}
\end{lyxcode}
Example: get the list of the y coordinates of all atoms
\begin{lyxcode}
y\_list=mmcif\_dict{[}'\_atom\_site.Cartn\_y'{]}
\end{lyxcode}
\subsubsection*{Can I access the header information?}
Thanks to Christian Rother you can access some information from the
PDB header. Note however that many PDB files contain headers with
incomplete or erroneous information. Many of the errors have been
fixed in the equivalent mmCIF files. \emph{Hence, if you are interested
in the header information, it is a good idea to extract information
from mmCIF files using the} \texttt{\emph{MMCIF2Dict}} \emph{tool
described above, instead of parsing the PDB header. }
Now that is clarified, let's return to parsing the PDB header. The
structure object has an attribute called \texttt{header} which is
a python dictionary that maps header records to their values.
Example:
\begin{lyxcode}
resolution=structure.header{[}'resolution'{]}
keywords=structure.header{[}'keywords'{]}
\end{lyxcode}
The available keys are \texttt{name, head, deposition\_\-date, release\_\-date,
structure\_\-method, resolution, structure\_\-reference} (maps to
a list of references), \texttt{journal\_\-reference, author} and
\texttt{compound} (maps to a dictionary with various information about
the crystallized compound).
The dictionary can also be created without creating a \texttt{Structure}
object, ie. directly from the PDB file:
\begin{lyxcode}
file=open(filename,'r')
header\_dict=parse\_pdb\_header(file)
file.close()
\end{lyxcode}
\subsubsection*{Can I use Bio.PDB with NMR structures (ie. with more than one model)?}
Sure. Many PDB parsers assume that there is only one model, making
them all but useless for NMR structures. The design of the \texttt{Structure}
object makes it easy to handle PDB files with more than one model
(see section \ref{sub:The-Structure-object}).
\subsubsection*{How do I download structures from the PDB?}
This can be done using the \texttt{PDBList} object, using the \texttt{retrieve\_pdb\_file}
method. The argument for this method is the PDB identifier of the
structure.
\begin{lyxcode}
pdbl=PDBList()
pdbl.retrieve\_pdb\_file('1FAT')
\end{lyxcode}
The \texttt{PDBList} class can also be used as a command-line tool:
\begin{lyxcode}
python~PDBList.py~1fat
\end{lyxcode}
The downloaded file will be called \texttt{pdb1fat.ent} and stored
in the current working directory. Note that the \texttt{retrieve\_pdb\_file}
method also has an optional argument \texttt{pdir} that specifies
a specific directory in which to store the downloaded PDB files.
The \texttt{retrieve\_pdb\_file} method also has some options to specify
the compression format used for the download, and the program used
for local decompression (default \texttt{.Z} format and \texttt{gunzip}).
In addition, the PDB ftp site can be specified upon creation of the
\texttt{PDBList} object. By default, the RCSB PDB server (\url{ftp://ftp.rcsb.org/pub/pdb/data/structures/divided/pdb/})
is used. See the API documentation for more details. Thanks again
to Kristian Rother for donating this module.
\subsubsection*{How do I download the entire PDB?}
The following commands will store all PDB files in the \texttt{/data/pdb}
directory:
\begin{lyxcode}
python~PDBList.py~all~/data/pdb~
python~PDBList.py~all~/data/pdb~-d~
\end{lyxcode}
\noindent The API method for this is called \texttt{download\_entire\_pdb}.
Adding the \texttt{-d} option will store all files in the same directory.
Otherwise, they are sorted into PDB-style subdirectories according
to their PDB ID's. Depending on the traffic, a complete download will
take 2-4 days.
\subsubsection*{How do I keep a local copy of the PDB up-to-date?}
This can also be done using the \texttt{PDBList} object. One simply
creates a \texttt{PDBList} object (specifying the directory where
the local copy of the PDB is present) and calls the \texttt{update\_pdb}
method:
\begin{lyxcode}
pl=PDBList(pdb='/data/pdb')
pl.update\_pdb()
\end{lyxcode}
One can of course make a weekly \texttt{cronjob} out of this to keep
the local copy automatically up-to-date. The PDB ftp site can also
be specified (see API documentation).
\texttt{PDBList} has some additional methods that can be of use. The
\texttt{get\_all\_obsolete} method can be used to get a list of all
obsolete PDB entries. The \texttt{changed\_this\_week} method can
be used to obtain the entries that were added, modified or obsoleted
during the current week. For more info on the possibilities of \texttt{PDBList},
see the API documentation.
\subsubsection*{What about all those buggy PDB files?}
It is well known that many PDB files contain semantic errors (I'm
not talking about the structures themselves know, but their representation
in PDB files). Bio.PDB tries to handle this in two ways. The PDBParser
object can behave in two ways: a restrictive way and a permissive
way (THIS IS NOW THE DEFAULT). The restrictive way used to be the
default, but people seemed to think that Bio.PDB 'crashed' due to
a bug (hah!), so I changed it. If you ever encounter a real bug, please
tell me immediately!
Example:
\begin{lyxcode}
\#~Permissive~parser
parser=PDBParser(PERMISSIVE=1)
parser=PDBParser()~\#~The~same~(default)
\#~Strict~parser
strict\_parser=PDBParser(PERMISSIVE=0)
\end{lyxcode}
In the permissive state (DEFAULT), PDB files that obviously contain
errors are 'corrected' (ie. some residues or atoms are left out).
These errors include:
\begin{itemize}
\item Multiple residues with the same identifier
\item Multiple atoms with the same identifier (taking into account the altloc
identifier)
\end{itemize}
These errors indicate real problems in the PDB file (for details see
the Bioinformatics article). In the restrictive state, PDB files with
errors cause an exception to occur. This is useful to find errors
in PDB files.
Some errors however are automatically corrected. Normally each disordered
atom should have a non-blanc altloc identifier. However, there are
many structures that do not follow this convention, and have a blank
and a non-blank identifier for two disordered positions of the same
atom. This is automatically interpreted in the right way.
Sometimes a structure contains a list of residues belonging to chain
A, followed by residues belonging to chain B, and again followed by
residues belonging to chain A, i.e. the chains are 'broken'. This
is also correctly interpreted.
\subsubsection*{Can I write PDB files?}
Use the PDBIO class for this. It's easy to write out specific parts
of a structure too, of course.
Example: saving a structure
\begin{lyxcode}
io=PDBIO()
io.set\_structure(s)
io.save('out.pdb')
\end{lyxcode}
If you want to write out a part of the structure, make use of the
\texttt{Select} class (also in \texttt{PDBIO}). Select has four methods:
\begin{lyxcode}
accept\_model(model)
accept\_chain(chain)
accept\_residue(residue)
accept\_atom(atom)
\end{lyxcode}
By default, every method returns 1 (which means the model/\-chain/\-residue/\-atom
is included in the output). By subclassing \texttt{Select} and returning
0 when appropriate you can exclude models, chains, etc. from the output.
Cumbersome maybe, but very powerful. The following code only writes
out glycine residues:
\begin{lyxcode}
class~GlySelect(Select):
~~~~def~accept\_residue(self,~residue):
~~~~~~~~if~residue.get\_name()=='GLY':
~~~~~~~~~~~~return~1
~~~~~~~~else:
~~~~~~~~~~~~return~0
io=PDBIO()
io.set\_structure(s)
io.save('gly\_only.pdb',~GlySelect())
\end{lyxcode}
If this is all too complicated for you, the \texttt{Dice} module contains
a handy \texttt{extract} function that writes out all residues in
a chain between a start and end residue.
\subsubsection*{Can I write mmCIF files?}
No, and I also don't have plans to add that functionality soon (or
ever - I don't need it at all, and it's a lot of work, plus no-one
has ever asked for it). People who want to add this can contact me.
\subsection{The Structure object\label{sub:The-Structure-object}}
\subsubsection*{What's the overall layout of a Structure object?}
The \texttt{Structure} object follows the so-called \texttt{SMCRA}
(Structure/\-Model/\-Chain/\-Residue/\-Atom) architecture :
\begin{itemize}
\item A structure consists of models
\item A model consists of chains
\item A chain consists of residues
\item A residue consists of atoms
\end{itemize}
This is the way many structural biologists/bioinformaticians think
about structure, and provides a simple but efficient way to deal with
structure. Additional stuff is essentially added when needed. A UML
diagram of the \texttt{Structure} object (forget about the \texttt{Disordered}
classes for now) is shown in Fig. \ref{cap:SMCRA}.
%
\begin{figure}[tbh]
\begin{center}\includegraphics[%
width=100mm,
keepaspectratio]{images/smcra.png}\end{center}
\caption{\label{cap:SMCRA}UML diagram of SMCRA architecture of the \texttt{Structure}
object. Full lines with diamonds denote aggregation, full lines with
arrows denote referencing, full lines with triangles denote inheritance
and dashed lines with triangles denote interface realization. }
\end{figure}
\subsubsection*{How do I navigate through a Structure object?}
The following code iterates through all atoms of a structure:
\begin{lyxcode}
p=PDBParser()
structure=p.get\_structure('X',~'pdb1fat.ent')
for~model~in~structure:
~~for~chain~in~model:~
~~~~for~residue~in~chain:
~~~~~~for~atom~in~residue:
~~~~~~~~print~atom
\end{lyxcode}
There are also some shortcuts:
\begin{lyxcode}
\#~Iterate~over~all~atoms~in~a~structure
for~atom~in~structure.get\_atoms():
~~~~print~atom
\#~Iterate~over~all~residues~in~a~model
for~residue~in~model.get\_residues():
~~~~print~residue
\end{lyxcode}
Structures, models, chains, residues and atoms are called \texttt{Entities}
in Biopython. You can always get a parent \texttt{Entity} from a child
\texttt{Entity}, eg.:
\begin{lyxcode}
residue=atom.get\_parent()
chain=residue.get\_parent()
\end{lyxcode}
You can also test wether an \texttt{Entity} has a certain child using
the \texttt{has\_id} method.
\subsubsection*{Can I do that a bit more conveniently?}
You can do things like:
\begin{lyxcode}
atoms=structure.get\_atoms()
residue=structure.get\_residues()
atoms=chain.get\_atoms()
\end{lyxcode}
You can also use the \texttt{Selection.unfold\_entities} function:
\begin{lyxcode}
\#~Get~all~residues~from~a~structure
res\_list=Selection.unfold\_entities(structure,~'R')
\#~Get~all~atoms~from~a~chain
atom\_list=Selection.unfold\_entities(chain,~'A')
\end{lyxcode}
Obviously, \texttt{A=atom, R=residue, C=chain, M=model, S=structure}.
You can use this to go up in the hierarchy, eg.\ to get a list of
(unique) \texttt{Residue} or \texttt{Chain} parents from a list of
\texttt{Atoms}:
\begin{lyxcode}
residue\_list=Selection.unfold\_entities(atom\_list,~'R')
chain\_list=Selection.unfold\_entities(atom\_list,~'C')
\end{lyxcode}
For more info, see the API documentation.
\subsubsection*{How do I extract a specific \texttt{Atom/\-Residue/\-Chain/\-Model}
from a Structure?}
Easy. Here are some examples:
\begin{lyxcode}
model=structure{[}0{]}
chain=model{[}'A'{]}
residue=chain{[}100{]}
atom=residue{[}'CA'{]}
\end{lyxcode}
Note that you can use a shortcut:
\begin{lyxcode}
atom=structure{[}0{]}{[}'A'{]}{[}100{]}{[}'CA'{]}
\end{lyxcode}
\subsubsection*{What is a model id?}
The model id is an integer which denotes the rank of the model in
the PDB/mmCIF file. The model is starts at 0. Crystal structures generally
have only one model (with id 0), while NMR files usually have several
models.
\subsubsection*{What is a chain id?}
The chain id is specified in the PDB/mmCIF file, and is a single character
(typically a letter).
\subsubsection*{What is a residue id?}
This is a bit more complicated, due to the clumsy PDB format. A residue
id is a tuple with three elements:
\begin{itemize}
\item The \textbf{hetero-flag}: this is \texttt{'H\_'} plus the name of
the hetero-residue (eg. \texttt{'H\_GLC'} in the case of a glucose
molecule), or \texttt{'W'} in the case of a water molecule.
\item The \textbf{sequence identifier} in the chain, eg. 100
\item The \textbf{insertion code}, eg. 'A'. The insertion code is sometimes
used to preserve a certain desirable residue numbering scheme. A Ser
80 insertion mutant (inserted e.g. between a Thr 80 and an Asn 81
residue) could e.g. have sequence identifiers and insertion codes
as follows: Thr 80 A, Ser 80 B, Asn 81. In this way the residue numbering
scheme stays in tune with that of the wild type structure.
\end{itemize}
The id of the above glucose residue would thus be \texttt{('H\_GLC',
100, 'A')}. If the hetero-flag and insertion code are blanc, the sequence
identifier alone can be used:
\begin{lyxcode}
\#~Full~id
residue=chain{[}('~',~100,~'~'){]}
\#~Shortcut~id
residue=chain{[}100{]}
\end{lyxcode}
The reason for the hetero-flag is that many, many PDB files use the
same sequence identifier for an amino acid and a hetero-residue or
a water, which would create obvious problems if the hetero-flag was
not used.
\subsubsection*{What is an atom id?}
The atom id is simply the atom name (eg. \texttt{'CA'}). In practice,
the atom name is created by stripping all spaces from the atom name
in the PDB file.
However, in PDB files, a space can be part of an atom name. Often,
calcium atoms are called \texttt{'CA..'} in order to distinguish them
from C$\alpha$ atoms (which are called \texttt{'.CA.'}). In cases
were stripping the spaces would create problems (ie. two atoms called
\texttt{'CA'} in the same residue) the spaces are kept.
\subsubsection*{How is disorder handled?}
This is one of the strong points of Bio.PDB. It can handle both disordered
atoms and point mutations (ie. a Gly and an Ala residue in the same
position).
Disorder should be dealt with from two points of view: the atom and
the residue points of view. In general, I have tried to encapsulate
all the complexity that arises from disorder. If you just want to
loop over all C$\alpha$ atoms, you do not care that some residues
have a disordered side chain. On the other hand it should also be
possible to represent disorder completely in the data structure. Therefore,
disordered atoms or residues are stored in special objects that behave
as if there is no disorder. This is done by only representing a subset
of the disordered atoms or residues. Which subset is picked (e.g.
which of the two disordered OG side chain atom positions of a Ser
residue is used) can be specified by the user.
\textbf{Disordered atom positions} are represented by ordinary \texttt{Atom}
objects, but all \texttt{Atom} objects that represent the same physical
atom are stored in a \texttt{Disordered\-Atom} object (see Fig. \ref{cap:SMCRA}).
Each \texttt{Atom} object in a \texttt{Disordered\-Atom} object can
be uniquely indexed using its altloc specifier. The \texttt{Disordered\-Atom}
object forwards all uncaught method calls to the selected Atom object,
by default the one that represents the atom with the highest
occupancy. The user can of course change the selected \texttt{Atom}
object, making use of its altloc specifier. In this way atom disorder
is represented correctly without much additional complexity. In other
words, if you are not interested in atom disorder, you will not be
bothered by it.
Each disordered atom has a characteristic altloc identifier. You can
specify that a \texttt{Disordered\-Atom} object should behave like
the \texttt{Atom} object associated with a specific altloc identifier:
\begin{lyxcode}
atom.disordered\_select('A')~\#~select~altloc~A~atom
atom.disordered\_select('B')~\#~select~altloc~B~atom~
\end{lyxcode}
A special case arises when disorder is due to \textbf{point mutations},
i.e. when two or more point mutants of a polypeptide are present in
the crystal. An example of this can be found in PDB structure 1EN2.
Since these residues belong to a different residue type (e.g. let's
say Ser 60 and Cys 60) they should not be stored in a single \texttt{Residue}
object as in the common case. In this case, each residue is represented
by one \texttt{Residue} object, and both \texttt{Residue} objects
are stored in a single \texttt{Disordered\-Residue} object (see Fig.
\ref{cap:SMCRA}).
The \texttt{Dis\-ordered\-Residue} object forwards all un\-caught
methods to the selected \texttt{Residue} object (by default the last
\texttt{Residue} object added), and thus behaves like an ordinary
residue. Each \texttt{Residue} object in a \texttt{Disordered\-Residue}
object can be uniquely identified by its residue name. In the above
example, residue Ser 60 would have id 'SER' in the \texttt{Disordered\-Residue}
object, while residue Cys 60 would have id 'CYS'. The user can select
the active \texttt{Residue} object in a \texttt{Disordered\-Residue}
object via this id.
Example: suppose that a chain has a point mutation at position 10,
consisting of a Ser and a Cys residue. Make sure that residue 10 of
this chain behaves as the Cys residue.
\begin{lyxcode}
residue=chain{[}10{]}
residue.disordered\_select('CYS')
\end{lyxcode}
In addition, you can get a list of all \texttt{Atom} objects (ie.
all \texttt{DisorderedAtom} objects are 'unpacked' to their individual
\texttt{Atom} objects) using the \texttt{get\_unpacked\_list} method
of a \texttt{(Disordered)\-Residue} object.
\subsubsection*{Can I sort residues in a chain somehow?}
Yes, kinda, but I'm waiting for a request for this feature to finish
it :-).
\subsubsection*{How are ligands and solvent handled?}
See 'What is a residue id?'.
\subsubsection*{What about B factors?}
Well, yes! Bio.PDB supports isotropic and anisotropic B factors, and
also deals with standard deviations of anisotropic B factor if present
(see \ref{sub:Analysis}).
\subsubsection*{What about standard deviation of atomic positions?}
Yup, supported. See section \ref{sub:Analysis}.
\subsubsection*{I think the SMCRA data structure is not flexible/\-sexy/\-whatever
enough...}
Sure, sure. Everybody is always coming up with (mostly vaporware or
partly implemented) data structures that handle all possible situations
and are extensible in all thinkable (and unthinkable) ways. The prosaic
truth however is that 99.9\% of people using (and I mean really using!)
crystal structures think in terms of models, chains, residues and
atoms. The philosophy of Bio.PDB is to provide a reasonably fast,
clean, simple, but complete data structure to access structure data.
The proof of the pudding is in the eating.
Moreover, it is quite easy to build more specialised data structures
on top of the \texttt{Structure} class (eg. there's a \texttt{Polypeptide}
class). On the other hand, the \texttt{Structure} object is built
using a Parser/\-Consumer approach (called \texttt{PDBParser/\-MMCIFParser}
and \texttt{Structure\-Builder}, respectively). One can easily re-use
the PDB/mmCIF parsers by implementing a specialised \texttt{Structure\-Builder}
class. It is of course also trivial to add support for new file formats
by writing new parsers.
\subsection{\label{sub:Analysis}Analysis}
\subsubsection*{How do I extract information from an \texttt{Atom} object?}
Using the following methods:
\begin{lyxcode}
a.get\_name()~\#~atom~name~(spaces~stripped,~e.g.~'CA')
a.get\_id()~\#~id~(equals~atom~name)
a.get\_coord()~\#~atomic~coordinates
a.get\_vector()~\#~atomic~coordinates~as~Vector~object
a.get\_bfactor()~\#~isotropic~B~factor
a.get\_occupancy()~\#~occupancy
a.get\_altloc()~\#~alternative~location~specifier
a.get\_sigatm()~\#~std.~dev.~of~atomic~parameters
a.get\_siguij()~\#~std.~dev.~of~anisotropic~B~factor
a.get\_anisou()~\#~anisotropic~B~factor
a.get\_fullname()~\#~atom~name~(with~spaces,~e.g.~'.CA.')
\end{lyxcode}
\subsubsection*{How do I extract information from a \texttt{Residue} object?}
Using the following methods:
\begin{lyxcode}
r.get\_resname()~\#~return~the~residue~name~(eg.~'GLY')
r.is\_disordered()~\#~1~if~the~residue~has~disordered~atoms
r.get\_segid()~\#~return~the~SEGID
r.has\_id(name)~\#~test~if~a~residue~has~a~certain~atom
\end{lyxcode}
\subsubsection*{How do I measure distances?}
That's simple: the minus operator for atoms has been overloaded to
return the distance between two atoms.
Example:
\begin{lyxcode}
\#~Get~some~atoms
ca1=residue1{[}'CA'{]}
ca2=residue2{[}'CA'{]}
\#~Simply~subtract~the~atoms~to~get~their~distance
distance=ca1-ca2
\end{lyxcode}
\subsubsection*{How do I measure angles?}
This can easily be done via the vector representation of the atomic
coordinates, and the \texttt{calc\_angle} function from the \texttt{Vector}
module:
\begin{lyxcode}
vector1=atom1.get\_vector()
vector2=atom2.get\_vector()
vector3=atom3.get\_vector()
angle=calc\_angle(vector1,~vector2,~vector3)
\end{lyxcode}
\subsubsection*{How do I measure torsion angles?}
Again, this can easily be done via the vector representation of the
atomic coordinates, this time using the \texttt{calc\_dihedral} function
from the \texttt{Vector} module:
\begin{lyxcode}
vector1=atom1.get\_vector()
vector2=atom2.get\_vector()
vector3=atom3.get\_vector()
vector4=atom4.get\_vector()
angle=calc\_dihedral(vector1,~vector2,~vector3,~vector4)
\end{lyxcode}
\subsubsection*{How do I determine atom-atom contacts?}
Use \texttt{NeighborSearch}. This uses a KD tree data structure coded
in C++ behind the screens, so it's pretty darn fast (see \texttt{Bio.KDTree}).
\subsubsection*{How do I extract polypeptides from a \texttt{Structure} object?}
Use \texttt{PolypeptideBuilder}. You can use the resulting \texttt{Polypeptide}
object to get the sequence as a \texttt{Seq} object or to get a list
of C$\alpha$ atoms as well. Polypeptides can be built using a C-N
or a C$\alpha$-C$\alpha$ distance criterion.
Example:
\begin{lyxcode}
\#~Using~C-N~
ppb=PPBuilder()
for~pp~in~ppb.build\_peptides(structure):~
~~~~print~pp.get\_sequence()
\#~Using~CA-CA
ppb=CaPPBuilder()
for~pp~in~ppb.build\_peptides(structure):~
~~~~print~pp.get\_sequence()
\end{lyxcode}
Note that in the above case only model 0 of the structure is considered
by \texttt{PolypeptideBuilder}. However, it is possible to use \texttt{PolypeptideBuilder}
to build \texttt{Polypeptide} objects from \texttt{Model} and \texttt{Chain}
objects as well.
\subsubsection*{How do I get the sequence of a structure?}
The first thing to do is to extract all polypeptides from the structure
(see previous entry). The sequence of each polypeptide can then easily
be obtained from the \texttt{Polypeptide} objects. The sequence is
represented as a Biopython \texttt{Seq} object, and its alphabet is
defined by a \texttt{ProteinAlphabet} object.
Example:
\begin{lyxcode}
>\,{}>\,{}>~seq=polypeptide.get\_sequence()
>\,{}>\,{}>~print~seq
Seq('SNVVE...',~<class~Bio.Alphabet.ProteinAlphabet>)
\end{lyxcode}
\subsubsection*{How do I determine secondary structure?}
For this functionality, you need to install DSSP (and obtain a license
for it - free for academic use, see \url{http://www.cmbi.kun.nl/gv/dssp/}).
Then use the \texttt{DSSP} class, which maps \texttt{Residue} objects
to their secondary structure (and accessible surface area). The DSSP
codes are listed in Table \ref{cap:DSSP-codes}. Note that DSSP (the
program, and thus by consequence the class) cannot handle multiple
models!
%
\begin{table}
\subsubsection*{\begin{tabular}{|c|c|}
\hline
Code&
Secondary structure\tabularnewline
\hline
\hline
H&
$\alpha$-helix\tabularnewline
\hline
B&
Isolated $\beta$-bridge residue\tabularnewline
\hline
E&
Strand \tabularnewline
\hline
G&
3-10 helix \tabularnewline
\hline
I&
$\Pi$-helix \tabularnewline
\hline
T&
Turn\tabularnewline
\hline
S&
Bend \tabularnewline
\hline
-&
Other\tabularnewline
\hline
\end{tabular}}
\caption{\label{cap:DSSP-codes}DSSP codes in Bio.PDB.}
\end{table}
\subsubsection*{How do I calculate the accessible surface area of a residue?}
Use the \texttt{DSSP} class (see also previous entry). But see also
next entry.
\subsubsection*{How do I calculate residue depth?}
Residue depth is the average distance of a residue's atoms from the
solvent accessible surface. It's a fairly new and very powerful parameterization
of solvent accessibility. For this functionality, you need to install
Michel Sanner's MSMS program (\url{http://www.scripps.edu/pub/olson-web/people/sanner/html/msms_home.html}).
Then use the \texttt{ResidueDepth} class. This class behaves as a
dictionary which maps \texttt{Residue} objects to corresponding (residue
depth, C$\alpha$ depth) tuples. The C$\alpha$ depth is the distance
of a residue's C$\alpha$ atom to the solvent accessible surface.
Example:
\begin{lyxcode}
model=structure{[}0{]}
rd=ResidueDepth(model,~pdb\_file)
residue\_depth,~ca\_depth=rd{[}some\_residue{]}
\end{lyxcode}
You can also get access to the molecular surface itself (via the \texttt{get\_surface}
function), in the form of a Numeric python array with the surface points.
\subsubsection*{How do I calculate Half Sphere Exposure?}
Half Sphere Exposure (HSE) is a new, 2D measure of solvent exposure.
Basically, it counts the number of C$\alpha$ atoms around a residue
in the direction of its side chain, and in the opposite direction
(within a radius of 13 ). Despite its simplicity, it outperforms
many other measures of solvent exposure. An article describing this
novel 2D measure has been submitted.
HSE comes in two flavors: HSE$\alpha$ and HSE$\beta$. The former
only uses the C$\alpha$ atom positions, while the latter uses the
C$\alpha$ and C$\beta$ atom positions. The HSE measure is calculated
by the \texttt{HSExposure} class, which can also calculate the contact
number. The latter class has methods which return dictionaries that
map a \texttt{Residue} object to its corresponding HSE$\alpha$, HSE$\beta$
and contact number values.
Example:
\begin{lyxcode}
model=structure{[}0{]}
hse=HSExposure()
\#~Calculate~HSEalpha
exp\_ca=hse.calc\_hs\_exposure(model,~option='CA3')
\#~Calculate~HSEbeta
exp\_cb=hse.calc\_hs\_exposure(model,~option='CB')
\#~Calculate~classical~coordination~number~exp\_fs=hse.calc\_fs\_exposure(model)
\#~Print~HSEalpha~for~a~residue
print~exp\_ca{[}some\_residue{]}
\end{lyxcode}
\subsubsection*{How do I map the residues of two related structures onto each other?}
First, create an alignment file in FASTA format, then use the \texttt{StructureAlignment}
class. This class can also be used for alignments with more than two
structures.
\subsubsection*{How do I test if a Residue object is an amino acid?}
Use \texttt{is\_aa(residue)}.
\subsubsection*{Can I do vector operations on atomic coordinates?}
\texttt{Atom} objects return a \texttt{Vector} object representation
of the coordinates with the \texttt{get\_vector} method. \texttt{Vector}
implements the full set of 3D vector operations, matrix multiplication
(left and right) and some advanced rotation-related operations as
well. See also next question.
\subsubsection*{How do I put a virtual C$\beta$ on a Gly residue?}
OK, I admit, this example is only present to show off the possibilities
of Bio.PDB's \texttt{Vector} module (though this code is actually
used in the \texttt{HSExposure} module, which contains a novel way
to parametrize residue exposure - publication underway). Suppose that
you would like to find the position of a Gly residue's C$\beta$ atom,
if it had one. How would you do that? Well, rotating the N atom of
the Gly residue along the C$\alpha$-C bond over -120 degrees roughly
puts it in the position of a virtual C$\beta$ atom. Here's how to
do it, making use of the \texttt{rotaxis} method (which can be used
to construct a rotation around a certain axis) of the \texttt{Vector}
module:
\begin{lyxcode}
\#~get~atom~coordinates~as~vectors
n=residue{[}'N'{]}.get\_vector()~
c=residue{[}'C'{]}.get\_vector()~
ca=residue{[}'CA'{]}.get\_vector()
\#~center~at~origin
n=n-ca~
c=c-ca~
\#~find~rotation~matrix~that~rotates~n~
\#~-120~degrees~along~the~ca-c~vector
rot=rotaxis(-pi{*}120.0/180.0,~c)
\#~apply~rotation~to~ca-n~vector
cb\_at\_origin=n.left\_multiply(rot)
\#~put~on~top~of~ca~atom
cb=cb\_at\_origin+ca
\end{lyxcode}
This example shows that it's possible to do some quite nontrivial
vector operations on atomic data, which can be quite useful. In addition
to all the usual vector operations (cross (use \texttt{{*}{*}}), and
dot (use \texttt{{*}}) product, angle, norm, etc.) and the above mentioned
\texttt{rotaxis} function, the \texttt{Vector} module also has methods
to rotate (\texttt{rotmat}) or reflect (\texttt{refmat}) one vector
on top of another.
\subsection{Manipulating the structure}
\subsubsection*{How do I superimpose two structures?}
Surprisingly, this is done using the \texttt{Superimposer} object.
This object calculates the rotation and translation matrix that rotates
two lists of atoms on top of each other in such a way that their RMSD
is minimized. Of course, the two lists need to contain the same amount
of atoms. The \texttt{Superimposer} object can also apply the rotation/translation
to a list of atoms. The rotation and translation are stored as a tuple
in the \texttt{rotran} attribute of the \texttt{Superimposer} object
(note that the rotation is right multiplying!). The RMSD is stored
in the \texttt{rmsd} attribute.
The algorithm used by \texttt{Superimposer} comes from \textit{Matrix
computations, 2nd ed. Golub, G. \& Van Loan (1989)} and makes use
of singular value decomposition (this is implemented in the general
\texttt{Bio.\-SVDSuperimposer} module).
Example:
\begin{lyxcode}
sup=Superimposer()
\#~Specify~the~atom~lists
\#~'fixed'~and~'moving'~are~lists~of~Atom~objects
\#~The~moving~atoms~will~be~put~on~the~fixed~atoms
sup.set\_atoms(fixed,~moving)
\#~Print~rotation/translation/rmsd
print~sup.rotran
print~sup.rms~
\#~Apply~rotation/translation~to~the~moving~atoms
sup.apply(moving)
\end{lyxcode}
\subsubsection*{How do I superimpose two structures based on their active sites?}
Pretty easily. Use the active site atoms to calculate the rotation/translation
matrices (see above), and apply these to the whole molecule.
\subsubsection*{Can I manipulate the atomic coordinates?}
Yes, using the \texttt{transform} method of the \texttt{Atom} object,
or directly using the \texttt{set\_coord} method.
\section{Other Structural Bioinformatics modules}
\subsubsection*{Bio.SCOP}
Info coming soon.
\subsubsection*{Bio.FSSP}
Info coming soon.
\section{You haven't answered my question yet!}
Woah! It's late and I'm tired, and a glass of excellent \textit{Pedro
Ximenez} sherry is waiting for me. Just drop me a mail, and I'll answer
you in the morning (with a bit of luck...).
\section{Contributors}
The main author/maintainer of Bio.PDB is yours truly. Kristian Rother
donated code to interact with the PDB database, and to parse the PDB
header. Indraneel Majumdar sent in some bug reports and assisted in
coding the \texttt{Polypeptide} module. Many thanks to Brad Chapman,
Jeffrey Chang, Andrew Dalke and Iddo Friedberg for suggestions, comments,
help and/or biting criticism :-).
\section{Can I contribute?}
Yes, yes, yes! Just send me an e-mail (thamelry@binf.ku.dk) if you
have something useful to contribute! Eternal fame awaits!
\section{Biopython License Agreement}
Permission to use, copy, modify, and distribute this software and
its documentation with or without modifications and for any purpose
and without fee is hereby granted, provided that any copyright notices
appear in all copies and that both those copyright notices and this
permission notice appear in supporting documentation, and that the
names of the contributors or copyright holders not be used in advertising
or publicity pertaining to distribution of the software without specific
prior permission.
THE CONTRIBUTORS AND COPYRIGHT HOLDERS OF THIS SOFTWARE DISCLAIM ALL
WARRANTIES WITH REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES
OF MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL THE CONTRIBUTORS
OR COPYRIGHT HOLDERS BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA
OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE
OF THIS SOFTWARE.
\end{document}
|