iGC—an integrated analysis package of gene expression and copy number alteration

Background With the advancement in high-throughput technologies, researchers can simultaneously investigate gene expression and copy number alteration (CNA) data from individual patients at a lower cost. Traditional analysis methods analyze each type of data individually and integrate their results using Venn diagrams. Challenges arise, however, when the results are irreproducible and inconsistent across multiple platforms. To address these issues, one possible approach is to concurrently analyze both gene expression profiling and CNAs in the same individual. Results We have developed an open-source R/Bioconductor package (iGC). Multiple input formats are supported and users can define their own criteria for identifying differentially expressed genes driven by CNAs. The analysis of two real microarray datasets demonstrated that the CNA-driven genes identified by the iGC package showed significantly higher Pearson correlation coefficients with their gene expression levels and copy numbers than those genes located in a genomic region with CNA. Compared with the Venn diagram approach, the iGC package showed better performance. Conclusion The iGC package is effective and useful for identifying CNA-driven genes. By simultaneously considering both comparative genomic and transcriptomic data, it can provide better understanding of biological and medical questions. The iGC package’s source code and manual are freely available at https://www.bioconductor.org/packages/release/bioc/html/iGC.html. Electronic supplementary material The online version of this article (doi:10.1186/s12859-016-1438-2) contains supplementary material, which is available to authorized users.


Index 12
create_gene_cna Load and map CNA gain/loss onto human gene location by genome reference

Description
The function reads through in all sample CNA data given by the sample description sample_desc and returns a joint CNA gain/loss table based on gene regions across samples.
gain_threshold CNA expression above this will be considered as gain region. By default log 2 2.5− 1 loss_threshold CNA expression below this will be considered as loss region. By default log 2 1.5− 1 read_fun Custom reader function, see its own section for more detail.
progress Whether to display a progress bar. By default TRUE.
progress_width The text width of the shown progress bar. By default is 48 chars wide.
parallel Enable parallelism by plyr. One has to specify a parallel engine beforehand. See example for more information.
... Arguments passed to the custom reader function specified in read_fun.

Details
A gene is considered to have CNA gain if the overlapped CNA record expression is higher than the given threshold. Similarly, a gene is considered CNA loss if the overlapped CNA record is lower than the given threshold. If multiple CNA records map onto the same gene region with both gain and loss, the majority wins. If none of the records map to the gene, NA is given.
By default it assumes the data to be of TCGA level 3 file format. For other data formats (e.g. raw data or other experiments from GEO), one should implement a custom reader function that accepts the filepath as the first argument. See section Custom reader function for full specification.
Currently the package ships a custom genome reference hg19, hg19DBNM, for gene region look up. Each gene's region is defined by the widest splicing form it has in NCBI curated records. The defined region includes intron regions. This limitation may change in the future.
Value data.table of CNA gain/loss on each gene region for all samples, whose rows represent regions of genes and columns represent sample names. First column GENE contains the corresponding gene names.

Custom reader function
Custom reader function is given by read_fun = your_reader_fun. It takes the filepath to CNA data as the first argument and returns a data.table with at least the following four columns: Chromosome, Start, End, and Segment_Mean of type character, integer, integer and numeric respectively.
Rest arguments of create_gene_cna(...) will be passed to this reader function.
Note: all string-like columns should NOT be of type factor. Remember to set stringsAsFactors = FALSE.

See Also
read.table and fread for custom reader function implementation; create_sample_desc for creating sample description. If the gene information already exists in the data, try direct_gene_cna to skip the genome reference lookup.
progress Whether to display a progress bar. By default TRUE.
progress_width The text width of the shown progress bar. By default is 48 chars wide.
... Arguments passed to the custom reader function specified in read_fun.

Details
By default it assumes the data to be of TCGA level 3 file format. However, nearly all real world data fail to have the same format as TCGA. In this case, one needs to tell the function how to parse the data by implementing a custom reader function that accepts the filepath as the first argument. See Detail section for full specification. The function naively concatenates all return expression as if all gene expressions are stated in the same gene order as columns in a new data.table.
Value data.table of all samples gene expression, whose rows are gene expression and columns are sample names. First column GENE contains the corresponding gene names.

Custom reader function
Custom reader function is given by read_fun = your_reader_fun. It takes the filepath as the first argument and return a data.table with the first two columns being GENE and Expression of type character and double.
The output joint gene expression table has first column GENE store the gene name, which are are determined by the first sample being evaluated.
Rest arguments of create_gene_exp(...) will be passed to this reader function.
Note: all string-like columns should NOT be of type factor. Remember to set stringsAsFactors = FALSE.

Note
The function assumes row order for all samples' gene expressions are the same.

See Also
read.table and fread for custom reader function implementation; create_sample_desc for creating sample description.
sample_names character vector of distinct sample names. Samples will be referenced by the given name through out the analysis process. They should be valid R data.table column names.
cna_filepaths character vector of filepaths to CNA data.
ge_filepaths character vector of filepaths to gene expression data.
sample_root path to the root of sample data. If given, this path will be appended before all given filepaths.
Value data.table of sample description having the following columns in order: Sample, CNA_filepath, and GE_filepath. Each row contains a sample's unique name and the corresponding filepaths to CNA and gene expression data.

Note
One could convert the relative file paths into absolute paths by passing the root folder path to sample_root.
If for some special reasons, for example gene expression of all samples have been collected or the CNA records for each gene exist, but do not have the file paths to either CNA or gene expression data, pass it with empty character vector of correct length, such as rep( , num_samples).

## End(Not run) direct_gene_cna
Load the existed CNA gain/loss based on gene location.

Description
This function aims to complement create_gene_cna. Instead of mapping CNA records onto genes by genome reference, it reads the existed column containing the gene each CNA lies on. Two functions share the same interface but they have different requirement for the read_fun implementation.
gain_threshold CNA expression above this will be considered as gain region. By default log 2 2.5− 1 loss_threshold CNA expression below this will be considered as loss region. By default log 2 1.5− 1 read_fun Custom reader function, see its own section for more detail.
progress Whether to display a progress bar. By default TRUE.
progress_width The text width of the shown progress bar. By default is 48 chars wide.
parallel Enable parallelism by plyr. One has to specify a parallel engine beforehand. See example for more information.
... Arguments passed to the custom reader function specified in read_fun.
Value data.table of CNA gain/loss on each gene region for all samples, whose rows represent regions of genes and columns are sample names. First column GENE contains the corresponding gene names.

Custom reader function
Similar to that of create_gene_cna, the reader function takes the filepath as the first argument. It will return a data.table with at least two columns: GENE and Segment_Mean of type character and numeric respectively.

Description
The function finds CNA-driven differentially expressed gene and returns the corresponding p-value, false discovery rate, and associated statistics. The result includes three tables which collects information for gain-, loss-, and both-driven genes.  progress Whether to display a progress bar. By default TRUE.
progress_width The text width of the shown progress bar. By default is 48 chars wide.
parallel Enable parallelism by plyr. One has to specify a parallel engine beforehand. See example for more information.

Details
The gene is considered CNA-gain if the proportion of the sample exhibiting gain exceeds the threshold gain_prop, that is, number of samples having gain_loss = 1. Reversely, the gene is considered CNA-loss if %samples that gain_loss = -1 is below a given threshold loss_prop.
When performing the t-test, sample grouping depends on the analysis scenario being either CNAgain or CNA-loss driven. In CNA-gain driven scenario, two groups, CNA-gain and the other samples, are made. In CNA-loss driven scenario, group CNA-loss and the others are made. Genes that appear in both scenarios will be collected into a third table and excluded from their original tables.
See the vignette for usage of this function by a thorough example.

Value
List of three data.table objects for CNA-driven scenarios: gain, loss, and both, which can be accessed by names: 'gain_driven', 'loss_driven' and 'both'.

Description
The human genome reference used here is RefSeq transcripts in version hg19 from UCSC Genome Browser. The transcripts with NM marker ID, which are protein-codeing, were selected to be our reference database and provided as hg19DBNM.rda.

Usage hg19DBNM
Format A data frame with 39997 rows and 7 variables: Marker.ID RefSeq name with its corrsponding gene symbol

Details
This reference provides region information, including chromosome number, starting position, ending position, strand and gene symbols, for converting copy number alteration data into human genes.
Value data.

create_sample_desc
The create_sample_desc function is provided for creating a sample description table containing all required inputs.

create_gene_exp function
The create_gene_exp function is used to rearrange the input gene expression files into a gene expression list of entire samples.

create_gene_cna function
The create_gene_cna function maps CNA data to human genes and then defines the mapped human genes as CN gain or loss based on the CN threshold, whose default values are set as 2.5 for gain and 1.5 for loss. These mapped genes will be assigned values in +1, -1 or 0, where +1 stands for CNA-gain, -1 stands for CNA-loss and 0 stands for neutral.

find_cna_driven_gene function
The find_cna_driven_gene function identifies CNA-driven differentially expressed genes. The input mapped genes remain for further analyses if its ratio of the number of CN changed samples, CNAgain (G) or CNA-loss (L), to the number of total samples is larger than a given threshold. Here the default setting is that only genes showing CNAs in at least 20 statistical tests, T-test and Wilcoxon rank sum test, are performed in the GE level by classifying the samples as G and L plus Nertral (N) groups or L and G plus N groups, depending on the CN of the interested gene increases or decreases.