cfdnakit vignette
Introduction
This package provides basic functions for analyzing next-generation sequencing data of cell-free DNA (cfDNA). The package focuses on extracting the length of cfDNA sample and visualization of genome-wide enrichment of short-fragmented cfDNA. The aberration of fragmentation profile, denoted modified ctDNA estimation score (CES), allows quantification of circulating tumor DNA (ctDNA). We recommend using this package to analysis shallow whole-genome sequencing data (~0.3X or more). This package was complemented by Bioconductor packages e.g. QDNAseq, Rsamtools and GenomicRanges which could further expand the functionality of this package in the future.
Installation
Install via the Bioconductor repository
Install the latest version via github
To install this package is via this github repository, please follow the instruction below.
Install prerequisites packages
if(! "devtools" %in% rownames(installed.packages()))
install.packages("devtools")
if(! "BiocManager" %in% rownames(installed.packages()))
install.packages("BiocManager")Install cfdnakit package
The installation should work fine without non-zero exit status. Try load cfdnakit package into current R session
Prepare input BAM
A coordination-sorted BAM file of cfDNA from any liquid biopsy source (e.g. blood-plasma, CSF, urine) should be applicable. You should examine if the sequencing coverage reach coverage threshold before the analysis. Please make sure that the sequencing reads are mapped onto the following version of reference genome. cfdnakit supports the human reference genome GRCh37(hg19) and GRCh38(hg38). Preliminary fragment-length analysis can be performed using mouse cfDNA when mapping onto the mouse reference genome GRCm38(mm10).
Read the BAM file with read_bamfile
Let’s read sequence alignment file (.bam) using function read_bamfile. A BAM index file (.bai) is necessary for efficiently reading the file. If it doesn’t exist in the same directory, this function will automatically create one. This function will split sequencing reads into equal-size non-overlapping windows. Possible size of bin are 100, 500, and 1000 KB. A path to the input file is given to the function read_bamfile. A SampleBam object will be created as the result.
For demonstration, we read an example sequence file ex.plasma.bam.
library(cfdnakit)
sample_bamfile <- system.file("extdata",
"ex.plasma.bam",
package = "cfdnakit")
plasma_SampleBam <- read_bamfile(sample_bamfile,
apply_blacklist = FALSE)
#> The BAM index file (.bai) is missing. Creating an index file
#> Bam index file created.
#> Reading bamfileBy default, running read_bamfile will split reads into 1000 KB non-overlapping bins based-on the human reference genome GRCh37.
Reading the file should take a while depending on the size of your BAM. We recommend to save the result as RDS file using saveRDS function.
Analyse the Fragment Length Distribution
Fragment-length distribution of the sample can be visualized with function plot_fragment_dist. In the top-right legend, the modal length (bp) is shown in the parenthesis behind each sample name. The x-axis is the length of cfDNA; y-axis is the proportion of cfDNA fragment having a specific length.
This document will demonstate by using SampleBam object (example_patientcfDNA_SampleBam.RDS). Now we load this file into R session.
example_RDS <- "example_patientcfDNA_SampleBam.RDS"
example_RDS_file <-
system.file("extdata",example_RDS,
package = "cfdnakit")
sample_bins <- readRDS(example_RDS_file)In general, plasma cfDNA should show non-random fragmentation pattern. The modal length of cfDNA is 167 bp and 10-bp periodic peaks. However, tumor-derived cfDNA are observed to be shorter (<150 bp) than cfDNA of non-tumor origin. Here, we compare the fragment length distribution of patient’s cfDNA with healthy individual. We derived a healthy cfDNA from Snyder et al. (2016) and create a RData file “BH01_chr15.RDS”. This file can be loaded in R environment with readRDS function.
control_rds<-"BH01_CHR15.SampleBam.rds"
control_RDS_file <-
system.file("extdata",control_rds,
package = "cfdnakit")
control_bins <-
readRDS(control_RDS_file)To provide visual comparison of cell-free DNA fragmentation, cfdnakit provide a function that allows plotting multiple distribution plots. We create a list of SampleBam files and plot their distribution together with function plot_fragment_dist. Each element in the list must be given a distinct sample name (e.g. Healthy.cfDNA).
comparing_list <- list("Healthy.cfDNA"=control_bins,
"Patient.1"=sample_bins)
plot_fragment_dist(comparing_list)
Optional Parameters
maximum_length = Maximum length of cfDNA (Default 550)
minimum_length = Minimum length of cfDNA (Default 20)
Quantification of Short Fragmented CfDNA
We can extract genome-wide fragment-length profile. First we define the short fragment as fragments having size (by default) between 100 - 150 base and 151 - 250 as the long fragment. Providing a SampleBam object to the function get_fragment_profile will return a SampleFragment object.
This SampleFragment contains a dataframe named sample_profile. This table contains information about the BAM file and the ratio number of Short/Long fragments. The table below describes important variables.
sample_profile$sample_profile
#> Sample.ID Total.Fragments Read.Pairs.in.range Read.Pairs.in.range..corrected.
#> 1 Patient.1 68094 65421 62091.41
#> N.Short N.Long Mode Median Mean Mad S.L.Ratio S.L.Ratio_corrected
#> 1 23503 41918 167 159 163.23 20.76 0.56 0.56
#> Bin.Size.KB.
#> 1 1000| Variable | Description |
|---|---|
| Total.Fragments | Number of DNA fragments (not the number of reads) |
| Read.Pairs.in.range | Number of fragments within the defined range (short and long) |
| Mode | Fragment length of the majority (bp) |
| Median | Median Fragment length (bp) |
| Mean | Average Fragment length (bp) |
| Mad | Fragment length median absolute deviation (MAD) |
| N.Short | Number of short fragments (n) |
| N.Long | Number of long fragments (n) |
| S.L.Ratio | N.Short/N.Long; Ratio of short fragment over long fragment |
| S.L.Ratio_corrected | Ratio of short fragment over long fragment after GC-bias correction |
| Bin.Size.KB. | Size of genomic bin of the analysis (KB) |
Plot Genome-wide Short-fragmented Ratio
We can plot genome-wide short-fragment ratio with the function plot_sl_ratio. Given short-fragment profile, short-fragment ratio per bin infer contribution of circulating tumor DNA (ctDNA) into cfDNA. The enrichment of short-fragment cfDNA in a large genomic region could infer the copy-number aberration status. The higher short-fragment ratio indicate amplification event whereas deletion would have relatively lower short-fragment ratio.
## For this demenstration, we load a real patient-derived cfDNA profile.
patient.SampleFragment.file <-
system.file("extdata",
"example_patientcfDNA_SampleFragment.RDS",
package = "cfdnakit")
patient.SampleFragment <- readRDS(patient.SampleFragment.file)
plot_sl_ratio(patient.SampleFragment)
The range of S.L.Ratio is broad. It could be as high as 2 or more in sample with high contribution of ctDNA. Adjusting the plot is possible with parameter ylim or with other ggplot functions.
Save SampleFragment as RDS file
Up to this step, we rather save the SampleFragment object as RDS for later use and for creation of a Panel-of-Normal.
Create Panel of Normal (PoN)
Creating list of PoN files
Making a Panel-of-Normal is necessary for downstream analysis as we want to compare fragment profile between a cfDNA from patient with pooled of healthy individuals. First, we create a text file where each line is a full path to rds file of the SampleFragment object.
Example content of PoN list file (let’s name it Pon-list.txt)
### Example content of Pon-list.txt
Path.to/SampleFragment_healthyplasma-01.rds
Path.to/SampleFragment_healthyplasma-02.rds
Path.to/SampleFragment_healthyplasma-03.rds
Path.to/SampleFragment_healthyplasma-04.rds
Path.to/SampleFragment_healthyplasma-05.rdsInferring CNV from short fragment cfDNA
Normalizing Short-fragmented Ratio
The contribution of short-fragmented cfDNA into genomic regions could infer copy-number aberration harbored in the tumor cells. For this demonstration, we load ready-to-use PoN profile (ex.PoN.rds).
PoN_rdsfile <- system.file("extdata",
"ex.PoN.rds",
package = "cfdnakit")
PoN.profiles <- readRDS(PoN_rdsfile)cfdnakit transforms S.L.Ratio of each genomic windows with the distribution of samples in PoN by function get_zscore_profile.
Circular Binary Segmentation (CBS)
cfdnakit implement the circular binary segmentation algorithm. Function segmentByCBS of package PSCBS were used through function segmentByPSCB. We can visualize both transformed S.L.Ratio and the segmentation by function plot_transformed_sl.
sample_zscore_segment <- segmentByPSCB(sample_zscore)
plot_transformed_sl(sample_zscore,sample_zscore_segment)
Optional parameters ylim = Y-axis of plot (Default c(-20,20))
Estimating tumor fraction (TF) and CNV calling
Base on hypothesis that short-fragmented cfDNA originate from tumor cells. Enrichment of short-fragmented cfDNA in a large genomic segments correlates with the absolute copy in the tumor origin. cfdnakit heuristically search for fittest solution per absolute ploidy of 2, 3 and 4. The distance of each solution is calculated similarly to ACE (Poell JB et.al 2019) method.
Function call_cnv returns 3 solutions of the CNV calling and TF estimation. The analysis perform only on autosomal chromosomes. We can plot the distance matrix of all solution where solutions with * are optimal solutions (minimum distance) per absolute genome ploidy.
Function call_cnv returns 3 solutions of the CNV calling and TF estimation. The analysis perform only on autosomal chromosomes. We can plot the distance matrix of all solution where solutions with * are optimal solutions (minimum distance) per absolute genome ploidy.
The available solution can be obtained with function get_solution_table. The table shows the best solution (minimum distance) per absolute genome ploidy and rank them by the distance.
Copy-number Abnormality Score
As the result of copy-number solution fitting, the tumor fraction (tf ) indicates the estimated quantity of ctDNA from the amplitude of signals. cfdnakit also implements the ctDNA estimation score (CES) score to quantify the tumor-derived cell-free DNA from the segmentation and the SLRatio.
This score is robust to coverage bias and noisy fragmented signals. Briefly, the Gaussian noise does not affect the score because the z-scores of segments, instead of the z-score of bins, are considered. Second, the average segment length is used as a penalty for sample quality. The signal of a bad quality sample does not strongly affect the score whereas a true highly unstable genome would overcome this penalty.
Session info
Output of sessionInfo on the system on which this document was compiled:
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.1 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] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] DNAcopy_1.81.0 cfdnakit_1.5.0 BiocStyle_2.35.0
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.6 xfun_0.49 bslib_0.8.0
#> [4] ggplot2_3.5.1 Biobase_2.67.0 vctrs_0.6.5
#> [7] tools_4.4.2 bitops_1.0-9 generics_0.1.3
#> [10] stats4_4.4.2 parallel_4.4.2 tibble_3.2.1
#> [13] fansi_1.0.6 pkgconfig_2.0.3 R.oo_1.27.0
#> [16] PSCBS_0.67.0 S4Vectors_0.45.0 lifecycle_1.0.4
#> [19] GenomeInfoDbData_1.2.13 R.cache_0.16.0 farver_2.1.2
#> [22] compiler_4.4.2 Rsamtools_2.23.0 Biostrings_2.75.1
#> [25] munsell_0.5.1 codetools_0.2-20 GenomeInfoDb_1.43.0
#> [28] htmltools_0.5.8.1 sys_3.4.3 buildtools_1.0.0
#> [31] sass_0.4.9 yaml_2.3.10 pillar_1.9.0
#> [34] crayon_1.5.3 jquerylib_0.1.4 R.utils_2.12.3
#> [37] BiocParallel_1.41.0 cachem_1.1.0 parallelly_1.39.0
#> [40] tidyselect_1.2.1 digest_0.6.37 future_1.34.0
#> [43] dplyr_1.1.4 listenv_0.9.1 labeling_0.4.3
#> [46] maketools_1.3.1 fastmap_1.2.0 grid_4.4.2
#> [49] colorspace_2.1-1 cli_3.6.3 magrittr_2.0.3
#> [52] utf8_1.2.4 withr_3.0.2 scales_1.3.0
#> [55] UCSC.utils_1.3.0 aroma.light_3.37.0 rmarkdown_2.29
#> [58] XVector_0.47.0 httr_1.4.7 matrixStats_1.4.1
#> [61] globals_0.16.3 R.methodsS3_1.8.2 evaluate_1.0.1
#> [64] knitr_1.49 GenomicRanges_1.59.0 IRanges_2.41.0
#> [67] rlang_1.1.4 glue_1.8.0 BiocManager_1.30.25
#> [70] BiocGenerics_0.53.1 jsonlite_1.8.9 R6_2.5.1
#> [73] zlibbioc_1.52.0- Introduction
- Prepare input BAM
- Read the BAM file with read_bamfile
- Analyse the Fragment Length Distribution
- Quantification of Short Fragmented CfDNA
- Plot Genome-wide Short-fragmented Ratio
- Save SampleFragment as RDS file
- Create Panel of Normal (PoN)
- Inferring CNV from short fragment cfDNA
- Copy-number Abnormality Score
- Session info