Identifying CNV events

The expression values should be normalized using TPM or FPKM.

library(CaSpER)
data (yale_meningioma)
data <- yale_meningioma$data
loh <-  yale_meningioma$loh
loh.name.mapping <-  yale_meningioma$loh.name.mapping
control.sample.ids <-  yale_meningioma$control.sample.ids
genoMat <-  yale_meningioma$genoMat
samps <-  yale_meningioma$samps
kable(yale_meningioma$data[1:5, 1:5])
MN-1171 MN-1237 MN-1137 MN-82152 MN-60835
ENSG00000223972 4.864083 3.543738 2.413422 5.058923 4.002694
ENSG00000227232 10.179492 9.386185 8.895256 9.054355 10.149439
ENSG00000243485 1.087406 1.127410 1.846561 1.945983 3.089998
ENSG00000237613 3.321832 2.520711 1.846561 4.019347 4.002694
ENSG00000268020 0.000000 1.127410 0.000000 0.000000 0.000000

Normalized gene expression matrix rows represents genes either with ensembl gene id or hgnc gene symbol. Columns represent samples or cells.

Cytoband information can be downloaded from UCSC.

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/cytoBand.txt.gz
data(hg19_cytoband)
kable(cytoband[1:5, ])
V1 V2 V3 V4
1 1.00e+00 125000001 p
1 1.25e+08 249250622 q
2 1.00e+00 93300001 p
2 9.33e+07 243199374 q
3 1.00e+00 91000001 p

Centromere information can be downloaded from UCSC.

#curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/cytoBand.txt.gz" | gunzip -c | grep acen | head

Annotation data which contains the positions of each gene along each chromosome in the genome can be generated using the following functions:

annotation <- generateAnnotation(id_type="ensembl_gene_id", genes=rownames(yale_meningioma$data), ishg19=T, centromere)
data <- data[match( annotation$Gene,rownames(data)), ]
kable(annotation[1:5, ])
Gene GeneSymbol EntrezID Chr start end band cytoband isCentromer Position new_positions
266 ENSG00000223972 DDX11L1 100287596 1 11869 14412 p36.33 1p no 13140.5 1
302 ENSG00000227232 WASH7P 100287171 1 14363 29806 p36.33 1p no 22084.5 2
457 ENSG00000243485 MIR1302-11 100422919 1 29554 31109 p36.33 1p no 30331.5 3
415 ENSG00000237613 FAM138A 654835 1 34554 36081 p36.33 1p no 35317.5 4
514 ENSG00000268020 OR4G4P NA 1 52473 54936 p36.33 1p no 53704.5 5

BAFExtract output files are read using the following functions:

loh <- readBAFExtractOutput ( path="./meningioma_baf\\", sequencing.type="bulk")
names(loh) <- gsub(".snp", "", names(loh))

The main functions you will need to use to identify CNV events are CreateCasperObject() and runCaSpER(casper_object).

1. Creating casper object

The casper object is required for performing CNV analysis on single-cell and bulk RNA-Seq. Casper object is created using following functions

object <- CreateCasperObject(raw.data=data,loh.name.mapping=loh.name.mapping, sequencing.type="bulk", 
  cnv.scale=3, loh.scale=3,
  annotation=annotation, method="iterative", loh=loh, 
  control.sample.ids=control.sample.ids, cytoband=cytoband)

Above CaSpER object stores all information associated with the dataset, including expression data, smoothed expression data in n (default:3) different scales, original baf signal, smoothed baf signal, annotations, control sample ids.

loh.name.mapping: contains the sample id mappings between expression signal and baf signal. For bulk sequencing data there will be one to one mapping since baf signal is generated for all samples seperately.

kable(yale_meningioma$loh.name.mapping[1:5, ])
loh.name sample.name
MN-60835 MN-60835
MN-60 MN-60
MN-5 MN-5
MN-1137 MN-1137
MN-1161 MN-1161

2. Pairwise comparison of scales from BAF and expression signals

After creating the casper object we perform HMM on all scales of expression signal and integrate HMM segment states with Allele Frequency Shift Information. CaSpER algorithm outputs CNV calls using all the pairwise comparisons of expression signal and BAF signal decompositions.

final.objects <- runCaSpER(object, removeCentromere=T, cytoband=cytoband, method="iterative")

3. Harmonization and Summarization of CNV calls from multiple scales and from multiple pairwise comparison of BAF and Expression Signals

The pairwise comparison and assignment of CNV calls generates a large number of per-scale information that must be summarized such that each position of the genome is assigned a final call about its CNV status, i.e., deletion/amplification/neutral. We use a consistency-based approach for harmonizing the pairwise comparisons: The events are put together and we assign the final CNV for a gene or large-scale event if the CNV calls are consistent among at least a certain number of pairwise scale comparisons. CaSpER harmonizes and summarizes the CNV calls by dividing them into large-scale, gene-based, and segment-based CNV calls

  1. Large-Scale CNV Summarization.

We assign a large-scale CNV call to every chromosome arm for each of the N×M pairwise scale comparisons. Next, for each chromosome arm, we ask whether the large-scale CNV call is consistent among at least \(\gamma\) of the N×M large-scale CNV calls. N denotes the index for the highest smoothing scale for expression signal. M denotes the index for the highest smoothing scale for baf signal. thr represents minimum percentage, 75% (at least 7 out of 9 scales), of consistent CNV calls (Out of N×M comparisons of expression scales and BAF scales) while assigning a final CNV (amp/del/neutral) call to a segment/gene/chromosome arm.

finalChrMat <- extractLargeScaleEvents (final.objects, thr=0.75)
  1. Segment based CNV Summarization.

The segments-based summarization aims at generating a final set of CNV calls for a final set of segments that are computed by comparison of scales. We first compare the segments from different expression scales and generate the consistent set of segments. For each segment in the final set, if there are more than \(\gamma\) (default=6) consistent CNV calls among N×M CNV calls, we assign the consistent CNV call to segment. When there is no consistency among the calls, we assign a neutral CNV state to segment.

gamma <- 7
all.segments <- do.call(rbind, lapply(final.objects, function(x) x@segments))
segment.summary <- extractSegmentSummary (final.objects)
loss <- segment.summary$all.summary.loss
gain <- segment.summary$all.summary.gain
loss.final <- loss[loss$count>=gamma, ]
gain.final <- gain[gain$count>=gamma, ]

Calculate large scale event accuracy using snp genotyping array as gold standard.

data("yale_meningioma")
genoMat <-  yale_meningioma$genoMat
common <- intersect(order.sampleNames, intersect(rownames(finalChrMat), rownames(genoMat)))
finalChrMat <- finalChrMat[match(common, rownames(finalChrMat)), ]
genoMat <- genoMat[match(common, rownames(genoMat)), ]
calcROC(chrMat=finalChrMat, chrMat2=genoMat)
  1. Gene based CNV Summarization.

Similar to the large-scale summarization, we generate a matrix where rows are the samples (cells) and columns are the genes. The matrix entry of 0 corresponds to no alteration, 1 corresponds to amplification and -1 corresponds to deletion. If an alteration is consistent in more than \(\gamma\) scale comparisons (out of N×M comparisons), we report that alteration event for that sample.

all.summary<- rbind(loss.final, gain.final)
colnames(all.summary) [2:4] <- c("Chromosome", "Start",   "End")
rna <-  GRanges(seqnames = Rle(gsub("q", "", gsub("p", "", all.summary$Chromosome))), 
    IRanges(all.summary$Start, all.summary$End))   
ann.gr <- makeGRangesFromDataFrame(final.objects[[1]]@annotation.filt, keep.extra.columns = TRUE, seqnames.field="Chr")
hits <- findOverlaps(geno.rna, ann.gr)
genes <- splitByOverlap(ann.gr, geno.rna, "GeneSymbol")
genes.ann <- lapply(genes, function(x) x[!(x=="")])
all.genes <- unique(final.objects[[1]]@annotation.filt[,2])
all.samples <- unique(as.character(final.objects[[1]]@segments$ID))
rna.matrix <- gene.matrix(seg=all.summary, all.genes=all.genes, all.samples=all.samples, genes.ann=genes.ann)

Calculate gene based TPR and FPR using genotyping array as gold standard

segments <- yale_meningioma$segments
segments$type<- segments$Type_Corrected
geno.gr <-  GRanges(seqnames = Rle(gsub("q", "", gsub("p", "", segments$chr))), IRanges(segments$start, segments$end))   
hits <- findOverlaps(geno.gr, ann.gr)
genes <- splitByOverlap(ann.gr, geno.gr, "GeneSymbol")
genes.ann <- lapply(genes, function(x) x[!(x=="")])
gt.matrix <- gene.matrix(seg=segments, all.samples=all.samples, all.genes=all.genes, genes.ann=genes.ann)
calcROC(chrMat=rna.matrix, chrMat2=gt.matrix)

4.Visualization

In order to comprehensively visualize the CNV events, CaSpER first generates visuals that combine all the samples so that user can view the events jointly over all samples and co-occurrence and mutually exclusive patterns can be visually inspected. To characterize these patterns, CaSpER performs clustering of samples and clustering of detected CNVs and generates clustering plots. We summarize these visualizations below.

  1. plotHeatmap: Visualization of the genomewide gene expression signal plot at different smoothing scales CaSpER outputs the expression signal at different scales. In these plots, each row is a sample and the columns are the chromosomes. These can be useful for comparison of outputs from different scales using the panoramic inspection of the expression signal.
obj <- final.objects[[9]]
plotHeatmap(object=obj, fileName="heatmap.png",cnv.scale= 3, cluster_cols = F, cluster_rows = T, show_rownames = T, only_soi = T)

  1. plotLargeScaleEvent: Visualization of the large-scale CNV events among all the samples/cells. Large scale event summarization is useful for summarizing the detected large-scale CNV events (deletions and amplifications) over multiple samples. This plot summarizes the large scale CNVs and may reveal the patterns that may otherwise be missed when data is visualized at smaller scales.
plotLargeScaleEvent (object=obj, fileName="large.scale.events.png")

  1. plotGEAndGT: Plot large scale events called from genotyping and RNA-Seq (can be used only with small sample size)
plotGEAndGT (chrMat=finalChrMat, genoMat=genoMat, fileName="RNASeqAndGT.png")

  1. plotBAFAllSamples: Visualization of BAF shift signal for all samples together. The inspection of BAF shift signal is useful especially when compared to the expression signal to analyze the CNV and LOH events. The BAF shift plots show the BAF shift signal such that each row is the genomewide BAF shift signal profile.
plotBAFAllSamples (loh = obj@loh.median.filtered.data,  fileName="LOHAllSamples.png")

  1. plotBAFOneSample: Visualization of BAF shift signal in different scales for one sample. This option plots the BAF shift signal for one sample at different scales. Similar to the multiscale smoothing of expression signal, this information enables panoramic assessment and identification of CNV and LOH events.
plotBAFOneSample (object, fileName="LOHPlotsAllScales.pdf")

  1. plotBAFInSeperatePages: Visualization of BAF deviation for each sample in separate pages. This option creates the BAF shift signal for each sample separately. This way the user can visualize each sample by itself.
plotBAFInSeperatePages (loh=obj@loh.median.filtered.data, folderName="LOHPlots")

  1. plotGEAndBAFOneSample: Gene expression and BAF signal for one sample in one plot. This option generates the visualization of the gene expression signal and BAF shift signal together so that the user can jointly assess them in one page.
plotGEAndBAFOneSample (object=obj, cnv.scale=3, loh.scale=3, sample= "MN-5")