Skip to content

Commit ddfae4c

Browse files
trayc7trayc7
committed
added biopy
1 parent 914fbec commit ddfae4c

2 files changed

Lines changed: 379 additions & 6 deletions

File tree

‎_episodes/06-biopython.md‎

Lines changed: 379 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,379 @@
1+
---
2+
title: Introduction to Biopython
3+
teaching: 1
4+
exercises: 0
5+
questions:
6+
- "What does Biopython do?"
7+
- "How does Biopython handle sequences?"
8+
- "How can I access sequences and data from Genbank?"
9+
objectives:
10+
- "Learn about the `Seq` and `SeqRecord` objects."
11+
- "Read in sequences from FASTA files."
12+
- "Download a sequence record directly from Genbank using the NCBI E-utilities."
13+
keypoints:
14+
- "Biopython is a very useful toolbox for working with sequence data."
15+
---
16+
17+
18+
# Biopython Background
19+
20+
Biopython is a freely available package for working with molecular biological data.
21+
In this lesson, we will just cover some basics of working with Biopython.
22+
The developers of this package have written a comprehensive [tutorial and cookbook](http://biopython.org/DIST/docs/tutorial/Tutorial.html).
23+
24+
The tutorial we are working with today was written by [Dr. Iddo Friedberg](https://iddo-friedberg.net/) and [Dr. Stuart Brown](https://scholar.google.com/citations?user=ig0QSzAAAAAJ&hl=en)
25+
26+
## What can Biopython do?
27+
28+
The [documentation page for Biopython](http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc3) provides a list of the many
29+
different tools in the package:
30+
31+
- The ability to parse bioinformatics files into Python utilizable data structures, including support for the following formats:
32+
- Blast output – both from standalone and WWW Blast
33+
- Clustalw
34+
- FASTA
35+
- GenBank
36+
- PubMed and Medline
37+
- ExPASy files, like Enzyme and Prosite
38+
- SCOP, including ‘dom’ and ‘lin’ files
39+
- UniGene
40+
- SwissProt
41+
- Files in the supported formats can be iterated over record by record or indexed and accessed via a Dictionary interface.
42+
- Code to deal with popular on-line bioinformatics destinations such as:
43+
- NCBI – Blast, Entrez and PubMed services
44+
- ExPASy – Swiss-Prot and Prosite entries, as well as Prosite searches
45+
- Interfaces to common bioinformatics programs such as:
46+
- Standalone Blast from NCBI
47+
- Clustalw alignment program
48+
- EMBOSS command line tools
49+
- A standard sequence class that deals with sequences, ids on sequences, and sequence features.
50+
- Tools for performing common operations on sequences, such as translation, transcription and weight calculations.
51+
- Code to perform classification of data using k Nearest Neighbors, Naive Bayes or Support Vector Machines.
52+
- Code for dealing with alignments, including a standard way to create and deal with substitution matrices.
53+
- Code making it easy to split up parallelizable tasks into separate processes.
54+
- GUI-based programs to do basic sequence manipulations, translations, BLASTing, etc.
55+
- Extensive documentation and help with using the modules, including this file, on-line wiki documentation, the web site, and the mailing list.
56+
- Integration with BioSQL, a sequence database schema also supported by the BioPerl and BioJava projects.
57+
58+
# Getting Started
59+
60+
## Download Example Files
61+
62+
This lesson will use the example files in the
63+
[`biopython`](https://github.com/EEOB-BioData/BCB546X-Fall2019/tree/master/course-files/biopython)
64+
folder of the course
65+
files in the `BCB546X-Fall2019` repository.
66+
Download these files and make sure they are in the same directory
67+
where you are creating your Jupyter notebook.
68+
69+
70+
71+
## Install Biopython and Create a Jupyter Notebook
72+
73+
The easiest way to install the Biopython tools is to use `pip`. From your terminal, you simply need to execute the following:
74+
75+
```
76+
$ pip install biopython
77+
```
78+
79+
Now create a new Jupyter notebook for this lesson.
80+
81+
82+
# Working with Biopython
83+
84+
## The `Seq` Object
85+
86+
The `Seq` object class is simple and fundamental for a lot of
87+
Biopython work. A Seq object can contain DNA, RNA, or protein.
88+
It contains a string (the sequence) and a defined alphabet for that string.
89+
The alphabets are actually defined objects such as `IUPACAmbiguousDNA` or
90+
`IUPACProtein`. A Seq object with a DNA alphabet has some different methods than one with an Amino Acid alphabet.
91+
92+
First, import the `Seq` and alphabet from Biopython
93+
94+
~~~
95+
from Bio.Seq import Seq
96+
from Bio.Alphabet import IUPAC
97+
~~~
98+
{: .python}
99+
100+
Now we can create a `Seq` object:
101+
102+
~~~
103+
my_seq = Seq("AGTACACTGGT", IUPAC.unambiguous_dna)
104+
my_seq
105+
~~~
106+
{: .python}
107+
108+
~~~
109+
Seq('AGTACACTGGT', IUPACUnambiguousDNA())
110+
~~~
111+
{: .output}
112+
113+
We can create protein sequence by specifying the alphabet:
114+
~~~
115+
my_prot = Seq("AGTACACTGGT", IUPAC.protein)
116+
my_prot
117+
~~~
118+
{: .python}
119+
120+
~~~
121+
Seq('AGTACACTGGT', IUPACProtein())
122+
~~~
123+
{: .output}
124+
125+
The nice thing about the sequence object is
126+
that it can be treated
127+
just like a Python string object.
128+
129+
~~~
130+
print(my_seq[:3])
131+
~~~
132+
{: .python}
133+
~~~
134+
AGT
135+
~~~
136+
{: .output}
137+
138+
`Seq` objects also have string methods like `.count()`
139+
~~~
140+
my_seq.count('AC')
141+
~~~
142+
{: .python}
143+
~~~
144+
2
145+
~~~
146+
{: .output}
147+
148+
And you can use functions that act on strings like `len()`
149+
150+
~~~
151+
len(my_seq)
152+
~~~
153+
{: .python}
154+
~~~
155+
11
156+
~~~
157+
{: .output}
158+
159+
160+
`Seq` objects also have special methods. For example, you can get the reverse complement of a sequence:
161+
162+
~~~
163+
my_seq = Seq("GATCGATGGGCCTATATAGGATCGAAAATCGC", IUPAC.unambiguous_dna)
164+
print(my_seq.reverse_complement())
165+
~~~
166+
{: .python}
167+
~~~
168+
GCGATTTTCGATCCTATATAGGCCCATCGATC
169+
~~~
170+
{: .output}
171+
172+
Just like strings in Python, the `Seq` object is immutable,
173+
meaning you cannot change it. If you try to change one of the sites in this
174+
sequence, you will get an error. If you want an editable sequence object, you
175+
will need to create a `MutableSeq` object.
176+
177+
~~~
178+
from Bio.Seq import MutableSeq
179+
mutable_seq = MutableSeq("GCCATTGTAATGGGCCGCTGAAAGGGTGCCCGA", IUPAC.unambiguous_dna)
180+
~~~
181+
{: .python}
182+
183+
Now you can try changing the nucleotide at index 3 to `'G'`.
184+
185+
## The `SeqRecord` Object
186+
187+
Biopython's `SeqRecord` is a complex object that contains a
188+
`Seq` object as well as other fields for attributes of that
189+
sequence (i.e., metadata). These attributes are also called
190+
"annotation fields":
191+
192+
- `.seq`
193+
- The sequence itself, typically a Seq object.
194+
- `.id`
195+
- The primary ID used to identify the sequence – a string. In most cases this is something like an accession number.
196+
- `.name`
197+
- A “common” name/id for the sequence – a string. In some cases this will be the same as the accession number, but it could also be a clone name. I think of this as being analogous to the LOCUS id in a GenBank record.
198+
- `.description`
199+
- A human readable description or expressive name for the sequence – a string.
200+
- `.letter_annotations`
201+
- Holds per-letter-annotations using a (restricted) dictionary of additional information about the letters in the sequence. The keys are the name of the information, and the information is contained in the value as a Python sequence (i.e. a list, tuple or string) with the same length as the sequence itself. This is often used for quality scores (e.g. Section 20.1.6) or secondary structure information (e.g. from Stockholm/PFAM alignment files).
202+
- `.annotations`
203+
- A dictionary of additional information about the sequence. The keys are the name of the information, and the information is contained in the value. This allows the addition of more “unstructured” information to the sequence.
204+
- `.features`
205+
- A list of SeqFeature objects with more structured information about the features on a sequence (e.g. position of genes on a genome, or domains on a protein sequence)
206+
- `.dbxrefs`
207+
- A list of database cross-references as strings.
208+
209+
You can create a `SeqRecord` by giving the constructor a `Seq` object:
210+
~~~
211+
from Bio.SeqRecord import SeqRecord
212+
simple_seq = Seq("GATC")
213+
simple_seq_r = SeqRecord(simple_seq)
214+
~~~
215+
{: .python}
216+
217+
And you can provide attributes:
218+
219+
~~~
220+
simple_seq_r.id = "AC12345"
221+
simple_seq_r.description = "This sequence is pretend."
222+
print(simple_seq_r)
223+
~~~
224+
{: .python}
225+
~~~
226+
ID: AC12345
227+
Name: <unknown name>
228+
Description: This sequence is pretend.
229+
Number of features: 0
230+
Seq('GATC')
231+
~~~
232+
{: .output}
233+
234+
## Reading Sequences from FASTA files
235+
236+
`SeqIO` enables reading in sequences from FASTA files and storing the data in a `SeqRecord`. Addtionally `SeqIO` provides tools for writing sequence data to a file.
237+
238+
We will read in the example file [`NC_005816.fna`](https://github.com/EEOB-BioData/BCB546X-Fall2019/blob/master/course-files/biopython/NC_005816.fna) using `SeqIO`.
239+
240+
~~~
241+
from Bio import SeqIO
242+
record = SeqIO.read("NC_005816.fna", "fasta")
243+
record
244+
~~~
245+
{: .python}
246+
<!-- ~~~
247+
SeqRecord(seq=Seq('TGTAACGAACGGTGCAATAGTGATCCACACCCAACGCCTGAAATCAGATCCAGG...CTG', SingleLetterAlphabet()), id='gi|45478711|ref|NC_005816.1|', name='gi|45478711|ref|NC_005816.1|', description='gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence', dbxrefs=[])
248+
~~~
249+
{: .output}
250+
-->
251+
252+
> ## Find out more about this sequence
253+
>
254+
> Use string methods and the `SeqRecord` attributes to get the length of the sequence
255+
> and the species name.
256+
>
257+
<!-- > > ## Solution
258+
> >
259+
> > Get the length using `len()`.
260+
> > ~~~
261+
> > len(record.seq)
262+
> > ~~~
263+
> > {: .python}
264+
> > ~~~
265+
> > 9609
266+
> > ~~~
267+
> > {: .output}
268+
> >
269+
> > The species name is given in the description of this FASTA file
270+
> > ~~~
271+
> > record.description
272+
> > ~~~
273+
> > {: .python}
274+
> > ~~~
275+
> > 'gi|45478711|ref|NC_005816.1| Yersinia pestis biovar Microtus str. 91001 plasmid pPCP1, complete sequence'
276+
> > ~~~
277+
> > {: .output}
278+
> >
279+
> {: .solution} -->
280+
{: .challenge}
281+
282+
Using `SeqIO` we can read in several sequences from a file and store
283+
them in a list of `SeqRecord` objects from a file. The file
284+
[`example.fasta`](https://github.com/EEOB-BioData/BCB546X-Fall2019/blob/master/course-files/biopython/example.fasta) looks like this:
285+
286+
```
287+
>EAS54_6_R1_2_1_413_324
288+
CCCTTCTTGTCTTCAGCGTTTCTCC
289+
>EAS54_6_R1_2_1_540_792
290+
TTGGCAGGCCAAGGCCGATGGATCA
291+
>EAS54_6_R1_2_1_443_348
292+
GTTGCTTCTGGCGTGGGTGGGGGGG
293+
```
294+
295+
With Biopython, we can use the `SeqIO.parse()` function to obtain the three sequences in this file
296+
297+
~~~
298+
handle = open("example.fasta", "r")
299+
seq_list = list(SeqIO.parse(handle, "fasta"))
300+
handle.close()
301+
print(seq_list[0].seq)
302+
~~~
303+
{: .python}
304+
~~~
305+
CCCTTCTTGTCTTCAGCGTTTCTCC
306+
~~~
307+
{: .output}
308+
309+
In the example above, we open the file and assign it to the variable `handle`
310+
which acts as a pointer to the file contents.
311+
312+
## Direct Access to GenBank
313+
314+
BioPython has modules that can directly access databases over the Internet using
315+
the `Entrez` module. This uses the NCBI Efetch service,
316+
which works on many NCBI databases including protein and PubMed
317+
literature citations.
318+
With a few tweaks, this code could be used to download a list of
319+
GenBank ID’s and save them as FASTA or GenBank files.
320+
321+
Before using the online NCBI resources, it is important to be aware of the user
322+
requirements. If you abuse their system (whether on purpose or on accident),
323+
they will block your access for some time.
324+
You can find the requirements in the
325+
[NCBI E-utilities guide](https://www.ncbi.nlm.nih.gov/books/NBK25497/#chapter2.Usage_Guidelines_and_Requiremen).
326+
327+
First, you are required to provide NCBI with
328+
your identity so that you can be contacted
329+
if there is a problem. This also limits abuse of this system so that
330+
their servers aren't overwhelmed.
331+
If you are identified as someone who is excessively using the E-utilities,
332+
NCBI will contact you before you are blocked.
333+
334+
The quote below from the
335+
[NCBI guide](https://www.ncbi.nlm.nih.gov/books/NBK25497/#chapter2.Usage_Guidelines_and_Requiremen)
336+
gives you a sense of what constitutes appropriate usage of
337+
the E-utility servers:
338+
339+
>In order not to overload the E-utility servers, NCBI recommends that users post no more than three URL requests per second and limit large jobs to either weekends or between 9:00 PM and 5:00 AM Eastern time during weekdays.
340+
341+
342+
Enter _your own email address_ in place of `<enter your email address>`:
343+
~~~
344+
from Bio import Entrez
345+
Entrez.email = "<enter your email address>"
346+
~~~
347+
{: .python}
348+
349+
Now we can fetch a Genbank record:
350+
~~~
351+
handle = Entrez.efetch(db="nucleotide", id="DQ137224", rettype="gb", retmode="text")
352+
record = SeqIO.read(handle, "genbank")
353+
print(record)
354+
~~~
355+
{: .python}
356+
~~~
357+
ID: DQ137224.1
358+
Name: DQ137224
359+
Description: Megadyptes antipodes voucher JD64A cytochrome b (cytb) gene, partial cds; mitochondrial
360+
Number of features: 3
361+
/molecule_type=DNA
362+
/topology=linear
363+
/data_file_division=VRT
364+
/date=26-JUL-2016
365+
/accessions=['DQ137224']
366+
/sequence_version=1
367+
/keywords=['']
368+
/source=mitochondrion Megadyptes antipodes (Yellow-eyed penguin)
369+
/organism=Megadyptes antipodes
370+
/taxonomy=['Eukaryota', 'Metazoa', 'Chordata', 'Craniata', 'Vertebrata', 'Euteleostomi', 'Archelosauria', 'Archosauria', 'Dinosauria', 'Saurischia', 'Theropoda', 'Coelurosauria', 'Aves', 'Neognathae', 'Sphenisciformes', 'Spheniscidae', 'Megadyptes']
371+
/references=[Reference(title='Multiple gene evidence for expansion of extant penguins out of Antarctica due to global cooling', ...), Reference(title='Direct Submission', ...)]
372+
Seq('ACACAAATTCTAACTGGCCTCCTACTGGCCGCCCACTACACTGCAGACACAACC...AGC', IUPACAmbiguousDNA())
373+
~~~
374+
{: .output}
375+
376+
377+
378+
379+

‎_includes/syllabus.html‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -40,12 +40,6 @@ <h2 id="schedule">Schedule</h2>
4040
{% endfor %}
4141
{% assign hours = current | divided_by: 60 %}
4242
{% assign minutes = current | modulo: 60 %}
43-
<tr>
44-
{% if multiday %}<td></td>{% endif %}
45-
<td class="col-md-1"></td>
46-
<td class="col-md-3">BioPython</td>
47-
<td class="col-md-7"></td>
48-
</tr>
4943
</table>
5044

5145

0 commit comments

Comments
 (0)