CSC 15 DNA Lab 

The target for this lab is to write some procedures for string
matching.  We will apply these procedures to DNA sequence analysis.
Some of the procedures will also involve file I/O.

######### Prelim: Some operations on strings that may be useful:
# Given some string
s = "abcdefg"
# the expression s[a:b] will return the substring of s from 
# s[a] to to s[b-1], so s[2:5] will return "cde".  The length of the substring
# is always b-a.  The expression s[:3] is the same as s[0:3] and s[3:] is
# the same as s[3:len(s)].  For example:
s = s[0:2] + 'x' + s[3:]  # will change 'c' to 'x', s becomes "abxdefg"
# However, for large strings, it would be more efficient to change individual
# characters by first converting the string into an array of characters.
#########

                            PART I:

0. Write a procedure that generates a random DNA sequence of length n:

def randomDNA(n):

   Think about the algorithm first: you'll start with an empty string "",
   then run a loop n times.  Inside the loop, you need to generate a 
   randint(0,3) and depending on the value of this random number, you
   will add either an "A", "C", "G" or "T" to the end of your string.
   For example, to add the character "C" to the end of string S, use the
   statement S = S+"C".

   The function should return the string generated.

   To use the randint function, put  *from random import randint*  at the
   top of your program.  Then randint(a,b) generates a random number between
   a and b inclusively.



#######
1. Write a function count(X,S) that returns the number of times that the
character X appears in string S.  For example, count("A","GATCAATC") should
return 3.  Test it on your random DNA sequences.


#######
The following procedures may require procedures to transform strings to
arrays of chars and vice versa.  Transforming a string into an array of chars
will make it easier to modify the string, as long as you transform it back
into a string at the end.

def stringtoarray(S):
        A = []
        i = 0
        while i<len(S):
            A.append(S[i])
            i = i+1
        return A
#string to array: usage : A = stringtoarray("ACGT")

def arraytostring(A):
            S = ""
            for c in A:
                S = S+c
            return S
#arraytostring: usage: s = arraytostring(['A','C','T','G'])

# Alternatively, you can do
from string import join
#...
s = join(['A','C','T','G'],"") # this will also set s to "ACTG"


######
2. A DNA sequence mutates when some of its code (symbols) change randomly.
   Write a procedure that randomly mutates up to N of the symbols in the
   DNA sequence.  For example, mutate("ATTGCGA",2) may return
  "ATAGCGC".  

        Algorithm: first convert the DNA string into an array of chars
        so that it can be easily modified.  Then run a loop n times.
        inside the loop, generate a random array index, and modify the
        array content at the index.  Calling your randomDNA function
        with an argument of 1 should return a random letter: "A", "C",
        "G" or "T".  After the loop is finished, convert the array
        back to a string and return the string.

        Note that it's possible that a mutation could end up with the
        same letter as before.  That's why the function is specified
        to mutate "up to" n of the symbols.

        If you convert the string to an array first, remember to convert
        it back to a string before returning it.

    CHALLENGE:  Write another version of this function that mutates EXACTLY
    N codes.


def mutate(dna,n):
dnaprep.txt
######
3. Write a function "splice":  (READ THE SPECS CAREFULLY!)

  def splice(A,B,j)

  Where A is a DNA string, B is a DNA string and j is a number.

  The function needs to replace a segment of B with the contents of A.
  j will indicate the starting index in B where this replacement will
  take place.  For example, splice("AAA","ATGTTGC",2) should return the
  string "ATAAAGC": the contents of string B from index 2 (third code) 
  is replaced by string A.  Please note that you're not "inserting" A into
  B, but replacing a substring in B with something else.  The length of B
  should remain the same after the operation.

  Your function also needs to detect error conditions: it's possible that
  the length of A is larger than the length of B.  It's also possible that
  j+len(A) is larger than the length of B.  In either case, your function
  should return B unchanged.
   
  This time, you'll have to come up with the algorithm yourself.


######
4. Using your functions above, generate several dna sequences large and
   small.  mutate and splice the smaller sequences into the larger ones.
   Write the dna sequences to file.  Please do not create DNA sequences 
   with lengths longer than 2**20 (1 megabyte).


You can use the following function to write a string to a file:

# write a string S to a file:
def writefile(S, filename):
    fd = open(filename,"w")
    fd.write(S)
    fd.close()
# write to file (function does not return value, just writes to file)

and the following function to read from a file:

def readfile(filename):  # return a string containing the file contents
    fd = open(filename,"r")
    S = fd.read()
    fd.close()
    return S
# read from file and returns string read (could be very big string)

Given string s, writefile(s,"file.txt") will write string to named file
t = readfile("somefile.txt") will read string in named file and return it.


############################################################################

#   Part II:

After years of painstaking research geneticists have identified the
DNA sequence responsible for causing some computer science students to
doze off during class.

The Dozer Gene:

GTGTGCAGCGTCCCCGATGCTACTCCATTCCTTTCGCGTGTATCTTGAGAGCAATGATGGGATTTACAATACAGATACGTAAACGCAACGTATCAACAGGTAGGGGTTCTGTTTGGAGCCAACTAGTC


Four DNA samples were randomly taken from students at an unnamed university in
Nassau county.  Their files are found on the homepage as sd1.dna, sd2.dna,
sd3.dna, and sd4.dna.

You are to determine the probability that a student carries the dozer gene.
This amounts to finding if one string A is a substring of another string B,
WHILE FORGIVING A NUMBER OF ERRORS.  For example, "abcde" matches "abudg"
with 2 errors.  You need to write a function that takes two strings and
a number indicating the maximum number of errors allowed.  The function
needs to determine the best possible match.  It should compute two values:
the starting position in the longer sequence where the best match occurs,
and the number of errors in the best match.  If there is no match that
contains less than or equal to the maximum number of allowed errors, then
the starting position value should be -1.  For example, calling

 subsequence("defgh","abcefkhdfsfdefgms",3)   // 3 indicates max errors

should return the tuple (11,1): that is, the starting position of the
best match is 11 and the number of errors is 1, indicating the
substring "defgm", which occurs starting at index 11 (indices start
from 0) in the longer string.  Note that the substring "cefkh", which
starts at position 2, is also a match, but has 2 errors.  Your
function needs to report the best possible match.

One way to write the subsequence function is to modify the following
function, which only checks if a is a substring of b and report the
position in b where a is found.  It does not tolerate errors:

# returns starting position of string a in b, -1 if not a substring
def substringi(a,b):  # 
    answer = -1  # assume worst case for there-exists type of loop
    j = 0  # index b
    while j<=len(b)-len(a) and answer==-1: # there exists j in b loop
        bx = j  # assume a is found in b starting at b[j]
        i = 0   # i.e., assume best case for bx in a for-all type of loop
        while i<len(a) and bx==j: # for all i in a
            if a[i] != b[j+i]: bx = -1
            i +=1
        # while i
	answer = bx   # it's possible to just use one variable (answer)
        j+=1          # but then it will be harder to generalize algorithm.
    # while j
    return answer
## substringi

If the entire dozer gene is found embedded in a dna file, with no
errors, then that person is considered to have a 100% chance of being
a dozer.  In general, if there are n errors, then the chance that the
carrier is a dozer is (100 - 2**(n-1))%.  So 1 error (n=1) means
there's still a 99% chance, 2 errors means 98%, 3 errors means 96%, 4
means 92%, 5 means 84% (100-16), 6 means 68% and 7 means 46%.  If there are
more than 7 errors then the sample is considered to have no significant 
chance of causing the dozer effect.  So your experiment should tolerate up
to a maximum of 7 errors. (Please note this is not necessarily the method
used in real DNA matching - ask your biology teacher for that!)

Your program should produce output in the following FORMAT (exact values
will differ)

DNA SAMPLE 1:
A best match with 4 errors was found starting at position 2018.
There is a 92% chance that this student is a dozer


#########################################################################

While conducting their research on the dozer gene, our brilliant
geneticists noticed another phenomenon: many of the students who were
dozing off during class were also drooling.  There's sticky sliva found
all over the desks, from which DNA samples were extracted.  The gene that 
causes this behavior was also identified:

The Drooler Gene:

ACTAAGGTTGCGTCGATTACAAAAATGGCGCTTGTAGCTTTACAAGCTACGTTTTTGCCGTGCTTGTAGGATAGAGTCAGGCGGAGCGGGCCTAAGAGTGTGAATACATGGGCGAGTCGGGGGATTATCTTTAGGTAAGGCATGATACCAGTTTACTCCTTCTTCTGTCCGATGTGTGACCTAGTCGGTGGAAATGATTGTTTATCGAGACCCATTGTAGACTGAAGCGAAACTTTGGTGTTCTGAGGTCTCAAGT


One geneticist boldly proposed the following, controversial hypothesis:

The "D & D" Hypothesis:  

   "Every carrier of the drooler gene also carries the dozer gene"


You need to modify your experiments to either lend proof to the hypothesis,
or disprove it by finding a counter-example (i.e., a student who drools
while awake). 

Report your findings.

----------------------------------------

EXTRA!

The work of a former, famous gneticist has been found.  After years of
working alone, he had discovered a gene that causes some computer
science students to forget everything they've learned right after a
class is over.  Unfortuantely, he left the gene stored in binary
format, in a a file called "forgetful.dna".  Each sequence of 16 DNA
symbols is encoded as a 32-bit signed (two's complement) integer.  The
first 32-bit integer in the file is the length of the DNA sequence (in
number of symbols). The sequence is padded with 'A' symbols so the length on
file is a multiple 16.  A is represented by 00, C by 01, G by 10 and T by
11.  However, if the first symbol is G or T, it is interpreted as -2
and -1 respectively.

You are to extract the forgetful DNA sequence as a string of A,C,G,T symbols.
Only then will we be able to formulate a cure for this frustrating disease!

### EXTRA EXTRA!

  Could it be that the forgetful gene is related to the dozer/drooler genes?
Find out and get the next Nobel prize!