<jupyter_start><jupyter_text>Week 2 This week I will create different variations of Phylogenetic trees to analyze Covid-19 genome variations and evolution, focusing on genome distances from the concensus sequence.  I. Get and Align Covid-19 Genome DataFollowing what my Instructor put in his guide, I downloaded sequence data from the Covid-19 Galaxy Project and used Virulign (found here) to align the downloaded information for easier stand comparisions.<jupyter_code>import wget

url = 'https://covid19.galaxyproject.org/genomics/4-Variation/current_complete_ncov_genomes.fasta'
file = './data/current_complete_ncov_genomes.fasta'
wget.download(url, file)<jupyter_output><empty_output><jupyter_text>Then, using a locally downloaded executable of Virulign, I ran the following command to align the sequence information. The downloaded sequence data and the aligned csv file are available here. The Spike protein acts as a entrance to a human cell and begins the infection process. **Spike proteins are very important for analysis because if you can prevent the proteins from attatching and fusing to human cells, you can prevent people from becoming ill.** ./virulign ./virulign/references/SARS-Cov2-2/S.xml ./data/current_complete_ncov_genomes.fasta --exportAlphabet Nucleotides --exportKind PositionTable > ./position_table.csv The last thing I am going to do is to put the aligned data into an easy to use table using pandas<jupyter_code>import pandas as pd

position_table = pd.read_csv('./data/position_table.csv')
position_table<jupyter_output><empty_output><jupyter_text>II. Build A Distance Matrix From The Consensus SequenceAccording to Wikipedia, a consensus sequence is "the calculated order of most frequent residues, either nucleotide of amino acid, found at each position in a sequence alignment" and represents a protein binding site. First, lets pull out the sequence.<jupyter_code>concensus_seq = position_table.drop('seqid',axis=1).mode(axis=0).T[0]
concensus_seq
position_table = position_table.set_index('seqid')<jupyter_output><empty_output><jupyter_text>Now, lets see which sequences are the farthest from the consensus sequence. Lets pick n=10.<jupyter_code>distance_from_concensus_seq = position_table.apply(lambda row: sum(row != concensus_seq),axis=1)
distance_from_concensus_seq_sorted = distance_from_concensus_seq.sort_values(ascending=False)
subset_seqs = distance_from_concensus_seq_sorted[:10].index
subset_seqs<jupyter_output><empty_output><jupyter_text>Now, it's time to create that distance matrix using our 10 sequences based off of their distance from the consensus sequence.<jupyter_code>distances = {}

for i,seqid1 in enumerate(subset_seqs):
    distances[seqid1,seqid1]=0
    for j in range(i+1,len(subset_seqs)):
        seqid2 = subset_seqs[j]
        distances[seqid1,seqid2] = sum(position_table.loc[seqid1] != position_table.loc[seqid2])
        distances[seqid2,seqid1] = distances[seqid1,seqid2]
distances = pd.Series(distances).unstack()
distances<jupyter_output><empty_output><jupyter_text>III. Build The Phylogenetic Tree and Analyze For creating the Phylogenetic Tree, I will utilize BioPython. BioPython is an awesome resource for bioinformatics algorithm analysis, specifcally for Python! You can find more information about BioPython here. To begin, I will construct I distance matrix specifically in the BioPython format so that it can be later used in the tree. Then, we use BioPython to construct the tree! Super easy. Distance Tree<jupyter_code>from Bio.Phylo.TreeConstruction import DistanceMatrix, DistanceTreeConstructor
from Bio import Phylo
import numpy as np

### Distance Matrix Construction ###
matrix = np.tril(distances.values).tolist()
for i in range(len(matrix)):
    matrix[i] = matrix[i][:i+1]
dm = DistanceMatrix(list(distances.index), matrix)


### Tree Construction ###
dt_constructor = DistanceTreeConstructor()
dt_tree = dt_constructor.nj(dm)
dt_tree.ladderize()

Phylo.draw(dt_tree)<jupyter_output><empty_output><jupyter_text>Parsimony Tree<jupyter_code>from Bio import AlignIO
from Bio import SeqIO
from Bio import Seq
from Bio import Phylo
from Bio.Phylo.TreeConstruction import *

sequences = list(SeqIO.parse('./data/current_complete_ncov_genomes.fasta', 'fasta'))
sequences = sequences[:10]
maxlen = max(len(sequence.seq) for sequence in sequences)

# pad to ensure all sequences have same length
for sequence in sequences:
    if len(sequence.seq) != maxlen:
        sequence.seq = Seq.Seq(str(sequence.seq).ljust(maxlen, '.'))
        
assert all(len(sequence.seq) == maxlen for sequence in sequences)
        
padded_f = "./data/padded_ncov_genomes.fasta"
with open(padded_f, 'w') as pf:
    SeqIO.write(sequences, pf, 'fasta')

aln = AlignIO.read(padded_f, 'fasta')
scorer = ParsimonyScorer()
searcher = NNITreeSearcher(scorer)
pars_constructor = ParsimonyTreeConstructor(searcher)
pars_tree = pars_constructor.build_tree(aln)

Phylo.draw(pars_tree)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08271281948078084, 'fraction_numerical': 0.007848661702555847, 'mean_word_length': 6.110157367668097, 'pattern_counts': {'":': 0, '<': 29, '<?xml version=': 0, '>': 30, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '271745', 'license': 'MIT', 'n_tokens_mistral': 1547, 'n_tokens_neox': 1480, 'n_words': 472, 'path': 'students/nabishop/week_2.ipynb', 'repo_name': 'anderson-github-classroom/csc-448-project', 'size': '4969'}