Last data update: 2014.03.03

R: TITAN: Subclonal copy number and LOH prediction whole genome...
TitanCNA-packageR Documentation

TITAN: Subclonal copy number and LOH prediction whole genome sequencing of tumours

Description

TITAN is a software tool for inferring subclonal copy number alterations (CNA) and loss of heterozygosity (LOH). The algorithm also infers clonal group cluster membership for each event and the tumour proportion, or cellular prevalence, for each event.

Details

Package: TitanCNA
Type: Package
Version: 1.9.0
Date: 2016-02-17
License: see LICENSE

example("TitanCNA-package") for quick tour of functionality and visualization

vignette("TitanCNA") for detailed example

Author(s)

Gavin Ha, Sohrab P Shah Maintainer: Gavin Ha <gavinha@broadinstitute.org>

References

Ha, G., Roth, A., Khattra, J., Ho, J., Yap, D., Prentice, L. M., Melnyk, N., McPherson, A., Bashashati, A., Laks, E., Biele, J., Ding, J., Le, A., Rosner, J., Shumansky, K., Marra, M. A., Huntsman, D. G., McAlpine, J. N., Aparicio, S. A. J. R., and Shah, S. P. (2014). TITAN: Inference of copy number architectures in clonal cell populations from tumour whole genome sequence data. Genome Research, 24: 1881-1893. (PMID: 25060187)

Examples

message('Running TITAN ...')
#### LOAD DATA ####
infile <- system.file("extdata", "test_alleleCounts_chr2.txt", package = "TitanCNA")
data <- loadAlleleCounts(infile)

#### LOAD PARAMETERS ####
message('titan: Loading default parameters')
numClusters <- 2
params <- loadDefaultParameters(copyNumber = 5, 
                                numberClonalClusters = numClusters, skew = 0.1)

#### READ COPY NUMBER FROM HMMCOPY FILE ####
message('titan: Correcting GC content and mappability biases...')
tumWig <- system.file("extdata", "test_tum_chr2.wig", package = "TitanCNA")
normWig <- system.file("extdata", "test_norm_chr2.wig", package = "TitanCNA")
gc <- system.file("extdata", "gc_chr2.wig", package = "TitanCNA")
map <- system.file("extdata", "map_chr2.wig", package = "TitanCNA")
cnData <- correctReadDepth(tumWig, normWig, gc, map)
logR <- getPositionOverlap(data$chr, data$posn, cnData)
data$logR <- log(2^logR) #transform to natural log

#### FILTER DATA FOR DEPTH, MAPPABILITY, NA, etc ####
data <- filterData(data, c(1:22,"X"), minDepth = 10, maxDepth = 200, map = NULL)

#### EM (FWD-BACK) TO TRAIN PARAMETERS ####
#### Can use parallelization packages ####
K <- length(params$genotypeParams$alphaKHyper)
params$genotypeParams$alphaKHyper <- rep(500, K)
params$ploidyParams$phi_0 <- 1.5 
convergeParams <- runEMclonalCN(data, gParams = params$genotypeParams, 
                                nParams = params$normalParams, 
                                pParams = params$ploidyParams, 
                                sParams = params$cellPrevParams, 
                                maxiter = 3, maxiterUpdate = 500, 
                                txnExpLen = 1e9, txnZstrength = 1e9, 
                                useOutlierState = FALSE, 
                                normalEstimateMethod = "map", 
                                estimateS = TRUE, estimatePloidy = TRUE)                                
#### COMPUTE OPTIMAL STATE PATH USING VITERBI ####
optimalPath <- viterbiClonalCN(data, convergeParams)

#### FORMAT RESULTS ####
results <- outputTitanResults(data, convergeParams, optimalPath, 
                              filename = NULL, posteriorProbs = FALSE,
                              subcloneProfiles = TRUE)
                              
#### REMOVE EMPTY CLUSTERS ####
corrResults <- removeEmptyClusters(convergeParams, results, proportionThreshold = 0.001,
																		proportionThresholdClonal = 0.3)
convergeParams <- corrResults$convergeParams
results <- corrResults$results

#### PLOT RESULTS ####
norm <- tail(convergeParams$n, 1)
ploidy <- tail(convergeParams$phi, 1)

par(mfrow=c(4, 1))    
plotCNlogRByChr(results, chr = 2, ploidy = ploidy, normal = norm, geneAnnot = NULL, 
                ylim = c(-2, 2), cex = 0.5, xlab = "", main = "Chr 2")
plotAllelicRatio(results, chr = 2, geneAnnot = NULL, ylim = c(0, 1), cex = 0.5, 
                xlab = "", main = "Chr 2")
plotClonalFrequency(results, chr = 2, normal = norm, geneAnnot = NULL, 
                    ylim = c(0, 1), cex = 0.5, xlab = "", main = "Chr 2")
plotSubcloneProfiles(results, chr = 2, cex = 2, main = "Chr 2")

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(TitanCNA)
Loading required package: foreach
Loading required package: IRanges
Loading required package: BiocGenerics
Loading required package: parallel

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:parallel':

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following objects are masked from 'package:stats':

    IQR, mad, xtabs

The following objects are masked from 'package:base':

    Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    as.data.frame, cbind, colnames, do.call, duplicated, eval, evalq,
    get, grep, grepl, intersect, is.unsorted, lapply, lengths, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, rank,
    rbind, rownames, sapply, setdiff, sort, table, tapply, union,
    unique, unsplit

Loading required package: S4Vectors
Loading required package: stats4

Attaching package: 'S4Vectors'

The following objects are masked from 'package:base':

    colMeans, colSums, expand.grid, rowMeans, rowSums

Loading required package: GenomicRanges
Loading required package: GenomeInfoDb
Loading required package: Rsamtools
Loading required package: Biostrings
Loading required package: XVector
> png(filename="/home/ddbj/snapshot/RGM3/R_BC/result/TitanCNA/TitanCNA-package.Rd_%03d_medium.png", width=480, height=480)
> ### Name: TitanCNA-package
> ### Title: TITAN: Subclonal copy number and LOH prediction whole genome
> ###   sequencing of tumours
> ### Aliases: TitanCNA-package TitanCNA
> ### Keywords: IO manip package
> 
> ### ** Examples
> 
> message('Running TITAN ...')
Running TITAN ...
> #### LOAD DATA ####
> infile <- system.file("extdata", "test_alleleCounts_chr2.txt", package = "TitanCNA")
> data <- loadAlleleCounts(infile)
titan: Loading data /home/ddbj/local/lib64/R/library/TitanCNA/extdata/test_alleleCounts_chr2.txt
> 
> #### LOAD PARAMETERS ####
> message('titan: Loading default parameters')
titan: Loading default parameters
> numClusters <- 2
> params <- loadDefaultParameters(copyNumber = 5, 
+                                 numberClonalClusters = numClusters, skew = 0.1)
> 
> #### READ COPY NUMBER FROM HMMCOPY FILE ####
> message('titan: Correcting GC content and mappability biases...')
titan: Correcting GC content and mappability biases...
> tumWig <- system.file("extdata", "test_tum_chr2.wig", package = "TitanCNA")
> normWig <- system.file("extdata", "test_norm_chr2.wig", package = "TitanCNA")
> gc <- system.file("extdata", "gc_chr2.wig", package = "TitanCNA")
> map <- system.file("extdata", "map_chr2.wig", package = "TitanCNA")
> cnData <- correctReadDepth(tumWig, normWig, gc, map)
Reading GC and mappability files
Slurping: /home/ddbj/local/lib64/R/library/TitanCNA/extdata/gc_chr2.wig
Parsing: fixedStep chrom=2 start=1 step=1000 span=1000
Sorting by decreasing chromosome size
Slurping: /home/ddbj/local/lib64/R/library/TitanCNA/extdata/map_chr2.wig
Parsing: fixedStep chrom=2 start=1 step=1000 span=1000
Sorting by decreasing chromosome size
Loading tumour file:/home/ddbj/local/lib64/R/library/TitanCNA/extdata/test_tum_chr2.wig
Slurping: /home/ddbj/local/lib64/R/library/TitanCNA/extdata/test_tum_chr2.wig
Parsing: fixedStep chrom=2 start=1 step=1000 span=1000
Sorting by decreasing chromosome size
Loading normal file:/home/ddbj/local/lib64/R/library/TitanCNA/extdata/test_norm_chr2.wig
Slurping: /home/ddbj/local/lib64/R/library/TitanCNA/extdata/test_norm_chr2.wig
Parsing: fixedStep chrom=2 start=1 step=1000 span=1000
Sorting by decreasing chromosome size
Correcting Tumour
Applying filter on data...
Correcting for GC bias...
Correcting for mappability bias...
Correcting Normal
Applying filter on data...
Correcting for GC bias...
Correcting for mappability bias...
Normalizing Tumour by Normal
> logR <- getPositionOverlap(data$chr, data$posn, cnData)
> data$logR <- log(2^logR) #transform to natural log
> 
> #### FILTER DATA FOR DEPTH, MAPPABILITY, NA, etc ####
> data <- filterData(data, c(1:22,"X"), minDepth = 10, maxDepth = 200, map = NULL)
> 
> #### EM (FWD-BACK) TO TRAIN PARAMETERS ####
> #### Can use parallelization packages ####
> K <- length(params$genotypeParams$alphaKHyper)
> params$genotypeParams$alphaKHyper <- rep(500, K)
> params$ploidyParams$phi_0 <- 1.5 
> convergeParams <- runEMclonalCN(data, gParams = params$genotypeParams, 
+                                 nParams = params$normalParams, 
+                                 pParams = params$ploidyParams, 
+                                 sParams = params$cellPrevParams, 
+                                 maxiter = 3, maxiterUpdate = 500, 
+                                 txnExpLen = 1e9, txnZstrength = 1e9, 
+                                 useOutlierState = FALSE, 
+                                 normalEstimateMethod = "map", 
+                                 estimateS = TRUE, estimatePloidy = TRUE)                                
titan: Running HMM...
fwdBack: Iteration 1 chr: 1 
Coord Descent[1]=-31433 n=0.5 s=[0.001,0.200] phi=1.5
Coord Descent[2]=-27411 n=0.5232 s=[0.00022,0.37516] phi=1.463
Coord Descent[3]=-26067 n=0.4002 s=[0.00019,0.44657] phi=1.516
Coord Descent[4]=-25747 n=0.3577 s=[0.00032,0.49378] phi=1.54
Coord Descent[5]=-25662 n=0.3349 s=[0.00064,0.51602] phi=1.553
Coord Descent[6]=-25643 n=0.3259 s=[0.0013,0.5272] phi=1.56
Coord Descent[7]=-25638 n=0.3218 s=[0.0022,0.5327] phi=1.563
Coord Descent[8]=-25636 n=0.3194 s=[0.0035,0.5357] phi=1.565
Coord Descent[9]=-25636 n=0.3178 s=[0.0048,0.5376] phi=1.567
Coord Descent[10]=-25635 n=0.3166 s=[0.0061,0.5388] phi=1.568
Coord Descent[11]=-25635 n=0.3155 s=[0.0073,0.5398] phi=1.569
Using Coordinate Descent iteration 11 with Fval=-25635 and n=0.3155 (map), s=[0.0073,0.5398], phi=1.569
fwdBack: loglik=-31916.2307
fwdBack: priorN=0.2591
fwdBack: priorS=-2.7370
fwdBack: priorVar=-807.8681
fwdBack: priorPhi=-0.8158
fwdBack: priorPiG=23.6845
fwdBack: priorPiZ=0.0189
fwdBack: EM iteration 1 complete loglik=-32703.6890
fwdBack: Elapsed time for iteration 1: 0.0260m
fwdBack: Iteration 2 chr: 1 
Coord Descent[1]=-24741 n=0.3155 s=[0.0073,0.5398] phi=1.569
Coord Descent[2]=-24569 n=0.2738 s=[0.0031,0.5561] phi=1.571
Coord Descent[3]=-24568 n=0.2703 s=[0.0065,0.5578] phi=1.57
Coord Descent[4]=-24567 n=0.2674 s=[0.0087,0.5577] phi=1.57
Coord Descent[5]=-24567 n=0.2661 s=[0.01,0.56] phi=1.57
Coord Descent[6]=-24566 n=0.2651 s=[0.011,0.558] phi=1.57
Coord Descent[7]=-24566 n=0.2644 s=[0.012,0.559] phi=1.571
Coord Descent[8]=-24566 n=0.2639 s=[0.013,0.559] phi=1.571
Coord Descent[9]=-24566 n=0.2633 s=[0.013,0.559] phi=1.571
Coord Descent[10]=-24566 n=0.2628 s=[0.014,0.560] phi=1.571
Coord Descent[11]=-24566 n=0.2623 s=[0.014,0.560] phi=1.572
Using Coordinate Descent iteration 11 with Fval=-24566 and n=0.2623 (map), s=[0.014,0.560], phi=1.572
fwdBack: loglik=-24456.0513
fwdBack: priorN=0.1494
fwdBack: priorS=-2.0717
fwdBack: priorVar=-930.7287
fwdBack: priorPhi=-0.8046
fwdBack: priorPiG=23.6845
fwdBack: priorPiZ=-0.1211
fwdBack: EM iteration 2 complete loglik=-25365.9435
fwdBack: Elapsed time for iteration 2: 0.0244m
fwdBack: Iteration 3 chr: 1 
Coord Descent[1]=-24590 n=0.2623 s=[0.014,0.560] phi=1.572
Coord Descent[2]=-24586 n=0.2557 s=[0.014,0.562] phi=1.572
Coord Descent[3]=-24586 n=0.2551 s=[0.015,0.562] phi=1.572
Coord Descent[4]=-24585 n=0.2544 s=[0.016,0.562] phi=1.572
Coord Descent[5]=-24585 n=0.2539 s=[0.016,0.562] phi=1.573
Coord Descent[6]=-24585 n=0.2535 s=[0.017,0.563] phi=1.573
Coord Descent[7]=-24585 n=0.2531 s=[0.017,0.563] phi=1.573
Coord Descent[8]=-24585 n=0.2527 s=[0.018,0.563] phi=1.573
Coord Descent[9]=-24585 n=0.2523 s=[0.018,0.563] phi=1.574
Coord Descent[10]=-24585 n=0.252 s=[0.019,0.563] phi=1.574
Coord Descent[11]=-24585 n=0.2516 s=[0.019,0.564] phi=1.574
Using Coordinate Descent iteration 11 with Fval=-24585 and n=0.2516 (map), s=[0.019,0.564], phi=1.574
fwdBack: loglik=-24132.7787
fwdBack: priorN=0.1221
fwdBack: priorS=-1.7853
fwdBack: priorVar=-935.8662
fwdBack: priorPhi=-0.7964
fwdBack: priorPiG=23.6845
fwdBack: priorPiZ=-0.1343
fwdBack: EM iteration 3 complete loglik=-25047.5543
fwdBack: Elapsed time for iteration 3: 0.0243m
fwdBack: Total elapsed time: 0.0907m
Warning message:
executing %dopar% sequentially: no parallel backend registered 
> #### COMPUTE OPTIMAL STATE PATH USING VITERBI ####
> optimalPath <- viterbiClonalCN(data, convergeParams)
> 
> #### FORMAT RESULTS ####
> results <- outputTitanResults(data, convergeParams, optimalPath, 
+                               filename = NULL, posteriorProbs = FALSE,
+                               subcloneProfiles = TRUE)
>                               
> #### REMOVE EMPTY CLUSTERS ####
> corrResults <- removeEmptyClusters(convergeParams, results, proportionThreshold = 0.001,
+ 																		proportionThresholdClonal = 0.3)
> convergeParams <- corrResults$convergeParams
> results <- corrResults$results
> 
> #### PLOT RESULTS ####
> norm <- tail(convergeParams$n, 1)
> ploidy <- tail(convergeParams$phi, 1)
> 
> par(mfrow=c(4, 1))    
> plotCNlogRByChr(results, chr = 2, ploidy = ploidy, normal = norm, geneAnnot = NULL, 
+                 ylim = c(-2, 2), cex = 0.5, xlab = "", main = "Chr 2")
> plotAllelicRatio(results, chr = 2, geneAnnot = NULL, ylim = c(0, 1), cex = 0.5, 
+                 xlab = "", main = "Chr 2")
> plotClonalFrequency(results, chr = 2, normal = norm, geneAnnot = NULL, 
+                     ylim = c(0, 1), cex = 0.5, xlab = "", main = "Chr 2")
> plotSubcloneProfiles(results, chr = 2, cex = 2, main = "Chr 2")
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>