From this excellent paper on biological variability in RNA-Seq experiments (bold highlights are mine):
“Biological variability has important implications for the design, analysis and interpretation of RNA-sequencing experiments. […] If only a few biological replicates are available, it will be impossible to estimate the level of biological variability in expression for each gene in a study. Supplementary Table 1 summarizes a large number of published RNA-sequencing studies over the past three years. In every case, except for the two studies we analyzed here, conclusions were based on a small number (n ≤ 2) of biological replicates. One goal of RNA-sequencing studies may be simply to identify and catalog expression of new or alternative transcripts. However, all of these studies make broader biological statements on the basis of a very small set of biological replicates.
Our analysis has two important implications for studies performed with a small number of biological replicates. First, significant results in these studies may be due to biological variation and may not be reproducible; and second, it is impossible to know whether expression patterns are specific to the individuals in the study or are a characteristic of the study populations. These ideas are now widely accepted for DNA microarray experiments, where a large number of biological replicates are now required to justify scientific conclusions. Our analysis suggests that as biological variability is a fundamental characteristic of gene expression, sequencing experiments should be subject to similar requirements.”
If you are doing RNA-Seq, be very vigilant in your experimental design and find a way to incorporate more replicates, even at the expense of testing fewer comparisons. It’s better to test one comparison (tissue X vs. Y, for example) with 5 or more replicates than to test three comparisons (Tissue X vs. Y, Y vs. Z, and X vx Z) with only 2 replicates for each tissue type.
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.
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.”
Browse to one of the output CSV files and select it. It will open in a new window like in Fig 3 below:
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:
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:
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.
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.
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.
A recent collaboration with a clinician here at UI hospital and clinics introduced me to the idea of mutational signatures in cancer. Characterizing mutational signatures is made possible by the falling cost and increasing accuracy of whole-genome sequencing methods. Tumors are sequenced across the entire genome and the catalog of somatic mutations (i.e, SNPs) is used to compute the mutational signatures of a tumor’s genome.
The idea is that the collection of somatic mutations found in a tumor are the result of a variety of defective DNA-repair or DNA-replication machinery combined with the action of known or unknown mutagens and environmental exposures. The processes operate over time and leave a “footprint” in the tumor DNA that can be examined. These sum of all of the mutational processes operating within a tumor cell is a distinct mutational “signature” that differs by tumor types.
For example, in lung cancer, the bulk of somatic mutations are C>A transversions resulting from chronic exposure to tobacco smoke. In melanoma, the predominant mutation type is C>T and CC>TT at dipyrimidines, a mutation type associated with UV-light exposure. And in colorectal cancer, defective DNA mismatch repair contributes the majority of the mutations.
A recent paper in Nature has formalized this notion of mutational signatures in tumors and provided a mathematical framework (written in MatLab) for assessing how many and which signatures are operational within an uncharacterized tumor type (generally there between 2 and 6 processes).
In the paper, the authors analyzed almost 5 million somatic cancer SNPs and identified 21 unique signatures of mutational processes through a mathematical process of deconvolution, followed by experimental validation. A curated catalog of the most current signatures based on available sequence data can be found at the COSMIC database.
In part 2 of this post, I’ll go into more detail on the mutational signatures and link to some python code I’ve written to help get flat-file lists of SNPs into the correct form for easy input into the MatLab framework.