<jupyter_start><jupyter_text>🔬 Sequence Comparison of DNA using `BioPython` 🦠 `Covid-19`, `SARS`, `MERS`, and `Ebola`   Analysis Techniques:* Compare their DNA sequence and Protein (Amino Acid) sequence* GC Content* Freq of Each Amino Acids* Find similarity between them  * Alignment  * hamming distance* 3D structure of each| DNA Sequence     | Datasource                                                    ||:-----------------|:--------------------------------------------------------------|| Latest Sequence  | https://www.ncbi.nlm.nih.gov/genbank/sars-cov-2-seqs/         || Wuhan-Hu-1       | https://www.ncbi.nlm.nih.gov/nuccore/MN908947.3?report=fasta  || Covid19          | https://www.ncbi.nlm.nih.gov/nuccore/NC_045512.2?report=fasta || SARS             | https://www.ncbi.nlm.nih.gov/nuccore/NC_004718.3?report=fasta || MERS             | https://www.ncbi.nlm.nih.gov/nuccore/NC_019843.3?report=fasta || EBOLA            | https://www.ncbi.nlm.nih.gov/nuccore/NC_002549.1?report=fasta | 1. Analysis Techniques<jupyter_code>import warnings
warnings.filterwarnings('ignore')
import pandas as pd 
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
from Bio.Seq import Seq

# Create our sequence 
seq1 = Seq('ACTCGA')
seq2 = Seq('AC')<jupyter_output><empty_output><jupyter_text>GC Contents In DNA  * `GC-content` (or guanine-cytosine content) is the **percentage of nitrogenous bases** in a DNA or RNA molecule that are either guanine (`G`) or cytosine (`C`) Usefulness  * In polymerase chain reaction (PCR) experiments, the GC-content of short oligonucleotides known as primers is often used to predict their **annealing temperature** to the template DNA.    * A `high` GC-content level indicates a relatively higher melting temperature.    * DNA with `low` GC-content is less stable than DNA with high GC-content  > Question: which sequence is more stable when heat is applied?<jupyter_code>from Bio.SeqUtils import GC
# Check GC (guanine-cytosine) percentage in sequence

print(f"{GC(seq1)}% \t({seq1})")
print(f"{GC(seq2)}% \t({seq2})")<jupyter_output>50.0% 	(ACTCGA)
50.0% 	(AC)<jupyter_text>Sequence Alignment* `Global alignment` finds the best concordance/agreement between all characters in two sequences* `Local Alignment` finds just the subsequences that align the best<jupyter_code>from Bio import pairwise2
from Bio.pairwise2 import format_alignment

print('seq1 =', seq1, '\nseq2 =', seq2, '\n\n')

# Global alignment
alignments = pairwise2.align.globalxx(seq1, seq2)

print(f'Alignments found: {len(alignments)}')
print(*alignments)
# Print nicely
print(format_alignment(*alignments[0]))
# 2nd alignment
print(format_alignment(*alignments[1]))
# To see all possible alignments
for a in alignments:
    print(format_alignment(*a), '\n')
# Get the number of possible sequence alignments
alignment_score = pairwise2.align.globalxx(seq1,seq2,one_alignment_only=True,score_only=True)
alignment_score<jupyter_output><empty_output><jupyter_text>Sequence Similarity* Fraction of nucleotides that is the same/ total number of nucleotides * 100%<jupyter_code>alignment_score/len(seq1)*100<jupyter_output><empty_output><jupyter_text>Hamming Distance: `How Many Subsitutions are Required to Match Two Sequences?`* Hamming distance between two strings of equal length is the number of positions at which the corresponding symbols are different.* In other words, it measures the minimum number of substitutions required to change one string into the other, or the minimum number of errors that could have transformed one string into the other* It is used for error detection or error correction* It is used to quantify the similarity of DNA sequences   Edit Distance   * Is a way of quantifying how dissimilar two strings (e.g., words) are to one another by counting the minimum number of operations required to transform one string into the other (e.g. Levenshtein distance)<jupyter_code>def hamming_distance(lhs, rhs):
    return len([(x,y) for x,y in zip(lhs,rhs) if x != y])
hamming_distance('TT', 'ACCTA')
def hammer_time(s1, s2, verbose=True):
    """Take two nucleotide sequences s1 and s2, and display 
       the possible alignments and hamming distance.
    """
    if verbose:
        print('s1 =', s1, '\ns2 =', s2, '\n\n')
        print('Hamming Distance:', hamming_distance(s1, s2), '\n(min substitutions for sequences to match)')
        print('\nAlignment Options:\n\n')
    alignments = pairwise2.align.globalxx(s1, s2)
    for a in alignments:
        print(format_alignment(*a), '\n')
    

s1 = 'ACTCGAA'
s2 = 'ACGA'
hammer_time(s1, s2)<jupyter_output>s1 = ACTCGAA 
s2 = ACGA 


Hamming Distance: 2 
(min substitutions for sequences to match)

Alignment Options:


ACTCGAA
|  || |
A--CG-A
  Score=4
 

ACTCGAA
||  | |
AC--G-A
  Score=4
 

ACTCGAA
|  ||| 
A--CGA-
  Score=4
 

ACTCGAA
||  || 
AC--GA-
  Score=4<jupyter_text>Dot Plot* A dot plot is a graphical method that allows the **comparison of two biological sequences** and identify regions of **close similarity** between them.* Simplest explanation: put a dot wherever sequences are identical UsefulnessDot plots can also be used to visually inspect sequences for    - Direct or inverted repeats    - Regions with low sequence complexity    - Similar regions    - Repeated sequences    - Sequence rearrangements    - RNA structures    - Gene orderAcknowledgement: https://stackoverflow.com/questions/40822400/how-to-create-a-dotplot-of-two-dna-sequence-in-python<jupyter_code>def delta(x,y):
    return 0 if x == y else 1

def M(seq1,seq2,i,j,k):
    return sum(delta(x,y) for x,y in zip(seq1[i:i+k],seq2[j:j+k]))

def makeMatrix(seq1,seq2,k):
    n = len(seq1)
    m = len(seq2)
    return [[M(seq1,seq2,i,j,k) for j in range(m-k+1)] for i in range(n-k+1)]

def plotMatrix(M,t, seq1, seq2, nonblank = chr(0x25A0), blank = ' '):
    print(' |' + seq2)
    print('-'*(2 + len(seq2)))
    for label,row in zip(seq1,M):
        line = ''.join(nonblank if s < t else blank for s in row)
        print(label + '|' + line)

def dotplot(seq1,seq2,k = 1,t = 1):
    M = makeMatrix(seq1,seq2,k)
    plotMatrix(M, t, seq1,seq2) #experiment with character choice
# The dot plot: put a dot where the two sequences are identical

s1 = 'ACTCGA'
s2 = 'AC'

dotplot(s1, s2)
#  Identical proteins will show a diagonal line.
s1 = 'ACCTAG'
s2 = 'ACCTAG'

dotplot(s1, s2)

print('\n\n')
hammer_time(s1, s2, verbose=False)<jupyter_output>|ACCTAG
--------
A|■   ■ 
C| ■■   
C| ■■   
T|   ■  
A|■   ■ 
G|     ■



ACCTAG
||||||
ACCTAG
  Score=6<jupyter_text>🔬 2. Comparative Analysis of Virus DNA 🦠 `Covid-19`, `SARS`, `MERS`, `Ebola`  * Covid19(`SARS-CoV2`) is a novel coronavirus identified as the cause of coronavirus disease 2019 (COVID-19) that began in Wuhan, China in late 2019 and spread worldwide.* MERS(`MERS-CoV`) was identified in 2012 as the cause of Middle East respiratory syndrome (MERS).* SARS(`SARS-CoV`) was identified in 2002 as the cause of an outbreak of severe acute respiratory syndrome (SARS). `fasta` DNA Sequence Files* Covid19 : https://www.rcsb.org/3d-view/6LU7* SARS: https://www.ncbi.nlm.nih.gov/nuccore/NC_004718.3?report=fasta* MERS: https://www.ncbi.nlm.nih.gov/nuccore/NC_019843.3?report=fasta* EBOLA:https://www.rcsb.org/structure/6HS4<jupyter_code>import pandas as pd
import numpy as np
from Bio import SeqIO

covid = SeqIO.read("../data/01_COVID_MN908947.3.fasta","fasta")
mers = SeqIO.read("../data/02_MERS_NC_019843.3.fasta","fasta")
sars = SeqIO.read("../data/03_SARS_rcsb_pdb_5XES.fasta","fasta")
ebola = SeqIO.read("../data/04_EBOLA_rcsb_pdb_6HS4.fasta","fasta")

# Convert imports to BioPython sequences
covid_seq = covid.seq
mers_seq = mers.seq
sars_seq = sars.seq
ebola_seq = ebola.seq

# Create dataframe
df = pd.DataFrame({'name': ['COVID19', 'MERS', 'SARS', 'EBOLA'],
                   'seq': [covid_seq, mers_seq, sars_seq, ebola_seq]})
df<jupyter_output><empty_output><jupyter_text>Length of Each Genome<jupyter_code>df['len'] = df.seq.apply(lambda x: len(x))

df[['name', 'len']].sort_values('len', ascending=False) \
                   .style.bar(color='#cde8F6', vmin=0, width=100, align='left')<jupyter_output><empty_output><jupyter_text>* `MERS`, `COVID` and `SARS` all have about the same genome length (30,000 base pairs)  Which of them is more heat stable?<jupyter_code># Check the GC content
df['gc_content'] = df.seq.apply(lambda x: GC(x))

df[['name', 'gc_content']].sort_values('gc_content', ascending=False) \
                          .style.bar(color='#cde8F6', vmin=0)<jupyter_output><empty_output><jupyter_text>* `MERS` is the most stable with a GC of `41.24` followed by Ebola  Translate RNA into proteins How many proteins are in each dna sequence?<jupyter_code># Translate the RNA into Proteins
df['proteins'] = df.seq.apply(lambda s: len(s.translate()))

df[['name', 'proteins']].sort_values('proteins', ascending=False) \
                        .style.bar(color='#cde8F6', vmin=0)<jupyter_output><empty_output><jupyter_text>How Many Amino Acids are Created?<jupyter_code>from Bio.SeqUtils.ProtParam import ProteinAnalysis
from collections import Counter
# Method 1
covid_analysed = ProteinAnalysis(str(covid_protein))
mers_analysed = ProteinAnalysis(str(mers_protein))
sars_analysed = ProteinAnalysis(str(sars_protein))
ebola_analysed = ProteinAnalysis(str(ebola_protein))
# Check for the Frequence of AA
covid_analysed.count_amino_acids()
# Method 2
from collections import Counter

# Find the Amino Acid Frequency
df['aa_freq'] = df.seq.apply(lambda s: Counter(s.translate()))
df<jupyter_output><empty_output><jupyter_text>Most Common Amino Acid<jupyter_code># For Covid
df[df.name=='COVID19'].aa_freq.values[0].most_common(10)
# Plot the Amino Acids of COVID-19

aa = df[df.name=='COVID19'].aa_freq.values[0]
plt.bar(aa.keys(), aa.values())
# All viruses -- same chart (not stacked)

for virus in df.name:
    aa = df[df.name==virus].aa_freq.values[0]
    plt.bar(aa.keys(), aa.values())
plt.show()<jupyter_output><empty_output><jupyter_text>Dot Plots of Opening Sequences<jupyter_code># COVID and MERS
dotplot(covid_seq[0:10],mers_seq[0:10])
# COVID and SARS
n = 10
dotplot(covid_seq[0:n],sars_seq[0:n])
# Plotting function to illustrate deeper matches

def dotplotx(seq1, seq2, n):
    seq1=seq1[0:n]
    seq2=seq2[0:n]
    plt.imshow(np.array(makeMatrix(seq1,seq2,1)))
    # on x-axis list all sequences of seq 2
    xt=plt.xticks(np.arange(len(list(seq2))),list(seq2))
    # on y-axis list all sequences of seq 1
    yt=plt.yticks(np.arange(len(list(seq1))),list(seq1))
    plt.show()
    
dotplotx(covid_seq, sars_seq, n=100)<jupyter_output><empty_output><jupyter_text>Notice the large diagonal line for the second half of the first 100 nucleotides - indicating these are the same for `COVID19` and `SARS`<jupyter_code>dotplotx(covid_seq, ebola_seq, n=100)<jupyter_output><empty_output><jupyter_text>No corresponding matches for `EBOLA` and `COVID`  Calculate Pairwise Alignment for the First 100 Nucleotides<jupyter_code>def pairwise_alignment(s1, s2, n):
    if n == 'full': n = min(len(s1), len(s2))
    alignment = pairwise2.align.globalxx(s1[0:n], s2[0:n], one_alignment_only=True, score_only=True)
    print(f'Pairwise alignment: {alignment:.0f}/{n} ({(alignment/n)*100:0.1f}%)')

    
# SARS and COVID 
pairwise_alignment(covid_seq, sars_seq, n=100)
pairwise_alignment(covid_seq, sars_seq, n=10000)
pairwise_alignment(covid_seq, sars_seq, n=len(sars_seq))<jupyter_output>Pairwise alignment: 24654/29751 (82.9%)<jupyter_text>* `82.9`% of the COVID19 genome is exactly the same as SARS<jupyter_code>pairwise_alignment(covid_seq, mers_seq, n='full')
pairwise_alignment(covid_seq, ebola_seq, n='full')<jupyter_output>Pairwise alignment: 12377/18959 (65.3%)<jupyter_text>* `COVID19` and `SARS` have a `82.9`% similarity. Both are of the same genus and belong to `Sars_Cov`.* `COVID19` and `EBOLA` have a `65.3`% similarity since they are from a different family of virus  Example of the Opening Sequence of `COVID19` and `SARS`Sequencing found similar structure from `40:100` so lets use our functions to visualise it.<jupyter_code>s1 = covid_seq[40:100]
s2 = sars_seq[40:100]

print('Similarity matrix (look for diagonal)')
dotplotx(s1, s2, n=100)

print('Possible alignment pathways: \n\n')
hammer_time(s1, s2, verbose=False)<jupyter_output>Similarity matrix (look for diagonal)
{'fraction_non_alphanumeric': 0.13253012048192772, 'fraction_numerical': 0.03145467777148864, 'mean_word_length': 4.578709968425801, 'pattern_counts': {'":': 0, '<': 64, '<?xml version=': 0, '>': 64, 'https://': 11, 'lorem ipsum': 0, 'www.': 10, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '11513', 'license': 'MIT', 'n_tokens_mistral': 4721, 'n_tokens_neox': 4290, 'n_words': 1176, 'path': 'notebooks/02-SequenceComparison.ipynb', 'repo_name': 'lukexyz/DNA-Toolbox', 'size': '12367'}