Introduction

To evaluate the aneuploidy and prevalence of clonal or quasiclonal tumors, we developed a novel tool to characterize the mosaic tumor genome on the basis of one major assumption: the genome of a heterogeneous multi-cell tumor biopsy can be sliced into a chain of segments that are characterized by homogeneous somatic copy number alternations (SCNAs) and B allele frequencies (BAFs). The model, termed BubbleTree, utilizes both SCNAs and the interconnected BAFs as markers of tumor clones to extract tumor clonality estimates. BubbleTree is an intuitive and powerful approach to jointly identify ASCN, tumor purity and (sub)clonality, which aims to improve upon current methods to characterize the tumor karyotypes and ultimately better inform cancer diagnosis, prognosis and treatment decisions.

Quickstart to Using BubbleTree

To perform a BubbleTree analysis, data pertaining to the position and B allele frequency of heterozygous snps in the tumor sample and segmented copy number information including the position, number of markers/segment and log2 copy number ratio between tumor and normal samples must first be obtained. The section “Preparing Data for BubbleTree” below demonstrates just one way of generating the proper input files from next generation sequencing (NGS) workflows. Example data in the desired format is provided as part of this package as GRanges objects and can be loaded as follows.

library(BubbleTree)
data(hetero.gr) #loads sequence variants
data(cnv.gr) #loads copy number variation data

The sequence variation and copy number data is then combined into a data frame object while simutaneously adding cytoband annotation. This function also calculates Heterozygosity-Deviation Scores (HDS) for each variant and summarizes mean HDS and HDS Standard Deviation/ CNV segment.

rbd=getRBD(hetero.gr,cnv.gr)
## Segments with high SD:
## [1] seg.id     hds.median hds.sd     num.mark   seg.mean   chr       
## [7] start      end        cyto.band 
## <0 rows> (or 0-length row.names)

Finally a BubbleTree diagram is created using the rbd data frame.

plotBubbles(rbd)

Preparing Data for BubbleTree

BubbleTree was developed using both whole exome sequencing (WES) and whole genome sequening (WGS) NGS data from paired tumor/normal biopsies, but this model should also be applicable to array comparative genomic hybridization (aCGH) and single nucleotide polymorphism (SNP) array data.

There many methods for generating and processing sequencing data in preparation for BubbleTree analysis. In the following section we provide example workflows starting from WES NGS which can be adapted as needed to alternate inputs.

Preparing Sequence Variation Data

The primary BubbleTree requirement for sequence variant information is a GRanges object containing the B alelle frequencies and genomic positions of variants known to be heterozygous in the paired normal sample.

Mapped reads from tumor and normal tissue can be processed with mutation caller software such as VARSCAN or MUTECT. In this example, we use a hypothetical vcf file output which contains mutation calls from both normal and tumor samples.

The B-allele frequency data is extracted using the Bioconductor package VariantAnnotation and converted from string to numeric format.

library(VariantAnnotation)
vc=readVcf("MyN_TSamples.vcf",genome="hg19")
fq=geno(vc)$FREQ
freq <- data.frame(fq)
freq[,] <- as.numeric(gsub("%", "", as.matrix(freq[,])))/100
colnames(freq)=paste(colnames(freq),"freq",sep=".")

As a QC criteria, we use only variants that meet a certain read depth. Depth information can also be extracted from the vcf file and combined with B-allele frequency data and snp position data as follows.

dp=geno(vc)$DP
colnames(dp)=paste(colnames(dp),"dp",sep=".")
#combine all with chr and position info
snp.dat=data.frame("CHROM"=as.vector(seqnames(vc)),
                   "POS"=start(vc),freq,dp)

We use a small function to select heterozygous mutations defined as having a B allele frequency between 0.4 and 0.6, and use the function to screen heterozygous calls in the normal sample with a read depth >=15.

is.hetero <- function(x, a=0.4, b=0.6) {
  (x - a)  *  (b - x) >= 0
} 

snp.ss=subset(snp.dat, ! CHROM %in% c("chrX", "chrY") & normal.dp >= 15 &  is.hetero(normal.freq, 0.4, 0.6))

Finally, a GRanges object containing the tumor B-allele frequency of mutations shown to be heterozygous in normal is created. BubbleTree functions expect tumor allele frequencies to be in a Metadata column labeled “freq”.

library(GRanges)
snp.gr <- GRanges(snp.ss$CHROM, IRanges(snp.ss$POS, snp.ss$POS), mcols=snp.ss[,"tumor.freq"])
names(elementMetadata(snp.gr))[grep(".freq",names(elementMetadata(snp.gr)))]<-"freq"

Preparing Copy Number Variation Data

BubbleTree requires segmented copy number information for sample analysis. Specifically, the genomic position, the number of markers and the mean log2 tumor/normal copy number ratio for each segment.

Alternative ways to generate this data exist, but for our WES NGS example we use the ExomeCNV.R package workflow as documented at this link; https://secure.genome.ucla.edu/index.php/ExomeCNV_User_Guide

The file demo.eCNV created in the above referenced workflow contains log2 ratios for each exon (prior to segmentation). We smooth and perform segmentation using the Bioconductor package DNAcopy.

library(DNAcopy)
#create a CNA object
CNA.object <- CNA(demo.eCNV$logR, demo.eCNV$chr,
                  demo.eCNV$probe_end, data.type = "logratio", sampleid = "test")
#smooth
smoothed.CNA.object <- smooth.CNA(CNA.object)
#segment
seg=segment(smoothed.CNA.object)

Next, a GRanges object is created containing the required information with a restriction that each segment contain at least 10 markers. The Metadata columns containing the number of markers/segment and the mean log2 ratio/segment must be named “num.mark” and “seg.mean” respectively.

library(GenomicRanges)
min.num <- 10
cnv.gr <- with(subset(seg$output, num.mark >= min.num & ! chrom %in% c("chrX", "chrY")) , GRanges(chrom, IRanges(loc.start, loc.end), mcols=cbind(num.mark, seg.mean)))

Main Bubbletree Functions

BubbleTree model and diagram

BubbleTree is a model based on three valid assumptions: 1) the paired normal specimen expresses the common diploid state, 2) variant sites are bi-allelic, and 3) genome segments (rather than the whole genome) with homogeneous copy number ratio and BAFs, exist in the profiled tumor genome. The first two assumptions generally hold, whereas the last homogeneity assumption can also be satisfied even in the case of a complex tumor clonal structure.

As the three assumptions are all generally plausible, we therefore developed a model for the BubbleTree diagram. For one homogenous genomic segment (x:y;p), we have,

Expected copy number, (CN)=2×(1-p)+(x+y)×p

Copy Ratio, R=(CN)/2=(1-p)+(x+y)/2×p (1)

B allele frequency, BAF=(y×p+1×(1-p))/((x+y)×p+2×(1-p))

and the homozygous-deviation score (HDS),

HDS= ∣BAF-0.5∣=(p×∣y-x∣)/(2×[(x+y)×p+2×(1-p)]) (2)

Based on equations (1) and (2), we are able to calculate an R score (copy ratio) and HDS for a segment (x:y; p). For example, (0:1; 0.75) will provide 0.625 and 0.3 for the R scores and HDS, respectively.

These calculations are performed by the getRBD() function

library(BubbleTree)
data(hetero.gr) 
data(cnv.gr)
rbd=getRBD(snp.gr=hetero.gr,cnv.gr=cnv.gr)
## Segments with high SD:
## [1] seg.id     hds.median hds.sd     num.mark   seg.mean   chr       
## [7] start      end        cyto.band 
## <0 rows> (or 0-length row.names)
head(rbd)
##   seg.id hds.median     hds.sd num.mark seg.mean  chr     start       end
## 1      1    0.12270 0.07103365      553   0.5339 chr1    762098   3440694
## 2      2    0.13615 0.05958364     5351   0.4673 chr1   3447666  49118903
## 3      3    0.29650 0.06771401     4272  -0.6061 chr1  49128694 121310027
## 4      8    0.25810 0.04683305     1188   1.0343 chr1 150418684 155954980
## 5      9    0.02815 0.04344431     6934   0.0513 chr1 155979184 247835485
## 6     10    0.28095 0.05170589       65  -0.6147 chr1 247875098 249210700
##   cyto.band
## 1    p36.33
## 2    p36.32
## 3       p33
## 4     q21.3
## 5       q22
## 6       q44

A plot of R score and HDS at various ploidy states forms the branches of a BubbleTree plot which can be generated as follows. Normally, this function is called internally by the plotBubbles() function.

drawBranches()