Skip to content

Commit 6e82d20

Browse files
trayc7trayc7
committed
added Biopython
1 parent ddfae4c commit 6e82d20

2 files changed

Lines changed: 174 additions & 7 deletions

File tree

‎_episodes/05-seaborn-viz.md‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -60,7 +60,7 @@ surveys_complete = pd.read_csv('surveys_complete.csv', index_col=0)
6060

6161
## A Simple Scatterplot
6262

63-
Let's start with a basic scatterplot. We'll simply plot the weight on the x-axis and the hind foot length on the y-axis. This uses the seaborn function `.lmplot()`. This function can take a Pandas DataFrame directly. It also will fit a regression line, by default. Since we may not want to visualize these data with a regression line, we will use the `fit_reg=False` argument.
63+
Let's start with a basic scatterplot. We'll simply plot the weight on the horizontal axis and the hind foot length on the vertical axis. This uses the seaborn function `.lmplot()`. This function can take a Pandas DataFrame directly. It also will fit a regression line, by default. Since we may not want to visualize these data with a regression line, we will use the `fit_reg=False` argument.
6464

6565
~~~
6666
sns.lmplot("weight", "hindfoot_length", data=surveys_complete, fit_reg=False)
@@ -162,7 +162,7 @@ my_fig.set_axis_labels('Weight (g)', 'Hindfoot Length (mm)')
162162
> Remember how to access subsets of a DataFrame based on conditional criteria?
163163
> Plot the scatter plot above for only plot number `12` and color by `sex`. (Make the markers larger circles.)
164164
>
165-
<!-- > > ## Solution
165+
> > ## Solution
166166
> >
167167
> > ~~~
168168
> > my_fig = sns.lmplot("weight", "hindfoot_length", data=surveys_complete[surveys_complete.plot_id == 12],
@@ -173,7 +173,7 @@ my_fig.set_axis_labels('Weight (g)', 'Hindfoot Length (mm)')
173173
> > {: .python}
174174
> >
175175
> > ![png](../fig/05-seaborn-scatter-9.png)
176-
> {: .solution} -->
176+
> {: .solution}
177177
{: .challenge}
178178
179179
@@ -205,7 +205,7 @@ ax.set(xlabel='Species ID', ylabel='Hindfoot Length (mm)')
205205
206206
Note that for the `.boxplot()` function, we can simply set the x and y values by giving the column names from our DataFrame. Then we must use the `data=surveys_complete` argument to indicate our DataFrame object.
207207
208-
Now let's use a violin plot to visualize the weight data. For these data, we would also like the weight to be on the log10 scale. We will make this plot horizontal, so we'll set the x-axis to be on the log scale.
208+
Now let's use a violin plot to visualize the weight data. For these data, we would also like the weight to be on the log10 scale. We will make this plot horizontal, so we'll set the horizontal axis to be on the log scale.
209209
An easy way to do this is to create a new graph variable from our `.violinplot()` function and then call the `.set_xscale()` that is a member method of the graph variable.
210210
211211
~~~
@@ -224,7 +224,7 @@ ax.set(ylabel='Species ID', xlabel='Weight (g)')
224224
> measurements for different sexed animals from a single
225225
> species, [*Onychomys leucogaster* (OL)](https://en.wikipedia.org/wiki/Northern_grasshopper_mouse), one of the coolest rodent species:
226226
>
227-
<!-- > > ## Solution
227+
> > ## Solution
228228
> >
229229
> > ~~~
230230
> > fig, ax = plt.subplots(figsize=plot_dims)
@@ -233,7 +233,7 @@ sns.violinplot(x = 'sex', y = 'weight', data=surveys_complete[surveys_complete.s
233233
> > {: .python}
234234
> >
235235
> > ![png](../fig/05-seaborn-violinplot-2.png)
236-
> {: .solution} -->
236+
> {: .solution}
237237
{: .challenge}
238238
239239
# Histograms
@@ -248,7 +248,7 @@ sns.distplot(surveys_complete['weight'], color='c')
248248
249249
![png](../fig/05-seaborn-distplot-1.png)
250250
251-
By default, the `.distplot()` function plots the histogram as a density plot with the kernel density estimate overlaid as a darker line on the graph. This may not be necessary and we may instead want the y-axis to represent the counts in each bin. We can further adjust the appearance of the histogram to make the bars more apparent and increase the number of bins.
251+
By default, the `.distplot()` function plots the histogram as a density plot with the kernel density estimate overlaid as a darker line on the graph. This may not be necessary and we may instead want the vertical axis to represent the counts in each bin. We can further adjust the appearance of the histogram to make the bars more apparent and increase the number of bins.
252252
253253
~~~
254254
fig, ax = plt.subplots(figsize=plot_dims)

‎_episodes/06-biopython.md‎

Lines changed: 167 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -373,7 +373,174 @@ Seq('ACACAAATTCTAACTGGCCTCCTACTGGCCGCCCACTACACTGCAGACACAACC...AGC', IUPACAmbiguo
373373
~~~
374374
{: .output}
375375
376+
## BLAST
376377
378+
BioPython makes it easy to work with NCBI's BLAST. To run
379+
blast over the internet, we can use `qblast()`.
380+
For this we must import the `NCBIWWW` module:
377381
382+
~~~
383+
from Bio.Blast import NCBIWWW
384+
~~~
385+
{: .python}
386+
387+
You can call the `help()` function on `NCBIWWW.qblast` to inspect how this
388+
function works. This will return all of the parameters of `qblast` so that
389+
you can understand how to specify your query correctly.
390+
~~~
391+
help(NCBIWWW.qblast)
392+
~~~
393+
{: .python}
394+
395+
396+
Next we can read in a sequence that is stored in a FASTA file:
397+
~~~
398+
query = SeqIO.read("test.fasta", format="fasta")
399+
~~~
400+
{: .python}
401+
402+
The variable we created called `query` is a sequence stored in a `SeqRecord`
403+
object.
404+
405+
<!-- ~~~
406+
query.description
407+
~~~
408+
{: .python}
409+
-->
410+
411+
To run a BLAST search on the sequence from our FASTA file, we simply
412+
have to specify the search program (`blastn`) and the database (`nt`).
413+
The last argument is the `Seq` object stored in our `query`.
414+
~~~
415+
result_handle = NCBIWWW.qblast("blastn", "nt", query.seq)
416+
~~~
417+
{: .python}
418+
419+
Note that this might not work for everyone in class. It is possible
420+
for NCBI to throttle non-interactive users.
421+
422+
Once we have the results of our BLAST, we can store them in an XML file.
423+
~~~
424+
blast_file = open("my_blast.xml", "w")
425+
blast_file.write(result_handle.read())
426+
~~~
427+
{: .python}
428+
429+
Once we have stored the results, it's best to close all of our open
430+
file handles:
431+
~~~
432+
blast_file.close()
433+
result_handle.close()
434+
~~~
435+
{: .python}
436+
437+
438+
We created an XML file containing our BLAST results. This is now easy to
439+
parse using the `NCBIXML` tools:
440+
~~~
441+
from Bio.Blast import NCBIXML
442+
handle = open("my_blast.xml")
443+
blast_record = NCBIXML.read(handle)
444+
~~~
445+
{: .python}
446+
447+
448+
Now that we have read in the file, we can
449+
print each of the hits:
450+
~~~
451+
for hit in blast_record.descriptions:
452+
print(hit.title)
453+
print(hit.e)
454+
~~~
455+
{: .python}
456+
457+
~~~
458+
gi|1105484513|ref|XM_002284686.3| PREDICTED: Vitis vinifera cold-regulated 413 plasma membrane protein 2 (LOC100248690), mRNA
459+
0.0
460+
gi|1420088022|gb|MG722853.1| Vitis vinifera cold-regulated 413 inner membrane protein 2 mRNA, complete cds
461+
0.0
462+
gi|123704572|emb|AM483681.1| Vitis vinifera, whole genome shotgun sequence, contig VV78X045699.9, clone ENTAV 115
463+
0.0
464+
465+
...
466+
467+
gi|1027107741|ref|XM_008238505.2| PREDICTED: Prunus mume cold-regulated 413 plasma membrane protein 1-like (LOC103335494), mRNA
468+
1.04564e-118
469+
~~~
470+
{: .output}
471+
472+
We can also view the alignments for each of the BLAST hits:
473+
~~~
474+
for hit in blast_record.alignments:
475+
for hsp in hit.hsps:
476+
print(hit.title)
477+
print(hsp.expect)
478+
print(hsp.query[0:75] + '...')
479+
print(hsp.match[0:75] + '...')
480+
print(hsp.sbjct[0:75] + '...')
481+
~~~
482+
{: .python}
483+
484+
~~~
485+
gi|1105484513|ref|XM_002284686.3| PREDICTED: Vitis vinifera cold-regulated 413 plasma membrane protein 2 (LOC100248690), mRNA
486+
0.0
487+
TACTCTACAGTCTCTGACTTTGTAAGCTTCGCGCTTCTTCTCCTTTTTCTCTCTGGGGAAAGATTTTCCCTTTCT...
488+
|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||...
489+
TACTCTACAGTCTCTGACTTTGTAAGCTTCGCGCTTCTTCTCCTTTTTCTCTCTGGGGAAAGATTTTCCCTTTCT...
490+
491+
...
492+
493+
gi|1027107741|ref|XM_008238505.2| PREDICTED: Prunus mume cold-regulated 413 plasma membrane protein 1-like (LOC103335494), mRNA
494+
1.04564e-118
495+
ATTGAAGATGGGGAAAAAGGGTTACTTGGCGATGAGGACTGACACTGATACTACTGATTTGATCAGTTCTGATCT...
496+
|||| |||||| ||| || | |||||| |||| |||||| | ||| | | ||| |||||| || |||||...
497+
ATTGGAGATGGCAAAACAGAGCTACTTGAAAATGATGACTGAACCAGATGCAAATGAATTGATCCACTCCGATCT...
498+
~~~
499+
{: .output}
378500

501+
Often a BLAST search will return many matches for a single query, as is the
502+
case for this example. This is why it is best to save them in an XML file.
503+
Using `NCBIXML.parse()` enables us to evaluate each of the BLAST records.
504+
We can specify a threshold so that we can easily inspect the closest
505+
matches.
506+
~~~
507+
E_VALUE_THRESH = 1e-150
508+
for record in NCBIXML.parse(open("my_blast.xml")):
509+
for align in record.alignments:
510+
for hsp in align.hsps:
511+
if hsp.expect < E_VALUE_THRESH:
512+
print("MATCH: %s " % align.title[:60])
513+
print(hsp.expect)
514+
~~~
515+
{: .python}
516+
517+
~~~
518+
MATCH: gi|1105484513|ref|XM_002284686.3| PREDICTED: Vitis vinifera
519+
0.0
520+
MATCH: gi|1420088022|gb|MG722853.1| Vitis vinifera cold-regulated 4
521+
0.0
522+
MATCH: gi|123704572|emb|AM483681.1| Vitis vinifera, whole genome sh
523+
0.0
524+
MATCH: gi|1217007653|ref|XM_021787586.1| PREDICTED: Hevea brasilien
525+
9.7761e-151
526+
MATCH: gi|1217007651|ref|XM_021787585.1| PREDICTED: Hevea brasilien
527+
9.7761e-151
528+
~~~
529+
{: .output}
530+
531+
532+
Our BLAST search has matched our sequence with _Vitis vinifera_.
533+
Let's check to see if it got it right:
534+
~~~
535+
query.description
536+
~~~
537+
{: .python}
538+
539+
540+
<!--
541+
~~~
542+
543+
~~~
544+
{: .python}
545+
-->
379546

0 commit comments

Comments
 (0)