{'fraction_non_alphanumeric': 0.033543143972591824, 'fraction_numerical': 0.0070512309776113455, 'mean_word_length': 4.629625476564252, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Machine Learning for Genomic Data

This report explores the application of machine learning techniques on short timeseries gene expression data. Although standard machine learning algorithms work well on longer time-series', they often fail to find meaningful insights from fewer timepoints. In this report, we explore model-based clustering techniques. We combine popular unsupervised learning techniques like K-Means, Gaussian Mixture Models, Bayesian Networks, Hidden Markov Models with the well-known Expectation Maximization algorithm. K-Means and Gaussian Mixture Models are fairly standard, while Hidden Markov Model and Bayesian Networks clustering are more novel ideas that suit time-series gene expression data.


Introduction
Machine learning is a field of artificial intelligence that uses statistical techniques to give a computer the ability to learn from data, without explicitly programming it to do so.
Genomics is a branch of molecular biology concerned with the structure, function, evolution, and mapping of an organism's complete set of DNA.
In this project, we apply machine learning techniques to viral genomic data recovered from human cells, and attempt to derive useful insights about the EV-A71 virus, also known as hand-foot-and-mouth disease (HFMD).

Machine Learning
With the large volume of data available nowadays, it is impossible to apply statistical techniques and derive useful insights without making use of a computer.
Considered one of the biggest innovations since the microchip, machine learning enables humans to build models and capture patterns within data. It allows one to define simple rules to improve the performance of these models. The applications of machine learning are manifold -from spam filters in e-mail to online banking, it is omnipresent. Machine learning's presence is growing, and its use cases continue increasing.
A machine learning model can be defined by its training data and the set of hypotheses functions.
Machine learning is broadly divided into the following 3 categories:-• Supervised Learning • Unsupervised Learning

• Reinforcement Learning
In supervised learning, the training data are labelled, i.e. each data point is attached with its true value, or ground truth. Typical problems include classification (binary/multi-class) and regression. E.g. If given a training dataset consisting of pictures of cats and dogs, our model accurately predicts whether a new image is that of a cat or a dog. This is a classification problem. Another simple example is of a dataset with people's heights and weights. Given a new person's height, we should be able to accurately predict weight. This is a regression problem.

INTRODUCTION
In unsupervised learning, unlike supervised learning, the training data are unlabelled. The number of labels may or may not be specified beforehand. Unsupervised learning is commonly used as an initial exploratory technique to draw useful patterns from the dataset. E.g. Given many different books, classifying them into different categories, or clusters, can help us identify the genres they belong to.
Reinforcement learning is different from both supervised and unsupervised learning. It focuses on maximising reward potential from the data by balancing the trade-off between exploration and exploitation. E.g. In recommendation engines, the items showed to a particular user may be ordered based on how much they liked previous products or how different they are from previous products. This is done to maximise the number of clicks, i.e. the reward.
Since the data that we have in this project are unlabelled, our focus is on unsupervised learning.

Genomics
Genomes contain the complete set of genes present in an organism, and genomics describe how this genetic information is stored and interpreted in the cell.

Key Definitions and Terms
• Nucleobases: Components of nucleotides, that are the building blocks of DNA -Adenine (A), Thymine (T), Cytosine (C), Guanine (G). Additionally, Uracil (U) is used instead of Thymine (T) for RNA.
• Codons: A sequence of three DNA or RNA nucleotides that corresponds with a specific amino acid or stop signal during protein synthesis. For eukaryotes, there are a total of 64 such codons, which map to 20 amino acids (and stop signals) through the genetic code.
• Amino acid: Organic molecules that make up the building block of proteins. There are 20 canonical amino acids.
• DNA: Deoxyribonucleic acid, a double helix molecule, that contains the genetic material for all organisms on Earth (including viruses).
• RNA: Ribonucleic acid, a single helix, polymeric molecule that is essential in various biological roles in coding, decoding, regulation and expression of genes. Biologically active RNAs, including transport, ribosomal and small nuclear RNA (tRNA, rRNA, snRNAs) fold into unique structures guided by

INTRODUCTION
complementary pairing between nucleotide bases. Cellular organisms use messenger RNA (mRNA) nucleotides to direct synthesis of specific proteins.
• Protein: Sequences of amino acids that carry out the majority of cellular functions such as motility, DNA regulation, and replication.
• Transcription: The process in which a DNA sequence of a gene is rewritten, or transcribed, to make an RNA molecule. (In eukaryotes, the RNA polymerase makes a strand of mRNA from DNA).
• Translation: The sequence of the mRNA is decoded, or translated, to specify an amino acid sequence (i.e. protein).
• Gene: A region of DNA (deoxyribonucleic acid) coding either for the messenger RNA encoding the amino acid sequence in a polypeptide chain or for a functional RNA molecule.
• Chromosome: Cellular structures that contain genes. A chromosome comprises of a single DNA molecule that may be either circular or linear.
• Virus: A small infectious agent that replicates only inside the living cells of other organisms.
• Gene Expression: The generation of a functional gene product from the information encoded by a gene, through the processes of transcription and translation. Gene products are often proteins. 1 The central dogma of biology states that: DNA is transcribed by RNA polymerases into mRNA (messenger RNA), which is read by ribosomes to generate protein through translation. [9]

INTRODUCTION
In this project, our focus is on translation. More formally, we will be studying the translational interaction between the EV-A71 virus and its host cell, by observing gene expression values at different timepoints.

Virus-Host Interaction
When a host gets infected with a virus, the virus hijacks the host cell's translation machinery, which affects protein production.
In the context of EV-A71, the viral replication process begins when a virus infects its host by attaching to the host cell and penetrating the cell wall or membrane. The viral genome then hijacks the host cell's translation mechanism, forcing it to replicate the viral genome by producing viral proteins to make new capsids (the protein shell of a virus that protects it). The viral particles are then assembled into virions. These virions burst out of the host cell during a process called lysis, killing the host cell. Some viruses take a portion of the host's membrane during the lysis process to form an envelope around the capsid. [1] After infection with any virus, overall translation of all genes decreases. In Figure 2, Following viral replication, the new viruses may go on to infect new hosts. Many viruses such as influenza, chicken pox, AIDS, the common cold, and rabies cause diseases in humans.

Hand, Foot, and Mouth Disease
Hand, foot, and mouth disease (HFMD) is a common infection caused by a group of viruses. The one we will be examining in this report is the enterovirus A71, or the EV-A71 virus. HFMD mostly affects small children, occasionally causing symptoms to develop in adults. It typically begins with a fever, which is followed by rashes and

INTRODUCTION
bumps on other parts of the body. Currently, there is no treatment that specifically targets the disease, nor is there a vaccine that is approved for use in Singapore [3]. Viruses depend on host cells for propagation and dissemination. Host cell mechanisms that allow viral entry, facilitate viral replication, and enable viral egress, are targeted for exploitation by viral pathogens [2].
Guanosine-5'-triphosphate (GTP) is a purine nucleoside triphosphate. It is one of the building blocks needed for the synthesis of RNA during the transcription process, and is involved in energy transfer within the cell. GTPases are a class of proteins that bind and hydrolyze GTP. GTPases include Rab, Ras, Rac, Ran, Rho, Arf, and Sar proteins.
Our focus is on the proteins Rab and Rho. Emerging evidence demonstrates that viruses have evolved numerous strategies to modulate Rab proteins' functions [2].
In polioviruses, the 2C protein possesses weak GTPase activity. The 2C protein in EV-A71 viruses plays an important role in viral replication [27]. Our hypothesis is that a certain member protein of GTPase is important for establishing the infection

INTRODUCTION
in EV-A71 viruses, and could be linked to Rho GTPases. We focus on the translation level changes in attempt to find information supporting our conjecture.

Objective
When a virus infects a cell, most of the processes in the cell get shut down and some proteins do not get generated. There are groups of:-1. Genes that escape this shut down by the virus 2. Genes are affected by the shutdown 3. Genes that are not influenced by the virus [1] is the most important, as these are the ones that the virus needs for its own production.
Thus, our objective is to identify these groups from the set of genes that we have available. We do this by using unsupervised learning techniques on gene expression data to generate clusters.

Data Generation through Genome Sequencing Techniques
RNA sequencing refers to techniques used to determine the sequence of RNA molecules. It includes high-throughput shotgun sequencing of cDNA molecules obtained by reverse transcription from RNA, and next-generation sequencing technologies to sequence the RNA molecules within a biological sample in an effort to determine the primary sequence and relative abundance of each RNA molecule.
For this project, we primarily use gene expression data that has been generated using the following techniques:-• RNASeq: RNASeq uses next-generation sequencing (NGS) to reveal the presence and total quantity of RNA in a biological sample at a given moment.
• Ribosome Profiling (RPF): RPF uses specialised messenger RNA (mRNA) sequencing to determine which mRNAs are being actively translated. Unlike RNASeq, which sequences all of the mRNA of a given sequence present in a sample, RPF only targets mRNA sequences that are being actively translated.

Methodology
We posit that host factors regulated in viral translation may be important for establishing viral infection. The infection and analysis was carried out as follows -

Clustering Time Series Gene Expression Values
Clustering techniques are essential in the data mining process, as they reveal natural structures and identify interesting patterns in the underlying data. For gene expression data observed at different time-points, clustering can help identify biologically relevant groups of genes and samples.
Through the principle of guilt-by-association [23], if we find a known cluster that contains some uncharacteristic genes, we can hypothesise that these genes are significant.
As mentioned in [11], there is no one-size-fits-all solution to clustering, or even a consensus of what an ideal clustering should be like. Clusters may have arbitrary shapes and sizes, and each clustering criterion imposes a certain structure on the data. If the data happen to conform to the requirements of a particular criterion, the true clusters are recovered; however, very rarely is this the case. Each algorithm imposes its own set of biases on the clusters it determines. Most (reasonable) clustering algorithms would yield similar results on synthetically constructed datasets; however, in practice, they can give widely differing results on real-world,noisy, gene expression data.
For temporal gene expression data, the most commonly used clustering algorithms [12,16,24]

Self Organising Maps
For all our model based clustering techniques, we determine partitions by using the expectation-maximisation (EM) algorithm (described in 2.2.1) for maximum likelihood estimates of model parameters .
In this report, we will use hierarchical clustering for exploratory analysis, and apply all the model based clustering techniques to the different datasets we are provided with.

Expectation Maximisation (EM) Algorithm
The expectation maximisation (EM) algorithm is a general technique for finding maximum likelihood solutions for probabilistic models having latent variables [5]. EM is used for all the model based techniques covered in this report.

Procedure
Consider a probabilistic model in which we collectively denote all the observed variables by X and the hidden variables by Z. The joint distribution is governed by a set of parameters, denoted θ. We try to maximise the likelihood function, that is given by If Z is continuous, we replace the summation with integration as appropriate. The maximum likelihood estimator is Direct optimisation of p(X|θ) is difficult due to the presence of latent variables, so we attempt to optimise the expected value of the log-likelihood of the data given our model. We are now ready to define the general EM algorithm in algorithm 1.

Algorithm 1 Expectation Maximisation
Initialisation: For l = 0, make an initial estimate of the parameters θ. E Step (l th iteration): Evaluate p(Z|X,θ), i.e find the conditional probability distribution of the latent variables given our current model θ (l) and the data X. M Step (l th iteration): Update θ as where End: Terminate when one of the following criteria is met 1. Iterate till l reaches the maximum number of iterations specified 2. Till the parameter estimates stop changing 3. Till the likelihood function stop changing

Limitations
Convergence of the EM algorithm is guaranteed by the Majorisation-Minimisation technique [26]. However, the EM Algorithm typically converges to the local optima. Thus, it needs to be run multiple times in order to find the best fitting model and its parameters.

K-Means
K-Means is an algorithm to identify groups, or clusters, of data points in a multidimensional space [5].
We group data points based on their feature similarities according to a given distance metric (Euclidean Distance in our case), by minimising a given cost function.
Given a set of data points D = {x t } n t=1 and fixed number of clusters K ≤ n, the algorithm works as follows [26]-
Step l ∈ N, Assignment (E) Step: Assign each point x t to its closest mean Step l ∈ N, Update (M) Step: Recompute µ j 's as means of the assigned points: End: Terminate when cluster assignments j t do not change.
For a given partition this algorithm minimises the overall squared error: where D 1 , D 2 , . . . , D n are the clusters and µ 1 , µ 2 , . . . , µ n are cluster means. The two steps of updating the clusterings and re-estimating the centers correspond to the E step and M step of the EM Algorithm respectively.

Equation (2.2.2) is also known as the cost function.
The main problem of k-means is its dependency on the initially chosen centers. If the centers are chosen arbitrarily, then the algorithm converges very slowly or results in poor clustering. However, this can be solved by using the K-Means++ [4] algorithm, which selects centers that are far apart from each other in a probabilistic manner.
Through K-Means++, we choose the initial K centers as follows [26]-1. Take a random data point µ 1 chosen uniformly from dataset D.

Calculate new center
denotes the shortest distance from the data point x to its closest, already chosen center) 3. Repeat the previous step till we have K centers.

Limitations
Due to its ease of implementation and termination in a finite number of steps, K-Means is one of the most popular (and usually one of the first) algorithms applied for unsupervised learning problems.
However, due to the nature of the distance function chosen (Euclidean), K-means doesn't capture temporal dependencies in the data. Thus, it is not entirely suitable for our time series dataset. Moreover, k-means does hard clustering, i.e. a point can belong to only one cluster, while in real life data, a point can belong to more than one cluster.

Gaussian Mixture Models
Gaussian mixture models (GMMs) are mixtures, or linear combinations, of models with Gaussian distributions. They are widely used in data mining, statistical analysis and pattern recognition problems [5]. A Gaussian mixture model is represented by where the equation above represents a mixture of K Gaussians, N (x|µ k , Σ k ) is a component of the mixture model (each having its own mean µ k and covariance Σ k , and π k are the component weights, such that k π k = 1. To use GMMs for clustering, we assume that our data are generated i.i.d from a mixture of K Gaussians, and we apply EM to estimate the K means, covariances, and weights, and thus determine the clusterings. We achieve this by introducing a latent variable z in (5). This latent variable has a 1-of-K representation -thus z k ∈ {0, 1} and k z k = 1 (for each x), and find p(z) and p(x|z), to get Another quantity that is important is the conditional probability of z given x, also known as the responsibility that component k takes for producing the observation x, denoted by These formulas containing latent variable z help us formulate EM for GMM, given in 3.
The clusterings are indicated by the latent variable z.

Algorithm 3 EM for Gaussian Mixture Models
Init: Initialise the means and covariances µ k and Σ k and mixing co-efficients π k , and evaluate the initial value of the log-likelihood.
Step l, E step: Evaluate the responsibilities from (7) using the current parameter values.
Step l, M Step: Re-estimate the parameters using the current responsibilities.
and terminate if one of the following criteria are met-1. Paramater estimates have converged (stopped changing).

Log likelihood has converged.
3. Maximum umber of iterations l have been reached. otherwise return to Step 2.

Covariance Matrix
Having different settings of the covariance matrix for each component help in restricting the number of parameters that need to be estimated (thus optimising the runtime of the algorithm). However, some of the settings (diagonal, tied and spherical) restrict the contour of the clusters obtained, as seen in figure 5. In the diagonal setting, the elliptical clusters obtained are aligned with the co-ordinate axes, while the spherical clusters are circular, while all clusters have the same shape in the tied setting. Assuming we have k components, here are 4 different covariance settings, summarised in the table 1.

Limitations
As discussed for previous models, GMM converges to a local optima, so we need to run the model multiple times to determine the best fit. Moreover, GMM takes

Hidden Markov Models
In a first-order Markov chain describing a sequence of possible events, the probability of each event depends only on the state attained in the previous event.
A hidden Markov model can be viewed as a Markov chain of latent variables, with each observation conditioned on the state of the corresponding latent variable. The states of the Markov chain are not observed directly, i.e. the Markov chain  A single time slice of HMM corresponds to a mixture distribution, with component densities given by p(x|z). Therefore HMMs can also be interpreted as an extension of a mixture model in which the choice of mixture component for each observation is not selected independently but depends on the choice of component for the previous observation [5].
HMMs are widely used in speech recognition, natural language modelling, on-line handwriting recognition and for the analysis of biological sequences such as proteins and DNA [5]. Thus, we can use a Hidden Markov Model (HMM) to model temporal dependencies in the data.
An HMM contains 3 essential components -observed sequence of data x i , hidden states z i , transition probabilities p(z i |z i−1 ) and emission probabilities p(x i |z i ).
The hidden states in an HMM are analogous to latent variables in GMMs, responsible for generating the given observations. The probability distribution of a hidden state z n depends upon the state of the previous hidden state z n−1 . This is encoded by a transition probability matrix. For an HMM with M hidden states, the transition probability matrix has M (M − 1) independent parameters.
The emission probabilities are conditional distributions of the observed variables are p(x n |z n , φ), where φ is a set of parameters governing the distribution. The φ may be given by Gaussians if the elements of x are continuous variables, or by conditional probability tables if x is discrete [5,13].
There are 3 basic problems tackled by hidden Markov models [26,13] -1. Given an observation sequence x and a model λ, compute the probability p(x|λ). This is achieved using the Forward-Backward algorithm.
2. Given an observation sequence x, compute the most likely observed state sequence. This is done using the Viterbi algorithm.
3. Adjust the model parameters of λ to maximise p(x|λ) using the Baum-Welch algorithm.

HMM-based Clustering
We can apply the EM algorithm defined in Algorithm 1 to do clustering using HMMs.
The EM procedure for HMM based clustering is fairly complex, as we need to use a sub-EM procedure (Baum-Welch algorithm) in the M step, K times.
Starting from an initial collection of HMMs, EM iteratively finds cluster assignments that maximise the joint likelihood of the clustering.
Initialisation: Given sequences O i , where i ∈ [1,n], we randomly generate K HMMs λ 1 , λ 2 , ..., λ k and a partitioning C = C 1 , C 2 , ...C k , where the assignment to clusters is done to maximise the following objective function: where L(O i |λ k ) defines the likelihood function, i.e. the probability density for generating sequence O i by model λ k .
Assuming that we have a collection of K HMMs, our algorithm is as follows:

Limitations
The biggest limitation with HMM-based clustering is convergence to a local optima, so the algorithm needs to be run multiple times with different random initialisations to find the best fit. However, like GMM training, HMM based clustering takes significant computation time -especially due to the Baum-Welch EM procedure in the M step of the clustering.

Bayesian Networks
Bayesian networks are probabilistic graphical models that represent a set of variables and their conditional dependencies via directed acyclic graphs.
In [26], we studied the Chow-Liu algorithm to learn tree-structured graphical models. In [21], the authors combined Chow-Liu Trees and Expectation Maximisation to derive a model-based clustering method, called the CL Multinet Classifier.
If there are K classes, then K CL trees are learned. Each tree distribution P r k approximates the joint probability distribution of the attributes, given a specific class. The root attribute has no parents, while the other attributes can have at most ONE parent, i.e. if the attributes are expressed as r.vs {X 1 , X 2 , ..., X n }, then Π X i = {X j =i } has only one attribute, and Π Xroot = {∅}.
An inherent advantage of applying Bayesian networks case is that we can infer the structure behind each cluster too. For each cluster, we can see which timepoints affect the other.
Using Bayes Theorem, the posterior probability of a Chow Liu Tree can be expressed as P r(C = k|X 1 , X 2 , ..., X n ) = P r(C = k)P r(X 1 , X 2 , ..., X n |C = k) k P r(C = k )P r(X 1 , X 2 , ..., X n |C = k ) We observe that the denominator is constant with respect to the class, so we express it as 1 β . Then, the previous formula can be re-written as The edge weights of a Chow Liu tree are encoded by the mutual information, which is derived from the empirical probability distributions observed from the data.
I(X, Y ) = x∈X y∈Y P r(x, y) log P r(x, y) P r(X)P r(y) A maximum weighted spanning tree (either Prim's/Kruskal's) algorithm is then used to find the Chow-Liu tree for each k th class.
Unsupervised training of the CL Trees (henceforth referred to as CL Multinets) is carried out assuming that the data have been generated from a mixture of K Bayesian networks. This mixture can be described as where The α k are the mixing coefficients and the f k are the joint probability distributions of the CL Multinets. Given N observations, K clusters (CL Multinets), and initial partitions P 0 ={P 0 0 ,P 0 1 ,...,P 0 K }, the EM procedure for the m th iteration is detailed in algorithm 4.
C-Step: Update the partition P m ={P m 0 ,P m 1 ,...,P m K } by assigning each x r to the cluster that provides the max posterior probability. M-Step: For k = 1, ..., K, maximise the Classification Maximum Likelihood (CML) criteria and re-estimate the parameters θ m using the new partitions P m k .
where n k are the number of samples belonging to the k th cluster, and α k are the cluster weights.
The parameters for the second term are re-estimated (maximised) using the following formula: The first term (after some manipulation) contains the mutual information computed by the observed data (empirical distributions). Using the maximum weighted spanning tree algorithm, we maximise (15), resulting in a new tree distribution. Terminate: If any of the following criteria are reached: 1. Maximum number of iterations reached.

Limitations
Although Bayesian Networks are directed acyclic graphs and thus use fewer parameters than HMMs, unlike HMMs, we need to discretise the data in order to use the empirical frequencies to estimate joint probability distributions. Discretisation also results in loss of information, as the resultant dataset is dependent upon the number of bins chosen. Moreover, just like HMM-based clustering, there exist only general information theoretic criteria to determine the number of clusters K, which is not always be applicable [21].

Hierarchical Clustering
All the previous models that we covered can be classified under partitional clustering. The other category, hierarchical clustering, seeks to build a hierarchy of clusters using a given distance metric for datum and a similarity criteria for the clusters obtained. However, unlike partitional clustering, there exist no definitive statistical techniques determine the optimal number of clusters -instead, we need to examine the resulting hierarchical structure (dendogram) and visually determine the optimal number of clusters [18].
There are two approaches to hierarchical clustering -• Agglomerative: Each datum initially represents a cluster of its own. Clusters are then successively merged until the desired cluster structure is obtained.
• Divisive: All objects initially belong to one cluster. This main cluster is then divided into sub-clusters, which are successively sub-divided into their own sub-clusters, and this process continues until the desired cluster structure is obtained.
We set the distance metric (for within cluster distances) as Euclidean. The linkage criterion in hierarchical clustering specifies the difference between any two sets (clusters) obtained. It determines which clusters to merge by calculating the differences between two clusters as a function of the pairwise distances between observations. The different types of linkages are-• Single: the distance between two clusters is equal to the shortest distance from any member of one cluster to any member of the other cluster.
• Average: the distance between two clusters is equal to the average distance from any member of one cluster to any member of the other cluster.
• Complete: the distance between two clusters is equal to the longest distance from any member of one cluster to any member of the other cluster.
• Ward: the distance is defined as the error function of the unified cluster minus the error functions of the individual clusters, where the error function is the average distance of each data point in a cluster to the cluster centroid.
For this report, we use agglomerative hierarchical clustering not in model building, but in Section 4.1 for feature engineering , so as to obtain a rough estimate of the number clusters we may expect from the data.

Evaluation Metrics for Model Selection
Given data, model selection is the task of selecting the best statistical model from a set of candidate models. In our clustering analysis, the main conundrum is determining the number of clusters.
The estimation of the true number of classes has been recognised as "one of the most difficult problems in cluster analysis" [6]. The correct choice of K is often ambiguous -sometimes, choosing the best fitting model is not the wisest choice. Consider increasing K without any penalty -this will keep reducing the error but the model is less than ideal (e.g. In K-Means, if we have n clusters for n points, each point is its own center, reducing the cost function to 0). Moreover, given multiple candidate models having similar accuracy, the simplest model is most likely to be the best choice (Occam's razor).
Thus, the optimal choice of K balances accuracy and number of parameters required for the model. This is tackled in 2.

and 2.3.2.
A secondary conundrum is comparing two different clustering models, given that we obtain "ground truth" values for some of the clusters. This is handled in 2.3.3.

Information Criteria
When a statistical model is used to represent the process that generated the data, the representation will almost never be exact. Some information is thus lost by the model. Bayesian Information Criteria (BIC) and Akaike Information Criteria (AIC) [7] deal with estimating this lost information, and calculate the trade-off between the goodness of fit of the model and the complexity of the model.

BIC and AIC
BIC and AIC penalise the total number of parameters. Both are functions of the maximised log-likelihood of the model and the estimated number of parametershowever, BIC penalises model complexity more heavily. We define the AIC and BIC as follows- where L(θ) is the maximised log-likelihood of the model θ, P is the number of free parameters and N is the number of data points. For a range of K values, we calculate AIC/BIC, and choose the K for which this criterion is maximum. Thus, higher AIC or BIC values indicate better fitting models.

Elbow Method
The elbow method [17] is a popular method to determine the true number of clusters by plotting the cost function for a given clustering algorithm against different K, and observing the point where the graph makes a sharp elbow shape, as seen in Figure 9. More formally, this method looks at the percentage of variance explained as a function of the number of clusters. In essence, the elbow method observes the point of sharpest decrease in the cost function, i.e. the most improvement to the cost function, and estimates that to be the ideal number of clusters. However, the limitation of this method is that the elbow point cannot always be unambiguously identified -sometimes there is no elbow, or sometimes there are several elbows.

Misclassification Error (ME) Distance
ME distance is used to calculate the similarity between two different clusterings of the same data. Given the "ground truth" (i.e. true classes of each point), we can compare the clustering obtained by our technique with the ground truth to see how our model compares. ME distance thus represents the well known cost of classification, minimised over all permutations of the labels [19].
It is calculated as follows - Then the ME distance is defined as where Π k contains all permutations of the label K.
We want the ME Distance between our models and the ground truth to be close to 0.

Gene Ontology
Gene Ontology (GO) is a framework for the model of biology, unifying the representation of gene and gene product attributes across all species. We define concepts and classes used to describe gene functions, and observe the relationships between them.
The ontology covers 4 categories-1. Biological Process: operations or sets of molecular events with a defined beginning and end, pertinent to the functioning of integrated living units: cells, tissues, organs, and organisms.
2. Cellular Component: the parts of a cell or its extracellular environment.
3. Molecular Function: the elemental activities of a gene product at the molecular level.
4. KEGG Pathway: pathways representing molecular interaction, reaction and relation networks for various biological processes.
We use DAVID [15] to perform gene ontology analysis on the clusters we generate.
Each cluster obtained from our model(s) contains genes that share similar biological properties. We perform GO analysis on all these different clusters, and attempt to identify terms that are related to infection -like viral process, translation, DNA repair, etc. If we do find terms related to viral process, we use guilt-by-association to find other significant terms. Figure 10: Sample GO for a cluster p-value is the probability of seeing at least x number of genes out of the total n genes in the list annotated to a particular GO term, given the proportion of genes in the whole genome that are annotated to that GO term. In essence, the GO terms shared by the genes in the user's list are compared to the background distribution of the annotation. Thus, the smaller the p-value is, the more significant the particular GO term associated with the group of genes is (i.e. the less likely the observed annotation of the particular GO term to a group of genes occurs by chance) [15].
As described in Section 1.3, RNASeq represents the total amount of genes within the cell, while RPF represents the amount of genes being translated. TE measures the rate at which genes get translated. We first filter the dataset for cdReads ≥ 10 to remove noise; this removes genes with very small reads, likely generated due to experimental errors.

Data Preprocessing
Then, we transform the cdRPKM values to log 2 , which is a standard practice in bioinformatics.
We can calculate TE data in Figure 12 by dividing the raw RPF value by the raw RNASeq value, and transforming it to log 2 . Alternatively, we can subtract the log 2 RNASeq value from the log 2 RPF value to also obtain TE.

Initial Analysis
For the RNASeq and RPF datasets, we plot individual timepoint log 2 values in Figures 13 and 14. From the plots, we infer that all timepoints follow a Gaussian distribution. we posit that Gaussian Mixture Models and Hidden Markov Models with continuous emissions will be useful in clustering our time-series gene expression values.
For both datasets, we observe that the median keeps shifting to the left for later timepoints. This is due to the infection spreading-with the passage of time, more viral translation leads to an overall decrease in gene expression.
To standardise the medians for all timepoints, we again follow a standard practice in bioinformatics -we median-normalise each timepoint (subtract the median from every value) so that each timepoint has the same median (0).
After median-normalisation, the distribution of the maximum and minimum values are generally similar for all time points.

Results
We apply the models from Section 2.1 to the datasets described in Section 3, and attempt to identify those groups of genes that are critical for viral translation.
There are two tasks we need to perform -firstly, identify an optimal K, and secondly, decide the benchmark algorithm (the most optimal model) that we will compare all our other techniques with.
We evaluate the performance of the models through the methods described in Section 2.3 and interpret the significance of the results obtained.
We perform our analysis in a systematic way -1. A priori, we know that the datasets will not contain more than 20-30 clusters.
We restrict our analysis for 2 ≤ K ≤ 30.
2. We apply model-based clustering techniques on the complete datasets, and use the evaluation criteria to determine the optimal K.
3. If the results are inconclusive, we implement feature engineering in Section 4.1 and re-run the model-based clustering techniques on this modified data.
4. Attempt to determine K through BIC/AIC.

5.
If no optimal value is found, do clustering for specific K (as determined by biologists, as in Figure 15).
6. Interpret results through GO analysis and/or other analyses.

Feature Engineering
Since this is time-series data, we also quantify the change between successive timepoints.
From the RNASeq and RPF datasets, we find the amount of change between timepoints and add these as features to the original dataset. Formally known as log 2 fold change in bioinformatics, it is calculated as follows -• Calculate the log 2 of the gene expression value, cdRPKM, for each timepoint.
Let this be x i .
• Subtract x 0 (the value at timepoint 0) from each successive timepoint x i .
• Add this as a new attribute to the RNASeq and RPF data.
For TE data, since the values are already log ratios, we merely subtract each timepoint's value from the first timepoint.
Even with the new features, we might encounter problems due to noise. Fortunately, we have a biological approach of dealing with this.
When a cell is infected with a virus, the RNASeq does not change much (as it measures the total content of all genes) while those genes that help in virus production will get translated more, i.e. the RPF of these genes will increase. Since TE measures translational efficiency (T E = RP F RN ASeq ), the TE of such genes will also keep increasing across timepoints.
We filter out those genes that do not change much at the RNASeq level. We retain genes which have foldChange values for each timepoint between −0.5 and 0.5. More formally, if x i is the foldChange for 1 ≤ i ≤ 5, then −0.5 ≤ x i ≤ 0.5 for each x i (later on, we reduce this threshold further to [-0.3,0.3] to check whether the GO terms we obtained are indeed important or not).
After this filtering, we are left with 5576 genes that are highly translated.
As we are interested in genes that are highly translated between timepoints, we apply hierarchical clustering only on the foldChange features in the TE dataset. We plot heatmaps for different linkages, and check the general trends to find the number of clusters. In a heatmap, the values of each gene are plotted according to the colour key (bottom left).
The values in the heatmaps are sorted accorded to how they have been clustered.
In our case, we find that complete linkage performs the best. We can identify 4-6 distinct clusters from the heatmap, and cut the dendogram accordingly (represented by the green line). The result is shown in Figure 15. After applying hierarchical clustering and doing heatmap analysis on the filtered data, we infer that K lies between 4 and 6. We use these values of K for clustering for Step 5 of Procedure 4.

K-Means
We first apply K-Means on all the datasets to determine the optimal K.
We initialise the algorithm using K-Means++ for optimal clusterings, and apply the Elbow method and BIC to evaluate our results. Results for RNASeq are given in Figure 16.   From the plots, we observe that determining the optimal K is difficult. There is no clear 'elbow' point, and the BIC keeps increasing with K; thus K-Means is not penalising enough. (Using AIC is futile, as the penalisation would be even lower in this case). The elbow plots and GO analysis do not mutually agree on a K.
Despite the mixed results from the information criteria metrics, we went ahead and performed GO analysis as per the K determined in Table 2. However, the GO results were also found to be inconclusive -there were multiple clusters associated with the 'viral process' term, and the p-value was quite high (i.e. the terms were not enriched).
Thus, there is no reliable estimate for K. This is most likely due to 2 reasons -• K-means cannot capture temporal dependencies in data.
• Data is too noisy.
Both can be solved through feature engineering -we add the changes between the timepoints as features.
We filter the data as described in Section 3, and using the heatmap results, perform K-Means clustering for 4,5 and 6 clusters on the foldChange TE data.
We plot heatmaps to observe the general trends of clusters from this heuristic method, and identify those clusters which have an increasing trend (the clusters are plotted top to bottom in the heatmap and demarcated by black lines).  This heuristic method gives us GO results (specifically, terms related to GTPases from Section 1.2.4.1) that are relevant to our hypothesis. This is covered in greater detail in Section 4.4.

Gaussian Mixture Models
As seen in the K-Means approach, it is futile to do GO analysis on complete datasets (especially for the noisy RNASeq data, as we will most likely observe no significant trends by merely tracking overall changes in gene expression).
To test this, we still apply GMMs on the complete TE data (which are most likely to show certain trends due to genes being translated). As expected, we have no definitive results. The GO is unclear, and the BIC keeps increasing with K for all covariance settings, offering no clear peak. The plots are in Appendix A, Section A.2.
We decide to apply GMMs on the feature-engineered data from Section 4.1. We use the full covariance matrix and do BIC analysis for 2 ≤ K ≤ 20.  We see in Figure 19 that Cluster 1 has an increasing trend. The GO analysis of this cluster is similar to that obtained through K-Means -we again observe terms related to GTPases.
We still perform GMM-based clustering for K = 4, 5, 6 for evaluating GMMs performance against the benchmark algorithm in Section 4.3.

HMM-based Clustering
We speculated that HMM-based clustering might perform better than other algorithms on the complete and unfiltered datasets, as HMMs inherently take the time dependencies into account.
We assume 3 hidden states for simplicity, and apply HMM-based clustering on all the datasets. The plots are in Section A.4.1, while the results are summarised in Table 3.

Data
BIC RNASeq (5 points) 13* RPF (5 points) 8* TE (5 points) 16* We can deduce from BIC's inability to determine an optimal K that we may be estimating too many parameters and unnecessarily complicating the models (Occam's Razor).  We proceed to the next steps of the analysis. We apply HMM-based clustering on the feature-engineered data for K = 4, 5 (we skip 6 as HMMs take the longest time to train) and generate heatmaps in Figure 20. We then check the heatmaps for up-regulated clusters, and perform GO analysis on the Clusters 2 and 1 respectively from Figure 20.
The trend graphs of each cluster are given in Section A.4.2.
The GO results from HMMs are comparable to those from GMMs and K-Means, and we again observe GTPases in the gene ontology of the up-regulated clusters.

Bayesian Networks
A limitation of Bayesian networks is that the current implementation is only applicable to discrete random variables (or discrete features).
Although there is some literature that deals with hybrid Bayesian Networks that model data containing both discrete and continuous attributes [8], data with purely continuous variables are difficult to model using Bayesian Networks.
Thus, as a first step, we discretise our gene expression data [21].
According to [21], there is no conclusive way of determining K for CL Multinets based clustering, so we skip AIC/BIC analysis here. Instead, we perform clustering using K = 4, 5, 6 to determine the structure of the underlying Bayesian network for the complete datasets.
The heatmap results for the TE dataset showed some increasing trends (covered in Section 29; the GO analysis for the complete TE dataset included GTPases in more abundance than HMM-based clustering. Moreover, the underlying structure for each tree (determined by the clustering, for different K) unfailingly resulted in the same line graph (in Figure 21), with the timepoints in order (timepoint 1, followed by timepoint 2, and so on). This, along with the GO results, shows that Bayesian networks cluster even noisy data with a certain level of accuracy. The heatmaps for 6 clusters are given in Figure 22. The blue lines demarcate the clusters.  Although the heatmaps are not as uniform and do not display concrete trends (as those from K-Means clustering) the GO results are reliable and contain more GTPase-related terms, along with some new terms related to viral translation (e.g. G2/M transition of mitotic cell [10]).

Performance Evaluation
We compare algorithms against a benchmark for the feature engineered TE dataset (with the foldChange features).
Bayesian networks consider the inherent time-dependency between attributes, resulting in a graph with the timepoints in correct sequence each time.
The GO results are also significant as they capture those genes that are upregulated yet also distinct from each other despite having similar expression values (unlike K-Means, which clusters according to Euclidean distance between absolute attribute values).
We pick Bayesian networks as our benchmark algorithm, and treat the results obtained from it as the ground truth. We use the ME distance metric to compare the other techniques against Bayesian networks for K = 4, 5, 6. For ME Distance, the closer it is to 0, the better the clustering is. We calculate the accuracy as 100(1-d M E )%.
*compared to Bayesian networks

Biological Results
Existing literature on picornavirus (which the EV-A71 virus is a type of) suggests that the proteins in the virus interact with certain GTPases and disrupt normal cellular functions (like nuclear transport of other viral and nuclear proteins). [22].
Looking more closely at the collection of heatmaps obtained by all the algorithms, we identify all those clusters that are up-regulated and contain terms containing/related to GTPases in their GO analysis. To be more rigorous about the terms we observe, we also reduce the filter to From the GO analysis, we find plenty of evidence from our datasets to support our hypothesis that a certain class/member of GTPase is important for establishing this infection.

Conclusion
This thesis explored the application of unsupervised learning techniques on viral genomic data -more specifically, human cells infected with the EV-A71 virus (which causes HFMD). We attempted to determine host factors important to the propagation of this virus, and identified several key biological terms based on our clustering results.

Summary
As described in Section 1.3.1, we collected RNASeq and RPF gene expression data at different timepoints. In Figure 2, we saw that the space occupied by the virus increases with time. Initial analysis of the complete datasets (RNASeq, RPF, TE) were consistent with this trend. Next, we attempted to determine groups of genes crucial to viral translation. Initial clustering results contained some relevant GO terms (e.g. viral process, DNA repair) that appeared in multiple clusters; however, the trend graphs and heatmaps proved to be inconclusive, suggesting that our datasets were too noisy, containing many points insignificant to our analysis.
We then focused on those genes which were highly translated at the RPF level but not at the RNASeq level. Collating data from all our GO analyses, we identified the following terms -Positive regulation of GTPase activity, Activation of GTPase activity, Regulation of small GTPase mediated signal transduction, Positive regulation of transcription, DNA templated and Signal transduction.
These terms imply that some proteins in the class GTPase are definitely important to viral translation, possibly increasing the production of viral genes at key phases.
Following up on these results, we conduct more extensive literature reviews to verify the link between small GTpases and picornavirus infection.

Biological Experiments
We also design experiments to verify whether GTPases are relevant to the EV-A71 virus.
The first experiment tracks changes in GTPase activity throughout infection, focusing more on Rho GTPases. For GTPase protein-specific activity, there are also targeted studies to check for specificity of GAP/GEF, which are related to activation of GTPases.
The second experiment inhibits specific GTPases to check their effect on viral cells. 3 kinds of drugs targeting different GTPases are used -• Specific Rho drugs e.g. Rhosin or C3 transferase against RhoA/B/C, but not cdc42/Rac1 • Specific cdc42 drugs e.g. Casin, ML141 • Specific Rho kinase inhibitors e.g. Rhodblock series Thus, the appearance of terms related to viral process in the GO analysis for all methods suggests that our clustering was useful. The abundance of GTPase-related terms also support our hypothesis that they are indeed vital to the propagation of the virus inside the cell.
6 Future Work 6.1 Algorithms

Tuning
We can improve the performance of our current model-based clustering techniques by incorporating additional data or adjusting the hyperparameters for the prior distributions we use for initialisation. We can consider discretising the data for Hidden Markov Model based clustering [28]; conversely, we can consider continuous variables for Bayesian networks [14].

Different Algorithms
We can consider self-organizing maps (SOMs) for clustering gene expression data too [25]. SOMs impose partial structure on the clusters (in contrast to the rigid structure of hierarchical clustering, the strong prior hypotheses used in Bayesian clustering, and the non-structure of k-means clustering) and facilitate easy visualization and interpretation. SOMs have good computational properties and are easy to implement, reasonably fast, and scalable to large data sets too.
We can also consider soft clustering, or fuzzy clustering; where one gene may belong to more than one cluster. Existing fuzzy-clustering techniques can indeed be modified for short time-series gene expression data [20].