1 Introduction

Tumor viruses cause a significant proportion of human cancers by hijacking cell signaling and inducing proliferation. For example, HPV E6 and E7 proteins can promote tumorigenesis by inducing degradation of tumor suppressors p53 and pRb. Despite the importance of the dosages of the viral genes like these, existing tools cover viral integration, emphasizing the need for viral copy number analysis framework.

1.1 Motivation for submitting to Bioconductor

ELViS (Estimating Copy Number Levels of Viral Genomic Segments) is an R package that addresses this need through the analysis of copy numbers within viral genomes. ELViS utilizes base-resolution read depth data over viral genome to find copy number segments with two-dimensional segmentation approach. It also provides publish-ready figures, including histograms of read depths for sample filtering, coverage line plots over viral genome annotated with copy number change events and viral genes, and heatmaps showing overall pictures with multiple types of data with integrative clustering of samples. Taken together, our framework helps computational biologists and bioinformaticians to analyze viral copy number changes with ease.

2 Installation

To install this package, start R (version “4.5”) and enter:

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("ELViS")

3 Run Example

3.1 Data preparation

3.1.1 Virus Reference Genome File

The following reference files are included in extdata directory of ELViS.

  • HPV16.fa
  • HPV16.fa.fai
  • HPV16REF_PaVE.gff

HPV16.fa is a reference sequence file with the associated fasta index file HPV16.fa.fai HPV16REF_PaVE.gff are viral gene annotation file in GFF3 format. They were downloaded from PaVE database (https://pave.niaid.nih.gov/locus_viewer?seq_id=HPV16REF)

You can get your own files(.fa,.fai,.gff) for the virus of your interest.

But, you should check if this viral reference genome FASTA(.fa) file and gene annotation GFF3(.gff) file have the same set of sequence names.

For the FASTA file, we recommend you add viral genome files to the host genome file (i.e. human reference genome hg38) to simultaneously align your reads to both host and viral genome. Multiple viral genome can be added. The reason we recommend this is as such. If there is a significant portion of viral reads in the sequencing data, and you align reads to either host or viral genome, there is a higher chance of false alignment and time it takes to finish the read alignment step can be extended.

# linux bash
cat hg38.fa virus1.fa virus2.fa > hg38_with_virus.fa

Even if you use FASTA file containing both viral genomes and host genome, ELViS require a separate GFF file for each virus. Don’t merge GFF3 files together for using ELViS.

FASTA index file(.fai) can be created by running samtools faidx. For more detail, please refer to https://www.htslib.org/doc/samtools-faidx.html

# linux bash
samtools faidx hg38_with_virus.fa

3.1.2 BAM Files

The following BAM files are included in extdata directory of ELViS.

  • Control_100X_1102.bam
  • Control_100X_1102.bam.bai
  • Control_100X_1119.bam
  • Control_100X_1119.bam.bai

These are simulation data made using w-WESSIM-2. (https://github.com/GeorgetteTanner/w-Wessim2) Briefly, sequencing reads were simulated using this tool and aligned to HPV16.fa to generate the BAM file.

For your case, you should align your sequencing reads in FASTQ format to the reference genome prepared according to instructions in 2.2.1. As mentioned, we recommend aligning reads to a FASTQ file contains host genome and viral genomes.

However, you may have discovered our package after aligning all the reads to the host genome only already. No worries since there is a work around.

In that case, you can extract the unaligned reads and align that to viral genome FASTA file.

For example, say you have a BAM file seq_align_to_host.bam that is aligned to host genome only.

Extract unaligned reads using samtools as follows.

samtools view seq_align_to_host.bam --incl-flags 12 -b | # extract reads that are unaligned or those with unaligned mates
samtools sort -n | # sort reads by read name
samtools fastq -1 seq_r1.fastq -2 seq_r2.fastq # save as read1 and read2 files.

Align seq_r1.fastq and seq_r2.fastq to merged reference of viral genomes using sequence aligners of your choice, such as BWA MEM or Bowtie2.

3.2 Generate Raw Read Depth Matrix with Toy Examples

Load required libraries.

library(ELViS)
library(ggplot2)
library(glue)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(ComplexHeatmap)
#> Loading required package: grid
#> ========================================
#> ComplexHeatmap version 2.25.0
#> Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
#> Github page: https://github.com/jokergoo/ComplexHeatmap
#> Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
#> 
#> If you use it in published research, please cite either one:
#> - Gu, Z. Complex Heatmap Visualization. iMeta 2022.
#> - Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
#>     genomic data. Bioinformatics 2016.
#> 
#> 
#> The new InteractiveComplexHeatmap package can directly export static 
#> complex heatmaps into an interactive Shiny app with zero effort. Have a try!
#> 
#> This message can be suppressed by:
#>   suppressPackageStartupMessages(library(ComplexHeatmap))
#> ========================================
theme_set(theme_bw())

Prepare BAM file name vector.

analysis_dir = tempdir()
dir.create(analysis_dir,showWarnings = FALSE)

package_name = "ELViS"

# load toy example meta data
data(toy_example,package = package_name)

# get lust of bam file paths
ext_path = system.file("extdata",package = package_name)
bam_files = list.files(ext_path,full.names = TRUE,pattern = "bam$")

Generate base-resolution read depth matrix from a list of BAM files. Parallel package is used to read BAM files fast.

os_name = Sys.info()["sysname"]
if( os_name == "Windows" ){
    N_cores <- 1L
}else{
    N_cores <- 2L
}


# the name of the reference viral sequence the reads were aligned to
target_virus_name = "gi|333031|lcl|HPV16REF.1|"

# temporary file directory
tmpdir=tempdir()

# generate read depth matrix
system.time({
    mtrx_samtools_reticulate__example = 
        get_depth_matrix(
            bam_files = bam_files,coord_or_target_virus_name = target_virus_name
            ,mode = "samtools_basilisk"
            ,N_cores = N_cores
            ,min_mapq = 30
            ,tmpdir=tmpdir
            ,condaenv = "env_samtools"
            ,condaenv_samtools_version="auto"
        )
})
#> Warning in check_mode_os(mode): mode=samtools_basilisk is not supported for
#> Windows. Changing mode to Rsamtools...
#>    user  system elapsed 
#>    0.25    0.00    0.25

3.3 Filtering Out Low Depth Samples

Determine sample filtering threshold using histogram and filter out low depth samples

# loading precalculated depth matrix
data(mtrx_samtools_reticulate)

# threshold
th = 50
# histogram with adjustable thresholds for custom function
depth_hist(mtrx_samtools_reticulate,th=th,smry_fun=max)
#> Warning in scale_x_continuous(trans = "log10"): log-10 transformation
#> introduced infinite values.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

depth_hist(mtrx_samtools_reticulate,th=th,smry_fun=quantile,prob=0.75)
#> Warning in scale_x_continuous(trans = "log10"): log-10 transformation
#> introduced infinite values.
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.


# filtered matrix
base_resol_depth = filt_samples(mtrx_samtools_reticulate,th=th,smry_fun=max)
print(base_resol_depth[1:4,1:4])
#>      Control_100X_58.bam Control_100X_61.bam Control_100X_64.bam
#> [1,]                  55                  57                  60
#> [2,]                  56                  59                  62
#> [3,]                  57                  61                  62
#> [4,]                  57                  61                  63
#>      Control_100X_83.bam
#> [1,]                  49
#> [2,]                  49
#> [3,]                  52
#> [4,]                  53

3.4 Run ELViS using the Filtered Depth Matrix

Running ELViS using the filtered read depth matrix(base_resol_depth).

system.time({
    result = run_ELViS(
        X = base_resol_depth
        ,N_cores=N_cores
        ,reduced_output=TRUE
    )
    
})

ELViS_toy_run_result = result
use_data(ELViS_toy_run_result)

# 4min for 120 samples using 10 threads

3.5 Plotting Figures

Prepare plotting data

# ELViS run result
data(ELViS_toy_run_result)
result = ELViS_toy_run_result

# Directory where figures will be saved
figure_dir = glue("{analysis_dir}/figures")
dir.create(figure_dir)

# give the gff3 file of the virus of your interest. Sequence name or chromosome name should match with that in the reference genome FASTA file.
gff3_fn = system.file("extdata","HPV16REF_PaVE.gff",package = package_name)

Raw read depth profile line plots.

# Plotting raw depth profile
gg_lst_x = 
    plot_pileUp_multisample(
        result = result,
        X_raw = base_resol_depth,
        plot_target = "x",
        gff3 = gff3_fn,
        baseline=1,
        exclude_genes = c("E6*","E1^E4","E8^E2"),
    )
#> Import genomic features from the file as a GRanges object ...
#> Warning in .local(con, format, text, ...): gff-version directive indicates
#> version is 3.1.26, not 3
#> OK
#> Prepare the 'metadata' data frame ... OK
#> Make the TxDb object ...
#> Warning in .extract_transcripts_from_GRanges(tx_IDX, gr, mcols0$type, mcols0$ID, : the transcript names ("tx_name" column in the TxDb object) imported from the
#>   "Name" attribute are not unique
#> Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information
#> is not available for this TxDb object
#> OK

# Save to pdf file, set SKIP = FALSE if you want to save as pdf
SKIP = TRUE
if(!SKIP){
    pdf(glue("{figure_dir}/Raw_Depth_CNV_call.pdf"),height=4,width=6)
    gg_lst_x
    dev.off()
}

# an example of raw read depth line plot
print(gg_lst_x[[1]])

You can adjust baseline after examining depth profile plots.

# set the longest segment as a new baseline
new_baseline = get_new_baseline(result,mode="longest")

# Plotting raw depth profile with new baseline
gg_lst_x = 
    plot_pileUp_multisample(
        result = result,
        X_raw = base_resol_depth,
        plot_target = "x",
        gff3 = gff3_fn,
        baseline=new_baseline,
        exclude_genes = c("E6*","E1^E4","E8^E2"),
    )
# Save to pdf file, set SKIP = FALSE if you want to save as pdf
SKIP = TRUE
if(!SKIP){
    # Save to pdf file
    pdf("figures/Raw_Depth_new_baseline_CNV_call.pdf",height=4,width=6)
    gg_lst_x
    dev.off()
}

# an example of raw read depth line plot with new baseline
gg_lst_x[[1]]

Normalized read depth profile line plots.

# Plotting normalized depth profile
gg_lst_y = 
    plot_pileUp_multisample(
        result = result,
        X_raw = base_resol_depth,
        plot_target = "y",
        gff3 = gff3_fn,
        baseline=new_baseline,
        exclude_genes = c("E6*","E1^E4","E8^E2"),
    )

# Save to pdf file
SKIP = TRUE
if(!SKIP){
    pdf("figures/Normalized_Depth_CNV_call.pdf",height=4,width=6)
    gg_lst_y
    dev.off()
}

# an example of normalized read depth line plot with new baseline
gg_lst_y[[1]]

Robust Z-score profile line plots.

# Plotting robust Z-score profile
gg_lst_z = 
    plot_pileUp_multisample(
        result = result,
        X_raw = base_resol_depth,
        plot_target = "z",
        gff3 = gff3_fn,
        baseline=new_baseline,
        exclude_genes = c("E6*","E1^E4","E8^E2")
    )

SKIP = TRUE
if(!SKIP){
    # Save to pdf file
    pdf("figures/Robust-Z-score_CNV_call.pdf",height=4,width=6)
    gg_lst_z
    dev.off()
}

# an example of Z-score line plot with new baseline
gg_lst_z[[1]]

Generating heatmaps with integrative clustering.

Calculation of viral loads. - Get total aligned base using tools such as picard. Here we use randomly generated numbers instead.

data(total_aligned_base__host_and_virus)

viral_load = (10^6)*(apply(base_resol_depth,2,\(x) sum(x)) )/total_aligned_base__host_and_virus

# distribtuion of overall viral load
viral_load %>%log10 %>% hist

Generate heatmaps with integrative clustering using data transformed in various ways.

exclude_genes = c("E6*","E1^E4","E8^E2")
integ_ht_result = integrative_heatmap(
    X_raw = base_resol_depth,
    result = result,
    gff3_fn = gff3_fn,
    exclude_genes = exclude_genes,
    baseline=1,
    total_aligned_base__host_and_virus = total_aligned_base__host_and_virus
)
#> Import genomic features from the file as a GRanges object ...
#> Warning in .local(con, format, text, ...): gff-version directive indicates
#> version is 3.1.26, not 3
#> OK
#> Prepare the 'metadata' data frame ... OK
#> Make the TxDb object ...
#> Warning in .extract_transcripts_from_GRanges(tx_IDX, gr, mcols0$type, mcols0$ID, : the transcript names ("tx_name" column in the TxDb object) imported from the
#>   "Name" attribute are not unique
#> Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information
#> is not available for this TxDb object
#> OK
#> Import genomic features from the file as a GRanges object ...
#> Warning in .local(con, format, text, ...): gff-version directive indicates
#> version is 3.1.26, not 3
#> OK
#> Prepare the 'metadata' data frame ... OK
#> Make the TxDb object ...
#> Warning in .extract_transcripts_from_GRanges(tx_IDX, gr, mcols0$type, mcols0$ID, : the transcript names ("tx_name" column in the TxDb object) imported from the
#>   "Name" attribute are not unique
#> Warning in .extract_transcripts_from_GRanges(tx_IDX, gr, mcols0$type, mcols0$ID, : genome version information is not available for this TxDb object
#> OK
#> `use_raster` is automatically set to TRUE for a matrix with more than
#> 2000 rows. You can control `use_raster` argument by explicitly setting
#> TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> `use_raster` is automatically set to TRUE for a matrix with more than
#> 2000 rows. You can control `use_raster` argument by explicitly setting
#> TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> `use_raster` is automatically set to TRUE for a matrix with more than
#> 2000 rows. You can control `use_raster` argument by explicitly setting
#> TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> `use_raster` is automatically set to TRUE for a matrix with more than
#> 2000 rows. You can control `use_raster` argument by explicitly setting
#> TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> `use_raster` is automatically set to TRUE for a matrix with more than
#> 2000 rows. You can control `use_raster` argument by explicitly setting
#> TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> Warning: The input is a data frame-like object, convert it to a matrix.
#> `use_raster` is automatically set to TRUE for a matrix with more than
#> 2000 rows. You can control `use_raster` argument by explicitly setting
#> TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> Warning: The heatmap has not been initialized. You might have different results
#> if you repeatedly execute this function, e.g. when row_km/column_km was
#> set. It is more suggested to do as `ht = draw(ht); column_order(ht)`.

# top annotation
top_ant =
    HeatmapAnnotation(
        `Log2 Overall\nViral Load` = anno_points(log2(viral_load)),
        annotation_name_side = "left",annotation_name_rot=0)

Generate heatmap showing maximum number of intact copies - min copy of the overlapping copy segments - ratio relative to certain gene(gene_ref)

gene_ref="E7"

gene_cn = 
    gene_cn_heatmaps(
        X_raw = base_resol_depth,
        result = result,
        gff3_fn = gff3_fn,
        baseline = new_baseline,
        gene_ref = gene_ref,
        exclude_genes = exclude_genes
    )
#> Import genomic features from the file as a GRanges object ...
#> Warning in .local(con, format, text, ...): gff-version directive indicates
#> version is 3.1.26, not 3
#> OK
#> Prepare the 'metadata' data frame ... OK
#> Make the TxDb object ...
#> Warning in .extract_transcripts_from_GRanges(tx_IDX, gr, mcols0$type, mcols0$ID, : the transcript names ("tx_name" column in the TxDb object) imported from the
#>   "Name" attribute are not unique
#> Warning in .makeTxDb_normarg_chrominfo(chrominfo): genome version information
#> is not available for this TxDb object
#> OK

Generate final heatmap in a single panel.

draw(top_ant %v% integ_ht_result$Heatmap %v% gene_cn$Heatmaps$intact_gene_cn %v% gene_cn$Heatmaps$rel_dosage)

4 sessionInfo

sessionInfo()
#> R version 4.5.0 (2025-04-11 ucrt)
#> Platform: x86_64-w64-mingw32/x64
#> Running under: Windows Server 2022 x64 (build 20348)
#> 
#> Matrix products: default
#>   LAPACK version 3.12.1
#> 
#> locale:
#> [1] LC_COLLATE=C                          
#> [2] LC_CTYPE=English_United States.utf8   
#> [3] LC_MONETARY=English_United States.utf8
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.utf8    
#> 
#> time zone: America/New_York
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] grid      stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#> [1] ComplexHeatmap_2.25.0 dplyr_1.1.4           glue_1.8.0           
#> [4] ggplot2_3.5.2         ELViS_1.1.1           knitr_1.50           
#> [7] BiocStyle_2.37.0     
#> 
#> loaded via a namespace (and not attached):
#>   [1] DBI_1.2.3                   bitops_1.0-9               
#>   [3] httr2_1.1.2                 biomaRt_2.65.0             
#>   [5] rlang_1.1.6                 magrittr_2.0.3             
#>   [7] clue_0.3-66                 GetoptLong_1.0.5           
#>   [9] matrixStats_1.5.0           compiler_4.5.0             
#>  [11] RSQLite_2.4.1               GenomicFeatures_1.61.3     
#>  [13] dir.expiry_1.17.0           png_0.1-8                  
#>  [15] vctrs_0.6.5                 txdbmaker_1.5.4            
#>  [17] stringr_1.5.1               pkgconfig_2.0.3            
#>  [19] shape_1.4.6.1               crayon_1.5.3               
#>  [21] fastmap_1.2.0               magick_2.8.7               
#>  [23] dbplyr_2.5.0                XVector_0.49.0             
#>  [25] labeling_0.4.3              Rsamtools_2.25.0           
#>  [27] rmarkdown_2.29              UCSC.utils_1.5.0           
#>  [29] tinytex_0.57                bit_4.6.0                  
#>  [31] xfun_0.52                   cachem_1.1.0               
#>  [33] GenomeInfoDb_1.45.4         jsonlite_2.0.0             
#>  [35] progress_1.2.3              blob_1.2.4                 
#>  [37] DelayedArray_0.35.1         uuid_1.2-1                 
#>  [39] BiocParallel_1.43.3         prettyunits_1.2.0          
#>  [41] parallel_4.5.0              cluster_2.1.8.1            
#>  [43] R6_2.6.1                    bslib_0.9.0                
#>  [45] stringi_1.8.7               RColorBrewer_1.1-3         
#>  [47] reticulate_1.42.0           rtracklayer_1.69.0         
#>  [49] GenomicRanges_1.61.0        jquerylib_0.1.4            
#>  [51] Rcpp_1.0.14                 bookdown_0.43              
#>  [53] SummarizedExperiment_1.39.0 iterators_1.0.14           
#>  [55] zoo_1.8-14                  IRanges_2.43.0             
#>  [57] Matrix_1.7-3                igraph_2.1.4               
#>  [59] tidyselect_1.2.1            dichromat_2.0-0.1          
#>  [61] abind_1.4-8                 yaml_2.3.10                
#>  [63] doParallel_1.0.17           codetools_0.2-20           
#>  [65] curl_6.3.0                  lattice_0.22-7             
#>  [67] tibble_3.3.0                withr_3.0.2                
#>  [69] Biobase_2.69.0              KEGGREST_1.49.0            
#>  [71] evaluate_1.0.3              segclust2d_0.3.3           
#>  [73] BiocFileCache_2.99.5        xml2_1.3.8                 
#>  [75] circlize_0.4.16             Biostrings_2.77.1          
#>  [77] pillar_1.10.2               BiocManager_1.30.26        
#>  [79] filelock_1.0.3              MatrixGenerics_1.21.0      
#>  [81] foreach_1.5.2               stats4_4.5.0               
#>  [83] generics_0.1.4              RCurl_1.98-1.17            
#>  [85] hms_1.1.3                   S4Vectors_0.47.0           
#>  [87] scales_1.4.0                tools_4.5.0                
#>  [89] BiocIO_1.19.0               data.table_1.17.4          
#>  [91] GenomicAlignments_1.45.0    XML_3.99-0.18              
#>  [93] Cairo_1.6-2                 AnnotationDbi_1.71.0       
#>  [95] colorspace_2.1-1            patchwork_1.3.0            
#>  [97] basilisk_1.21.5             restfulr_0.0.15            
#>  [99] cli_3.6.5                   rappdirs_0.3.3             
#> [101] viridisLite_0.4.2           S4Arrays_1.9.1             
#> [103] gtable_0.3.6                sass_0.4.10                
#> [105] digest_0.6.37               BiocGenerics_0.55.0        
#> [107] SparseArray_1.9.0           rjson_0.2.23               
#> [109] farver_2.1.2                memoise_2.0.1              
#> [111] htmltools_0.5.8.1           lifecycle_1.0.4            
#> [113] httr_1.4.7                  GlobalOptions_0.1.2        
#> [115] bit64_4.6.0-1