<jupyter_start><jupyter_text>BMI 6018 Final Project <NAME> In this notebook, I am going to walk you through some of the functions and methods I wrote for DNA sequence manipulation. My eventual goal was to be able to take two sequences and perform sequence alignment on them, using a version of the Smith-Waterman algorithm for local alignment, or the Needleman-Wunsch algorithm for global alignment. However, implementing these algorithms ended up being far beyond the scope of my rudimentary Python programming skills, and they are well implemented in the BioPython package, so I ended up not trying to re-invent the wheel.  DNA manipulation<jupyter_code>import DNA
import mutation
dir(DNA)<jupyter_output><empty_output><jupyter_text>Having a function that can parse fasta files is very important if you are going to be analyzing sequence data. The parser will take a fasta file and parse it into a dictionary, with the id code as the key, and the sequence as the value.<jupyter_code>with open("ex_ref.fasta", "r") as fo:
    parsed_dna = DNA.parse_fasta(fo)
parsed_dna<jupyter_output><empty_output><jupyter_text>Let's take the value of one of our dictionary entries and set it as a string.<jupyter_code>dna1 = parsed_dna["1111750"]
dna1<jupyter_output><empty_output><jupyter_text>DNA has been built into a class, inheriting from the string class, with a wide variety of functions available to manipulate the DNA string, in addition to the normal string methods.<jupyter_code>dna1 = DNA.DNA(dna1)
help(DNA.DNA)<jupyter_output>Help on class DNA in module DNA:

class DNA(builtins.str)
 |  Class for representing DNA sequences as strings
 |  
 |  Method resolution order:
 |      DNA
 |      builtins.str
 |      builtins.object
 |  
 |  Methods defined here:
 |  
 |  __init__(self, seq)
 |      Create DNA object to string seq.
 |  
 |  __str__(self)
 |      Return str(self).
 |  
 |  codons(self, rf=0)
 |      Finds the list of codons in a sequence in a given reading frame.
 |      
 |      In a sequence there are 6 possible reading frames, 3 on forward strand, 3 on complement.
 |      Arguments: sequence for splitting in to codons
 |      rf: reading frame. 0 = 1st forward frame; 1 = 2nd forward frame; 2 = 3rd forward frame
 |          3 = 1st complement frame; 4 = 2nd complement frame; 5 = 3rd complement frame
 |      Returns: list of codons in sequence for specified reading frame
 |  
 |  complement(self)
 |      Finds the complement of the sequence by list comprehension using dict of complement values
 |  
 [...]<jupyter_text>The first method is a simple reversal. The two paired sequences of a DNA strand run anti-parallel, where the complementary sequences run in the opposite direction of each other. DNA will only get transcribed in the 5' to 3' direction, so it is important to be able to get the reversed complement strand.<jupyter_code>dna1.reverse()
comp_dna = DNA.DNA(dna1.complement())
comp_dna.reverse()<jupyter_output><empty_output><jupyter_text>When DNA is transcribed into RNA, the RNA forms 3 letter "codes" for amino acids that will make up a protein during translation. These codes, or codons, are important to know as the various mutations that change a sequence can make a drastic change to the shape or folding of a protein because of the chemical differences between the amino acids. Most of the time, the codons are view in RNA code, but due to most sequencing being done on DNA, a DNA codon translation table is still useful.<jupyter_code>dna1.codons()<jupyter_output><empty_output><jupyter_text>RNA is similar to DNA, but has Uracil in place of Thymine in the sequence.<jupyter_code>dna1.transcribe()<jupyter_output><empty_output><jupyter_text>GC content is important to know in molecular biology. Due to the three hydrogen bonds that form between a guanine and cytosine vs. the two hydrogen bonds that form between an adenine and thymine or adenine and uracil, sequences with higher GC content will have a higher melting point in PCR assays.<jupyter_code>dna1.gc_content()<jupyter_output><empty_output><jupyter_text>The DNA class can also translate the sequence into the amino acid codes. The translate function runs until it hits a stop codon at which point the translation stops. A sequence of DNA actually will have 6 potential reading frames for translating to protein sequences. I could not get my class defined six frames function to work, but did get one standalone in the module. The six frames function will translate our sequence to RNA, and then align the amino acid code in the middle of their codon, 3 on the forward strand, and 3 on the complement. The visualization is easier to see if the sequence is shorted a bit, so I have limited the sequence to 100 base pairs.<jupyter_code>dna1.translate()
DNA.six_frames(dna1.transcribe()[:100])<jupyter_output>R  A  A  D  G  *  V  T  R  R  N  V  P  Q  S  E  G  *  A  Q  E  T  E  S  N  T  A  Y  G  L  R  I 
  T  S  G  G  R  V  S  N  A  *  E  R  A  P  K  *  G  I  S  S  G  N  *  V  *  Y  R  I  W  S  S  D  * 
 D  E  R  R  T  G  E  *  R  V  G  T  C  P  K  V  R  D  K  L  R  K  L  S  L  I  P  H  M  V  F  G  L 
GACGAGCGGCGGACGGGUGAGUAACGCGUAGGAACGUGCCCCAAAGUGAGGGAUAAGCUCAGGAAACUGAGUCUAAUACCGCAUAUGGUCUUCGGAUUAA
CUGCUCGCCGCCUGCCCACUCAUUGCGCAUCCUUGCACGGGGUUUCACUCCCUAUUCGAGUCCUUUGACUCAGAUUAUGGCGUAUACCAGAAGCCUAAUU
 L  L  A  A  C  P  L  I  A  H  P  C  T  G  F  H  S  L  F  E  S  F  D  S  D  Y  G  V  Y  Q  K  P  N 
  C  S  P  P  A  H  S  L  R  I  L  A  R  G  F  T  P  Y  S  S  P  L  T  Q  I  M  A  Y  T  R  S  L  I 
   A  R  R  L  P  T  H  C  A  S  L  H  G  V  S  L  P  I  R  V  L  *  L  R  L  W  R  I  P  E  A  *<jupyter_text>The last program in the DNA module is the open reading frame protein sequence finder. An open reading frame is defined as a sequence that starts at a start codon and then will get translated until it hits a stop codon in that sequence. If the sequence ends before a stop codon is reached, the protein is not added to the list of possible proteins. The orf_protein function parses a fasta file and then returns the set of proteins that can be found in any open reading frame, including the complement sequence that is infered by the DNA sequence.<jupyter_code>with open("ex_ref.fasta", "r") as fo:
    firstNlines=fo.readlines()[0:5]
    possible_proteins = DNA.orf_protein(firstNlines)
possible_proteins<jupyter_output><empty_output><jupyter_text>RNA and mutations<jupyter_code>help(mutation)<jupyter_output>Help on module mutation:

NAME
    mutation - This module is for various mutations of sequence data.

DESCRIPTION
    Functions include random RNA generator, mutation programs that modify the RNA strings
    and mutations finder programs that attempt to identify the mutation.
    
    Copyright <NAME>, 2017

FUNCTIONS
    data_gen(number)
        Given an input number, generates an RNA string, mutates the string, and
        then pickles the data for later use.
    
    deletion(rna)
        Function to take an RNA string and randomly choose an index for the string and then
        randomly delete 1-9 nucleotides sequentially and return the new RNA string
    
    find_codon_diff(protein, mutated_protein)
        Function to find the effect the SNP had on the protein.
        First compares the length of the proteins to see if the SNP caused a premature
        stop codon, otherwise, returns the new codon.
    
    find_deletion(rna, deleted_rna)
        Function to find the de[...]<jupyter_text>For the next portion of my project, I wanted to mess around with RNA and various mutations of RNA sequences. A goal was to create a program that given two sequences, could detect the mutation that occured between them. In order to test the various mutation functions, I created a random RNA generator. I did not have the time or the willpower to build these into an RNA class to mirror the DNA class. Maybe over christmas break...The random RNA function takes a length input, begins the RNA string with a start codon, ends with an arbitrarily chosen stop codon, and then fills in the middle with random codons, adjusting the end length of the string based on the length % 3, so there won't be any nucleotides that can't translate to an amino acid code.<jupyter_code>rna1 = mutation.random_rna(60)
rna1<jupyter_output><empty_output><jupyter_text>The first type of mutation function I built was a SNP mutator. A SNP, or single nucleotide polymorphism, is where a single nucleotide changes to a different nucleotide. These can be transitions (A to G, C to T/U or vice versa, aka purine to purine, pyrimidine to pyrimidine) or transversions (purine to pyrimidine, or vice versa.) The code gives equal weight to all possible mutations, although in nature, transitions occur more frequently due to the more similar structure between the nucleotides, but I wasn't sure how to implement the mutation with skewed odds. A while loop ensures that the replaced nucleotide isn't replaced with itself, looping back through until the randomly chosen nucleotide has a different value. A SNP may be silent, due to the redundancy in the amino acid table, may be a missense mutation, where an amino acid will change values, or a nonsense mutation, where a premature stop codon would be generated, truncating the protein.<jupyter_code>mutation.snp(rna1)<jupyter_output><empty_output><jupyter_text>I was finding it hard to eyeball the difference between the original and mutated string, so I built a SNP finder function.<jupyter_code>mutated_rna1 = mutation.snp(rna1)
mutation.find_snp(rna1, mutated_rna1)
protein1 = DNA.translate(rna1)
mutated_protein = DNA.translate(mutated_rna1)
mutation.find_codon_diff(protein1, mutated_protein)<jupyter_output><empty_output><jupyter_text>The next mutations were insertions and deletions. Depending on the amount of nucleotides that are inserted or deleted, as well as where they are inserted or deleted, drastic effects can happen to the protein after processing. Frameshift mutations can generate premature stop codons, and change many amino acid residues, greatly affecting the folding and function of a protein. Functions were also built to detect the insertions and deletions.<jupyter_code>mutated_rna2 = mutation.insertion(rna1)
print(rna1)
print(mutated_rna2)
mutation.find_insertion(rna1, mutated_rna2)
mutated_rna3 = mutation.deletion(rna1)
print(rna1)
print(mutated_rna3)
mutation.find_deletion(rna1, mutated_rna3)<jupyter_output><empty_output><jupyter_text>I then built a mutation randomizer, which will take an rna string and mutate it with any of the previously used functions. Again, the odds for the random choice were the same for all mutation types, as I couldn't find a way to skew the odds. SNPs should be much more common than deletion or duplications. Similarly, a mutation detection program was written to hopefully identify the mutation.<jupyter_code>mutated_rna4 = mutation.mutator(rna1)
print(rna1)
print(mutated_rna4)
mutation.find_mutation(rna1, mutated_rna4)<jupyter_output><empty_output><jupyter_text>Finally, I created a data set generator, data_gen(X) where an int is called and then the program will generate X RNA string, mutate them, find the mutation, translate both the original and mutated RNA, and then pickle the whole result into the gzipped "rna_mutation_data.pickle.gz" file.<jupyter_code>mutation.data_gen(15)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.04778698642214777, 'fraction_numerical': 0.00811144418973726, 'mean_word_length': 4.097977528089888, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 68, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '230310', 'license': 'MIT', 'n_tokens_mistral': 3563, 'n_tokens_neox': 3290, 'n_words': 1603, 'path': 'DNA Manipulation/Ben_Clyde_Project_Walkthrough.ipynb', 'repo_name': 'beniclyde/beniclyde.github.io', 'size': '11353'}