<jupyter_start><jupyter_text>Goals In this colab, you will:*   Learn how to clean up and preprocess genome data*   Convert genomic data into a feature matrix*   Build a logistic regression model predicting the country a SARS-CoV-2 lineage came from based on its genome<jupyter_code>#@title Set up the notebook
!pip install Biopython
from Bio import SeqIO
import numpy as np
import pandas as pd
import tqdm
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import confusion_matrix
from collections import Counter
from sklearn import model_selection, linear_model

import gdown
data_path = 'https://drive.google.com/uc?id=1f1CtRwSohB7uaAypn8iA4oqdXlD_xXL1'
cov2_sequences = 'SARS_CoV_2_sequences_global.fasta'
gdown.download(data_path, cov2_sequences, True)<jupyter_output>Collecting Biopython
[?25l  Downloading https://files.pythonhosted.org/packages/76/02/8b606c4aa92ff61b5eda71d23b499ab1de57d5e818be33f77b01a6f435a8/biopython-1.78-cp36-cp36m-manylinux1_x86_64.whl (2.3MB)
[K     |████████████████████████████████| 2.3MB 6.5MB/s 
[?25hRequirement already satisfied: numpy in /usr/local/lib/python3.6/dist-packages (from Biopython) (1.19.4)
Installing collected packages: Biopython
Successfully installed Biopython-1.78<jupyter_text>Data preprocessing  Read in and examine the data. We are going to read in a set of SARS-CoV-2 genomes from around the world. Note that sequence 0 is the "reference sequence"-- one of the original sequences from Wuhan. These global sequences come from the [NCBI database](https://www.ncbi.nlm.nih.gov/labs/virus/vssi//virus?SeqType_s=Nucleotide&VirusLineage_ss=SARS-CoV-2,%20taxid:2697049&SLen_i=29000%20TO%2031000&Completeness_s=complete&HostLineage_ss=Homo%20sapiens%20(human),%20taxid:9606).  You can examine the different sequences using the form below.<jupyter_code>sequences = [r for r in SeqIO.parse(cov2_sequences, 'fasta')]
sequence_num =  23#@param {type:"integer"}
print(sequences[sequence_num])<jupyter_output>ID: MN975262
Name: MN975262
Description: MN975262 |Severe acute respiratory syndrome coronavirus 2 isolate 2019-nCoV_HKU-SZ-005b_2020| complete genome|China
Number of features: 0
Seq('ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGT...---')<jupyter_text>**Exercise: How many sequences are there?**Note: Sequences have been uploaded/stored in a variable called ```sequences```.<jupyter_code>len(sequences)<jupyter_output><empty_output><jupyter_text>**Exercise: How different are the 1st (non-reference) and 10th SARS-CoV-2 sequences?**<jupyter_code>sequence_first = np.array(sequences[0]) 
sequence_tenth = np.array(sequences[9])

# Compute length, number of different bases, 
# and number of the same bases.

#### Fill in code to compute the length of sars-cov-2 (ignore -'s if possible!) ########
length_first = sum(sequence_first!='-') # ignore -'s
n_bases_different = sum(sequence_first!=sequence_tenth)
n_bases_same = sum(sequence_first==sequence_tenth)

print("1. Length of SARS-CoV-2 genome: ", length_first)
print("2. # of bases that differ from RTG13: ", n_bases_different)
print("3. Percent similarity: %", 100*n_bases_same/length_first)<jupyter_output>1. Length of SARS-CoV-2 genome:  29903
2. # of bases that differ from RTG13:  7
3. Percent similarity: % 99.9765909774939<jupyter_text>**Exercise (BONUS):  Make a histogram of the number of mutations each SARS-CoV-2 sequence has compared to the reference genome.**<jupyter_code>reference = np.array(sequences[0])
## Note: This can take a minute or two to run. Be patient.

#### FILL IN #####
mutations_per_seq = np.array([sum(np.array(s)!=reference) for s in sequences])
plt.hist(mutations_per_seq)
plt.xlabel('# mutations')
plt.ylabel('# sequences')
plt.show()<jupyter_output><empty_output><jupyter_text>Interestingly, it looks like there are a couple sequences with a LOT of mutations! We can investigate these sequences a little more.**Examine some of these sequences with high number of mutations by selecting the minimum  of mutations from the form below. What do you notice about the sequences? Discuss with your instructor and peers.**<jupyter_code>min_number_of_mutations  =  320#@param {type:"integer"}
idx = np.random.choice(np.where(mutations_per_seq>min_number_of_mutations)[0])
print("Sequence %i has > %.0f mutations! \n" % (idx, min_number_of_mutations))
print(sequences[idx], '\n')
print("The sequence is composed of: ")
Counter(np.array(sequences[idx]))<jupyter_output>Sequence 108 has > 320 mutations! 

ID: MT374101
Name: MT374101
Description: MT374101 |Severe acute respiratory syndrome coronavirus 2 isolate SARS-CoV-2/human/TWN/CGMH-CGU-02/2020| complete genome|Taiwan
Number of features: 0
Seq('ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGT...AAA') 

The sequence is composed of:<jupyter_text>Missing Data It is hard to see, but some of the sequences have `N` in them. Run the cell below for an example **Exercise: Calculate the number of sequences that have an ```N``` in them.****What do you think ```N``` means?**<jupyter_code>n_sequences_with_N = 0
for i in sequences:
  for j in i:
    if j=='N':
      n_sequences_with_N+=1
      break
print('%i sequences have at least 1 "N"!' % n_sequences_with_N)<jupyter_output>326 sequences have at least 1 "N"!<jupyter_text>`N` is not a nucleic acid- it just stands for "missing", or "low quality". "Missing" is different than ```_``` or a deletion. At the locations with ```N```, the sequencing machine had low quality data here, so it was unable to determine what base was at that location. We should remember this when we extact our features. Stay tuned for more on sequencing machines and how sequences are built in the bonus colab of this project!  Feature ExtractionWe are going to build a model that predicts the country a SARS-CoV-2 virus came from based on its genome.  **Exercise: Recall the structure of machine learning models.** **In general what two categories of data do we need to build a supervised machine learning model? What will we use for each category?**<jupyter_code>_1_  =  'data with labeled features (genomes)' #@param {type:"string"}
_2_  =  'data with labels of the country where the genome came from' #@param {type:"string"}

print('1. We need a set of FEATURES (X).\n',
      '  Our features will be the genomes of the different sequences.')
print('2. We need LABElS (Y).\n',
      '  Our labels will be the country that each sequence came from.')<jupyter_output>1. We need a set of FEATURES (X).
   Our features will be the genomes of the different sequences.
2. We need LABElS (Y).
   Our labels will be the country that each sequence came from.<jupyter_text>**Question: How will we turn our features into a numeric matrix?**  Extract Features (X) Remember that our input must be a *numeric* matrix/table.We are going to create a matrix where our features are the presence/absence of a specific mutation (given by ``````_``````).  Our columns will be 1_A, 1_T, 3_G, 4_A, etc.| Sequence ID | 1_A | 1_T | 3_G | 4_A  | ...||-------------|-----|-----|-----|------|----||Sequence 1   |  1  |  0  |   1 |    0 |  0 ||Sequence 2   |  0  |  0  |   1 |    0 |  0 ||Sequence 3   |  1  |  1  |   0 |    0 |  0 ||Sequence 4   |  0  |  0  |   0 |    1 |  1 ||Sequence 5   |  1  |  0  |   0 |    0 |  1 |<jupyter_code># Note: This can take a couple minutes to run! 
# but we can monitor our progress using the tqdm library
mutation_df = pd.DataFrame()
n_bases_in_seq = len(sequences[0])

# Iterate though all positions in this sequence.
for location in tqdm.tqdm(range(n_bases_in_seq)): # tqdm is a nice library that prints our progress.
  bases_at_location = np.array([s[location] for s in sequences])
  # If there are no mutations at this position, move on.
  if len(set(bases_at_location))==1: continue # If
  for base in ['A', 'T', 'G', 'C', '-']:
    feature_values = (bases_at_location==base)
    
    
    # Set the values of any base that equals 'N' to np.nan.
    feature_values[bases_at_location=='N'
                   ] = np.nan
    
    # Convert from T/F to 0/1.
    feature_values  = feature_values*1
    
    # Make the column name look like <location>_<base> (1_A, 2_G, 3_A, etc.)
    column_name = str(location)+"-"+base
    mutation_df[column_name] = feature_values

# Print the size of the feature matrix/table.
n_rows = np.shape(mutation_df)[0]
n_columns = np.shape(mutation_df)[1]
print("Size of matrix: %i rows x %i columns" %(n_rows, n_columns))

# Check what the matrix looks like:
mutation_df.head()<jupyter_output>100%|██████████| 29903/29903 [01:21<00:00, 367.16it/s]<jupyter_text>Extract Label (Y) We are going to use the region of the world that each sample came from as the **label**. ![alt text](https://upload.wikimedia.org/wikipedia/commons/3/3d/Flag-map_of_the_world_%282017%29.png)First, let's see how many samples we have from different countries/regions.<jupyter_code>#@title ###**Exercise: Explore the different number of samples that come from each country/region.**
country = "Pakistan" #@param dict_keys(['China', 'Kazakhstan', 'India', 'Sri Lanka', 'Taiwan', 'Hong Kong', 'Viet Nam', 'Thailand', 'Nepal', 'Israel', 'South Korea', 'Iran', 'Pakistan', 'Turkey', 'Australia', 'USA']
countries = [(s.description).split('|')[-1] for s in sequences]
print("There are %i sequences from %s." %
      (Counter(countries)[country], country))<jupyter_output>There are 2 sequences from Pakistan.<jupyter_text>Since some countries/regions only have a couple samples, we are going to use the geographical **region** of the world as our labels. Since we have a large number of samples from Asia, North America, and Oceania, we will filter ours sequences to just these geographical regions. We will convert our countries/regions to geographical regions using the code below. **Exercise: Convert each country/region to it's geographical region of the world.** **Use the code below to create a dictionary of ``````:`````` where ```region``` is either ```'Oceania'```, ```'North America'```, or ```'Asia'```, and convert each country/region to geographical region.**<jupyter_code>countries_to_regions_dict = {
         'Australia': 'Oceania',
         'China': 'Asia',
         'Hong Kong': 'Asia',
         'India': 'Asia',
         'Nepal': 'Asia',
         'South Korea': 'Asia',
         'Sri Lanka': 'Asia',
         'Taiwan': 'Asia',
         'Thailand': 'Asia',
         'USA': 'North America',
         'Viet Nam': 'Asia' 
}

regions = [countries_to_regions_dict[c] if c in 
           countries_to_regions_dict else 'NA' for c in countries]
mutation_df['label'] = regions<jupyter_output><empty_output><jupyter_text>**Now see how many samples there are from each geographical region of the world.**<jupyter_code>region = "Oceania" #@param ['Oceania', 'North America', 'Asia']
print("There are %i sequences from %s." %
      (Counter(regions)[region], region))<jupyter_output>There are 157 sequences from Oceania.<jupyter_text>Balancing the Data Recall that ML models work the best if we have *balanced* data- a dataset with equal numbers of samples with each label. Run the following code to remove duplicate samples from the dataset, and then balance the samples. **Exercise: Balance the data equally between samples from Asia, Oceania, and North America**<jupyter_code>balanced_df = mutation_df.copy()
balanced_df['label'] = regions
balanced_df = balanced_df[balanced_df.label!='NA']
balanced_df = balanced_df.drop_duplicates()
samples_north_america = balanced_df[balanced_df.label== 'North America']
samples_oceania = balanced_df[balanced_df.label== 'Oceania']
samples_asia = balanced_df[balanced_df.label== 'Asia']

# Number of samples we will use from each region.
n = min(len(samples_north_america),
        len(samples_oceania),
        len(samples_asia))

balanced_df = pd.concat([samples_north_america[:n],
                    samples_asia[:n],
                    samples_oceania[:n]])
print("Number of samples in each region: ", Counter(balanced_df['label']))<jupyter_output>Number of samples in each region:  Counter({'North America': 128, 'Asia': 128, 'Oceania': 128})<jupyter_text>Logistic Regression Model ***Congrats!***  We finally are done with preprocessing/cleaning our data! Although tedious, this is an important part of doing machine learning in biology. The data can be complex and messy, and if we don't do some cleaning up beforehand, our models will have poor performance.![](https://media.makeameme.org/created/we-did-it-3b3ac27d2a.jpg)Finally, run the code to set up a ```X``` feature matrix and a ```Y``` label list from our ```balanced_df```. You can explore the different values using the code below:<jupyter_code>X = balanced_df.drop('label', 1)
Y = balanced_df.label
data = "X (features)" #@param ['X (features)', 'Y (label)']
start = 1 #@param {type:'integer'}
stop =  500#@param {type:'integer'}

if start>=stop:print("Start must be < stop!")
else:
  if data=='X (features)':
    print(X.iloc[start:stop])
  if data=='Y (label)':
    print(Y[start:stop])<jupyter_output>0-A  0-T  0-G  0-C  0--  ...  29902-A  29902-T  29902-G  29902-C  29902--
323    1    0    0    0    0  ...        0        0        0        0        1
324    0    0    1    0    0  ...        1        0        0        0        0
325    1    0    0    0    0  ...        0        0        0        0        1
326    1    0    0    0    0  ...        0        0        0        0        1
327    0    0    1    0    0  ...        0        0        0        0        1
..   ...  ...  ...  ...  ...  ...      ...      ...      ...      ...      ...
290    0    0    0    0    1  ...        0        0        0        0        1
291    0    0    0    0    1  ...        0        0        0        0        1
292    0    0    0    0    1  ...        0        0        0        0        1
293    0    0    0    0    1  ...        0        0        0        0        1
294    0    0    0    0    1  ...        0        0        0        0        1

[383 rows x 12680 columns]<jupyter_text>Training We will be using the logistic regression model we have learned about with one modification. We will use the "multinomial" class of logistic regression model.  This is used when there are more than 2 categories in the label set. In our case, we have ```Asia```, ```Oceania```, and ```North America``` as our possible labels. **Exercise: Train the model using the standard pipeline you have mastered!**<jupyter_code>lm = linear_model.LogisticRegression(
    multi_class="multinomial", max_iter=1000,
    fit_intercept=False, tol=0.001, solver='saga', random_state=42)

# Split into training/testing set. Use a training size of .8
X_train, X_test, Y_train, Y_test = model_selection.train_test_split(X, Y, train_size=0.8)

# Train/fit model.
#### FILL IN ########
lm.fit(X_train, Y_train)<jupyter_output><empty_output><jupyter_text>Testing/EvaluationIn addition to printing the accuracy of a model, we can also use a *confusion matrix* to see how well the model performed. **Exercise: Evaluate the model on the test set.**<jupyter_code># Predict on the test set.
Y_pred = lm.predict(X_test)

# Compute accuracy.
accuracy = 100*np.mean(Y_pred == Y_test)
print("Accuracy: %", accuracy)

# Compute confusion matrix.
confusion_mat = pd.DataFrame(confusion_matrix(Y_test, Y_pred))
confusion_mat.columns = [c + ' predicted' for c in lm.classes_]
confusion_mat.index = [c + ' true' for c in lm.classes_]

print(confusion_mat)<jupyter_output>Accuracy: % 94.8051948051948
                    Asia predicted  North America predicted  Oceania predicted
Asia true                       23                        0                  2
North America true               1                       23                  0
Oceania true                     1                        0                 27<jupyter_text>Wrapping up ![alt text](https://memecrunch.com/meme/60HHF/done-and-done/image.jpg) ***Great job!*** You built a pretty accurate model that uses genomic data to predict what country a SARS-CoV-2 sample comes from.<jupyter_code>#@title #### **Exercise: To wrap up, write a couple of sentences explaining how your model performed.**
Response = "the model performed well with an accuracy of 94%" #@param {type:"string"}<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11526844603340936, 'fraction_numerical': 0.03525858349257227, 'mean_word_length': 3.5765867418899857, 'pattern_counts': {'":': 0, '<': 61, '<?xml version=': 0, '>': 63, 'https://': 6, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '159248', 'license': 'MIT', 'n_tokens_mistral': 5638, 'n_tokens_neox': 5093, 'n_words': 1773, 'path': 'From_Genomes_to_Machine_Learning.ipynb', 'repo_name': 'HowardHuang1/SARS-CoV-2', 'size': '16223'}