--- title: "ELViS Vignette" author: - name: Hyo Young Choi affiliation: Department of Preventive Medicine, Division of Biostatistics, University of Tennessee Health Science Center, Memphis, TN - name: Jin-Young Lee affiliation: Center for Cancer Research, University of Tennessee Health Science Center, Memphis, TN - name: Jeremiah R. Holt affiliation: Department of Medicine, College of Medicine, University of Tennessee Health Science Center, Memphis, TN - name: Katherine A. Hoadley affiliation: Lineberger Comprehensive Cancer Center, University of North Carolina, Chapel Hill, NC - name: D. Neil Hayes affiliation: Department of Preventive Medicine, Division of Biostatistics, University of Tennessee Health Science Center, Memphis, TN output: BiocStyle::html_document: toc_float: true BiocStyle::pdf_document: default vignette: | %\VignetteIndexEntry{Authoring R Markdown vignettes} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} library(knitr) knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # 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. ## 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. ------------------------------------------------------------------------ # Installation To install this package, start R (version "4.5") and enter: ```{r , echo=TRUE, eval=FALSE} if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("ELViS") ``` # Run Example ## Data preparation ### 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 ``` ### BAM Files \label{section-2-2-1} 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](#section-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. ## Generate Raw Read Depth Matrix with Toy Examples Load required libraries. ```{r setup} library(ELViS) library(ggplot2) library(glue) library(dplyr) library(ComplexHeatmap) theme_set(theme_bw()) ``` Prepare BAM file name vector. ```{r} 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. ```{r} 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_reticulate" ,N_cores = N_cores ,min_mapq = 30 ,tmpdir=tmpdir ,condaenv = "env_samtools" ,condaenv_samtools_version="auto" #,condaenv_samtools_version="1.22" ) }) ``` ```{r eval=FALSE, eval=FALSE, include=FALSE} # Actual depth matrix of 120 samples had been generated using the follpwing code SKIP=1 if(!SKIP){ # 120 bams # 5.5 sec system.time({ mtrx_Rsamtools = get_depth_matrix( bam_files = bam_files,coord_or_target_virus_name = target_virus_name ,mode = "Rsamtools" ,N_cores = N_cores ,max_depth = 1e5 ,min_mapq = 30 ,tmpdir="tmpdir" ) }) } # 4.7 sec system.time({ mtrx_samtools_reticulate = get_depth_matrix( bam_files = bam_files,coord_or_target_virus_name = target_virus_name ,mode = "samtools_reticulate" ,N_cores = N_cores ,min_mapq = 30 ,tmpdir="tmpdir" ,condaenv = "env_samtools" ,condaenv_samtools_version="auto" #,condaenv_samtools_version="1.22" ) }) SKIP=1 if(!SKIP){ #2.3 sec system.time({ mtrx_samtools = get_depth_matrix( bam_files = bam_files,coord_or_target_virus_name = target_virus_name ,mode = "samtools_custom" ,N_cores = N_cores ,min_mapq = 30 ,tmpdir="tmpdir" ,modules="module/file/example/1.2" # name of modulefile the user want to use ) }) } use_data(mtrx_samtools_reticulate) ``` ## Filtering Out Low Depth Samples Determine sample filtering threshold using histogram and filter out low depth samples ```{r,fig.width=5,fig.height=3} # 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) depth_hist(mtrx_samtools_reticulate,th=th,smry_fun=quantile,prob=0.75) # filtered matrix base_resol_depth = filt_samples(mtrx_samtools_reticulate,th=th,smry_fun=max) print(base_resol_depth[1:4,1:4]) ``` ## Run ELViS using the Filtered Depth Matrix Running ELViS using the filtered read depth matrix(`base_resol_depth`). ```{r , echo=TRUE, eval=FALSE} 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 ``` ## Plotting Figures Prepare plotting data ```{r,fig.width=7,fig.height=5} # 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. ```{r,fig.width=7,fig.height=5} # 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"), ) # 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. ```{r,fig.width=7,fig.height=5} # 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. ```{r,fig.width=7,fig.height=5} # 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. ```{r,fig.width=7,fig.height=5} # 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. ```{r eval=FALSE, include=FALSE,echo = FALSE} set.seed(54374373); total_aligned_base__host_and_virus = c( sample( (4:6)*(10^8),80,replace = TRUE), sample( (7:10)*(10^8),20,replace = TRUE), sample( (1:3)*(10^8),20,replace = TRUE) ) use_data(total_aligned_base__host_and_virus,overwrite = TRUE) ``` ```{r,fig.width=5,fig.height=4} 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. ```{r} 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 ) # 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`) ```{r} 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 ) ``` Generate final heatmap in a single panel. ```{r,fig.width=7,fig.height=8} draw(top_ant %v% integ_ht_result$Heatmap %v% gene_cn$Heatmaps$intact_gene_cn %v% gene_cn$Heatmaps$rel_dosage) ``` ------------------------------------------------------------------------ # sessionInfo ```{r} sessionInfo() ```