239 Discrete Mathematics for the Life Sciences

1 239 Discrete Mathematics for the Life SciencesLior PachterNotes by Qiaochu YuanSpring 20131 IntroductionThere are 20 (actually 21) amino acids which are coded from a 4-letter alphabet. Thisis the largest possible size of a comma-free code on 4 letters, so it was tempting toconclude that amino acids were coded via a comma-free code. However, this is nottrue. In fact, amino acids are coded by multiple 3-letter words (codons); usually 4,but sometimes 3,2, or 1 (as in tryptophan). This genetic code is very nearly universal(there are small variations in a handful of species).Codons which code for the same amino acid are called synonymous. Synonymouscodons are usually the same in the first 2 letters and differ only in the third; thisseems like it might be useful because mutations in the third letter don t affect theresulting amino explorer ( ) is a tool for describingfrequencies of codons in various genomes.

2 It returns, among other things, the CodonAdaptation Index, which roughly speaking describes how close the codons are to beingthe most common synonymous mystery of the two straight lines is the following. If we measure how frequentdifferent nucleotides are in different positions in codons, this correlates linearly verywell with GC-content (percentage of bases that are either guanine or cytosine).How do we fit lines to points in general? One way is the least-squares approxima-tion. (It is easy to do such things with R.) Least-squares can be done by minimizingthe sum of the squares of the vertical or horizontal distances (if we want to fit lines topoints in a plane). Another is using principal component analysis, which minimizesperpendicular distances to the data points instead.

3 One way to interpret the PCAline is that, after projecting all the points onto the line, it maximizes the remainingvariance in the human genome consists of billion sites and there are billion can be combined into an enormous matrix on which it would then be possibleto perform PCA. Part of this matrix has been collected (for Europeans) and the twolargest principal components have been plotted; the result agrees very well with thegeography of not only determine what amino acids are made but how many. Thiswas tested using a synthetic library of GFP (green fluorescent protein) genes, whosecodons were randomly replaced by synonyms and then let loose to observe how of-ten the corresponding proteins were produced. It turns out that replacing codons1with synonyms determines how RNA molecules (which are single-stranded) fold in onthemselves.

4 Different folding structures have different free energies, and this correlatesto fluorescence (hence how much of the protein is created).Nucleotides can be modified by a process called methylation, which adds a methylgroup to either cytosine or adenine. It s known that methylation affects transcrip-tion and therefore translation rates. Methylated cytosines also occasionally The genome as a random objectSuppose we have a collection ofNsequences (of either nucleotides or genes). Wesuspect that they may have a subsequence in common of lengthw, and the problemis to determine how likely this is. We don t demand that the subsequences are exactlythe same, but only that they are similar. This is relevant to binding site , we will think of a genome as randomly determined by proba-bilitiespiof a nucleotide appearing in a particular position, as well as probabili-tiesqj,A,qj,C,qj,G,qj,T(the position weight matrix or position-specific scoring matrix)wherejindexes positions in the motif.

5 These probabilities determine how likely nu-cleotides are to appear when the genome switches to expressing the motif. This is avaguely reasonable model because, although mutations occur randomly, binding sitesneed to continue to exist. (This isn t a great model but a lot of people use it.)Weblogo is an online resource for visualizing binding sites. There are columnsof nucleotides for each binding site, and the height of a given letter, sayA, isqi,ARiwhereRi= log24 (Hi+ei) and whereHi= qi,Alog2qi,Ais the Shannon entropyandeiis a correction for small sample want an algorithm that will learn thepiand theqj, , which we do as follows(assuming that we knoww). We first choose a random collection of starting positionsak(wherekruns over theNsequences). Next, we choose a random sequencezandlearn thepiandqj, by counting starting from theakin every sequence other these learned probabilities, we then compute the probability of every possiblestarting positionazin sequencezand pick a random start site inzaccording to thatdistribution.

6 We then repeat. This algorithm converges to the correct probabilitiesand is a special case of the Metropolis-Hastings not known, the problem of determiningwis the problem of model selec-2tion. It is necessary to penalize for longerwbecause of the existence of additionalparameters (which leads to overfitting), but none of the standard method for doingthis appear to work in this drawback of this algorithm is that it assumes that motifs appear one in eachsequence. This is usually false. We can address this by adding new parametersdetermining how often motifs appear in each variant of the algorithm involves cycling through the sequences rather thanchoosing them at random.(Lecture on Metropolis-Hastings and Gibbs sampling not given)3 SignificanceHow do you tell if a genome feature is significant?

7 Ak-meris a sequence ofkbase pairs. We want to compute probabilities like theprobabilities of seeing copies of somek-mer in a row given that we are looking at arandom genome in such a way that preservesk-mer frequencies. This is importantbecausek-mer frequencies are not uniform in real genomes.(In a plot wherexdescribes a givenk-mer frequency andydescribes how manyk-mers appear at a given frequency, we expect something like a normal distributionand don t see that; instead the distribution for 11-mers in the human genome has twopeaks to the left of the naive mean and has heavy tails.)One way to model the genome is using Markov chains. Let be a finite alphabet(here it will just be{0,1}). Let be a set of parameters (here 00, 01, 10, 11)between 0 and 1; these describe transition probabilities.

8 Let`be a length (here 4).Our Markov model will bepijkl= ij jk klz(1)wherez= i,j, ij jk kl.(2)The distributions we get in this way satisfy various relations. For example,p0010=p0100=p1001. Less obviously, we have quadratic relations likep20011=p0001p0111, andnine additional relations likep0111p1010=p0101p1110which turn out to generate the3ideal of all polynomial relations and which describe all probability distributions weget in this way. This is a special case of the Hammersley-Clifford above is an example of a log-linear model. In general a log-linear model takesthe form( 1,.. d)7 ( a1,.. am) aj(3)where aj=d i=1 aiji(4)andaijis an arbitrary matrix. Markov models are log-linear, and log-linear modelshave the property that maximum likelihood estimation is easy, which is nice.

9 (Thisis made precise by Birch s theorem, which will be on the homework.)Here is an algorithm for generating random sequences from a given sequence in away that preserves 2-mer counts. The idea will be to randomly swap subwords in away that preserves 2-mers. Ifsis a sequence andka constant, then we denote byXk(s) the set of sequences of the same length asswith the samek-mer frequenciesass, and the goal is to sample fromXk(s) more or less will do this by choosing four locationsa < b < c < din the sequence andtesting whether it is possible to swap the substring betweenaandbwith the substringbetweencanddwithout alteringk-mer also denote byYk(s) the subset ofXk(s) with the same firstk 1-mer. Wesee thatXk(s) = s Yk(s )(5)wheres runs over all cyclic permutations ofs.

10 We will sample fromYk(s) forsimplicity. Define a graphGk(s) whose vertices are the elements ofYk(s) and whoseedges are given by shuffles (if the shuffle is invalid, the edge is a loop). The firstquestion is whether this graph is connected. The second question is what the mixingtime of a random walk is on this (s)is connected. Its diameter is at most length(s) 2k+ : letj(u,v) be the first position from the left where the sequencesuandvdisagree and let (u,v) =n k+ 2 j(u,v). Then (u,v) = 0 ifu= (s) is regular of degree(n k+14). We want to perform a random walk on thisgraph and determine the mixing time . IfTis the corresponding transition matrix,the answer depends on the second largest eigenvalue , for fixedkand alphabet size, there is a polynomialp(n) such thatsuch that for any sequencesof lengthnand anym, if we performmp(n) steps, thenthe final distribution is less than 2 mfrom uniform in total variation alternative approach is the following.