This notebook is associated with the following accompanied reading:

- Yandell, Mark, and Daniel Ence. 2012. “A Beginner’s Guide to Eukaryotic Genome Annotation.” Nature Reviews Genetics 13 (5): 329–42. https://doi.org/10.1038/nrg3174.

- The Python Data Science Handbook
**Chapter 2**https://jakevdp.github.io/PythonDataScienceHandbook/

By the end of this module you should be able to:

- Find and use scientific and numerical functions in
`numpy`

. - Generate arrays of data and compute values on them.
- Understand the difference between lists and numpy arrays.
- Calculate N50 contig size.

This notebook should be completed while reading Chapter 2 of the Data Science Handbook that was assigned for your reading. The numpy library is a *third party* Python library, meaning that it is not distributed by default with every Python installation. It can be easily installed however, and provides a huge suite of tools for scientific computing. I think that the assigned Chapter introduces numpy very well, so this notebook will mostly consist of exercises to test your comprehension of the reading.

In [1]:

```
# start by importing numpy
import numpy as np
```

There are many ways to create a numpy array. Numpy has several built-in functions for generating arrays that are composed entirely of one value, or a range of values using `.zeros()`

, `.ones()`

, or `.arange()`

; or, we can also generate an array by passing in a list like in the last two examples below.

In [2]:

```
# create an array with ten items in it that are all zeros
np.zeros(10)
```

Out[2]:

In [3]:

```
# create an array with ten items that are all zeros as integers
np.zeros(10, dtype=int)
```

Out[3]:

In [4]:

```
# create an array with a range of values from 0-10
np.arange(0, 10)
```

Out[4]:

In [5]:

```
# create an array from a list
np.array([0, 3, 4, 10, 2, 2, 2, 2])
```

Out[5]:

In [6]:

```
# create an array that is made of 0 and 1 alternating
np.array([0, 1] * 5)
```

Out[6]:

When you create a numpy array from a list it tries to infer the datatype from contents of the list. Above, when we created a list of all int elements it created an int array. However, when we pass it a list below where some elements are ints and some are strings, it converts everything to strings. **This is because numpy works most efficiently by storing all data in an array as a single datatype**. You can create arrays with a mixed datatype but you lose much of the efficiency of numpy when you do so.

In [7]:

```
# mixed type lists will be converted to a single dtype array
np.array([0, 1, "apple", "orange"])
```

Out[7]:

Numpy arrays can be indexed just like Python list objects to select particular elements from them. In addition to the one dimension in which lists can be indexed, however, arrays can be indexed in multiple dimensions to select both rows and columns, and they can apply functions over these indices as well. If you need a refresher on how to `index`

and `slice`

objects in Python look back at our reading from last week (link), or google "python indexing".

In [8]:

```
# create a 2-dimensional array
np.zeros((4, 4))
```

Out[8]:

In [9]:

```
# create a 3-dimensional array
np.zeros((5, 3, 3))
```

Out[9]:

In [10]:

```
# index the third element in a 1-d array
np.arange(100)[3]
```

Out[10]:

In [11]:

```
# slice the third through tenth elements in a 1-d array
np.arange(100)[3:10]
```

Out[11]:

In [14]:

```
# one way
arr = np.zeros((3, 5), dtype=int)
arr[:] = 35
arr
```

Out[14]:

In [17]:

```
# another way
arr = np.empty((3, 5), dtype=int)
arr[:] = 35
arr
```

Out[17]:

In [19]:

```
# another way
arr = np.empty((3, 5), dtype=int)
arr.fill(35)
arr
```

Out[19]:

In [21]:

```
# another way
arr = np.full((3, 5), 35, dtype=int)
arr
```

Out[21]:

In [29]:

```
np.random.randint(0, 10, (3, 5))
```

Out[29]:

In [39]:

```
arr = np.arange(30)
```

In [40]:

```
arr[24:27]
```

Out[40]:

In [41]:

```
arr = np.arange(24).reshape((12, 2))
```

In [47]:

```
arr[7:, 0]
```

Out[47]:

When learning about lists and dictionaries we learned that a convenient way to use these objects is to first create an empty list or dictionary and then fill it with values as you iterate over an object and then append items to the list, or add new key/val pairs to the dictionary. In numpy you have to do things a bit differently. You start by creating an array that is the full size that you plan to work with, initialized with some null value like zeros, and then you update the values with new data.

The difference between these two approaches may seem subtle, but it can lead to huge speed improvements when done correctly in numpy. This is because numpy essentially reserves space in your computers memory for the size of the `array`

, whereas a `list`

that changes size as you extend it will need to keep updating the object in memory. As explained in your reading, this allows numpy to perform calculations using compiled functions in faster languages like C and fortran; it's a workaround that allows us to write pretty Python code but still benefit from super fast speed.

These details only begin to matter when you do pretty high level computing, but it's worth learning the motivation for why we use numpy, and why the code looks the way that it does.

In [48]:

```
# common with lists: start with an empty list and fill it as you iterate
empty = []
for i in range(100):
empty.append(i)
```

In [49]:

```
# filling an empty array doesn't work like with lists
empty = np.array([])
for i in range(100):
# ... a function does not exist to extend the size of arrays.
# ... only to generate new arrays of a given size
pass
```

In [50]:

```
# instead, you create the full sized array with null values and update them by indexing
empty = np.zeros(100)
for i in range(100):
empty[i] = i
```

In [51]:

```
# run this first so that we all start from the same random seed
import numpy as np
np.random.seed(12345)
```

In [52]:

```
# 1. comment: generate 50 random integers between 1K-50K
contig_sizes = np.random.randint(1000, 50000, size=50)
contig_sizes
```

Out[52]:

In [54]:

```
# 2. comment: sort integer array by size
contig_sizes.sort()
contig_sizes
```

Out[54]:

In [55]:

```
# 3. comment: reverse the order of the array
contig_sizes = contig_sizes[::-1]
contig_sizes
```

Out[55]:

In [57]:

```
# 4. comment: get the total sum (e.g., assembly size)
total_len = contig_sizes.sum()
total_len
```

Out[57]:

In [58]:

```
# 5. comment: get half of the total size
half_total_len = total_len / 2
half_total_len
```

Out[58]:

In [59]:

```
# 6. comment: make an array of zeros
contig_sum_lens = np.zeros(contig_sizes.size, dtype=int)
contig_sum_lens
```

Out[59]:

In [60]:

```
# 7. comment: interate over the items in the empty array and
# fill each index with the sum of all up to that point.
for i in range(contig_sizes.size):
contig_sum_lens[i] = contig_sizes[i:].sum()
contig_sum_lens
```

Out[60]:

In [61]:

```
# 8. comment: calculate which values are > half the total length
which_contigs_longer_than_half = contig_sum_lens > half_total_len
which_contigs_longer_than_half
```

Out[61]:

In [62]:

```
# 9. comment: keep the values that answered True above
contigs_longer_than_half = contig_sizes[which_contigs_longer_than_half]
contigs_longer_than_half
```

Out[62]:

In [63]:

```
# 10. comment: keep the lowest value from above, which is the N50
n50 = contigs_longer_than_half.min()
n50
```

Out[63]:

In [ ]:

```
```

In [66]:

```
def calculate_N50_from_contig_sizes(contig_sizes):
# 2. comment: sort integer array by size
contig_sizes.sort()
contig_sizes
# 3. comment: reverse the order of the array
contig_sizes = contig_sizes[::-1]
contig_sizes
# 4. comment: get the total sum (e.g., assembly size)
total_len = contig_sizes.sum()
total_len
# 5. comment: get half of the total size
half_total_len = total_len / 2
half_total_len
# 6. comment: make an array of zeros
contig_sum_lens = np.zeros(contig_sizes.size, dtype=int)
contig_sum_lens
# 7. comment: interate over the items in the empty array and
# fill each index with the sum of all up to that point.
for i in range(contig_sizes.size):
contig_sum_lens[i] = contig_sizes[i:].sum()
contig_sum_lens
# 8. comment: calculate which values are > half the total length
which_contigs_longer_than_half = contig_sum_lens > half_total_len
which_contigs_longer_than_half
# 9. comment: keep the values that answered True above
contigs_longer_than_half = contig_sizes[which_contigs_longer_than_half]
contigs_longer_than_half
# 10. comment: keep the lowest value from above, which is the N50
n50 = contigs_longer_than_half.min()
return n50
```

In [68]:

```
# run this first so that we all start from the same random seed
import numpy as np
np.random.seed(12345)
contig_sizes = np.random.randint(1000, 50000, size=50)
# run the new function
calculate_N50_from_contig_sizes(contig_sizes)
```

Out[68]: