SAAS-CNV: A Joint Segmentation Approach on Aggregated and Allele Specific Signals for the Identification of Somatic Copy Number Alterations with Next-Generation Sequencing Data
[['Profound somatic copy number alternations ( SCNAs ) are present in many types of tumors [1–3] , where they affect a larger fraction of the genome than other types of somatic variations [3 , 4] .', 'The roles of SCNAs in promoting oncogenesis and tumor progression are under intensive study .', 'In contrast to germline copy number variations ( CNVs ) , which are sparsely distributed along the genome and of small to moderate size , tumor SCNAs are large in size and have a much wider range of magnitudes in copy number .', 'Accurate detection and characterization of genome-wide SCNA profile are further complicated by aneuploidy and heterogeneity of tumor cells and contamination of normal cells [5] .', 'Array comparative genomic hybridization ( array-CGH ) [6] and single nucleotide polymorphism ( SNP ) array [7] were widely used in surveying genome-wide SCNAs in the past decade .', 'More recently , next-generation sequencing ( NGS ) technology provides unprecedented resolution to comprehensively characterize SCNAs at rapidly decreasing cost [8 , 9] .', 'A long list of tools has been developed and successfully applied in analyzing CNV data harvested from NGS; see [2 , 10 , 11] for reviews .', 'The available methods mainly fall into four categories:', '1 ) read depth ( RD ) ,', '2 ) pair-end mapping ( PEM ) ,', '3 ) split read ( SR ) and', '4 ) Assembling ( AS ) .', 'While they take advantage of complementary information , featuring different perspectives of CNVs , each of them encounters important limitations , due to the complications from the high dynamics of cancer genome .', 'With whole genome sequencing ( WGS ) data , the latter three methods are good at detecting CNVs of small to moderate sizes , locating CNV break points at higher resolution , and discovering copy neutral rearrangements , but they are less capable in characterizing large-size and wide-range copy number changes at the genome-wide scale .', 'Moreover , their applicability in whole exome sequencing ( WES ) data is limited by the fact that short reads from WES are concentrated in interspersed genomic regions .', 'Therefore , the RD-based methods are used more widely in the study of tumor CNVs with both WGS and WES data [10] .', 'Basically , the majority of RD-based methods , such as CNV-seq [12] , SegSeq [9] , ExomeCNV [13] and PatternCNV [14] , follows a “bottom-up” procedure: short-reads mapping , normalization of read depth , copy number estimation in a local region ( usually in a window of certain size or in an exonic region ) and segmentation to merge regions with the same copy number status [10] .', 'This strategy is sensitive in the detection of germline CNVs , but for tumor CNVs ( i . e . SCNAs ) , it is difficult for the local inference to correctly decide the baseline of ploidy and accurately discern weak signals of copy number change in the presence of aneuploidy and normal cell contamination , so that the genome-wide inference drawn from the later segmentation step is inclined to accumulate false positive findings from the earlier local inference step .', 'Further , these methods only consider the aggregated depth of sequencing reads carrying paternal and maternal alleles , aimed at estimating the total copy number , but largely ignore allele specific read depth , while the latter contains critical information of copy number change , copy-neutral loss of heterozygosity ( CN-LOH ) , and , importantly , genome ploidy .', 'In analogy to Illumina SNP array [7] , the allele specific signal from NGS can be quantified at heterozygote sites ( e . g . SNP loci ) and converted to the so-called B allele frequency ( BAF ) , which is defined as the proportion of the reads carrying “B” allele ( i . e . non-reference allele ) among all the mappable reads at that site .', 'At an SCNA-free locus , the expectation of BAF is close to 1/2 , indicating equal amount of paternal and maternal alleles .', 'With an SCNA event , the expectation of BAF might deviate from 1/2 , reflecting differential copy number changes in the two alleles .', 'The two data dimensions , total read depth and allele-specific signal ( i . e . , BAF ) , carry complementary and consensual information of SCNA events .', 'Incorporating both quantities could theoretically improve power and accuracy in identifying SCNA segment boundaries and characterizing alteration types ( i . e . , gains , losses , and CN-LOHs ) .', 'Control-FREEC [15] is a representative method that incorporates both types of information .', 'However , it uses them separately rather than jointly in the segmentation and SCNA calling , nor does it make good use of data from match normal sample to generate reliable allele-specific signal .', 'Herein , we propose a tumor SCNA analysis method , SAAS-CNV , based on a joint segmentation algorithm [16] , to accommodate both total read depth and BAF .', 'Using both data dimensions simultaneously , the method first takes a “top-down” strategy to partition the genome into segments with different alterations by joint segmentation , and subsequently , determines their alteration types .', 'The method was designed for paired normal-tumor settings , in which the spatial variability from non-uniform distribution of short reads along the genome can be alleviated [13] .', 'It is able to accommodate WGS and WES data as well as SNP array data .', 'In this paper , we first characterized the type-I-error properties of SAAS-CNV using a carefully designed “null” dataset containing no SCNAs ( H0 ) , created by pairing sequencing replicates of a single HapMap sample as normal-tumor pairs [17] ( Materials and Methods , Dataset I ) .', 'Further , we demonstrated its ability to analyze a real-world WGS study of 88 hepatocellular carcinoma ( HCC ) samples ( Ha ) , some of which exhibit highly complicated SCNA profiles [18 , 19] ( Materials and Methods , Dataset II ) .', 'In both settings ( Ho and Ha ) , we compared the performance of SAAS-CNV with other existing methods .', 'We also explored its potential extension in the estimation of ploidy and purity of tumor cells .'], ['We used GATK variant analysis pipeline [21–23] to process raw fastq files in this study [17] ( more details in S1B Text ) .', 'As a standard output of the pipeline , single nucleotide variants ( SNVs ) and small insertion and deletions ( Indels ) , called from mapped high-quality reads , were saved in the variant call format ( VCF ) files [24] .', 'At each variant locus , we extracted the genotype and the depth of reads carrying reference allele and alternative allele , all informative for SCNA inference .', 'Fig 1 overviews the workflow of our algorithm in the paired normal-tumor setting .', 'For both datasets , a standard implementation of NGS analysis pipeline following the GATK best practices for variant detection [17 , 21 , 22] was applied to the raw FASTQ files to generate recalibrated and deduplicated high-quality BAM files [27] , as well as VCF files containing detected SNVs and Indels .', 'Information retrieved from biallelic heterozygous sites were processed by our SCNA analysis pipeline ( Fig 1 ) .', 'Illumina SNP array data was processed to generate log2ratio and log2mBAF signals in a similar way for downstream analysis .', 'Detailed data analysis steps are described in S1B Text .'], ['First , we present a visualization of the processed signals and the results from SAAS-CNV .', 'Fig 2 demonstrates a typical example taken from the analysis of Dataset II .', 'The 2-dimensional profiles ( log2ratio and log2mBAF ) from SNP array ( Fig 2A and 2B ) and WGS ( Fig 2C and 2D ) were processed for genome-wide SCNA detection .', 'The signals in the two dimensions are altered by copy number gains and losses , showing consistent patterns of changes , for both platforms .', 'The alteration status assigned to each segment is highly consistent between the two platforms ( Fig 2A and 2C ) .', 'We projected the medians of log2ratio and log2mBAF of each segment onto a 2-D space for SNP array and WGS respectively , confirming the consistency of the inference drawn from the two platforms ( Fig 2B and 2D ) .', 'Segments with the same alteration status are clustered together .', 'Segments with less difference in log2ratio dimension may be distinguished in log2mBAF dimension , suggesting the valuable information added from BAF for SCNA inference .', 'It is noticed that the original baselines in both dimensions deviate substantially from zero ( Fig 2B and 2D ) , and simply using zero as baseline for SCNA inference would result in many false positive calls , underscoring the necessity of baseline adjustment .', 'We processed the HKU HCC WGS data , consisting of 88 pairs of tumor and adjacent normal tissues , with the standard implementation of GATK pipeline [17] .', 'The average read depths of 86 pairs are moderate , ranging from 24 . 8x to 71 . 5x , with two pairs sequenced at higher depth ( >120x ) .', 'The average effective read depths measured at heterozygous sites range from 18 . 2x to 99 . 8x , composing 68%~90% of available read depth ( S1 Fig ) .', 'We applied SAAS-CNV , CNVnorm and Control-FREEC on the 88 pairs of normal-HCC samples .', 'Further , we synthesized WES data by retrieving reads located in exonic regions from WGS data .', 'The synthesized WES data was used to test the performance of SAAS-CNV and ExomeCNV .', 'We also applied SAAS-CNV and GAP on the SNP array data from the same samples .', 'As stated above , the GAP results on SNP array data were used as “truth” in benchmarking other SCNA methods .', 'Some basic metrics from these analyses are summarized in Fig 6 .', 'WGS provides millions of data points ( loci ) per sample while SNP array and WES provides hundreds of thousands of data points ( loci ) ( Fig 6A ) .', 'Since SAAS-CNV takes heterozygous loci as input , it utilizes less number of loci than competing methods .', 'For example , only tens of thousands of loci were used by SAAS-SNV in WES data .', 'On SNP array data , SAAS-CNV and GAP produced comparable number of segments per sample ( Fig 6B ) , as well as comparable segment size in terms of locus number per segment ( Fig 6C ) and physical length ( Fig 6D ) .', 'On WES and WGS data , the competing methods tended to chop the genome into smaller segments ( Fig 6B ) than SAAS-CNV ( Fig 6C and 6D ) .', 'We highlighted the potential usage of BAF information in improving the estimation of absolute copy number and tumor purity on the HKU HCC sample PT116 in Dataset II ( Materials and Methods ) .', 'This sample exhibits complicated SCNA profile ( Fig 9 ) .', 'We applied SAAS-CNV and CNAnorm together to analyze this sample .', 'In CNAnorm analysis , the ratio of tumor versus normal read depth was calculated per 1kb window along the genome ( gray dots in Fig 10A ) .', 'A smoothing approach was then employed by CNAnorm [5] to reduce the random error variability ( fitted black curve in Fig 10A ) .', 'The distribution of smoothed ratio signal clearly showed seven major peaks corresponding to different copy numbers ( Fig 10B ) .', 'CNAnorm attempted to fit a Gaussian mixture model along with Akaike’s information criterion ( AIC ) onto the distribution in order to identify the number of components ( peaks ) and their locations .', 'However , CNAnorm was only able to identify five out of the seven peaks ( Fig 10C ) .', 'It then searched different configurations , where plausible copy numbers and the five identified peak centers were aligned in different ways , and chose the most likely configuration that resulted in the best correspondance .', 'Here , the correspondance was measured by the goodness-of-fit ( R2 value ) of the linear regression model of peak centers on estimated copy numbers ( Fig 10C ) .', 'In this case , CNAnorm assigned the most common component ( the highest peak in Fig 10A ) to copy number 2 ( i . e . tumor-normal ratio 1 ) .', 'Finally , CNAnorm estimated tumor purity ( ρ ) based on the relationship:', 'μiμCN=2−1= ( CNi2−1 ) ⋅ρ', '( 5 )', 'where μi is the identified peak center associated with copy number CNi and μCN = 2 corresponds to peak center associated with copy number 2 .', 'With the estimated tumor purity , the absolute copy number of each segment was obtained ( S13A Fig ) .', 'By a visual check of Fig 10B and 10C , it is obvious that CNAnorm failed to assign copy number correctly .', 'Without accounting for the allelic imbalance pattern information ( Fig 9 ) , CNAnorm was not able to correctly infer absolute copy number in such a complicated cancer genome .', 'For example , Chromosomes 3 , 9 , 11 , 12 , 14 and 15 were estimated to be diploid ( S13A Fig ) , contradictory to the strong pattern of allelic imbalance manifested in log2mBAF dimension ( Fig 9 ) ; Chromosomes 2 , 6 and 18 were estimated to be approximately triploid , also contradictory to the log2mBAF pattern showing allelic balance ( Fig 9 ) .', 'These observations motivated us to manually correct the inference from CNAnorm .', 'We adjusted the correspondence between peak centers and copy numbers and as a result of the correction , the R2 was substantially improved ( Fig 10D ) .', 'The tumor purity was re-estimated with the corrected correspondence between μi and CNi in Eq ( 5 ) , and revised from 64 . 92% to 80 . 88% .', 'The re-estimated copy number profile was improved in terms of better fit to integer copy numbers ( horizontal black segments overlap more with the horizontal gray lines in S13B Fig ) .', 'The correction also lead to strikingly better fit of the theoretical mBAF to the observed mBAF values ( Fig 11 ) .', 'The theoretical mBAF was calculated based on the estimated copy number and tumor purity ( detailed in S1K Text ) .', 'Interestingly , this tumor genome likely underwent a whole-genome doubling event [31] and was estimated to be tetraploid , based on which the baseline of “normal” copy number was anchored .', 'In summary , this case study illustrated the joint inference from BAF along with total copy number could result in more reliable estimation of tumor ploidy and purity .'], ['We have developed a joint segmentation and inference approach for SCNA analysis using both total and allele-specific sequencing depth and investigated its performance with a “null” dataset and a real-world large-scale dataset .', 'Compared with existing methods: ExomeCNV , PatternCNV , CNAnorm and Control-FREEC , our method demonstrated improved accuracy , scalability to large cancer sequencing studies , and the potential in improving the estimation of tumor ploidy and purity .', 'Our approach exhibits flexibility and applicability in a wide range of platforms , including both deep sequencing and SNP array .', 'These good properties can facilitate integrative cancer genomics study using multiple platforms .', 'An R package called saasCNV , which implements our proposed appoach , is avaiable at https://zhangz05 . u . hpc . mssm . edu/saasCNV/index . htm .', 'In contrast with germline CNV detection , characterization of SCNA in cancer genome gives rise to particular challenges , including ambiguous baseline due to aneuploidy and diluted signal pattern due to heterogeneity and normal cell contamination .', 'For this regard , allelic specific information adds valuable input .', 'Our method takes advantage of this information in both deep sequencing and SNP array data and achieves better capability in the correct identification of copy number baseline .', 'We have also demonstrated that incorporating BAF could substantially improve the inference of tumor ploidy and purity .', 'An important future work is the development of an integrated statistical model based on segment-level total sequencing depth and BAF to infer absolute copy number and tumor purity simultaneously .', 'SAAS-CNV features good efficiency and scalability by adopting the “top-down” strategy , which reduces the number of statistical inferences from millions to hundreds , and by taking advantage of condensed information from VCF file , which is about 1% of the size of BAM file .', 'On Dataset II , we showed WES of moderate sequencing depth can still provide comparable specificity and sensitivity as WGS in large SCNA detection .', 'At the same time , we acknowledge that quantifying SCNA at heterozygous loci has limitations in that not all information from mapped reads is fully utilized .', 'While this limitation has less influence in WGS , where SNPs and indels are densely distributed on the genome , it affects the resolution in a greater degree on WES data .', 'A possible improvement is to leverage BAMs to boost the SNR of log2ratio signal and the resolution by averaging over read depths within the windows surrounding or nearby heterozygous sites .', 'However , it is quite time-comsumig to manipulate large BAM files .', 'An alternative option is to utilize genomic VCF ( GVCF ) file , which is in the similar format as VCF , but record both variant sites and non-variant blocks .', 'Inspired by the normalization procedure of Illumina SNP array data [7] , the data processing step in our method can be further improved by incorparating multiple samples simultaneously , especially when manipulating GVCF makes it computationally feasible .', 'The matched tumor-normal design provides several desired features for the identification of SCNA:', '1 ) it is biologically sensible to take matched normal genome as reference to define somatic alterations in tumor genome;', '2 ) it is helpful to reduce the bias induced from the spatially non-uniform distribution of short reads across the genome , due to variability in GC content , exon capture efficiency , mappability of complex regions and so forth ( also see S1L Text ) ;', '3 ) it is able to improve the normality of signals ( S2 Fig ) , which is desirable for joint segmentation step;', '4 ) it can help alleviate allelic signal bias commonly observed in SNP array data ( S10 Fig ) [30] .', 'Lastly , our method can be used in accompany with paired-end mapping or split read methods , for example CREST [28] , to refine the resolution of break points up to base pair level and verify other types of genomic rearrangements , such as intra-chromosomal and inter-chromosomal translocations , which are commonly associated with SCNAs [32] .']]
['Introduction', 'Methods and Materials', 'Results', 'Discussion']
['Cancer genomes exhibit profound somatic copy number alterations ( SCNAs ) .', 'Studying tumor SCNAs using massively parallel sequencing provides unprecedented resolution and meanwhile gives rise to new challenges in data analysis , complicated by tumor aneuploidy and heterogeneity as well as normal cell contamination .', 'While the majority of read depth based methods utilize total sequencing depth alone for SCNA inference , the allele specific signals are undervalued .', 'We proposed a joint segmentation and inference approach using both signals to meet some of the challenges .', 'Our method consists of four major steps:', '1 ) extracting read depth supporting reference and alternative alleles at each SNP/Indel locus and comparing the total read depth and alternative allele proportion between tumor and matched normal sample;', '2 ) performing joint segmentation on the two signal dimensions;', '3 ) correcting the copy number baseline from which the SCNA state is determined;', '4 ) calling SCNA state for each segment based on both signal dimensions .', 'The method is applicable to whole exome/genome sequencing ( WES/WGS ) as well as SNP array data in a tumor-control study .', 'We applied the method to a dataset containing no SCNAs to test the specificity , created by pairing sequencing replicates of a single HapMap sample as normal/tumor pairs , as well as a large-scale WGS dataset consisting of 88 liver tumors along with adjacent normal tissues .', 'Compared with representative methods , our method demonstrated improved accuracy , scalability to large cancer studies , capability in handling both sequencing and SNP array data , and the potential to improve the estimation of tumor ploidy and purity .']
['Somatic copy number alterations ( SCNAs ) are essential in oncogensis and progression of a variety of cancers .', 'Accurate identification and quatification of SCNAs are fundamental in the effort of cataloging different variants in cancer genome .', 'This task has its own challenges due to complex nature of tumor SCNA profile and is further complicated by the heterogeneity of the cells collected from a tumor tissue and the contamination from adjacent normal cells , making it difficult for the methods well tailored for the detection of germline copy number variation ( CNV ) to fit in tumor SCNA detection .', 'Next generation sequencing provides an opportunity to comprehensively characterize SCNA at unprecedent resolution .', 'While total read depth information is commonly used in SCNA detection methods , the allele-specific read depth is less often considered , leading to sub-optimal solution .', 'By incorparating both pieces of information , we developed a segmentation-based pipeline to address aforementioned issues in SCNA detection .', 'This tool is applicable on both deep sequencing data as well as SNP array data and enables accurate and efficient characterization of genome-wide SCNA profile to facilitate large-scale cancer studies .']
[]