Mutation Archives - Michael's Bioinformatics Blog

January 27, 2016

Hands-on with cancer mutational signatures, part 2

In this second part of the “Hands On” series, I want to address how to create the input for the MatLab mutational signature framework from the output of my python code to prepare the SNP data for analysis.

First, creating a Matlab .mat file for input to the program. The code is expecting an input file that contains a set of mutational catalogues and metadata information about the cancer type and the mutational types and subtypes represented in the data.

Fig 1. The required data types within one .mat file to run the framework.

As you can see from Fig 1., you need to provide a 96 by X matrix, where X is the number of samples in your mutational catalogue. You also need an X by 1 cell array specifying sample names, a 96 by 1 cell array specifying the subtypes (ACA, ACC, ACG, etc…) and a 96 by 1 cell array specifying the types (C>A, C>A, C>A, etc…). These must correspond to the same order as specified in the “originalGenomes” matrix or the results won’t make sense.

My code outputs .csv files for all of these needed inputs. For example, when you run my python code on your SNP list, you will get a “subtypes.csv”, “types.csv”, “samples.csv”, and “samples_by_counts.csv” matrix (i.e., originalGenomes) corresponding to the above cell arrays in Fig 1.

Now, the challenge is to get those CSV files into MatLab. You should have downloaded and installed MatLab on your PC. Open MatLab and select “Import Data.”

Fig 2. Select the “Import Data” button.

Browse to one of the output CSV files and select it. It will open in a new window like in Fig 3 below:

Fig 3. The data import window from MatLab.

Be sure to select the correct data type in the “imported data” section. Also, select only row 1 for import (row 2 is empty). Once you’re finished, click Import Selection. It will create a 1×96 cell called “types.” It looks like this:

Fig 4. The new imported cell data "types." — Fig 4. The new imported cell data “types.”

We’re almost done, but we have to switch the cell to be 96×1 rather than 1×96. To do this, just double-click it and select “transpose” in the variable editor. Now you should repeat this process for the other CSV input files, being sure to select “matrix” as the data type for the “samples_by_counts” file. Pay special attention to make sure the dimensions and data types are correct.

Once you have everything in place you should be ready do run the mutational analysis framework from the paper. To do this, open the “example2.m” Matlab script included with the download. In the “Define parameters” section, change the file paths to the .mat file you just created:

Fig 5. Define your parameters for the signature analysis.

Here you can see in Fig 5, I’ve specified 100 iterations per core, a number of possible signatures from 1-8, and the relative path to the input and output files. The authors say that ~1000 iterations is necessary for accuracy, but I’ve found little difference in the predictions between 50-500 iterations. I would perform as many iterations as possible, given constraints from time and computational power.

Note also that you may need to change the part of the code that deals with setting up a parallel computing pool. Since MatLab 2014, I believe, the “matlabpool” processing command has been deprecated. Substituting the “parpool” command seems to work fine for me (Mac OS X 10.10, 8 core Macbook Pro Retina) as follows:

if ( parpool('local') == 0 )

parpool open; % opens the default matlabpool, if it is not already opened

end

This post is getting lengthy, so I will stop here and post one more part later about how to compare the signatures you calculate with the COSMIC database using the cosine similarity metric.

November 12, 2015January 27, 2016

Hands-on with cancer mutational signatures

In my last post, I wrote about the biological context of mutational signatures in cancer and how a recent series of papers has addressed this notion by creating an algorithm for extracting signatures from a catalogue of tumor SNPs.

In this follow-up post, I wanted to offer practical advice, based on my experience, about how to prepare data for mutational signature analysis and how to interpret the results.

First, in order to analyze your SNPs for mutational signatures, one needs to derive the trimer context of each SNP (i.e., the upstream base and downstream base flanking the SNP) for technical reasons described in the paper linked above.

In order to do this, a reference genome must be queried with the positions of the SNPs in order to find those bases adjacent to them. One could either query a remote database, like Entrez Nucleotide, or download an entire 40GB reference genome and search it locally.

In my code, I opted for the former option: querying the NCBI/Entrez Nucleotide database using HTTP requests. The advantage of this approach is that I can reuse the same code to query multiple reference genomes (e.g., hg38 vs. hg19), depending on the changing needs of future projects.

The relevant section of code is as follows:

def get_record(chrom, pos):               
    '''Fetch record, and minus/plus one base from Entrez'''
    try:
        handle = Entrez.efetch(db="nucleotide", 
                    id=hg19_chrom[chrom], 
                    rettype="fasta", 
                    strand=1, 
                    seq_start=int(pos) - 1, 
                    seq_stop=int(pos) + 1)
                       
        record = SeqIO.read(handle, "fasta")
        handle.close()
    except BaseException:
        return None
    
    return record

You can see from the code that I am using a dictionary (‘hg19_chrom’) to translate between chromosome numbers and their Entrez Nucleotide IDs in the eFetch request.

The disadvantage of this approach is that Entrez HTTP tools limits the user to 3 queries per second (in fact this limitation is hard-coded into Biopython). Even with my mediocre coding skills, this turns out to be the rate-limiting step. Thus, this code would have to be run pretty much overnight for any sizable number of SNPs (~50,000 SNPs would take ~4.6 hrs). However, it’s easy to let the script run overnight, so this wasn’t a deal breaker for me.

In the ~~third and final post~~ next two posts on this topic I will address how to create the MatLab .mat file from the output of this script and how to compare the signatures generated by the MatLab framework to the COSMIC reference database in a non-biased way.

November 20, 2013

From inhibitors to activators: rethinking drug action

Most small molecule drugs are designed to inhibit their target protein from carrying out its cellular function. Drug discovery typically focuses on disrupting biochemical systems in cells in order to induce apoptosis (cell death) or reduce the activity of an overactive pathway.

Interestingly, it is becoming apparent that there may be a novel way of looking at the problem. Instead of trying to muck up a cell’s function, new therapeutic approaches may seek to enhance the functioning of healthy biochemical pathways or systems that are under-activated owing to genetic mutation.

One example comes from the search for effective small molecule drugs against Parkinson’s disease (PD). Several studies have recently shown that malfunctioning lysosomes are involved in the progression of PD by failing to clear waste and allowing the accumulation of misfolded proteins. The lysosomal function is reduced in PD patients owing to mutations in a protein critical for proper functioning called PARK9.

Researchers are now looking for molecules that can stimulate the lysosomal autophagy pathway by interacting with PARK9 or other proteins. By increasing the abnormally lowered activity of the pathway, it is hoped that increased clearance of PD-related plaques may be achieved.

Interestingly, a compound from traditional Chinese medicine (TCM) has been found to be activating towards autophagic activity and is now in development in the biotech industry.

More information can be found here:

http://www.alzforum.org/new/detail.asp?id=3172