R Under development (unstable) (2024-02-16 r85931 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 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("PSCBS") PSCBS v0.67.0 successfully loaded. See ?PSCBS for help. > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Load SNP microarray data > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > data <- PSCBS::exampleData("paired.chr01") > str(data) 'data.frame': 73346 obs. of 6 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : int 1145994 2224111 2319424 2543484 2926730 2941694 3084986 3155127 3292731 3695086 ... $ CT : num 1.625 1.071 1.406 1.18 0.856 ... $ betaT : num 0.757 0.771 0.834 0.778 0.229 ... $ CN : num 2.36 2.13 2.59 1.93 1.71 ... $ betaN : num 0.827 0.875 0.887 0.884 0.103 ... > > # Drop single-locus outliers > dataS <- dropSegmentationOutliers(data) > > # Run light-weight tests > # Use only every 5th data point > dataS <- dataS[seq(from=1, to=nrow(data), by=5),] > # Number of segments (for assertion) > nSegs <- 3L > # Number of bootstrap samples (see below) > B <- 100L > > str(dataS) 'data.frame': 14670 obs. of 6 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : int 1145994 2941694 3710825 4240737 4276892 4464544 4714611 5095111 5034491 5158416 ... $ CT : num 1.63 1.4 1.41 1.17 1.16 ... $ betaT : num 0.7574 0.169 0.2357 0.2604 0.0576 ... $ CN : num 2.36 2.13 2.26 2.01 2.32 ... $ betaN : num 0.8274 0.5072 0.1671 0.1609 0.0421 ... > R.oo::attachLocally(dataS) > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Calculate DH > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > muN <- aroma.light::callNaiveGenotypes(betaN, censorAt=c(0,1)) > # SNPs are identifies as those loci that have non-missing 'betaT' & 'muN' > isSnp <- (!is.na(betaT) & !is.na(muN)) > isHet <- isSnp & (muN == 1/2) > rho <- rep(NA_real_, length=length(muN)) > rho[isHet] <- 2*abs(betaT[isHet]-1/2) > > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Paired PSCBS segmentation using TCN and DH only > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > fit <- segmentByPairedPSCBS(CT, rho=rho, + chromosome=chromosome, x=x, + seed=0xBEEF, verbose=-10) Segmenting paired tumor-normal signals using Paired PSCBS... Setup up data... 'data.frame': 14670 obs. of 4 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : num 1145994 2941694 3710825 4240737 4276892 ... $ CT : num 1.63 1.4 1.41 1.17 1.16 ... $ rho : num NA 0.662 NA NA NA ... Setup up data...done Dropping loci for which TCNs are missing... Number of loci dropped: 12 Dropping loci for which TCNs are missing...done Ordering data along genome... 'data.frame': 14658 obs. of 4 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : num 554484 730720 782343 878522 916294 ... $ CT : num 1.88 1.8 1.59 1.64 1.53 ... $ rho : num NA NA NA NA NA ... Ordering data along genome...done Keeping only current chromosome for 'knownSegments'... Chromosome: 1 Known segments for this chromosome: [1] chromosome start end <0 rows> (or 0-length row.names) Keeping only current chromosome for 'knownSegments'...done alphaTCN: 0.009 alphaDH: 0.001 Number of loci: 14658 Random seed temporarily set (seed=c(48879), kind="L'Ecuyer-CMRG") Produced 2 seeds from this stream for future usage Identification of change points by total copy numbers... Segmenting by CBS... Chromosome: 1 Random seed temporarily set (seed=c(10407, 1066287653, -51199871, 161854402, -1995183193, 1503453565, -747102133), kind="L'Ecuyer-CMRG") Segmenting by CBS...done List of 4 $ data :'data.frame': 14658 obs. of 4 variables: ..$ chromosome: int [1:14658] 1 1 1 1 1 1 1 1 1 1 ... ..$ x : num [1:14658] 554484 730720 782343 878522 916294 ... ..$ y : num [1:14658] 1.88 1.8 1.59 1.64 1.53 ... ..$ index : int [1:14658] 1 2 3 4 5 6 7 8 9 10 ... $ output :'data.frame': 3 obs. of 6 variables: ..$ sampleName: chr [1:3] NA NA NA ..$ chromosome: int [1:3] 1 1 1 ..$ start : num [1:3] 5.54e+05 1.44e+08 1.85e+08 ..$ end : num [1:3] 1.44e+08 1.85e+08 2.47e+08 ..$ nbrOfLoci : int [1:3] 7599 2668 4391 ..$ mean : num [1:3] 1.39 2.07 2.63 $ segRows:'data.frame': 3 obs. of 2 variables: ..$ startRow: int [1:3] 1 7600 10268 ..$ endRow : int [1:3] 7599 10267 14658 $ params :List of 5 ..$ alpha : num 0.009 ..$ undo : num 0 ..$ joinSegments : logi TRUE ..$ knownSegments:'data.frame': 1 obs. of 3 variables: .. ..$ chromosome: int 1 .. ..$ start : num -Inf .. ..$ end : num Inf ..$ seed : int [1:7] 10407 1066287653 -51199871 161854402 -1995183193 1503453565 -747102133 - attr(*, "class")= chr [1:2] "CBS" "AbstractCBS" - attr(*, "processingTime")= 'proc_time' Named num [1:5] 0.57 0 0.58 NA NA ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ... - attr(*, "pkgDetails")= chr "DNAcopy v1.76.0" - attr(*, "randomSeed")= int [1:7] 10407 1066287653 -51199871 161854402 -1995183193 1503453565 -747102133 Identification of change points by total copy numbers...done Restructure TCN segmentation results... chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean 1 1 554484 143926517 7599 1.3859 2 1 143926517 185449813 2668 2.0704 3 1 185449813 247137334 4391 2.6341 Number of TCN segments: 3 Restructure TCN segmentation results...done Total CN segment #1 ([ 554484,1.43927e+08]) of 3... Number of TCN loci in segment: 7599 Locus data for TCN segment: 'data.frame': 7599 obs. of 5 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : num 554484 730720 782343 878522 916294 ... $ CT : num 1.88 1.8 1.59 1.64 1.53 ... $ rho : num NA NA NA NA NA ... $ index : int 1 2 3 4 5 6 7 8 9 10 ... Number of loci: 7599 Number of SNPs: 2111 (27.78%) Number of heterozygous SNPs: 2111 (100.00%) Chromosome: 1 Segmenting DH signals... Segmenting by CBS... Chromosome: 1 Random seed temporarily set (seed=c(10407, 1797822437, 388243314, 91406689, -519578635, -1381924756, 1089253019), kind="L'Ecuyer-CMRG") Segmenting by CBS...done List of 4 $ data :'data.frame': 7599 obs. of 4 variables: ..$ chromosome: int [1:7599] 1 1 1 1 1 1 1 1 1 1 ... ..$ x : num [1:7599] 554484 730720 782343 878522 916294 ... ..$ y : num [1:7599] NA NA NA NA NA ... ..$ index : int [1:7599] 1 2 3 4 5 6 7 8 9 10 ... $ output :'data.frame': 1 obs. of 6 variables: ..$ sampleName: chr NA ..$ chromosome: int 1 ..$ start : num 554484 ..$ end : num 1.44e+08 ..$ nbrOfLoci : int 2111 ..$ mean : num 0.524 $ segRows:'data.frame': 1 obs. of 2 variables: ..$ startRow: int 10 ..$ endRow : int 7594 $ params :List of 5 ..$ alpha : num 0.001 ..$ undo : num 0 ..$ joinSegments : logi TRUE ..$ knownSegments:'data.frame': 1 obs. of 3 variables: .. ..$ chromosome: int 1 .. ..$ start : num 554484 .. ..$ end : num 1.44e+08 ..$ seed : int [1:7] 10407 1797822437 388243314 91406689 -519578635 -1381924756 1089253019 - attr(*, "class")= chr [1:2] "CBS" "AbstractCBS" - attr(*, "processingTime")= 'proc_time' Named num [1:5] 0.07 0 0.06 NA NA ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ... - attr(*, "pkgDetails")= chr "DNAcopy v1.76.0" - attr(*, "randomSeed")= int [1:7] 10407 1797822437 388243314 91406689 -519578635 -1381924756 1089253019 DH segmentation (locally-indexed) rows: startRow endRow 1 10 7594 int [1:7599] 1 2 3 4 5 6 7 8 9 10 ... DH segmentation rows: startRow endRow 1 10 7594 Segmenting DH signals...done DH segmentation table: dhStart dhEnd dhNbrOfLoci dhMean 1 554484 143926517 2111 0.5237 startRow endRow 1 10 7594 Rows: [1] 1 TCN segmentation rows: startRow endRow 1 1 7599 TCN and DH segmentation rows: startRow endRow 1 1 7599 startRow endRow 1 10 7594 NULL TCN segmentation (expanded) rows: startRow endRow 1 1 7599 TCN and DH segmentation rows: startRow endRow 1 1 7599 2 7600 10267 3 10268 14658 startRow endRow 1 10 7594 startRow endRow 1 1 7599 Total CN segmentation table (expanded): chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs tcnNbrOfHets 1 1 554484 143926517 7599 1.3859 2111 2111 (TCN,DH) segmentation for one total CN segment: tcnId dhId chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean 1 2111 554484 143926517 2111 0.5237 Total CN segment #1 ([ 554484,1.43927e+08]) of 3...done Total CN segment #2 ([1.43927e+08,1.8545e+08]) of 3... Number of TCN loci in segment: 2668 Locus data for TCN segment: 'data.frame': 2668 obs. of 5 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : num 1.44e+08 1.44e+08 1.44e+08 1.44e+08 1.44e+08 ... $ CT : num 2.1 2.1 2.09 1.8 2.34 ... $ rho : num NA NA NA NA NA NA NA NA NA NA ... $ index : int 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... Number of loci: 2668 Number of SNPs: 774 (29.01%) Number of heterozygous SNPs: 774 (100.00%) Chromosome: 1 Segmenting DH signals... Segmenting by CBS... Chromosome: 1 Random seed temporarily set (seed=c(10407, 1797822437, 388243314, 91406689, -519578635, -1381924756, 1089253019), kind="L'Ecuyer-CMRG") Segmenting by CBS...done List of 4 $ data :'data.frame': 2668 obs. of 4 variables: ..$ chromosome: int [1:2668] 1 1 1 1 1 1 1 1 1 1 ... ..$ x : num [1:2668] 1.44e+08 1.44e+08 1.44e+08 1.44e+08 1.44e+08 ... ..$ y : num [1:2668] NA NA NA NA NA NA NA NA NA NA ... ..$ index : int [1:2668] 1 2 3 4 5 6 7 8 9 10 ... $ output :'data.frame': 1 obs. of 6 variables: ..$ sampleName: chr NA ..$ chromosome: int 1 ..$ start : num 1.44e+08 ..$ end : num 1.85e+08 ..$ nbrOfLoci : int 774 ..$ mean : num 0.154 $ segRows:'data.frame': 1 obs. of 2 variables: ..$ startRow: int 15 ..$ endRow : int 2664 $ params :List of 5 ..$ alpha : num 0.001 ..$ undo : num 0 ..$ joinSegments : logi TRUE ..$ knownSegments:'data.frame': 1 obs. of 3 variables: .. ..$ chromosome: int 1 .. ..$ start : num 1.44e+08 .. ..$ end : num 1.85e+08 ..$ seed : int [1:7] 10407 1797822437 388243314 91406689 -519578635 -1381924756 1089253019 - attr(*, "class")= chr [1:2] "CBS" "AbstractCBS" - attr(*, "processingTime")= 'proc_time' Named num [1:5] 0.03 0 0.03 NA NA ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ... - attr(*, "pkgDetails")= chr "DNAcopy v1.76.0" - attr(*, "randomSeed")= int [1:7] 10407 1797822437 388243314 91406689 -519578635 -1381924756 1089253019 DH segmentation (locally-indexed) rows: startRow endRow 1 15 2664 int [1:2668] 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... DH segmentation rows: startRow endRow 1 7614 10263 Segmenting DH signals...done DH segmentation table: dhStart dhEnd dhNbrOfLoci dhMean 1 143926517 185449813 774 0.1542 startRow endRow 1 7614 10263 Rows: [1] 2 TCN segmentation rows: startRow endRow 2 7600 10267 TCN and DH segmentation rows: startRow endRow 2 7600 10267 startRow endRow 1 7614 10263 startRow endRow 1 1 7599 TCN segmentation (expanded) rows: startRow endRow 1 1 7599 2 7600 10267 TCN and DH segmentation rows: startRow endRow 1 1 7599 2 7600 10267 3 10268 14658 startRow endRow 1 10 7594 2 7614 10263 startRow endRow 1 1 7599 2 7600 10267 Total CN segmentation table (expanded): chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs tcnNbrOfHets 2 1 143926517 185449813 2668 2.0704 774 774 (TCN,DH) segmentation for one total CN segment: tcnId dhId chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 2 2 1 1 143926517 185449813 2668 2.0704 774 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean 2 774 143926517 185449813 774 0.1542 Total CN segment #2 ([1.43927e+08,1.8545e+08]) of 3...done Total CN segment #3 ([1.8545e+08,2.47137e+08]) of 3... Number of TCN loci in segment: 4391 Locus data for TCN segment: 'data.frame': 4391 obs. of 5 variables: $ chromosome: int 1 1 1 1 1 1 1 1 1 1 ... $ x : num 1.85e+08 1.85e+08 1.85e+08 1.86e+08 1.86e+08 ... $ CT : num 2.93 2.15 2.82 2.93 2.46 ... $ rho : num NA 0.0308 NA 0.2533 NA ... $ index : int 10268 10269 10270 10271 10272 10273 10274 10275 10276 10277 ... Number of loci: 4391 Number of SNPs: 1311 (29.86%) Number of heterozygous SNPs: 1311 (100.00%) Chromosome: 1 Segmenting DH signals... Segmenting by CBS... Chromosome: 1 Random seed temporarily set (seed=c(10407, 1797822437, 388243314, 91406689, -519578635, -1381924756, 1089253019), kind="L'Ecuyer-CMRG") Segmenting by CBS...done List of 4 $ data :'data.frame': 4391 obs. of 4 variables: ..$ chromosome: int [1:4391] 1 1 1 1 1 1 1 1 1 1 ... ..$ x : num [1:4391] 1.85e+08 1.85e+08 1.85e+08 1.86e+08 1.86e+08 ... ..$ y : num [1:4391] NA 0.0308 NA 0.2533 NA ... ..$ index : int [1:4391] 1 2 3 4 5 6 7 8 9 10 ... $ output :'data.frame': 1 obs. of 6 variables: ..$ sampleName: chr NA ..$ chromosome: int 1 ..$ start : num 1.85e+08 ..$ end : num 2.47e+08 ..$ nbrOfLoci : int 1311 ..$ mean : num 0.251 $ segRows:'data.frame': 1 obs. of 2 variables: ..$ startRow: int 2 ..$ endRow : int 4388 $ params :List of 5 ..$ alpha : num 0.001 ..$ undo : num 0 ..$ joinSegments : logi TRUE ..$ knownSegments:'data.frame': 1 obs. of 3 variables: .. ..$ chromosome: int 1 .. ..$ start : num 1.85e+08 .. ..$ end : num 2.47e+08 ..$ seed : int [1:7] 10407 1797822437 388243314 91406689 -519578635 -1381924756 1089253019 - attr(*, "class")= chr [1:2] "CBS" "AbstractCBS" - attr(*, "processingTime")= 'proc_time' Named num [1:5] 0.03 0 0.03 NA NA ..- attr(*, "names")= chr [1:5] "user.self" "sys.self" "elapsed" "user.child" ... - attr(*, "pkgDetails")= chr "DNAcopy v1.76.0" - attr(*, "randomSeed")= int [1:7] 10407 1797822437 388243314 91406689 -519578635 -1381924756 1089253019 DH segmentation (locally-indexed) rows: startRow endRow 1 2 4388 int [1:4391] 10268 10269 10270 10271 10272 10273 10274 10275 10276 10277 ... DH segmentation rows: startRow endRow 1 10269 14655 Segmenting DH signals...done DH segmentation table: dhStart dhEnd dhNbrOfLoci dhMean 1 185449813 247137334 1311 0.2512 startRow endRow 1 10269 14655 Rows: [1] 3 TCN segmentation rows: startRow endRow 3 10268 14658 TCN and DH segmentation rows: startRow endRow 3 10268 14658 startRow endRow 1 10269 14655 startRow endRow 1 1 7599 2 7600 10267 TCN segmentation (expanded) rows: startRow endRow 1 1 7599 2 7600 10267 3 10268 14658 TCN and DH segmentation rows: startRow endRow 1 1 7599 2 7600 10267 3 10268 14658 startRow endRow 1 10 7594 2 7614 10263 3 10269 14655 startRow endRow 1 1 7599 2 7600 10267 3 10268 14658 Total CN segmentation table (expanded): chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs tcnNbrOfHets 3 1 185449813 247137334 4391 2.6341 1311 1311 (TCN,DH) segmentation for one total CN segment: tcnId dhId chromosome tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 3 3 1 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean 3 1311 185449813 247137334 1311 0.2512 Total CN segment #3 ([1.8545e+08,2.47137e+08]) of 3...done chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 2 1 2 1 143926517 185449813 2668 2.0704 774 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean 1 2111 554484 143926517 2111 0.5237 2 774 143926517 185449813 774 0.1542 3 1311 185449813 247137334 1311 0.2512 Calculating (C1,C2) per segment... Calculating (C1,C2) per segment...done Number of segments: 3 Segmenting paired tumor-normal signals using Paired PSCBS...done Post-segmenting TCNs... Number of segments: 3 Number of chromosomes: 1 [1] 1 Chromosome 1 ('chr01') of 1... Rows: [1] 1 2 3 Number of segments: 3 TCN segment #1 ('1') of 3... Nothing todo. Only one DH segmentation. Skipping. TCN segment #1 ('1') of 3...done TCN segment #2 ('2') of 3... Nothing todo. Only one DH segmentation. Skipping. TCN segment #2 ('2') of 3...done TCN segment #3 ('3') of 3... Nothing todo. Only one DH segmentation. Skipping. TCN segment #3 ('3') of 3...done Chromosome 1 ('chr01') of 1...done Update (C1,C2) per segment... Update (C1,C2) per segment...done Post-segmenting TCNs...done chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 2 1 2 1 143926517 185449813 2668 2.0704 774 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean c1Mean c2Mean 1 2111 554484 143926517 2111 0.5237 0.3300521 1.055848 2 774 143926517 185449813 774 0.1542 0.8755722 1.194828 3 1311 185449813 247137334 1311 0.2512 0.9862070 1.647893 chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 2 1 2 1 143926517 185449813 2668 2.0704 774 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean c1Mean c2Mean 1 2111 554484 143926517 2111 0.5237 0.3300521 1.055848 2 774 143926517 185449813 774 0.1542 0.8755722 1.194828 3 1311 185449813 247137334 1311 0.2512 0.9862070 1.647893 > print(fit) chromosome tcnId dhId start end tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 2 1 2 1 143926517 185449813 2668 2.0704 774 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhNbrOfLoci dhMean c1Mean c2Mean 1 2111 2111 0.5237 0.3300521 1.055848 2 774 774 0.1542 0.8755722 1.194828 3 1311 1311 0.2512 0.9862070 1.647893 > > # Plot results > plotTracks(fit) > > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Bootstrap segment level estimates > # (used by the AB caller, which, if skipped here, > # will do it automatically) > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > fit <- bootstrapTCNandDHByRegion(fit, B=B, verbose=-10) Resample (TCN,DH) signals and re-estimate summaries for segment & changepoint... Already done? tcn_2.5% tcn_5% tcn_95% tcn_97.5% dh_2.5% dh_5% dh_95% dh_97.5% FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE c1_2.5% c1_5% c1_95% c1_97.5% c2_2.5% c2_5% c2_95% c2_97.5% FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE Bootstrapping (TCN,DH,C1,C2) segment mean levels... Identifying heterozygous & homozygous SNPs and non-polymorphic loci... Number of loci: 14658 Number of SNPs: 4196 Number of non-SNPs: 10462 Identifying heterozygous & homozygous SNPs and non-polymorphic loci...done num [1:3, 1:100, 1:4] NA NA NA NA NA NA NA NA NA NA ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : NULL ..$ : chr [1:4] "tcn" "dh" "c1" "c2" Segment #1 (chr 1, tcnId=1, dhId=1) of 3... chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean c1Mean c2Mean 1 2111 554484 143926517 2111 0.5237 0.3300521 1.055848 Number of TCNs: 7599 Number of DHs: 2111 int [1:7599] 1 2 3 4 5 6 7 8 9 10 ... int [1:14658] 1 2 3 4 5 6 7 8 9 10 ... int [1:7599] 1 2 3 4 5 6 7 8 9 10 ... Identify loci used to bootstrap DH means... Heterozygous SNPs to resample for DH: int [1:2111] 10 12 24 28 31 33 34 39 46 48 ... Identify loci used to bootstrap DH means...done Identify loci used to bootstrap TCN means... SNPs: int [1:2111] 10 12 24 28 31 33 34 39 46 48 ... Non-polymorphic loci: int [1:5488] 1 2 3 4 5 6 7 8 9 11 ... Heterozygous SNPs to resample for TCN: int [1:2111] 10 12 24 28 31 33 34 39 46 48 ... Homozygous SNPs to resample for TCN: int(0) Non-polymorphic loci to resample for TCN: int [1:5488] 1 2 3 4 5 6 7 8 9 11 ... Heterozygous SNPs with non-DH to resample for TCN: int(0) Loci to resample for TCN: int [1:7599] 1 2 3 4 5 6 7 8 9 10 ... Identify loci used to bootstrap TCN means...done Number of (#hets, #homs, #nonSNPs): (2111,0,5488) Bootstrapping while preserving (#hets, #homs, #nonSNPs)... Number of bootstrap samples: 100 Bootstrapping while preserving (#hets, #homs, #nonSNPs)...done Segment #1 (chr 1, tcnId=1, dhId=1) of 3...done Segment #2 (chr 1, tcnId=2, dhId=1) of 3... chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 2 1 2 1 143926517 185449813 2668 2.0704 774 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean c1Mean c2Mean 2 774 143926517 185449813 774 0.1542 0.8755722 1.194828 Number of TCNs: 2668 Number of DHs: 774 int [1:2668] 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... int [1:14658] 1 2 3 4 5 6 7 8 9 10 ... int [1:2668] 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... Identify loci used to bootstrap DH means... Heterozygous SNPs to resample for DH: int [1:774] 7614 7616 7626 7627 7628 7635 7638 7639 7640 7642 ... Identify loci used to bootstrap DH means...done Identify loci used to bootstrap TCN means... SNPs: int [1:774] 7614 7616 7626 7627 7628 7635 7638 7639 7640 7642 ... Non-polymorphic loci: int [1:1894] 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... Heterozygous SNPs to resample for TCN: int [1:774] 7614 7616 7626 7627 7628 7635 7638 7639 7640 7642 ... Homozygous SNPs to resample for TCN: int(0) Non-polymorphic loci to resample for TCN: int [1:1894] 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... Heterozygous SNPs with non-DH to resample for TCN: int(0) Loci to resample for TCN: int [1:2668] 7600 7601 7602 7603 7604 7605 7606 7607 7608 7609 ... Identify loci used to bootstrap TCN means...done Number of (#hets, #homs, #nonSNPs): (774,0,1894) Bootstrapping while preserving (#hets, #homs, #nonSNPs)... Number of bootstrap samples: 100 Bootstrapping while preserving (#hets, #homs, #nonSNPs)...done Segment #2 (chr 1, tcnId=2, dhId=1) of 3...done Segment #3 (chr 1, tcnId=3, dhId=1) of 3... chromosome tcnId dhId tcnStart tcnEnd tcnNbrOfLoci tcnMean tcnNbrOfSNPs 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhStart dhEnd dhNbrOfLoci dhMean c1Mean c2Mean 3 1311 185449813 247137334 1311 0.2512 0.986207 1.647893 Number of TCNs: 4391 Number of DHs: 1311 int [1:4391] 10268 10269 10270 10271 10272 10273 10274 10275 10276 10277 ... int [1:14658] 1 2 3 4 5 6 7 8 9 10 ... int [1:4391] 10268 10269 10270 10271 10272 10273 10274 10275 10276 10277 ... Identify loci used to bootstrap DH means... Heterozygous SNPs to resample for DH: int [1:1311] 10269 10271 10281 10282 10283 10293 10295 10297 10300 10302 ... Identify loci used to bootstrap DH means...done Identify loci used to bootstrap TCN means... SNPs: int [1:1311] 10269 10271 10281 10282 10283 10293 10295 10297 10300 10302 ... Non-polymorphic loci: int [1:3080] 10268 10270 10272 10273 10274 10275 10276 10277 10278 10279 ... Heterozygous SNPs to resample for TCN: int [1:1311] 10269 10271 10281 10282 10283 10293 10295 10297 10300 10302 ... Homozygous SNPs to resample for TCN: int(0) Non-polymorphic loci to resample for TCN: int [1:3080] 10268 10270 10272 10273 10274 10275 10276 10277 10278 10279 ... Heterozygous SNPs with non-DH to resample for TCN: int(0) Loci to resample for TCN: int [1:4391] 10268 10269 10270 10271 10272 10273 10274 10275 10276 10277 ... Identify loci used to bootstrap TCN means...done Number of (#hets, #homs, #nonSNPs): (1311,0,3080) Bootstrapping while preserving (#hets, #homs, #nonSNPs)... Number of bootstrap samples: 100 Bootstrapping while preserving (#hets, #homs, #nonSNPs)...done Segment #3 (chr 1, tcnId=3, dhId=1) of 3...done Bootstrapped segment mean levels num [1:3, 1:100, 1:4] 1.38 2.08 2.63 1.38 2.07 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : NULL ..$ : chr [1:4] "tcn" "dh" "c1" "c2" Calculating (C1,C2) mean levels from (TCN,DH) mean levels... num [1:3, 1:100, 1:4] 1.38 2.08 2.63 1.38 2.07 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : NULL ..$ : chr [1:4] "tcn" "dh" "c1" "c2" Calculating (C1,C2) mean levels from (TCN,DH) mean levels...done Calculating polar (alpha,radius,manhattan) for change points... num [1:2, 1:100, 1:2] -0.5588 -0.0962 -0.5365 -0.1285 -0.5378 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : NULL ..$ : chr [1:2] "c1" "c2" Bootstrapped change points num [1:2, 1:100, 1:5] -2.89 -1.78 -2.87 -1.86 -2.88 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : NULL ..$ : chr [1:5] "alpha" "radius" "manhattan" "d1" ... Calculating polar (alpha,radius,manhattan) for change points...done Bootstrapping (TCN,DH,C1,C2) segment mean levels...done Summarizing bootstrapped segment ('tcn', 'dh', 'c1', 'c2') data... num [1:3, 1:4, 1:4] NA NA NA NA NA NA NA NA NA NA ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : chr [1:4] "2.5%" "5%" "95%" "97.5%" ..$ : chr [1:4] "tcn" "dh" "c1" "c2" Field #1 ('tcn') of 4... Segment #1 of 3... Segment #1 of 3...done Segment #2 of 3... Segment #2 of 3...done Segment #3 of 3... Segment #3 of 3...done Field #1 ('tcn') of 4...done Field #2 ('dh') of 4... Segment #1 of 3... Segment #1 of 3...done Segment #2 of 3... Segment #2 of 3...done Segment #3 of 3... Segment #3 of 3...done Field #2 ('dh') of 4...done Field #3 ('c1') of 4... Segment #1 of 3... Segment #1 of 3...done Segment #2 of 3... Segment #2 of 3...done Segment #3 of 3... Segment #3 of 3...done Field #3 ('c1') of 4...done Field #4 ('c2') of 4... Segment #1 of 3... Segment #1 of 3...done Segment #2 of 3... Segment #2 of 3...done Segment #3 of 3... Segment #3 of 3...done Field #4 ('c2') of 4...done Bootstrap statistics num [1:3, 1:4, 1:4] 1.38 2.06 2.62 1.38 2.06 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : chr [1:4] "2.5%" "5%" "95%" "97.5%" ..$ : chr [1:4] "tcn" "dh" "c1" "c2" Statistical sanity checks (iff B >= 100)... Available summaries: 2.5%, 5%, 95%, 97.5% Available quantiles: 0.025, 0.05, 0.95, 0.975 num [1:3, 1:4, 1:4] 1.38 2.06 2.62 1.38 2.06 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : chr [1:4] "2.5%" "5%" "95%" "97.5%" ..$ : chr [1:4] "tcn" "dh" "c1" "c2" Field #1 ('tcn') of 4... Seg 1. mean=1.3859, range=[1.38092,1.3949], n=7599 Seg 2. mean=2.0704, range=[2.05747,2.08326], n=2668 Seg 3. mean=2.6341, range=[2.62068,2.64694], n=4391 Field #1 ('tcn') of 4...done Field #2 ('dh') of 4... Seg 1. mean=0.5237, range=[0.51753,0.532002], n=2111 Seg 2. mean=0.1542, range=[0.144468,0.16453], n=774 Seg 3. mean=0.2512, range=[0.242575,0.258832], n=1311 Field #2 ('dh') of 4...done Field #3 ('c1') of 4... Seg 1. mean=0.330052, range=[0.323996,0.336038], n=2111 Seg 2. mean=0.875572, range=[0.86318,0.887699], n=774 Seg 3. mean=0.986207, range=[0.975123,0.998982], n=1311 Field #3 ('c1') of 4...done Field #4 ('c2') of 4... Seg 1. mean=1.05585, range=[1.05006,1.06231], n=2111 Seg 2. mean=1.19483, range=[1.18417,1.2081], n=774 Seg 3. mean=1.64789, range=[1.63403,1.66098], n=1311 Field #4 ('c2') of 4...done Statistical sanity checks (iff B >= 100)...done Summarizing bootstrapped segment ('tcn', 'dh', 'c1', 'c2') data...done Summarizing bootstrapped changepoint ('alpha', 'radius', 'manhattan', 'd1', 'd2') data... num [1:2, 1:4, 1:5] NA NA NA NA NA NA NA NA NA NA ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : chr [1:4] "2.5%" "5%" "95%" "97.5%" ..$ : chr [1:5] "alpha" "radius" "manhattan" "d1" ... Field #1 ('alpha') of 5... Changepoint #1 of 2... Changepoint #1 of 2...done Changepoint #2 of 2... Changepoint #2 of 2...done Field #1 ('alpha') of 5...done Field #2 ('radius') of 5... Changepoint #1 of 2... Changepoint #1 of 2...done Changepoint #2 of 2... Changepoint #2 of 2...done Field #2 ('radius') of 5...done Field #3 ('manhattan') of 5... Changepoint #1 of 2... Changepoint #1 of 2...done Changepoint #2 of 2... Changepoint #2 of 2...done Field #3 ('manhattan') of 5...done Field #4 ('d1') of 5... Changepoint #1 of 2... Changepoint #1 of 2...done Changepoint #2 of 2... Changepoint #2 of 2...done Field #4 ('d1') of 5...done Field #5 ('d2') of 5... Changepoint #1 of 2... Changepoint #1 of 2...done Changepoint #2 of 2... Changepoint #2 of 2...done Field #5 ('d2') of 5...done Bootstrap statistics num [1:2, 1:4, 1:5] -2.92 -1.86 -2.91 -1.85 -2.87 ... - attr(*, "dimnames")=List of 3 ..$ : NULL ..$ : chr [1:4] "2.5%" "5%" "95%" "97.5%" ..$ : chr [1:5] "alpha" "radius" "manhattan" "d1" ... Summarizing bootstrapped changepoint ('alpha', 'radius', 'manhattan', 'd1', 'd2') data...done Resample (TCN,DH) signals and re-estimate summaries for segment & changepoint...done > print(fit) chromosome tcnId dhId start end tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 2 1 2 1 143926517 185449813 2668 2.0704 774 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhNbrOfLoci dhMean c1Mean c2Mean 1 2111 2111 0.5237 0.3300521 1.055848 2 774 774 0.1542 0.8755722 1.194828 3 1311 1311 0.2512 0.9862070 1.647893 > plotTracks(fit) > > > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > # Calling segments in allelic balance (AB) and > # in loss-of-heterozygosity (LOH) > # NOTE: Ideally, this should be done on whole-genome data > # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - > fit <- callAB(fit, verbose=-10) Calling segments of allelic balance from one-sided DH bootstrap confidence intervals... delta (offset adjusting for bias in DH): 0.3466649145302 alpha (CI quantile; significance level): 0.05 Calling segments... Number of segments called allelic balance (AB): 2 (66.67%) of 3 Calling segments...done Calling segments of allelic balance from one-sided DH bootstrap confidence intervals...done > fit <- callLOH(fit, verbose=-10) Calling segments of allelic balance from one-sided DH bootstrap confidence intervals... delta (offset adjusting for bias in C1): 0.771236438183453 alpha (CI quantile; significance level): 0.05 Calling segments... Number of segments called low C1 (LowC1, "LOH_C1"): 1 (33.33%) of 3 Calling segments...done Calling segments of allelic balance from one-sided DH bootstrap confidence intervals...done > print(fit) chromosome tcnId dhId start end tcnNbrOfLoci tcnMean tcnNbrOfSNPs 1 1 1 1 554484 143926517 7599 1.3859 2111 2 1 2 1 143926517 185449813 2668 2.0704 774 3 1 3 1 185449813 247137334 4391 2.6341 1311 tcnNbrOfHets dhNbrOfLoci dhMean c1Mean c2Mean abCall lohCall 1 2111 2111 0.5237 0.3300521 1.055848 FALSE TRUE 2 774 774 0.1542 0.8755722 1.194828 TRUE FALSE 3 1311 1311 0.2512 0.9862070 1.647893 TRUE FALSE > plotTracks(fit) > > proc.time() user system elapsed 2.50 0.15 3.00