Notes of Bioinformatics V

Notes of Genomic Data Science and Clustering (Bioinformatics V)

Chapter 9: How Did Yeast Become a Wine Maker? (Clustering Algorithms)

In this chapter, we learned the clustering algorithms in bioinformatics study. 

Some concepts:

gene expression analysis.

Good clustering Problem.

k-Means Clustering Problem: optimization problem

- Farthest First Traversal Problem

- Squared error distortion

- Lloyd Algorithm for k-means clustering: two steps: centers to clusters and clusters to centers.

- limitations: make a "hard" assignment of each point to only one cluster.

Soft k-means clustering:

- Expectation Maximization Algorithm: starts with a random choice of Parameters. It then alternates between the E-step, in which we compute a responsibility matrix HiddenMatrix for Data given Parameters, and the M-step, in which we re-estimate Parameters using HiddenMatrix.

- Centers to Soft Clusters (E-step) AND Soft Clusters to Centers (M-step)

- SoftKMeans Problem

Hierarchical Clustering

- how to use distance matrix to partition genes into clusters

- UPGMA in disguise

The Python implementations of all the algorithms are available on my github: https://github.com/aprilchunyuzhao/BioinformaticsFromCoursera.

Notes of Bioinformatics III

My NOTES of  Bioinformatics III: Comparing Genes, Proteins, and Genomes

Chapter 5: How Do We Compare Biological Sequences? (Dynamic Programming)

For my PhD research, I developed a novel pairwise protein structure alignment algorithm (UniAlign: http://sacan.biomed.drexel.edu/unialign/). And part of it is to calculate the sequence similarity between two protein sequences by glocal alignment. So, I am quite familiar with dynamic programming.

What is new to me in this Chapter is that, it applies dynamic programming to the problem of finding the longest common subsequence. Finding a longest common subsequence of two strings is equivalent to finding an alignment of these strings maximizing the number of matches.

The matches in an alignment of two strings define a common subsequence of the two strings, or a sequence of symbols appearing in the same order (yet not necessarily consecutively) in both strings.

An introduction to Dynamic Programming: the Change problem

1. The trick behind dynamic programming is to solve the smaller problems once rather than billions of times.

2. To find a longest path in an arbitrary DAG, first we need to order the nodes of the DAG so that every node falls after all its predecessors. 

Computational Problems:1. Manhattan Tourist Problem2. Output LCS Problem3. Longest Path in a DAG problem

Some Concepts:

1. scoring matrices
2. global alignment

3. local alignment

- when compute the s(i,j), add zero-weight edges from (0,0) to every node, which means that the source node (0, 0) is a predecessor of every node (i,j) (free taxi ride).

4. alignment with affine gap penalties problem

The Changing Faces of Sequence Alignment:

(1). Edit Distance: the minimum number of edit operations(insertion, deletion or substitution) needed to transform one string into another.

Hints: do a global alignment with mismatch penalty = -1 and gaps = -1, match = 0.

(2). Fitting Alignment: find a region within the longer protein sequence v that has high similarity with all of the shorter sequence w. In other words, "fitting" w to v required finding a substring v' of v that maximize the global alignment score between v' and w among all substrings of v.

Hints: "free taxi rides" only for the second string.

(3). Overlap Alignment: refers to the global alignment of a suffix of v with a prefix of w. 

Hints: "free taxi rides" only for the first string.

These three applications of sequence alignment using dynamic programming are very nice!!

Affine Gap Penalties:

- Build Manhattan on three levels: lower, middle and upper.

Space-Efficient Sequence Alignment:

1. The memory required to store the dynamic programming matrix is substantial O(n*m).

2. divide-and-conquer algorithms: divide phase splits a problem instance into smaller instances and solves them; conquer phase stitches the smaller solutions into a solution to the original problem. 

3. The idea of space reduction techniques is that we don't need to store any backtracking pointers if we are willing to spend a little more time.

4. The Middle Node Problem!! 

- In general, at each new step before the final step, we double the number of middle nodes found while halving the runtime required to find middle nodes.

5. The Middle Edge Problem

- middle edge: an edge in an optimal alignment path starting at the middle node (more than one middle edge may exist for a given middle node).

Multiple Sequences Alignment

- todo

Chapter 6: Are There Fragile Regions in the Human Genome? (Combinatorial Algorithms)

Biology:

reversal: the most common form of genome rearrangement.

Random breakage model of chromosome evolution: the breakage points of rearrangements are selected randomly.

- exponential distribution: of synteny block lengths

Concept:

reversal distance: the minimum number of reversals required to transform P into Q.

identity permutation.

Breakpoint Theorem.

Rearrange in multi-chromosomal genomes: 

- reversals and translocations, as well as fusions and fissions.

- genome graph

- breakpoint graph

- Cycle Theorem: given genome P and Q, any 2-break applied to  P can increase Cycles(P,Q) by at most 1.

- 2-Break Distance Theorem: the 2-break distance between P and Q is equal to Blocks(P)-Cycles(Q).

Fragile Breakage Model: every mammalian genome is a mosaic of long solid regions, which are rarely affected by rearrangements, as well as short fragile regions that serve as rearrangement hotspots and that account only for a small fraction of the genome. For humans and mice, these fragile regions make up approximately 3% of the genome.

Synteny Block Construction:

- construct synteny blocks from shared k-mers

- construct synteny blocks as connected components in graph

Computational Problem:

1. Greedy Sorting by Reversal Problem

2. 2-Break Distance Problem

3. 2-Break Sorting Problem

4. Finding Shared K-mers Problem

All the python codes are available on my github: https://github.com/aprilchunyuzhao/BioinformaticsFromCoursera  

Notes of Bioinformatics II

My NOTES of  Bioinformatics II: Genome Sequencing offered by UCSD on Coursera.

Chapter 3: How Do We Assemble Genomes? (Graph Algorithms)

Difficulty: 

1. double stranded DNA: no way of knowing which strand a given read derives from

2. sequencing errors

3. some regions of the genome may not be covered by any reads

4. Repeats complicate genome assembly: approximately 50% of human genome is made of repeats, e.g., the 300 nucleotide-long Alu sequence is repeated over a million times.

Some Concepts:

reads: sequence much shorter DNA fragments.

genome assembly: put a genome back together from its reads using overlapping information!

overlap graph: form a node for each k-mer in Patterns and connect k-mers Pattern and Pattern' by a directed edge if Suffix(Pattern) is equal to Prefix(Pattern').

hamiltonian path: a path in a graph visit every node once.

k-universal string: contains every binary k-mer exactly once.

de Bruijn graph:assign each k-mer in Patterns as edges, prefix and suffix of each edge as the node, and then glue identical nodes together.

eulerian path: visit every edge exactly once.

NP-complete: can be verified in polynomial time, yet no fast solution. (non-deterministic polynomial time). All NP complete problems are equivalent to each other.

Computational Problems:

(1) String Composition Problem

(2) String Spelled by a Genome Path Problem

(3) Overlap Graph Problem

(4) de Bruijn Graph Problem

(5) de Bruijn Graph From k-mers Problem

Concepts:

Eulerian Cycle: a cycle that traverse each edge of a graph exactly once.

Euler's Theorem: every balanced, strongly connected direct graph is Eulerian.

read pairs: pairs of reads separated by a fixed distance d in the genome.

Paired de Bruijn Graph:

Genome Assembly Faces Real Sequencing Data:

1. Breaking reads into shorted k-mers: yet smaller k may result in a more tangled de Bruijn graph

2. Splitting the genome into contigs (long, continuous fragment of the genome)

- maximal non-branching path: a non-branching path that cannot be extended into a longer non-branching path.

Computational Problems:

(1) Eulerian Cycle Problem

- STACK!!

(2) Eulerian Path Problem

(3) String Reconstruction Problem: reduced to finding an Eulerian Path in the de Bruijn graph generated from read!

(4) k-Universal Circular String Problem: reconstruct a circular string given its k-mer composition.

(5) String Reconstruction from Read-Pairs Problem

(6) Contig Generation Problem

Chapter 4: How Do We Sequence Antibiotics? (Brute Force Algorithms)

A difficult problem in antibiotics research is that of sequencing newly discovered antibiotics, or determining the order of amino acids making up the antibiotic peptide.

Necessary Biology:

1. There are three different ways to divide a DNA string into codons for translation, one starting at each of the first three starting positions of the string (reading frames).

2. Protein translation is carried out by a molecular machine called ribosome.

3. non-ribosomal peptides (NRPs): synthesized not by the ribosome, but by a giant protein called NRP synthetase.

NRPs have pharmaceutical applications because they have been optimized by eons of evolution as "molecular bullets" that bacteria and fungi use to kill their enemies. If these enemies happen to be pathogens, then researchers are eager to borrow these bullets as antibacterial drugs.

4. mass spectrometer: an expensive molecular scale that shatters molecules into pieces and then weighs the resulting fragments, in daltons (Da); 1 Da is approximately equal to the mass of a single nuclear particle.

The mass spectrometer can break each molecule of Tyrocidine B1 into two linear fragments, and it analyzes samples that may contain billions of identical copies of the peptide, with each copy breaking in its own way.

Our goal is to use all of these different fragments to sequence the peptide.

5. subpeptides: for a cyclic peptide with n amino acids, there are n * (n-1) subpeptides.

6. Cyclospectrum(Peptide): theoretical spectrum; the collection of all masses of its subpeptides, as well as mass 0 and mass of the entire peptide.

Assumption: for a given linear peptide Peptide, the mass of any subpeptide is equal to the difference between the masses of two prefixes of Peptide!!

7. For a given cyclic peptide: its theoretical spectrum are those found by LinearSpectrum and those corresponding to subpeptides wrapping around the end of Peptide.

Again, our goal is to reconstruct unknown peptide from its experimental spectrum.

8. A branch-and-bound algorithm for cyclopeptide sequencing: "grow" candidate linear peptides whose theoretical spectra are "consistent" with the experimental spectrum.

- branching step: to increase the number of candidate solutions

- bounding step: to remove hopeless candidates

9. Real Spectra: mass spectrometers generate "noisy" spectra: false masses and missing masses. THUS, in our algorithm design, instead of finding exact match, we need a scoring function that can score how match between the given theoretical spectrum and the given experimental spectrum. 

10. Replace Peptides with Leaderboard: hold the N highest scoring candidate for further extension (including ties).

11. Since there are more than 100 NRPs, we MUST determine the amino acid composition of a peptide from its spectrum so that we may run LeaderboardCyclopeptideSequencing on this smaller alphabet of amino acids.

12. spectral convolution: take the positive differences of masses in the spectrum.

13. Epilogue: from simulated to real spectra!

Computational Problems:

(1) Protein Translation Problem:

(2) Generating Theoretical Spectrum Problem: LinearSpectrum and CyclicSpectrum

(3) Counting Peptides with Given Mass Problem: dynamic programming algorithm

(4) Cyclopeptide Sequencing Problem: still not efficient enough

Notes of Bioinformatics I

My NOTES of  Bioinformatics I: Find Hidden Messages in DNA offered by UCSD on Coursera.

Chapter 1: Where in the Genome Does DNA Replication Begin?

In this chapter, we learned how to find the replication region (oriC) in bacterial genomes (mostly single circular chromosome).

Biology:

1. Feature of oriC of bacterial genome: typically a few hundred nucleotides long.

2. The initiation of DNA replication is mediated by DnaA, which is a protein binds to a short segment within the oriC region called DnaA box. And the DnaA box is the 'hidden messages' that we are trying to look for from the DNA sequence. 

The computational translation for this problem is: find the most frequent k-mer in the given oriC sequence ( if it maximizes Count(text, pattern) among all k-mers).

Related Computational tasks:

1. PatternCount(text, pattern): count the occurrence of pattern in text

2. FrequentWord(text, k): find the all most frequent k-mers in text

Charging Station:

The running time of the brute force implementation of FrequentWord problem is O(|Text|^2). To improve the algorithm, the data structure frequency array was introduced.
1. Frequency Array: given a integer k, the length of the frequency array is 4^k (since 4^k different integers for all 4^k k-mers lexicographically), and the i-th element of the array hold the number of occurrences of i-th k-mer in lexicographic order in text. In other words, for each k-mer, we need to calculate a unique integer for that k-mer used as the index in the frequency array. 
Considering the DNA sequence, we need a function SymbolToNumber(symbol) to transform symbol A,C,G,T to 0,1,2,3.

• PatternToNumber(pattern): transform a k-mer pattern into an integer (recursion)

PatternToNumber(pattern) = 4 * PatternToNumber(pattern[:-1]) + SymbolToNumber(pattern[-1])

• NumberToPattern(index,k): transform an integer between 0 and 4^k-1 into a k-mer
• FasterFrequentWords(text, k): running time O(4^k + |text|^k + 4^k) => still impractical when k is large.

2. FindingFrequentWorksBySorting(text,k): sorting can bring same things together! 

After learning how to find the 'hidden messages' from the DNA sequences, we continue our journal based on the observation (statistically interesting) that some hidden messages are more surprising than the others.
First some biological concepts: 

• each DNA strand has a direction: read in the 5' -> 3' direction
• DnaA protein does not care which of the two strands it binds to when the DnaA protein binds to DnaA boxes and initiates the replication.

Followed by some statistical concepts:

• overlapping words paradox: different k-mers have different probabilities of occurring multiple times as a substring of a random string.
• Pr(N,A,Pattern,t) is a complex problem because this probability depends heavily on the particular choice of Pattern.

Before we conclude that we have found the DnaA box, we need to check whether there are other short regions in the bacterial genome exhibiting frequent occurrence of that k-mer.
clumps: appear close to each other in a small region of the genome; we defined a k-mer as a "clump" if it appears many times within a short interval of the genome; (L, t) - clump

Related Computational tasks:
3. ClumpFinding(Genome, k, t, L)

Furthermore, we continued to learn a new algorithm to find the oriC based on some features of DNA replication biology.
Important Biology:
1. DNA polymerase does not wait for the two parent strands to completely separate before initiating replication; instead, it starts copying while the strands are unraveling.
2. To start replications, a DNA polymerase needs a primer, a short complementary segment that binds to the parent strand and jump starts the DNA polymerase.
3. DNA polymerase is unidirectional: only traverse a template strand of DNA in the 3' -> 5' direction, which is opposite from the 5' -> 3' direction of DNA.
4. Asymmetry of Replication: replication on reverse half-strand progresses continuously; yet a DNA polymerase on a forward half-strand has no choice but to wait again until the replication fork has opened another 2000 nucleotides or so. -> Okazaki fragments from multiple primers
5. Deamination: cytosine (C) has a tendency to mutate into thymine (T) through a process called deamination. Deamination rates increase 100-fold when DNA is single-stranded, which leads to a decrease in cytosine (C) on the forward half-strand.

THEREFORE, we can use the skew diagram #G-#C to find the oriC <- if this difference starts increasing, then we guess that we are on the forward half-strand; if this difference starts decreasing, then we guess we are on the reverse half-strand.
The skew should achieve a minimum at the position where the reverse half-strand ends and the forward half-strand begins, which is exactly the location of oriC!
Related Computational tasks:
4. Minimum Skew Problem


The LAST computational task of this chapter is Approximate Pattern Matching Problem! And the following are some thoughts:
First, we define a d-neighborhood Neighbors(Pattern, d) of a k-mer Pattern, as the set of all k-mers that are close to Pattern. (Recursive)
Then, for each neighbor only, we update the frequency array.
The approximate pattern matching by sorting and pattern matching by sorting employ the same idea that sorting bring same things together so we can count their appearance by comparing with the neighbor element....Remember this!



Chapter 2: Which DNA Patterns Play the Role of Molecular Clocks? (Randomized Algorithms)
In Chapter 1, we learned one type of hidden message in genomes, however, there are other types of hidden messages that have nothing to do with genome replication, for example, regulatory DNA motifs responsible for gene expression, that we will learn in Chapter 2. 
BIOLOGY:
1. A transcription factor regulates a gene by binding to a specific short DNA interval called a regulatory motif, or transcription factor binding site, in the gene's upstream region, a 600-1000 nucleotide-long region preceding the start of the gene.
2. Regulatory motif are not completely conserved, meaning they can vary at some positions for different genes.
THUS, comes our new COMPUTATIONAL task: develop algorithms for motif finding: discover the "hidden message" shared by a collection of strings.

The DIFFERENCE between Frequent Words Problem and Motif Finding Problem is that: a DnaA box is a pattern that clumps, or appears frequently, within a relatively short interval of the genome; WHILE a regulatory motif is a pattern that appears at least once (with variation) in each of my different regions that are dispersed throughout the genome.

(k,d) motif: a k-mer appears in every string from Dna with at most d mismatches.
Observation: any (k, d) motif must be at most d mismatches apart from some k-mer appearing in one of the strings of Dna.
(1) A brute force algorithm for motif finding.
Limitation: as long as a single sequence does not contain the transcription factor binding sites, a (k,d) motif does not exist!

Another model: select a k-mer from each string and score these k-mers of how similar they are to each other (motif matrix).
- define the Score(Motifs)
- minimize this score
Motifs -> Score -> Count -> Profile -> Consensus String
- Entropy: the entropy of completely conserved column is 0, and 2 for equally-likely nucleotides column.
- Motif logo
d(Pattern, Dna): the sum of distances between Pattern and all strings in Dna
(2) Median String Problem
Limitation: since median string problem has to consider 4^k k-mers, so toooo slow for cases like k = 15.
(3) Greedy Motif Search
Observation: a k-mer has a higher probability when it is more similar to the consensus string of a profile.
Improvement: pesudocounts
(4) Randomized Motif Search
- Monte Carlo Algorithm: not guaranteed to return exact solutions, but they do quickly find approximate solutions.
- Las Vegas Algorithm
- continue to iterate for as long as the constructed motifs keeps improving
Advantage: Randomized Motif Search can find longer motifs.
(5) Gibbs Sampling
Limitation: local optimum