# iCNV Vignette

## Maintainer

Zilu Zhou zhouzilu@pennmedicine.upenn.edu

This is a demo for using the iCNV package in R. iCNV is a normalization and copy number variation detection procedure for multiple study designs: WES only, WGS only, SNP array only, or any combination of SNP and sequencing data. iCNV applies platform specific normalization, utilizes allele specific reads from sequencing and integrates matched NGS and SNP-array data by a Hidden Markov Model (HMM). Figure below shows the overall pipeline. This package specifically emphasizes on the steps within the parenthesis. Below is an example on calling copy number variation using whole-exome sequencing data and array SNPs of 46 modified samples from 1000 Genome Project. Only small portion of chromosome 22 are analyzed for illustration purposes. We will separately shows normalization for WES and SNP array. We will also introduce integrated calling procedure as well as single platform procedure.

                  NGS                                             |         Array
     BAM        BED(UCSC for WES or bed_generator.R for WGS 2.2)  |  SNP Intensity(in standard format)
      |----------------|                                          |           |
      |                |                                          |           |icnv_array_input(2.4)
      |SAMTools(2.3)   |CODEX(2.2)                                |           |
      |                |                                          |           |-----------|
    Variants BAF(vcf) PLR                                         |      Array LRR   Array BAF
      |                |                                          |           |           |
      |                |                                          |           |SVD(2.4)   |
      |bambaf_from_vcf |                                          |           |           |
      | (2.3)          |                                          |   Normalized LRR      |
      |                |                                          |           |           |
      -------------------------------------------------------------------------------------
                                                      |
                                                      |iCNV_detection(2.5,2.6)
                                                      |
                                                 CNV calling
                                                      |
                                                      |icnv_output_to_gb(2.5,2.6)
                                                      |
                                                Genome Browser file

We strongly recommend combining platform when both WES data and SNP array are available. However, for high quality WGS data, SNP information isn’t so necessary.

## 1. Website and online forum

iCNV website with usage and installation information: https://github.com/zhouzilu/iCNV

## 2. iCNV workflow

We could separate the basic iCNV workflow into 6 sections: 1. package installation; 2. .bam file normalization; 3. sequence variants BAF calling; 4. SNP array LRR normalization and BAF; 5. CNV detection using `iCNV_detection` function. 6. CNV detection using `iCNV_detection` function with single platform. We will illustrate them one by one in the following sessions.

### 2.1 Install iCNV.

Install `CODEX` first

    ## try http:// if https:// URLs are not supported
    if (!requireNamespace("BiocManager", quietly=TRUE))
        install.packages("BiocManager")
    BiocManager::install("CODEX")

Install `iCNV` from Bioconductor:

    ## try http:// if https:// URLs are not supported
    if (!requireNamespace("BiocManager", quietly=TRUE))
        install.packages("BiocManager")
    BiocManager::install("iCNV")

Install the current release from Github:

    install.packages("devtools")
    library(devtools)
    install_github("zhouzilu/iCNV")

### 2.2 .bam file normalization using CODEX

For researchers work on WGS data, you need to generate your own BED file for CODEX normalization. We made a function bed_generator to help you with this.

After you have BAM file and BED file available, follow code below to get normalized NGS PLR. Code modified from https://github.com/yuchaojiang/CODEX with permission. If you have any additional question, please refer to `CODEX`. iCNV only adopt normalized Poisson Likelihood Ratio from `CODEX`, representing NGS ‘intensity’. In this example, we utilize toy data from the 1000 Genomes Project for illustration.

    ### 2.2 code chunk number 1: bambedObj1
    library(CODEX)
    library(WES.1KG.WUGSC) # Load Toy data from the 1000 Genomes Project.
    library(iCNV)
    # Assume your work dir set to demo
    dir <- system.file("extdata", package = "WES.1KG.WUGSC") # 'PATH/TO/BAM'
    bamFile <- list.files(dir, pattern = '*.bam$')
    bamdir <- file.path(dir, bamFile)
    sampname <- as.matrix(read.table(file.path(dir, "sampname")))
    bedFile <- file.path(dir, "chr22_400_to_500.bed")
    chr <- 22
    bambedObj <- getbambed(bamdir = bamdir, bedFile = bedFile,
                       sampname = sampname, projectname = "icnv.demo.", chr)
    bamdir <- bambedObj$bamdir; sampname <- bambedObj$sampname
    ref <- bambedObj$ref; projectname <- bambedObj$projectname; chr <- bambedObj$chr

    ### 2.2 code chunk number 2: coverageObj1
    coverageObj <- getcoverage(bambedObj, mapqthres = 20)
    Y <- coverageObj$Y; readlength <- coverageObj$readlength

    ### 2.2 code chunk number 3: gcmapp1
    gc <- getgc(chr, ref)
    mapp <- getmapp(chr, ref)

    ### 2.2 code chunk number 4: qcObj1
    qcObj <- qc(Y, sampname, chr, ref, mapp, gc, cov_thresh = c(20, 4000),
                length_thresh = c(20, 2000), mapp_thresh = 0.9, gc_thresh = c(20, 80))
    Y_qc <- qcObj$Y_qc; sampname_qc <- qcObj$sampname_qc; gc_qc <- qcObj$gc_qc
    mapp_qc <- qcObj$mapp_qc; ref_qc <- qcObj$ref_qc; qcmat <- qcObj$qcmat

    ### 2.2 code chunk number 5: normObj1
    normObj <- normalize(Y_qc, gc_qc, K = 1:9)
    Yhat <- normObj$Yhat; AIC <- normObj$AIC; BIC <- normObj$BIC
    RSS <- normObj$RSS; K <- normObj$K
    choiceofK(AIC, BIC, RSS, K, filename = paste0(projectname, chr, "_choiceofK.pdf"))
    save(qcObj,normObj,sampname,file=paste0(projectname,'plrObj_', chr,".rda"))
    # This saving step is necessary for large dataset, such that you don't need
    # to recompute everything again if you would like to view your choice of K
    # plot. You could just load this as in the next chunk of code.

CODEX reports all three statistical metrics (AIC, BIC, percent of Variance explained) and uses BIC as the default method to determine the number of latent Poisson factors. Since false positives can be screened out through a closer examination of the post-segmentation data, whereas CNV signals removed in the normalization step cannot be recovered, CODEX opts for a more conservative normalization that, when in doubt, uses a smaller value of K. For example,

Here we pick `optK=2`

    load(paste0(projectname,'plrObj_',chr,".rda"))
    optK <- K[which.max(BIC)] # by default or customize given curves

    ### 2.2 code chunk number 6: plr
    Y_qc <- qcObj$Y_qc; sampname_qc <- qcObj$sampname_qc; gc_qc <- qcObj$gc_qc
    mapp_qc <- qcObj$mapp_qc; ref_qc <- qcObj$ref_qc; qcmat <- qcObj$qcmat
    Yhat <- normObj$Yhat; AIC <- normObj$AIC; BIC <- normObj$BIC
    RSS <- normObj$RSS; K <- normObj$K
    ref_qc <- qcObj$ref_qc # IRanges object for exon target
    sampname_qc <- qcObj$sampname_qc # sample names
    Y_norm <- normObj$Yhat[[optK]] # normalized read count under null (no CNV)
    plr <- log(pmax(Y_qc,0.0001)/pmax(Y_norm,0.0001)) # log transformed z-scores
    ngs_plr <- lapply(seq_len(ncol(plr)), function(i) plr[,i])
    ngs_plr.pos <- lapply(seq_len(ncol(plr)),function(x){return(cbind(start(ref_qc),end(ref_qc)))})
    save(sampname_qc,ngs_plr,ngs_plr.pos,file=paste0(projectname,'plrObj_',chr,'_',optK,'.rda'))

For detailed illustration of CODEX, please visit https://github.com/yuchaojiang/CODEX

### 2.3 sequence variants BAF calling

For sequencing data without sophisticated pipeline and SNVs call set in VCF format, we manually call SNVs from quality controlled BAM files by mpileup module in samtools, and calculate B allele frequency(BAF) on heterogeneous loci by dividing AD (Number of high-quality non-reference bases, FORMAT) from DP (Number of high-quality bases, FORMAT). Example code are:

    # Prerequest: samtools, bcftools and reference fasta file
    cd PATH/TO/BAM
    for i in *bam; do PATH/TO/SAMTOOLS/samtools mpileup -ugI -t AD -t DP -f
      PATH/TO/REF/human_hg37.fasta $i | PATH/TO/BAFTOOLS/bcftools
      call -cv -O z -o PATH/TO/OUTPUT/$i.vcf.gz; done
    # The reference genome depends on your BAM built
    for i in *gz; do PATH/TO/BAFTOOLS/bcftools view -g het -i 'FORMAT/DP>10' -O z -o $i.filter $i; done

We could further extract variants BAF info from vcf file by function bambaf_from_vcf.

    ### 2.3 code chunk number 1: bambaf_from_vcf
    projectname <- "icnv.demo."
    dir <- system.file("extdata", package="iCNV")
    bambaf_from_vcf(dir,'bam_vcf.list',chr=22,projectname)
    # ignore chr argement if you want to convert all chromosome 1-22
    load(paste0(projectname,'bambaf_22.rda'))

### 2.4 SNP array LRR normalization and BAF

Because signal intensity files varies from platform, we set a standard signal intensity file format for each individual. The format looks as follow:

    Name,Chr,POS,sample_1.Log R Ratio,sample_1.B Allele Freq
    rs1,22,15462739,-0.096390619874,0.0443964861333
    rs2,22,15520991,-0.154103130102,0.963218688965
    rs3,22,15780940,-0.110297381878,0.0457459762692
    rs4,22,15863717,-0.21270519495,0.957377910614
    rs5,22,16532045,-0.0330782271922,0.0300635993481

First row is the rsid (Name), followed by chromosome (Chr) and position (POS), then the log R ratio (sample1.Log R Ratio) and BAF (sample1.B Allele Freq). We could use function get_array_input to convert into iCNV input. For example, in demo:

    ### 2.4 code chunk number 1: icnv_array_input
    # PATH/TO/FOLDER where you stored the array data
    dir <- system.file("extdata", package="iCNV")#PATH/TO/FOLDER where you stored the array data
    chr <- 22
    projectname <- "icnv.demo."
    pattern <- paste0('*.csv.arrayicnv$')
    get_array_input(dir,pattern,chr,projectname)
    load(paste0(projectname,'array_lrrbaf_',chr,'.rda'))

For some of the SNP array LRR data, we need to apply SVD normalization to remove high dimension noisy and preserve low dimension signal. The best way to decide data sanity is by plotting out the data by `plot_intensity` function. Noisy data has the feature of local strip across samples. Conventional way for identifying elbow points as the optK can also apply here. Example code for remove high dimension noisy and plot:

    ### 2.4 code chunk number 2: svd (optional)
    library(iCNV)
    projectname <- "icnv.demo."
    load(paste0(projectname,'array_lrrbaf_',chr,'.rda'))
    # standardize LRR
    lrr.sd <- mapply(function(x){
      x=pmax(pmin(x,100),-100);(x-mean(x,na.rm=TRUE))/sd(x,na.rm = TRUE)
      },snp_lrr,SIMPLIFY = TRUE)
    # You may want to substitute all the NAs in lrr by rowMeans
    lrr.sd.dena <- apply(lrr.sd,2,function(x){x[is.na(x)]=mean(x,na.rm=TRUE);return(x)})

    # lrr is LRR matrix, with row as sample and column as positions
    lrr.svd <- svd (t(lrr.sd.dena))
    save(lrr.svd,file=paste0(projectname,'array_lrrbaf_svd_',chr,'.rda'))

    # Plot out the variance explained with K; You can store the plot in a pdf
    # pdf(paste0(projectname,'array_lrrbaf_svd_',chr,'.pdf'),height=8,width=10)
    plot(x=seq_len(10),y=(lrr.svd$d[1:10])^2/sum(lrr.svd$d^2),
         xlab='rank',ylab='D square',main='Variance explained',type='l')

    # dev.off() # Close the plot connection

There is no universal rule of picking the optK. Usually we prefer those “elbow points”. We pick `optK=2`

    ### 2.4 code chunk number 3: normalization (optional)
    library(iCNV)
    optK <- 5
    D <- diag(lrr.svd$d)
    D.lowrank <- diag(c(rep(0,optK),lrr.svd$d[-(seq_len(optK))]))
    lrr.denoise <- t(lrr.svd$u %*% D.lowrank %*% t(lrr.svd$v))
    snp_lrr <- lapply(seq_len(ncol(lrr.denoise)), function(i) lrr.denoise[,i])
    save(snp_lrr,snp_baf,snp_lrr.pos,snp_baf.pos,
        file=file.path(dir,paste0(projectname,'array_lrrbaf_',chr,'_svded.rda')))

### 2.5 Mutiple platform CNV detection using iCNV

At this step, we should already have PLR and variants BAF from sequencing, normalized LRR and BAF from SNP array. All the input should be in list form. This is mainly to accommodate the fact that length for `ngs_baf` and `ngs_baf.pos` is sample specific

    ### 2.5 code chunk number 1: iCNV_detection with genotype
    str(ngs_plr) # List of n vector, each one is the PLR for an exon
    str(ngs_plr.pos) # List of n matrix (p x 2), each one is the start and end location for an exon
    str(ngs_baf) # List of n vector, each one is the variants BAF from .bam
    str(ngs_baf.pos) # List of n vector, each one is the variants BAF position
    str(snp_lrr) # List of n vector, each one is the normalized LRR for a SNP
    str(snp_lrr.pos) # List of n vector, each one is a SNP position
    str(snp_baf) # List of n vector, each one is the BAF for a SNP
    str(snp_baf.pos) # List of n vector, each one is the SNP BAF position
    projectname <- 'icnv.demo.'
    icnv_res0 <- iCNV_detection(ngs_plr,snp_lrr,
                             ngs_baf,snp_baf,
                             ngs_plr.pos,snp_lrr.pos,
                             ngs_baf.pos,snp_baf.pos,
                             projname=projectname,CN=1,mu=c(-3,0,2),cap=TRUE,visual = 1)
    icnv.output <- output_list(icnv_res0,sampname_qc,CN=1)
    head(icnv.output)

The results plot looks are follow. Each row is a sample and each column is a hidden state. The color indicates hidden state Z-score (large positive number prefers amplification, low negative number prefers deletion). Black dots represent amplification detected, while white dots show deletion detected.

We could further convert icnv.output to Genome Browser input file:

    ### 2.5 code chunk number 2: covert result to Genome Browser input
    icnv.output <- output_list(icnv_res=icnv_res0,sampleid=sampname_qc, CN=0, min_size=10000)
    gb_input <- icnv_output_to_gb(22,icnv.output)
    write.table(gb_input,file=paste0(projectname,'icnv_res_gb_chr',chr,'.tab'),
                quote=FALSE,col.names=FALSE,row.names=FALSE)

The color code for CNV in Genome Browser is:

This is different from the color code for us plotting the data. We are just try to make it coordinate with DGV notation. We could also plot out the data in DGV format by setting `col='DGV'` in `plotHMMscore`.

We could also plot information in a single individual using function `plotindi`. Example shows below:

    ### 2.5 code chunk number 3: plotindi
    I <- 24 # ploting out the 24th individual
    plotindi(ngs_plr,snp_lrr,
             ngs_baf,snp_baf,
             ngs_plr.pos,snp_lrr.pos,
             ngs_baf.pos,snp_baf.pos,
             icnv_res0,I)

### 2.6 Single platform CNV detection using iCNV

At this step, we should already have PLR and variants BAF from sequencing OR normalized LRR and BAF from SNP array. Please make sure all the input are in list form.

NGS only CNV detection using iCNV. As we showed in the paper, we highly recommend using iCNV for WGS CNV detection.

    ### 2.6 code chunk number 1: NGS only iCNV_detection
    str(ngs_plr) # List of n vector, each one is the PLR for an exon
    str(ngs_plr.pos) # List of n matrix (p x 2), each one is the start and end location for an exon
    str(ngs_baf) # List of n vector, each one is the variants BAF from .bam
    str(ngs_baf.pos) # List of n vector, each one is the variants BAF position
    projectname <- 'icnv.demo.ngs.'
    icnv_res0_ngs <- iCNV_detection(ngs_plr=ngs_plr, ngs_baf = ngs_baf,
                                 ngs_plr.pos = ngs_plr.pos,ngs_baf.pos = ngs_baf.pos,
                                 projname=projectname,CN=0,mu=c(-3,0,2),cap=TRUE,visual = 2)
    icnv.output <- output_list(icnv_res0_ngs,sampname_qc,CN=0)
    head(icnv.output)
    plotHMMscore(icnv_res0_ngs,title=projectname)

SNP array only CNV detection using iCNV

    ### 2.6 code chunk number 2: array only iCNV_detection
    str(snp_lrr) # List of n vector, each one is the normalized LRR for a SNP
    str(snp_lrr.pos) # List of n vector, each one is a SNP position
    str(snp_baf) # List of n vector, each one is the BAF for a SNP
    str(snp_baf.pos) # List of n vector, each one is the SNP BAF position
    projectname <- 'icnv.demo.snp'
    icnv_res0_snp <- iCNV_detection(snp_lrr=snp_lrr, snp_baf = snp_baf,
                                 snp_lrr.pos = snp_lrr.pos,snp_baf.pos = snp_baf.pos,
                                 projname=projectname,CN=1,mu=c(-3,0,2),cap=TRUE,visual = 2)
    icnv.output <- output_list(icnv_res0_snp,sampname_qc,CN=1)
    head(icnv.output)
    plotHMMscore(icnv_res0_snp,title=projectname)

## 3. Session information

    sessionInfo()

    ## R version 4.5.2 (2025-10-31)
    ## Platform: x86_64-pc-linux-gnu
    ## Running under: Ubuntu 24.04.3 LTS
    ##
    ## Matrix products: default
    ## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
    ## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
    ##
    ## locale:
    ##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
    ##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C
    ##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
    ##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
    ##  [9] LC_ADDRESS=C               LC_TELEPHONE=C
    ## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
    ##
    ## time zone: Etc/UTC
    ## tzcode source: system (glibc)
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] iCNV_1.31.0                       WES.1KG.WUGSC_1.43.0
    ##  [3] CODEX_1.43.0                      BSgenome.Hsapiens.UCSC.hg19_1.4.3
    ##  [5] BSgenome_1.79.1                   rtracklayer_1.71.3
    ##  [7] BiocIO_1.21.0                     GenomeInfoDb_1.47.2
    ##  [9] Rsamtools_2.27.0                  Biostrings_2.79.4
    ## [11] XVector_0.51.0                    GenomicRanges_1.63.1
    ## [13] IRanges_2.45.0                    S4Vectors_0.49.0
    ## [15] Seqinfo_1.1.0                     BiocGenerics_0.57.0
    ## [17] generics_0.1.4
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] dotCall64_1.2               SummarizedExperiment_1.41.0
    ##  [3] spam_2.11-3                 rjson_0.2.23
    ##  [5] xfun_0.56                   bslib_0.10.0
    ##  [7] Biobase_2.71.0              lattice_0.22-7
    ##  [9] vctrs_0.7.1                 tools_4.5.2
    ## [11] bitops_1.0-9                curl_7.0.0
    ## [13] parallel_4.5.2              tibble_3.3.1
    ## [15] pkgconfig_2.0.3             Matrix_1.7-4
    ## [17] data.table_1.18.2.1         RColorBrewer_1.1-3
    ## [19] cigarillo_1.1.0             truncnorm_1.0-9
    ## [21] lifecycle_1.0.5             compiler_4.5.2
    ## [23] fields_17.1                 codetools_0.2-20
    ## [25] maps_3.4.3                  htmltools_0.5.9
    ## [27] sys_3.4.3                   buildtools_1.0.0
    ## [29] sass_0.4.10                 RCurl_1.98-1.17
    ## [31] yaml_2.3.12                 pillar_1.11.1
    ## [33] crayon_1.5.3                jquerylib_0.1.4
    ## [35] tidyr_1.3.2                 BiocParallel_1.45.0
    ## [37] cachem_1.1.0                DelayedArray_0.37.0
    ## [39] abind_1.4-8                 tidyselect_1.2.1
    ## [41] digest_0.6.39               purrr_1.2.1
    ## [43] dplyr_1.1.4                 restfulr_0.0.16
    ## [45] maketools_1.3.2             fastmap_1.2.0
    ## [47] grid_4.5.2                  cli_3.6.5
    ## [49] SparseArray_1.11.10         magrittr_2.0.4
    ## [51] S4Arrays_1.11.1             XML_3.99-0.20
    ## [53] withr_3.0.2                 UCSC.utils_1.7.1
    ## [55] rmarkdown_2.30              httr_1.4.7
    ## [57] matrixStats_1.5.0           evaluate_1.0.5
    ## [59] knitr_1.51                  viridisLite_0.4.2
    ## [61] rlang_1.1.7                 Rcpp_1.1.1
    ## [63] glue_1.8.0                  jsonlite_2.0.0
    ## [65] R6_2.6.1                    MatrixGenerics_1.23.0
    ## [67] GenomicAlignments_1.47.0