iExCN algorithm employs Bayesian statistics to identify disease genes with significant associations between copy-number alterations and gene expression changes. (GitHub)
The integrated expression and copy-number analysis pipeline, iExCN, was developed to identify genes whose RNA expression is likely to be driven by DNA copy-number alteration. The core rationale is that a biologically meaningful copy-number driver should show concordance between genomic dosage and transcript abundance. Genes gained or amplified at the DNA level are expected, in aggregate, to show increased expression relative to copy-number-neutral tumors, whereas genes affected by copy-number loss are expected to show decreased expression. iExCN formalizes this relationship using Bayesian estimation rather than relying on linear correlation or permutation-based comparisons.
This section is intended to provide practical guidance for readers who wish to reproduce or adapt the iExCN methodology in their own datasets. Although the original study applied iExCN to rhabdomyosarcoma tumor samples profiled by SNP array and exon expression array, the pipeline can be adapted to other cancer types and other molecular platforms, provided that two essential inputs are available from the same specimens: gene-level copy-number status and normalized gene-expression measurements.
The iExCN workflow begins with paired copy-number and expression data from the same tumor specimens. Copy-number calls are first reduced to gene-level alteration states, and expression data are normalized at the gene or transcript level. For each gene, samples are divided into copy-number-defined groups, typically copy-number neutral, copy-number gain, and copy-number loss. Expression values are then compared across these groups using Bayesian estimation. Instead of asking whether a linear correlation coefficient is statistically significant, iExCN asks whether the posterior distribution supports the expected direction of expression change. For copy-number gains, the relevant probability is the probability that expression in the gained group is greater than expression in the neutral group. For copy-number losses, the relevant probability is the probability that expression in the loss group is lower than expression in the neutral group.
The practical workflow shown in the accompanying schematic consists of four main computational steps. First, splitdata.R prepares and partitions the copy-number and expression matrices for downstream analysis. Second, paraestim_MCMC-Gibs.R processes each paired copy-number and expression block and launches Bayesian parameter estimation. Third, BESTmcmc2.R performs the Bayesian estimation using Markov chain Monte Carlo sampling. Fourth, BayesStatistics.R summarizes the Bayesian results and extracts candidate disease genes whose expression is strongly associated with copy-number status. This workflow corresponds conceptually to the three-stage pipeline described in the Methods: matrix generation, Bayesian estimation, and final statistical summarization.
The minimum input for iExCN is a pair of matched, gene-by-sample matrices. The first matrix contains gene-level copy-number status, and the second contains normalized gene-expression values. The two matrices must have the same genes in the same order and the same samples in the same order. This point is critical because the workflow assumes that columns correspond directly between the copy-number and expression files. If sample order is mismatched, iExCN will compare the copy-number status of one tumor with the expression value of another tumor, producing invalid results.
The copy-number matrix should be a tab-delimited text file, conventionally named CNV.txt, with one row per gene and one column per sample. The first column should contain gene identifiers. In the current workflow, Ensembl gene identifiers are preferred because the summary script maps Ensembl IDs to gene symbols using org.Hs.eg.db. Copy-number states are represented numerically. In the default implementation, 0 indicates copy-number neutral, 1 indicates copy-number gain, and -1 indicates copy-number loss. The script paraestim_MCMC-Gibs.R uses these values explicitly when assigning samples to neutral, gain, and loss groups. (GitHub)
A simplified example of the required copy-number matrix is:
GenesENS Sample_001 Sample_002 Sample_003 Sample_004
ENSG00000141510 0 1 0 -1
ENSG00000171862 0 0 1 0
ENSG00000139618 -1 0 0 1
The expression matrix should also be a tab-delimited file, conventionally named Expression.txt, with the same first-column gene identifiers and the same sample columns in the same order. Expression values should be normalized before running iExCN. In the study described here, exon-array expression data were normalized using Robust Multiarray Average normalization. For RNA-seq datasets, users should provide a comparable normalized continuous expression measure, such as log-transformed normalized counts, variance-stabilized counts, or log-transformed TPM values, after appropriate filtering and batch correction.
A simplified example of the expression matrix is:
GenesENS Sample_001 Sample_002 Sample_003 Sample_004
ENSG00000141510 7.41 8.36 7.58 6.22
ENSG00000171862 5.94 5.88 7.10 6.01
ENSG00000139618 4.32 5.49 5.61 6.27
Before running iExCN, users should confirm that gene identifiers are harmonized between copy-number and expression datasets. If multiple probes or transcripts map to the same gene, the investigator should decide whether to collapse these features to a single gene-level value or retain them as separate entries. The current splitdata.R script makes duplicate gene identifiers unique by appending suffixes to repeated entries, but this should not be treated as a substitute for thoughtful preprocessing. (GitHub)
iExCN does not require raw SNP array, sequencing, or expression files as direct input. Instead, raw molecular data must first be processed into gene-level copy-number and expression matrices. In the rhabdomyosarcoma study, copy-number data from different sources were processed using a consistent segmentation strategy to minimize platform-specific bias. SNP-array and molecular inversion probe array data were imported into Nexus BioDiscovery, corrected for GC content, segmented with the SNP-FASST algorithm, and assigned copy-number states representing gain, amplification, single-copy loss, or homozygous deletion. Recurrent copy-number alterations were additionally evaluated using GISTIC.
For practical implementation, the key requirement is that each gene in each sample be assigned a discrete copy-number state suitable for comparison. A typical implementation would map segmented copy-number data to genes, then recode the gene-level result into the three-state format required by iExCN: loss, neutral, and gain. If the investigator wishes to distinguish high-level amplification from low-level gain, or homozygous deletion from single-copy loss, the default iExCN scripts must be modified because the standard workflow uses only -1, 0, and 1 categories.
When preparing copy-number calls, users should define the thresholds used to assign gain and loss before beginning the iExCN analysis. These thresholds should be applied consistently across all samples. If multiple platforms or batches are combined, platform-specific artifacts should be evaluated carefully. The biological interpretation of iExCN depends on accurate copy-number classification; systematic overcalling or undercalling of gains and losses will directly affect the downstream Bayesian comparisons.
Expression values should be normalized, quality controlled, and transformed to a scale appropriate for comparing expression distributions across samples. In array-based datasets, robust multiarray normalization or a comparable platform-specific normalization method should be applied. For RNA-seq, raw counts should not be used directly. Instead, investigators should use normalized and transformed expression values that reduce mean-variance dependence and permit meaningful between-sample comparison.
Genes with very low or absent expression across nearly all samples should be filtered before running iExCN, because such genes are unlikely to produce stable Bayesian estimates and may lead to uninterpretable results. Similarly, expression outliers should be examined. iExCN is designed to be more flexible than linear correlation analysis, but it is still sensitive to poor-quality input data. If the dataset contains known technical batches, such as sequencing center, expression platform, RNA quality, or library preparation batch, batch effects should be evaluated before the copy-number-expression integration step.
The first practical step in the workflow is to place CNV.txt and Expression.txt in the working directory. These files should be tab-delimited and should contain the same genes and samples in matched order. The splitdata.R script reads these two files, adds a genes column, and divides the genome-wide matrices into smaller blocks for downstream processing. The block size is passed as a command-line argument. For example:
Rscript splitdata.R 500In this example, the script splits the full input matrices into files containing approximately 500 genes each. The output files are named sequentially, for example:
CNV-1.txt
Expression-1.txt
CNV-2.txt
Expression-2.txt
CNV-3.txt
Expression-3.txt
...
The block size can be adjusted depending on computational resources. Smaller blocks make job scheduling easier and reduce the cost of re-running failed jobs. Larger blocks reduce the total number of files but may be less convenient on a high-performance computing cluster. In genome-wide studies, splitting the data is strongly recommended because the MCMC-based Bayesian estimation is computationally intensive.
Before proceeding, users should verify that each CNV-i.txt file has a corresponding Expression-i.txt file and that the number of rows in each paired block is identical. A simple manual check of the first few lines is often sufficient to detect formatting problems:
head CNV-1.txt
head Expression-1.txtAfter the matrices are split, each paired copy-number and expression block is analyzed using paraestim_MCMC-Gibs.R. This script reads a copy-number block and the corresponding expression block, separates samples into neutral, gain, and loss groups for each gene, and calls BESTmcmc2.R to perform Bayesian estimation. The current script reads the CNV and expression filenames as command-line arguments and writes an output file prefixed with BayesEst2-. (GitHub)
A typical command for one block is:
Rscript paraestim_MCMC-Gibs.R CNV-1.txt Expression-1.txtFor multiple blocks, the analysis can be run in a loop:
for i in $(seq 1 100)
do
Rscript paraestim_MCMC-Gibs.R CNV-${i}.txt Expression-${i}.txt
doneOn a cluster, the same operation should usually be submitted as an array job so that blocks are processed in parallel. For example, the conceptual structure of an array job is:
Rscript paraestim_MCMC-Gibs.R CNV-${TASK_ID}.txt Expression-${TASK_ID}.txtwhere TASK_ID is replaced by the job-array index used by the local scheduler. The original implementation was designed for use on a high-performance Linux cluster. Because MCMC is run repeatedly across thousands of genes, wall time can be substantial, and users should anticipate the need for parallel execution.
For each gene, the script calculates median expression in the loss, neutral, and gain groups. It then performs Bayesian estimation comparing neutral versus gain and neutral versus loss expression groups. For gains, the posterior quantity of interest is the fraction of MCMC samples in which expression in the gain group exceeds expression in the neutral group. For losses, the posterior quantity of interest is the fraction of MCMC samples in which expression in the loss group is below expression in the neutral group. These values are reported as RateGain and RateLoss, respectively. (GitHub)
The output file from each block has the following columns:
GenesENS ExpLoss ExpNor ExpGain RateLoss RateGain
These columns can be interpreted as follows. GenesENS is the gene identifier. ExpLoss is the median expression among samples with copy-number loss. ExpNor is the median expression among copy-number-neutral samples. ExpGain is the median expression among samples with copy-number gain. RateLoss is the posterior probability that expression is lower in the loss group than in the neutral group. RateGain is the posterior probability that expression is higher in the gain group than in the neutral group.
The Bayesian estimation is performed by BESTmcmc2.R, which implements a two-group Bayesian comparison using MCMC. In the current script, the default settings include 100,000 saved MCMC steps, 1,000 burn-in steps, three chains, and optional parallel execution when sufficient computational cores are available. The function returns posterior samples for group means and variance-related parameters and also records convergence-related attributes such as Gelman diagnostic values and effective sample size. (GitHub)
The important practical distinction between iExCN and a conventional frequentist test is that iExCN does not simply ask whether two expression groups have different means. Instead, it estimates the posterior probability that the expression difference is in the biologically expected direction. For a copy-number gain, the relevant question is: what is the posterior probability that the gain group has higher expression than the neutral group? For a copy-number loss, the relevant question is: what is the posterior probability that the loss group has lower expression than the neutral group?
This direction-aware interpretation is central to the method. A gene with a statistically detectable expression difference in the wrong direction should not be considered a canonical copy-number-driven expression candidate. Conversely, a gene whose expression increases with gain and decreases with loss, especially with a large posterior probability and consistent median ordering, is a strong iExCN candidate.
After all matrix blocks have been analyzed, the block-level output files should be combined into a single genome-wide result file. Care should be taken to retain the header only once. For example:
head -n 1 BayesEst2-CNV-1.txt > BayesEst2-CNV-all.txt
for f in BayesEst2-CNV-*.txt
do
tail -n +2 "$f" >> BayesEst2-CNV-all.txt
doneAfter concatenation, users should confirm that the combined file has the expected number of rows and that no blocks failed silently. A simple check is to count the number of genes in the original CNV.txt file and compare it with the number of data rows in the combined Bayesian output.
The final step is to summarize the Bayesian results using BayesStatistics.R. This script reads the combined Bayesian output and applies a posterior probability threshold. The default candidate-selection logic requires complete data, ordered median expression across copy-number states, and posterior probabilities greater than the selected cutoff. In practical terms, a stringent iExCN-positive gene should satisfy the following pattern:
ExpLoss < ExpNor < ExpGain
RateLoss > cutoff
RateGain > cutoff
For the threshold used in the study, the cutoff should be set to 0.999, corresponding to a posterior probability greater than 99.9% in the expected direction. A typical command is:
Rscript BayesStatistics.R BayesEst2-CNV-all.txt 0.999The output file is named according to the selected cutoff, for example:
BayesEst_CancerRelatedGene_cutf0.999.xls
The summary script also maps Ensembl gene identifiers to gene symbols using org.Hs.eg.db, which is why consistent gene identifiers are important. (GitHub)
Users should interpret the cutoff as a Bayesian posterior probability threshold, not as a p-value or a false-discovery-rate-adjusted q-value. A value of RateGain = 0.999, for example, means that 99.9% of posterior samples support higher expression in the copy-number gain group than in the neutral group. It does not mean that the frequentist probability of a false positive is 0.1%.
An iExCN-positive gene is a gene whose expression is strongly associated with copy-number status in the expected direction. Such genes are candidates for dosage-sensitive disease biology. However, iExCN should be interpreted as a prioritization method, not as definitive proof that a gene is a driver. Copy-number regions often contain many genes, and some genes may show expression changes because they are passengers within a recurrently altered region. Biological interpretation should therefore integrate additional evidence, including recurrence of the copy-number alteration, magnitude of expression change, known cancer biology, pathway context, dependency-screen results, and independent validation.
A particularly strong candidate would meet several criteria. It would reside in a recurrently altered genomic region, show ordered expression across loss, neutral, and gain states, pass the stringent Bayesian posterior threshold, have sufficient numbers of samples in each copy-number group, and retain the copy-number-expression relationship in independent or stratified analyses. In contrast, genes altered in only one or two samples should be treated cautiously, even if the posterior probability appears high, because limited group size can make estimates unstable or non-generalizable.
For cancer datasets with known molecular subtypes, investigators should also evaluate whether iExCN results are confounded by subtype. In rhabdomyosarcoma, for example, fusion status, histology, and other molecular features may influence both copy-number landscapes and expression programs. A gene may appear copy-number associated in the full cohort but actually reflect subtype-specific biology. Therefore, when sample size permits, iExCN analysis should be repeated within relevant biological strata or adjusted through complementary analyses.
Several quality-control steps should be performed before accepting iExCN results.
First, confirm sample matching. The copy-number and expression matrices must contain the same samples in the same order. A sample manifest should be maintained and archived with the analysis. If sample IDs were renamed or harmonized across datasets, the mapping file should be retained.
Second, confirm copy-number coding. The current implementation assumes that 0 denotes neutral copy number, 1 denotes gain, and -1 denotes loss. Other coding systems must be recoded before analysis. For example, if the upstream copy-number pipeline uses 2 for amplification or -2 for homozygous deletion, those values must either be collapsed into 1 and -1 or the iExCN scripts must be modified to treat them as separate categories.
Third, confirm expression normalization. iExCN should be run on normalized expression values, not raw intensities or raw read counts. Expression data should be inspected for outlier samples, batch effects, and genes with uniformly absent expression.
Fourth, confirm sufficient group sizes. The default script reports NA when there are too few samples in the neutral, gain, or loss groups. This behavior is useful because genes with very sparse copy-number events are difficult to evaluate robustly. If a study is designed to evaluate only gains or only losses, the summary logic can be modified, but such changes should be described explicitly.
Fifth, visually inspect top candidates. For each high-priority gene, users should generate a plot of expression values stratified by copy-number state. This simple visualization often identifies artifacts, such as a result driven by a single outlier sample, a batch effect, or incorrect copy-number coding.
Sixth, monitor MCMC behavior. Genome-wide MCMC analyses can produce occasional unstable estimates. When possible, users should inspect convergence diagnostics for representative genes, especially top candidates and borderline cases. If MCMC convergence is poor, increasing the number of saved steps, adjusting burn-in, or re-running problematic genes may be necessary.
To make iExCN analyses reproducible, manuscripts using the method should report the following information.
The study should describe the number of samples with paired copy-number and expression data, the source and platform of each data type, and any inclusion or exclusion criteria. If multiple cohorts or platforms were combined, the normalization and harmonization strategy should be stated clearly.
The copy-number-processing section should describe how raw copy-number data were generated, segmented, corrected, and mapped to genes. It should also state the thresholds used to define gain and loss and whether amplifications and deletions were collapsed into broader gain and loss categories.
The expression-processing section should describe the normalization method, transformation, filtering criteria, and any batch correction. For array-based data, the probe-to-gene mapping strategy should be reported. For RNA-seq data, the count normalization and transformation method should be provided.
The iExCN implementation section should specify the version of the scripts used, the software environment, the R packages required, and whether analyses were run on a cluster or local machine. The block size used by splitdata.R, the MCMC settings, and the posterior probability cutoff should be reported. If the default cutoff of 0.999 was used, the manuscript should state that genes were required to show greater than 99.9% posterior support for the expected expression direction.
The results section should report the number of genes tested, the number of genes with sufficient copy-number group sizes, and the number of genes passing the final iExCN threshold. For top candidates, it is useful to report median expression in loss, neutral, and gain groups, as well as RateLoss and RateGain.
The most common implementation error is mismatched sample order between CNV.txt and Expression.txt. Because the scripts assume that corresponding columns represent the same tumor sample, this error can produce plausible-looking but biologically meaningless results. Before running the analysis, users should programmatically confirm that column names are identical and identically ordered.
A second common problem is inconsistent gene identifiers. If copy-number data are summarized by gene symbol but expression data are summarized by Ensembl ID, the matrices must be harmonized before running iExCN. Users should avoid mixing identifier systems unless a stable mapping has been applied.
A third pitfall is inappropriate copy-number state coding. Values from upstream tools may not match the -1, 0, and 1 system expected by iExCN. All copy-number states should be recoded explicitly, and the recoding rules should be documented.
A fourth issue is insufficient numbers of altered samples. Genes rarely affected by gains or losses may not be analyzable by the default pipeline. This limitation should not be treated as a failure of the method. Rather, it reflects the fact that copy-number-expression relationships require enough samples in each group to estimate expression differences reliably.
A fifth issue is overinterpretation. iExCN identifies genes whose expression is associated with copy-number status; it does not by itself establish that the gene is required for tumor growth or that it is the functional target of the copy-number event. Downstream biological validation remains necessary.
A practical genome-wide iExCN run can be organized as follows:
# 1. Prepare required files in the working directory:
# CNV.txt
# Expression.txt
# 2. Split matrices into blocks of 500 genes.
Rscript splitdata.R 500
# 3. Run Bayesian estimation for each block.
# Replace 100 with the actual number of generated CNV/Expression block pairs.
for i in $(seq 1 100)
do
Rscript paraestim_MCMC-Gibs.R CNV-${i}.txt Expression-${i}.txt
done
# 4. Combine block-level outputs.
head -n 1 BayesEst2-CNV-1.txt > BayesEst2-CNV-all.txt
for f in BayesEst2-CNV-*.txt
do
tail -n +2 "$f" >> BayesEst2-CNV-all.txt
done
# 5. Extract high-confidence iExCN genes using the 99.9% posterior cutoff.
Rscript BayesStatistics.R BayesEst2-CNV-all.txt 0.999This example is intended to illustrate the logic of the analysis. In production use, especially for large cancer cohorts, the Bayesian estimation step should be parallelized on a high-performance computing cluster. Each CNV-i.txt and Expression-i.txt pair can be treated as an independent job.
When applying iExCN to a new disease context, investigators should first define the biological question. If the goal is to identify genes activated by copy-number gain, then gain-focused criteria may be appropriate. If the goal is to identify genes suppressed by deletion, then loss-focused criteria may be more relevant. The default summary script is conservative because it prioritizes genes showing the expected ordered pattern across loss, neutral, and gain states. This is powerful for identifying bidirectional dosage-sensitive genes but may miss genes that are recurrently gained without recurrent losses, or recurrently lost without recurrent gains.
Therefore, users may adapt the final filtering step depending on the study design. For gain-focused analysis, a candidate might be defined by ExpGain > ExpNor and RateGain > 0.999. For loss-focused analysis, a candidate might be defined by ExpLoss < ExpNor and RateLoss > 0.999. If these modified criteria are used, they should be described clearly and should not be presented as the unchanged default iExCN selection rule.
Users should also consider whether covariates need to be addressed. The default implementation compares expression distributions by copy-number group without modeling additional covariates. In datasets where subtype, purity, ploidy, sex, platform, or batch strongly affects both copy-number status and expression, iExCN results should be interpreted alongside stratified or covariate-aware analyses.
After identifying candidate iExCN genes, validation should proceed at several levels. Computational validation can include re-running the analysis in an independent cohort, performing leave-one-cohort-out analysis, testing robustness to alternative copy-number thresholds, and comparing results across molecular subtypes. Visualization of top candidates is strongly recommended, with expression plotted by copy-number state.
Biological validation may include comparison with known cancer genes, pathway enrichment analysis, dependency data, CRISPR/Cas9 perturbation experiments, RNA interference, or pharmacologic studies. In the described rhabdomyosarcoma study, copy-number status in cell lines used for CRISPR/Cas9 screens was also evaluated using whole-exome sequencing data, allowing candidate genes to be considered in the context of both tumor copy-number-expression association and experimental model copy-number status.
In practical terms, iExCN converts paired gene-level copy-number and expression data into a Bayesian, direction-aware ranking of candidate dosage-sensitive genes. The method is especially useful when the relationship between copy number and expression may be non-linear, heterogeneous, or difficult to capture with simple correlation analysis. Successful implementation depends on careful preprocessing, exact matching of copy-number and expression matrices, appropriate copy-number recoding, sufficient altered-sample counts, and transparent reporting of the Bayesian threshold used for candidate selection. When these steps are followed, iExCN can serve not only as an analytical method for a specific study but also as a reusable framework for discovering copy-number-driven expression programs across cancer datasets.