In this series of posts, I want to describe some work I’ve been doing attempting to apply an artificial neural network model to the problem of classifying copy number variation in a clinically-important gene, CYP2D6. I will be presenting a poster-paper on this topic at the 2018 ISMB/ISCB meeting in Chicago.
Here is the abstract from the poster:
Pharmacogenomics is a rapidly developing field that aims to deliver on the promise of personalized medicine by guiding pharmacological intervention using an understanding of a patient’s individual genotype for drug metabolism. By avoiding ineffective or dangerous treatments, patient outcomes could be dramatically improved and hospital costs reduced. We have developed a sequencing-based pharmacogenomics screening panel that uses targeted capture to perform deep sequencing of 200+ critical drug metabolism genes. Assessing copy-number variation is a critical part of correctly interpreting genotypes in key drug metabolism genes, such as CYP2D6. Historically, this has been a time-consuming step in our clinical pipeline involving large amounts of expert analyst time and data visualization. This study presents a novel application of an artificial neural network (ANN) machine learning algorithm to learn the complex patterns in CNV data. The result is a trained network that can quickly and accurately classify copy number events according to known training categories in the CYP2D6 gene. We show that a simple, one hidden-layer network is sufficient to achieve the extremely high accuracy and low false-positive rate required in a high-throughput clinical setting.
Motivating this work is the fact that interpreting copy number variation data from targeted capture sequencing is difficult owing to several factors. First, the data are noisy owing to biases in capture efficiency and GC content. Second, the copy number variation events in a gene like CYP2D6 are complex and subtle, but have dramatic impacts on the functional status of the gene. Third, most copy number variation detection methods are optimized for whole genome sequencing with smooth and even coverage across each chromosome.
The result is that, although most of the variant calling and interpretation is automated, a person still has to sit down with the copy number variation copy-ratio plots and make a manual determination of genotype for the CYP2D6 gene (and other genes if they contain clinically-actionable CNVs). Below is one such copy-ratio plot that illustrates the problems we face:
You can see from the plot that one sample (red) has what appears to be a deletion (copy ratio ~ 0.5). We do look at these plots individually, but the problem with noise and complexity remains. This means that the analyst must be, in effect, a domain expert on each gene with a clinically-actionable copy number variation. This limits pipeline throughput and is impractical if the goal is to process hundreds or thousands of patient samples per month.
A role for machine learning?
It occurred to me that the copy ratio data fall into distinct patterns that a human (like myself) learns by eye and with experience. These patterns could be classifiable by machine-learning methods. I considered applying a convolution neural network (CNN) to the plot above.
However, I only had 175 samples in my test/train set drawn from real patients and CEPH/Coriell depositories. Therefore, I thought it would be best to start with simple models and go from there. To establish the ground truth CYP2D6 genotype, the samples had all be processed by our current pipeline methods, with many (but not all) confirmed by Taqman assay.
It began to dawn on me that I first needed a better representation of the training data. I liked the output of the CNV-kit method (below) better than the method from our pipeline (above):
The important part of the plot above are the gray dots, representing copy ratio values across 19 different bins or segments along the CYP2D6/2D7 gene and pseudogene. CNV-kit also had the advantage of being well-documented, fast to run, and available as open-source software on github.
Tidying the training set
The data plotted above looked like this after I did some wrangling and tidying in python:
With the data in “tidy” form (each column a variable, each row an observation, thank you Hadley Wickham), I was ready to train some machine learning models to see if I could classify common copy number variation events in CYP2D6 and bypass the time-consuming visualization step.
I had no idea if this would work, given that there are only 175 samples (total) to train on and some of the more rare copy number variation events only have a handful of examples. In part 2, I will talk more about how I tried a LASSO regression model, which, while performing well, failed to yield the high accuracy needed for an automated clinical pipeline. I then tried a simple one-layer neural network approach. I’ll talk about the surprisingly good performance of neural network approach and also future directions for this project.
I recently had to inspect some genomic alignments as part of a project. Usually, I am just working with BAM files and if inspection is needed, I just visualize the pileups to see what is going on.
In this case, I just wanted a quick answer to how the reads were aligning to the reference, and I didn’t want to go through the process of subsetting and copying the BAM files to my local machine.
The SAM file is the uncompressed record of the read alignments produced by an aligner method (STAR, TopHat, BWA, etc….). This file can get very large, and so is usually compressed into BAM (faster for machine parsing, but not human readable) and the SAM file is discarded.
In my case, I still had the SAM files around to inspect. If you find yourself needing to read a SAM file, here are three helpful reference tools to make the process less painful:
1) This page has an enormous amount of detail about SAM files including this helpful chart that enumerates all of the fields that you can expect to find specified within each alignment:
2) This post from the blog “zenfractal.com” contains a great exposition on CIGAR strings and how to decode them:
3) And finally, if you’re trying to decode the SAM bitwise flags, you can calculate them using this tool from the Broad Institute: