Welcome to the SVS Microarray-Based Copy Number Variation (CNV) Univariate Analysis Tutorial!
Updated: January 28, 2021
Packages: CNV Analysis, Power Seat
This tutorial covers a basic workflow for whole genome CNV analysis and association testing using the univariate segmentation process in SVS. The tutorial is built around the Affymetrix 500K array, but the workflows are generally applicable to most SNP microarray platforms as well as most aCGH platforms. There are some anomalies with Illumina microarray data where certain analysis steps may not directly apply. For more information on what to be aware of with Illumina data, see CNV Analysis Tips for Illumina Data.
To complete this tutorial you will need to download and unzip the following file, which includes several datasets.
NOTE: If you have downloaded this file for any of the other Copy Number Variation Analysis Tutorials, you do not need to download it again.
Files included in the above ZIP file:
- Sample CNV Data Project file – contains Phenotype, PCA Corrected Data and Segmentation Output spreadsheets.
- The dataset in this project file is based on the results of the CNV QC tutorial. That tutorial has 200 samples that were reviewed for various quality measurements, with 190 samples eventually selected for analysis. Here, we start with those 190 samples after applying a PCA correction to the Log Ratio Data based on the first 7 principal components. We also provide univariate segmentation results for the full 200 subjects with the 190 active subjects for further analysis.
Copy number analysis a broad term that encompasses many aspects and workflows for interrogating CNVs. Applications include CNV genome-wide association studies (GWAS), cytogenetics, de novo copy number detection between parents and offspring, or “paired analysis”, such as comparing normal cell DNA to tumor cell DNA. Each application has its own set of workflows and data types to work with, but there are many processes that overlap among them. The general workflow is as follows:
- Normalize raw data and generate log ratios. Array manufacturers generate raw intensities for each probe on an array. Many factors can influence the distribution of raw intensity values generated from each chip, and normalizing the raw values is necessary before making comparisons between subjects. The measurement that is most commonly used to determine copy number status is the “log ratio” (LR), also called “log-2 ratio” or “log-R ratio” (“R” is commonly used as a variable to represent probe intensity). The standard formula for LR calculation is log2(observed intensity/reference intensity). LR calculation requires using a reference panel to determine the “normal” or baseline intensity expected at each marker. Choosing a proper reference panel is critical as it can affect all subsequent analyses.
- Perform quality assurance. With copy number data, we have found that significant bias can be introduced from differences between experiments – plate-to-plate, machine-to-machine, and site-to-site – collectively referred to as batch effects. Additional sources of variation include sample extraction and preparation procedures, cell types, temperature fluctuations and even ozone levels within the lab. Beyond batch effects and experimental variability, the signal-to-noise properties of each sample affect the accuracy of CNV detection. Batch effects and noisy data confound results, with complications ranging from poorly defined copy number segments to false-positive and non-replicable association findings. Therefore, there are a number of quality assurance protocols one should consider to either filter problematic samples entirely or correct the LR data before performing further analyses.
- Detect copy number regions. The objective in this step is to determine regions in the genome where a given sample’s mean LR value differs from the reference. These regions are referred to as copy number variants (CNVs). A mean LR around zero indicates that a sample has the same number of copies as the reference. A LR segment mean above zero typically means there is a copy number gain, and a LR segment mean below zero means there is a copy number loss. The higher or lower the LR segment mean, the more or less copies there are respectively.There are numerous algorithms available to detect copy number regions. SVS employs two optimal segmentation algorithms. 1) The univariate method, which considers only one sample at a time, is ideal for detecting rare and/or large CNVs. 2) The multivariate method, which simultaneously considers all subjects, is ideal for detecting small, common CNVs. Some CNV detection algorithms also assign actual copy number calls (0,1,2,3…) for a segment based on thresholds of LR segment means, though there are some reasons why this is problematic (see CNV Analysis Tips for Illumina Data).
- Analyze the copy number variations. Here you look at the LR segment means (or copy number calls) to find differences between two or more groups (case/control, tumor/normal, parent/offspring, etc.). This can be done visually using advanced plotting techniques, or statistically using simple summary statistics or advanced association and regression analysis. We advocate a combination of all to ensure copy number differences are indeed real and not just data processing anomalies.
- Make sense of the findings. This can be done in a number of ways, (looking at genomic annotations, pathway analysis, gene ontologies, etc.) and in some respects, is outside the scope of this tutorial, though we will touch on a few sense-making techniques.
Considering this generalized workflow, you may begin copy number analysis at any of the steps using SVS depending on the data you have access to. For example, if your data provider has only given you access to pre-processed and called copy number data, then you can simply import the copy number calls and begin visualizing the data or performing association analysis. If you’re an Illumina customer and have access to GenomeStudio, you will most likely be able to export the normalized LRs, with which you can export and perform quality assurance or go right into copy number detection.
Whenever possible (and depending on how comfortable you are with upstream analysis), we recommend you start with the raw intensity data as there are a number of things to consider during each step of the process that can have profound ramifications for downstream analysis. This tutorial leads you through each step listed above with a focus on CNV-based GWAS using Affymetrix microarray data, but please keep in mind that no single tutorial can cover every aspect of copy number analysis. As we have stated in our blog post :
“Let’s face it: CNV analysis is not easy. While SNP-based GWAS follows a largely standardized process, there are few conventions for CNV analysis. CNV studies typically require extensive craftsmanship to fully understand the contents of the data and explain the findings. As much as we would like to automate the process, there are always factors that require the researcher to make informed decisions or educated guesses about how to proceed.”
2. CNAM Optimal Segmenting
CNAM employs an optimal segmenting algorithm which uses dynamic programming to detect CNV segment boundaries. There are two methods available in SVS: univariate and multivariate. The univariate method, which considers only one sample at a time, is ideal for detecting rare and/or large CNVs. The multivariate method, which considers all samples simultaneously, is ideal for detecting small, common CNVs. This tutorial leads you through univariate segmentation.
A. Performing CNAM Optimal Segmenting
IMPORTANT: The segmentation algorithm implemented in CNAM is an optimal algorithm, which produces high quality results, but comes at the expense of computation time. As of SVS v7.4, CNAM can take advantage of video graphics cards, or graphical processing units (GPUs), in addition to your computer’s CPUs. This can dramatically decrease computation time while producing the same high-quality results. Depending on computer configuration, the speed increase can be 20 times or more when utilizing the GPU for segmentation. (For more information see Video Graphics and Genomics: A Real Game Changer?)
With that said, segmentation on a whole genome dataset still takes some time. Therefore, rather than have you perform the segmentation yourself, we’ve provided the finished segmentation results in the project, outlining the steps below by which we achieved the results.
- Open PCA-Corrected Data (Center by Marker assumed) – Sheet 1 and choose Numeric >CNAM Optimal Segmenting.
- Set the parameters as they are in Figure 2-1.
The options in Figure 2-1 reflect the fact that the computer used in developing this tutorial has an OpenCL-compatible video graphics card (Nvidia Quadro NVS 295) and is a quad-core (4 threads) machine. When you run this on your own you will need to set these parameters based on your own hardware. For a more thorough description of each parameter in this window, click on the Help button in the lower left corner.
- Click Run to begin the segmentation process. Again, you do not need to actually complete the segmentation for this tutorial as the results are already provided in the project.
When the segmentation finishes, three spreadsheets and a run log are created. The Segmentation Covariates Every Column spreadsheet (Figure 2-2) contains the mean LR value of a given segment with redundant information displayed for every marker in the segment. This spreadsheet is useful for specific plots of the segmentation covariates. The Segmentation Covariates First Column spreadsheet contains the same information as the Segmentation Covariates Every Column spreadsheet, but it only includes columns that correspond to the first column of a segment identified in any of the subjects in the data. This spreadsheet is better for association testing as it is smaller and more accurately reflects the true multiple-testing burden of non-redundant data. The Segment List spreadsheet (Figure 2-3) contains more detailed information about each segment for each subject in a list format. The information corresponding to the 10 samples excluded during the QC process have been inactivated in each spreadsheet.
B. Filtering on Autosomal Segment Count
After segmentation is complete, it is often useful to examine the total number of segments for each subject in the data. An unusually large number of segments is often indicative of data quality problems such as wave effects that were not detected earlier by the log ratios themselves.
- Open the Segment List spreadsheet. Choose Numeric > CNAM Output Analysis > Count Number of Segments Per Sample.
- From the resulting Segment Counts spreadsheet plot the Segment Count column by right-clicking on the column label header and choosing Plot Histogram.
The plot should look like Figure 2-4 (with a bin count of 64).
There appear to be outliers on the right side of the distribution. Sorting the Segment Counts spreadsheet will let us know which samples these are.
- From the Segment Counts spreadsheet right-click the Segment Count column header and choose Sort Descending.
This produces a second segment counts spreadsheet, Segment Counts – Sheet 2. The sample with the highest segment count is S59. Let’s take a look at this sample’s LRs and Segment results to get an idea of what is going on. To plot the sample data, we must first transpose the PCA-corrected LR spreadsheet as well as the Segmentation output, as SVS requires data to be oriented in columns in order to be plotted in the genome browser view.
- Open PCA-Corrected Data (Center by Marker assumed) – Sheet 2 and select Edit > Transpose Spreadsheet.
- Leave the default options and click OK.
- Repeat the transpose steps for Segmentation Covariates Every Column – Sheet 1.
Now we start by plotting the corrected LR values.
- Open PCA-Corrected Data (Center by Marker assumed) Transposed – Sheet 1 and choose GenomeBrowse > Numeric Value Plot.
- Check the S59 box and click Plot.
The plot should look like Figure 2-5.
Now let’s plot the segmentation results on top of the log ratios to see why there were so many segments found.
- From the newly-created log ratio plot of S59 click on the first S59 node in the Plot Tree.
- On the Add tab of the Controls dialog, click Add Item. Click the Project button and select the Segmentation Covariates Every Column Transposed – Sheet 1 then choose S59 from the list and click Plot & Close.
This creates a second S59 under the S59 node. We will need to change the color to differentiate between the two sets and move the new node below the original so we can see the data points
- Right-click on the top S59 item and rename it S59 – Segment Means for the sake of clarity.
- Click on the new S59 – Segment Means item and drag the item below the first S59 item while holding the mouse button.
- To better display the segmentation covariates, select S59 – Segment Means and under the Display tab change Connector to Mid Step and leave the weight at 1. Change the color of the data points under the Style tab by clicking the blue square and selecting green then change the weight to 1.
Your plot should now look like Figure 2-6.
Zoom in for a closer look.
- Double click on 12 in the Domain View to zoom into Chromosome 12.
We can already see a little bit of a wave effect. Applying a small median smooth to the log ratios will make this even more apparent.
- Select the S59 item (for the corrected log ratios) under the Display tab under the Smoothing drop-down select Median Symmetric with a Window Radius of 2.
Even with a median smooth of 2, the wave effect is much more apparent. Zoom in on the y-axis and it becomes more visible (Figure 2-7).
- Click the top S59 node and under the Display tab under the Y-Range click the Manual button and enter -1 – 1
In the previous Quality Assurance tutorial, sample S59 had an absolute wave factor that was just below the IQR cutoff. This suggests we may have benefited by using a more stringent outlier threshold for wave factor for this dataset.
For this tutorial, we will use IQR on the Segment Counts to determine which samples warrant exclusion.
- Open Segment Counts – Sheet 2, and select Numeric > Statistics (per Column).
- Leave 1.5 as the IQR Multiplier, and click OK.
In this case the Upper Outlier Threshold is 228.5. As we can see in the Segment Counts – Sheet 2 spreadsheet, there are 17 samples that have segment counts above this threshold (Figure 2-8).
Let’s exclude these samples from our Phenotype spreadsheet.
- From Segment Counts – Sheet 2, inactivate rows 1 – 17 by clicking once on the first row label, holding the shift key and clicking on the seventeenth row label. This will create a new spreadsheet, Segment Counts – Sheet 3.
- Now got to Select > Apply Current Selection to Second Spreadsheet
- Apply filtered rows to the following spreadsheet: Final Sample Set Dataset – Sheet 1 and click OK. Click OK again to finish.
- From the Final Sample Set Dataset – Sheet 1 create a subset of the active data by choosing Select > Subset Active Data.
- Rename the resulting spreadsheet in the Project Navigator to Final Sample Set Filtered.
- Close all open windows (except the Project Navigator) before continuing.
C. Discretizing Copy Number Segment Covariates
At this point you might continue on to run association tests based on the LR segment covariates. However, it is sometimes useful to discretize the segment means as two state (0,1) or three state (-1,0,1) covariates based on defined thresholds. Discretizing the covariates has the following benefits:
- Approximate copy number calls (potential: deletion, neutral, duplication, not a deletion, or not a duplication) can be made based on thresholds denoting approximate transitions between copy number states.
- Discretizing can magnify small, statistically significant differences between cases and controls.
- Using discretized values reduces the influence of outliers (extremely small or large logR values) on a p-value.
To do this, we first need to determine the appropriate number of copy number classes by examining the segment mean histogram and then discretize the covariates based on the determined thresholds.
The segment means can be accessed from the Segment List spreadsheet.
- Open the Segment List spreadsheet, right-click on the Segment Mean column header (4) and choose Plot Histogram.
The histogram in Figure 2-9 appears (with a bin count of 128).
At first glance we can see the distribution is centered around zero as we’d expect, and we can see some additional peaks to the left and right. If you zoom in on the y-axis you can see the peaks more clearly.
- To zoom click on the Graph 1 node in the Graph Control Interface, and then under the Graph tab enter 800 in the Y-Max text box under the Graph tab and hit Enter.
- You may also wish to adjust the X axis.
You can now begin to see the peaks more clearly (Figure 2-10).
In this case it actually looks like there are four copy number classes. For the purpose of this tutorial we will make it simple and look for thresholds that distinguish a copy number loss, neutral, or gain, which at first glance appears to be around -0.15 and 0.15 (you can zoom in on the X-axis to see this more clearly).
- Beyond visual inspection of the histogram, a more exact approach would be to run CNAM segmenting on the Mean Value column, but that’s beyond the scope of this tutorial.
- The “peaks” or “clusters” of values seen in this histogram indicate a change in intensity relative to the reference intensity. Although we often refer to the lower peak as “losses” and the upper peak as “gains,” it is important to remember that this is relative, and doesn’t always correlate to actual copy numbers of more than 2 copies or less than 2 copies. For example, if there is a common indel where most of the reference samples have the deletion, a subject with two copies present would appear at the far right of the distribution. We might call it a “gain”, but in reality it is 2 copies, which would typically be “neutral.” The discretization approach described here works well for rare CNVs, but isn’t always appropriate for CNVs that are more common, as the discretization thresholds shift according to the frequency of the CNV in the reference population.
- If you observe a uni-modal distribution (only one hill) regardless of Bin Count and/or a large number of outliers in your data, discretizing is not recommended. In this case using the PCA corrected segments as covariates for association testing is recommended.
Now that appropriate thresholds have been determined, you can use these values to discretize the copy number segment covariates (continuous LR values) as three-state (-1,0,1) covariates. We’ll discretize the Segmentation Covariates First Column spreadsheet as we’ll be using this for association testing.
- Open the Segmentation Covariates First Column and select Numeric > CNAM Output Analysis > Discretize CN Segment Covariates with Counts.
You will be prompted to choose a Three State Model or either of two Two State Models.
- Choose Three State Model (-1,0,1).
- Keep -0.15 as the Lower Threshold Level and 0.15 as the Upper Threshold Level and click OK.
This will create two new spreadsheets: Three State Covariates, with -1 for potential losses, 0 for potential neutrals, and 1 for potential gains, and Copy Number State Counts – Mapped Sheet 1 with various statistics about each marker in the segment. We will be using the first spreadsheet as covariates for association testing.
- Close all open windows before continuing except the Project Navigator.
3. CNV Association Analysis
In practice there are several data sources upon which to base CNV association tests, including the Log Ratios, PCA corrected LRs, the mean LR segmentation covariates and the two-state or three-state discretized segment covariates. This tutorial focuses on using the discretized segment covariates for association testing, though the process is similar for using the others.
A. CNV Association Analysis
Before running association tests we first need to join the phenotype information with the discretized segment covariates. Again, the phenotype information in this tutorial is simulated for demonstration purposes.
- Open the Final Sample Set Filtered spreadsheet and select File > Join or Merge Spreadsheets.
- Select the Three State Covariates spreadsheet and click OK.
- Leave the default options in the Join or Merge Spreadsheet window and click OK.
The new spreadsheet, Final Sample Set Filtered + Three State Covariates – Sheet 1, will be used for CNV Association tests.
- From the resulting spreadsheet, left-click the Phenotype 1 column label header once to turn the column magenta, denoting it as the dependent variable.
- Select Numeric > Numeric Association Tests.
- For this tutorial set the same parameters to those in Figure 3-1 and click Run.
The resulting Association Tests spreadsheet contains four columns. To view the most significant markers, sort the Corr/Trend P column ascending.
- Right-click on the Corr/Trend P column and choose ‘Inactivate Missings’. Create a subset by choosing Select > Subset Active Data.
- In Association Tests – Active Subset right-click on the Corr/Trend P column (1) and select Sort Ascending.
SNP_A-2186409 is the most significant marker with a p-value of about 0.00075… (Figure 3-2). It does not meet a Bonferroni correction (column 4). Keep in mind that the Bonferroni correction is based on the total number of tests performed, which includes a lot of redundant tests due to slight inconsistencies in the start and endpoints of common CNV segments, as revealed in the visualization steps below.
- To plot these results, right-click on the Corr/Trend –log10P column header and select Plot Variable in GenomeBrowse.
The resulting plot looks like Figure 3-3. Notice the most significant peak on Chromosome 22.
It may be desirable to visualize the underlying data in conjunction with the association test results. The best way to do this is with a heat map of the segment covariates, sorted by phenotype status.
- Open the Final Sample Set Filtered spreadsheet. Follow the steps outlined above to merge this spreadsheet with the Segmentation Covariates Every Column spreadsheet. Leave the default values and click OK to merge.
- Return to the P-value plot.
- Go to File >Add click Project. Choose Final Sample Set Filtered + Segmentation Covariates Every Column – Sheet 1 and click Heat Map then Plot & Close.
- Click on the Heat Map node in the Plot Tree and under the Display tab change the all of the middle color options to white to better see the gains/losses.
- Click on the Group by tab and under Field: choose Phenotype from the list and click OK.
- Double-click on Chromosome 22 in the Domain View to zoom to that chromosome, and zoom to the region around the p-value peak by clicking and dragging on the x-axis of the p-value plot.
This concludes the association testing aspects of this CNV study. The remaining step outlines some additional visualization and plotting techniques for CNV data.
4. Data Visualization Techniques
Outlined below are additional data visualization techniques that will allow for further examination of your copy number data.
A. Heat Map of All Sample Segment Covariates
Open up the plot previously created. Deselect the pvalue plot so that only the heat map is visible.
- Uncheck the first Corr/Trend -log10 P node in the Plot Tree.
- From the chromosome selection drop-down choose All.
- Select the Heat Map node and change the color options on the Display tab similar to Figure 4-1.
The copy number variations become more apparent with this color scheme. Since the heat map is sorted by case/control status, color-consistent streaks on either the lower or upper half of the plot could signal an association. Though at first glance there are not any obvious differences between cases and control there are a couple regions of interest. Figure 4-1 highlights large CNVs for individual samples in green, and copy variations among all the samples in orange.
Let’s look at one of the common CNVs first.
- Enter 17:38,879,225-42,979,664 in the genomic region box on the top of the window.
The plot should now look like Figure 4-2 with three common CNVs shown. Notice that not every sample has the exact same starting and ending boundaries for each region. This is common when using the univariate segmentation method. The consequence is that the output spreadsheet containing the “first column” of each CNV segment will contain a substantial amount of redundant data. As it is likely that the same feature (probably a common indel region) is being detected in each of the subjects, the overlapping area is probably the most correct representation of the underlying biology. The multivariate segmentation method will determine the most likely endpoints by comparing data across all samples and will result in uniform endpoints for common CNV regions.
Let’s add the Database of Genomic Variants annotation track to see how this region has been catalogued and the Affymetrix 500K probe track to see how dense the genotyping was around these areas.
- Click on the Add button at the top left of the plot window.
- Under the Public Annotations sources check Affymetrix 500K na30, GHI and DGV Variants 2011-05-26, DGV and click Plot & Close.
- Enter 17:41,551,267 – 41,723,784 in the genomic region box on the top of the window.
The plot now looks like Figure 4-3.
Notice that these regions have been extensively cataloged by the Database of Genomic Variants. You can zoom into each region to explore them further.
Now let’s take a look at the individual samples with larger chromosomal aberrations.
B. Investigating Individual Samples
You can now plot the LRs and segmentation covariates together to investigate these samples of interest.
- From the heat map created in Figure 4-3 select 2 from the chromosome drop-down.
Notice the large gain in the left of the plot (Figure 4-4).
To see which sample this is, zoom in on the Y-axis and click on the blue streak.
- Click and drag on the Y-axis to highlight an area around the blue streak to zoom into.
- Click on the blue streak. Notice in the Console that this is sample S150.
Let’s add this sample’s LRs to the plot.
- Click the Add button, click Project select PCA-Corrected Data (Center by Marker assumed) Transposed – Sheet 1 and click S150 then click Plot & Close.
Now let’s plot the segmentation covariates on top of the LRs as before.
- Select the S150 graph node (top) in the Plot Tree. Under the Add tab select Add Item(s).
- Click Project and select Segmentation Covariates Every Column Transposed – Sheet 1 then click the box for S150 and click Plot & Close.
- Move the first S150 graph item below the second as before by clicking and dragging it above the first S150 graph item. Then change the color to green and add the connect as before.
The graph should look like Figure 4-5. Notice the shift in the LRs and segmentation covariates for that region.
Congratulations! You have now worked your way through an entire copy number analysis project, no easy task. We wish you the best of luck on your own study. As always, if you have any questions or need help with any portion of this tutorial or your own analysis, please give us a call. We’d be happy to help!