vignettes/RAIDS.Rmd
RAIDS.Rmd
Package: RAIDS
Authors: Pascal Belleau [cre, aut] (https://orcid.org/0000-0002-0802-1071), Astrid Deschênes
[aut] (https://orcid.org/0000-0001-7846-6749), David A. Tuveson
[aut] (https://orcid.org/0000-0002-8017-2712), Alexander
Krasnitz [aut]
Version: 1.3.1
Compiled date: 2024-09-10
License: Apache License (>= 2)
The RAIDS
package and the underlying RAIDS code are
distributed under
the https://opensource.org/licenses/Apache-2.0 license. You
are free to use and redistribute this software.
If you use the RAIDS package for a publication, we would ask you to cite the following:
Pascal Belleau, Astrid Deschênes, Nyasha Chambwe, David A. Tuveson, Alexander Krasnitz; Genetic Ancestry Inference from Cancer-Derived Molecular Data across Genomic and Transcriptomic Platforms. Cancer Res 1 January 2023; 83 (1): 49–58. https://doi.org/10.1158/0008-5472.CAN-22-0682
Multiple methods have been implemented to infer ancestry from germline DNA sequence (Price et al. 2006; Pritchard, Stephens, and Donnelly 2000; Alexander, Novembre, and Lange 2009). However, genotyping of DNA from matched normal specimens is not part of standard clinical practice and is not performed routinely outside academic clinical centers. In sum, matched germline DNA sequence is often missing for cancer-derived molecular data. In such cases, having the possibility to infer ancestry from tumor-derived data would be beneficial.
The RAIDS package implements an inference procedure that has been specifically developed to accurately infer genetic ancestry from cancer-derived sequences. The current version can handle cancer-derived sequences of:
The RAIDS package implements a data synthesis method that, for any given cancer-derived sequence profile, enables on the one hand, profile-specific inference parameter optimization and on the other hand, a profile-specific inference accuracy estimate.
To install this package from Bioconductor, start R (version 4.3 or later) and enter:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("RAIDS")
This is an overview of genetic ancestry inference from cancer-derived molecular data:
The main steps are:
Step 1. Format reference data from the population reference dataset (optional)
Step 2.1 Optimize ancestry inference parameters
Step 2.2 Infer ancestry for the subjects of the external study
These steps are described in detail in the following. Steps 2.1 and 2.2 can be run together using one wrapper function.
A wrapper function encapsulates multiple steps of the workflow.
In summary, the wrapper function generates the synthetic dataset and uses it to selected the optimal parameters before calling the genetic ancestry on the current profiles.
According to the type of input data (RNA or DNA), a specific wrapper function is available.
The wrapper function, called runExomeAncestry(), requires 4 files as input:
In addition, a data.frame containing the general information about the study is also required. The data.frame must contain those 3 columns:
For demonstration purpose, a small population reference GDS file (called ex1_good_small_1KG.gds) and a small population reference SNV Annotation GDS file (called ex1_good_small_1KG_Annot.gds) are included in this package. Beware that those two files should not be used to run a real ancestry inference.The results obtained with those files won’t be reliable.
The required population reference GDS file and population reference SNV Annotation GDS file should be stored in the same directory. In the example below, this directory is referred to as pathReference.
The Profile SNP file can be either in a VCF format or in a generic format.
The Profile SNP VCF file follows the VCF standard with at least those genotype fields: GT, AD and DP. The identifier of the genotype in the VCF file must correspond to the profile identifier Name.ID. The SNVs must be germline variants and should include the genotype of the wild-type homozygous at the selected positions in the reference. One file per profile is need and the VCF file must be gzipped.
Note that the name assigned to the Profile SNP VCF file has to correspond to the profile identifier Name.ID in the following analysis. For example, a SNP file called “Sample.01.vcf.gz” would be associated to the “Sample.01” profile.
A generic SNP file can replace the VCF file. The Profile SNP Generic file format is coma separated and the mandatory columns are:
Beware that the starting position in the population reference GDS File is zero (like BED files). The Profile SNP Generic file should also start at position zero.
Note that the name assigned to the Profile SNP Generic file has to correspond to the profile identifier Name.ID in the following analysis. For example, a SNP file called “Sample.01.generic.txt.gz” would be associated to the “Sample.01” profile.
The Profile PED RDS file must contain a data.frame describing all the profiles to be analyzed. These 5 mandatory columns:
Important: The row names of the data.frame must be the profiles Name.ID.
This file is referred to as the Profile PED RDS file (PED for pedigree). Alternatively, the PED information can be saved in another type of file (CVS, etc..) as long as the data.frame information can be regenerated in R (with read.csv() or else).
This example run an ancestry inference on an exome sample. Both population reference files are demonstration files and should not be used for a real ancestry inference. Beware that running an ancestry inference on real data will take longer to run.
#############################################################################
## Load required packages
#############################################################################
library(RAIDS)
library(gdsfmt)
## Path to the demo 1KG GDS file is located in this package
dataDir <- system.file("extdata", package="RAIDS")
#############################################################################
## Load the information about the profile
#############################################################################
data(demoPedigreeEx1)
head(demoPedigreeEx1)
## Name.ID Case.ID Sample.Type Diagnosis Source
## ex1 ex1 ex1 Primary Tumor Cancer Databank B
#############################################################################
## The population reference GDS file and SNV Annotation GDS file
## need to be located in the same directory.
## Note that the population reference GDS file used for this example is a
## simplified version and CANNOT be used for any real analysis
#############################################################################
pathReference <- file.path(dataDir, "tests")
fileGDS <- file.path(pathReference, "ex1_good_small_1KG.gds")
fileAnnotGDS <- file.path(pathReference, "ex1_good_small_1KG_Annot.gds")
#############################################################################
## A data frame containing general information about the study
## is also required. The data frame must have
## those 3 columns: "study.id", "study.desc", "study.platform"
#############################################################################
studyDF <- data.frame(study.id="MYDATA",
study.desc="Description",
study.platform="PLATFORM",
stringsAsFactors=FALSE)
#############################################################################
## The Sample SNP VCF files (one per sample) need
## to be all located in the same directory.
#############################################################################
pathGeno <- file.path(dataDir, "example", "snpPileup")
#############################################################################
## Fix RNG seed to ensure reproducible results
#############################################################################
set.seed(3043)
#############################################################################
## Select the profiles from the population reference GDS file for
## the synthetic data.
## Here we select 2 profiles from the simplified 1KG GDS for each
## subcontinental-level.
## Normally, we use 30 profile for each
## subcontinental-level but it is too big for the example.
## The 1KG files in this example only have 6 profiles for each
## subcontinental-level (for demo purpose only).
#############################################################################
gds1KG <- snpgdsOpen(fileGDS)
dataRef <- select1KGPop(gds1KG, nbProfiles=2L)
closefn.gds(gds1KG)
## GenomeInfoDb and BSgenome are required libraries to run this example
if (requireNamespace("GenomeInfoDb", quietly=TRUE) &&
requireNamespace("BSgenome.Hsapiens.UCSC.hg38", quietly=TRUE)) {
## Chromosome length information
## chr23 is chrX, chr24 is chrY and chrM is 25
chrInfo <- GenomeInfoDb::seqlengths(BSgenome.Hsapiens.UCSC.hg38::Hsapiens)[1:25]
###########################################################################
## The path where the Sample GDS files (one per sample)
## will be created needs to be specified.
###########################################################################
pathProfileGDS <- file.path(tempdir(), "exampleDNA", "out.tmp")
###########################################################################
## The path where the result files will be created needs to
## be specified
###########################################################################
pathOut <- file.path(tempdir(), "exampleDNA", "res.out")
## Example can only be run if the current directory is in writing mode
if (!dir.exists(file.path(tempdir(), "exampleDNA"))) {
dir.create(file.path(tempdir(), "exampleDNA"))
dir.create(pathProfileGDS)
dir.create(pathOut)
#########################################################################
## The wrapper function generates the synthetic dataset and uses it
## to selected the optimal parameters before calling the genetic
## ancestry on the current profiles.
## All important information, for each step, are saved in
## multiple output files.
## The 'genoSource' parameter has 2 options depending on how the
## SNP files have been generated:
## SNP VCF files have been generated:
## "VCF" or "generic" (other software)
##
#########################################################################
runExomeAncestry(pedStudy=demoPedigreeEx1, studyDF=studyDF,
pathProfileGDS=pathProfileGDS,
pathGeno=pathGeno,
pathOut=pathOut,
fileReferenceGDS=fileGDS,
fileReferenceAnnotGDS=fileAnnotGDS,
chrInfo=chrInfo,
syntheticRefDF=dataRef,
genoSource="VCF")
list.files(pathOut)
list.files(file.path(pathOut, demoPedigreeEx1$Name.ID[1]))
#######################################################################
## The file containing the ancestry inference (SuperPop column) and
## optimal number of PCA component (D column)
## optimal number of neighbours (K column)
#######################################################################
resAncestry <- read.csv(file.path(pathOut,
paste0(demoPedigreeEx1$Name.ID[1], ".Ancestry.csv")))
print(resAncestry)
## Remove temporary files created for this demo
unlink(pathProfileGDS, recursive=TRUE, force=TRUE)
unlink(pathOut, recursive=TRUE, force=TRUE)
unlink(file.path(tempdir(), "exampleDNA"), recursive=TRUE, force=TRUE)
}
}
## sample.id D K SuperPop
## 1 ex1 14 6 AFR
The runExomeAncestry() function generates 3 types of files in the pathOut directory.
In addition, a sub-directory (named using the profile ID) is also created.
The inferred ancestry is stored in the ancestry inference CSV file (“.Ancestry.csv” file) which also contains those columns:
The process is the same as for the DNA but use the wrapper function called runRNAAncestry(). Internally the data is process differently. It requires 4 files as input:
A data.frame containing the general information about the study is also required. The data.frame must contain those 3 columns:
For demonstration purpose, a small population reference GDS file (called ex1_good_small_1KG.gds) and a small population reference SNV Annotation GDS file (called ex1_good_small_1KG_Annot.gds) are included in this package. Beware that those two files should not be used to run a real ancestry inference.The results obtained with those files won’t be reliable.
The required population reference GDS file and population reference SNV Annotation GDS file should be stored in the same directory. In the example below, this directory is referred to as pathReference.
The Profile SNP file can be either in a VCF format or in a generic format.
The Profile SNP VCF file follows the VCF standard with at least those genotype fields: GT, AD and DP. The identifier of the genotype in the VCF file must correspond to the profile identifier Name.ID. The SNVs must be germline variants and should include the genotype of the wild-type homozygous at the selected positions in the reference. One file per profile is need and the VCF file must be gzipped.
Note that the name assigned to the Profile SNP VCF file has to correspond to the profile identifier Name.ID in the following analysis. For example, a SNP file called “Sample.01.vcf.gz” would be associated to the “Sample.01” profile.
A generic SNP file can replace the VCF file. The Profile SNP Generic file format is coma separated and the mandatory columns are:
Beware that the starting position in the population reference GDS File is zero (like BED files). The Profile SNP Generic file should also start at position zero.
Note that the name assigned to the Profile SNP Generic file has to correspond to the profile identifier Name.ID in the following analysis. For example, a SNP file called “Sample.01.generic.txt.gz” would be associated to the “Sample.01” profile.
The Profile PED RDS file must contain a data.frame describing all the profiles to be analyzed. These 5 mandatory columns:
Important: The row names of the data.frame must be the profiles Name.ID.
This file is referred to as the Profile PED RDS file (PED for pedigree). Alternatively, the PED information can be saved in another type of file (CVS, etc..) as long as the data.frame information can be regenerated in R (with read.csv() or else).
This example run an ancestry inference on an RNA sample. Both population reference files are demonstration files and should not be used for a real ancestry inference. Beware that running an ancestry inference on real data will take longer to run.
#############################################################################
## Load required packages
#############################################################################
library(RAIDS)
library(gdsfmt)
## Path to the demo 1KG GDS file is located in this package
dataDir <- system.file("extdata", package="RAIDS")
#############################################################################
## Load the information about the profile
#############################################################################
data(demoPedigreeEx1)
head(demoPedigreeEx1)
## Name.ID Case.ID Sample.Type Diagnosis Source
## ex1 ex1 ex1 Primary Tumor Cancer Databank B
#############################################################################
## The population reference GDS file and SNV Annotation GDS file
## need to be located in the same directory.
## Note that the population reference GDS file used for this example is a
## simplified version and CANNOT be used for any real analysis
#############################################################################
pathReference <- file.path(dataDir, "tests")
fileGDS <- file.path(pathReference, "ex1_good_small_1KG.gds")
fileAnnotGDS <- file.path(pathReference, "ex1_good_small_1KG_Annot.gds")
#############################################################################
## A data frame containing general information about the study
## is also required. The data frame must have
## those 3 columns: "study.id", "study.desc", "study.platform"
#############################################################################
studyDF <- data.frame(study.id="MYDATA",
study.desc="Description",
study.platform="PLATFORM",
stringsAsFactors=FALSE)
#############################################################################
## The Sample SNP VCF files (one per sample) need
## to be all located in the same directory.
#############################################################################
pathGeno <- file.path(dataDir, "example", "snpPileupRNA")
#############################################################################
## Fix RNG seed to ensure reproducible results
#############################################################################
set.seed(3043)
#############################################################################
## Select the profiles from the population reference GDS file for
## the synthetic data.
## Here we select 2 profiles from the simplified 1KG GDS for each
## subcontinental-level.
## Normally, we use 30 profile for each
## subcontinental-level but it is too big for the example.
## The 1KG files in this example only have 6 profiles for each
## subcontinental-level (for demo purpose only).
#############################################################################
gds1KG <- snpgdsOpen(fileGDS)
dataRef <- select1KGPop(gds1KG, nbProfiles=2L)
closefn.gds(gds1KG)
## GenomeInfoDb and BSgenome are required libraries to run this example
if (requireNamespace("GenomeInfoDb", quietly=TRUE) &&
requireNamespace("BSgenome.Hsapiens.UCSC.hg38", quietly=TRUE)) {
## Chromosome length information
## chr23 is chrX, chr24 is chrY and chrM is 25
chrInfo <- GenomeInfoDb::seqlengths(BSgenome.Hsapiens.UCSC.hg38::Hsapiens)[1:25]
#############################################################################
## The path where the Sample GDS files (one per sample)
## will be created needs to be specified.
#############################################################################
pathProfileGDS <- file.path(tempdir(), "exampleRNA", "outRNA.tmp")
#############################################################################
## The path where the result files will be created needs to
## be specified
#############################################################################
pathOut <- file.path(tempdir(), "exampleRNA", "resRNA.out")
## Example can only be run if the current directory is in writing mode
if (!dir.exists(file.path(tempdir(), "exampleRNA"))) {
dir.create(file.path(tempdir(), "exampleRNA"))
dir.create(pathProfileGDS)
dir.create(pathOut)
#########################################################################
## The wrapper function generates the synthetic dataset and uses it
## to selected the optimal parameters before calling the genetic
## ancestry on the current profiles.
## All important information, for each step, are saved in
## multiple output files.
## The 'genoSource' parameter has 2 options depending on how the
## SNP files have been generated:
## SNP VCF files have been generated:
## "VCF" or "generic" (other software)
#########################################################################
runRNAAncestry(pedStudy=demoPedigreeEx1, studyDF=studyDF,
pathProfileGDS=pathProfileGDS,
pathGeno=pathGeno,
pathOut=pathOut,
fileReferenceGDS=fileGDS,
fileReferenceAnnotGDS=fileAnnotGDS,
chrInfo=chrInfo,
syntheticRefDF=dataRef,
blockTypeID="GeneS.Ensembl.Hsapiens.v86",
genoSource="VCF")
list.files(pathOut)
list.files(file.path(pathOut, demoPedigreeEx1$Name.ID[1]))
#########################################################################
## The file containing the ancestry inference (SuperPop column) and
## optimal number of PCA component (D column)
## optimal number of neighbours (K column)
#########################################################################
resAncestry <- read.csv(file.path(pathOut,
paste0(demoPedigreeEx1$Name.ID[1], ".Ancestry.csv")))
print(resAncestry)
## Remove temporary files created for this demo
unlink(pathProfileGDS, recursive=TRUE, force=TRUE)
unlink(pathOut, recursive=TRUE, force=TRUE)
unlink(file.path(tempdir(), "example"), recursive=TRUE, force=TRUE)
}
}
## sample.id D K SuperPop
## 1 ex1 14 6 AFR
The runRNAAncestry() function generates 3 types of files in the pathOut directory.
In addition, a sub-directory (named using the profile ID) is also created.
The inferred ancestry is stored in the ancestry inference CSV file (“.Ancestry.csv” file) which also contains those columns:
A population reference dataset with known ancestry is required to infer ancestry. The population must be large enough to ensure ???
Three important reference files, containing formatted information about the reference dataset, are required:
The format of those files are described the Population reference dataset GDS files vignette.
The reference files associated to the Cancer Research associated paper are available. Note that these pre-processed files are for 1000 Genomes (1KG), in hg38. The files are available here:
https://labshare.cshl.edu/shares/krasnitzlab/aicsPaper
The size of the 1KG GDS file is 15GB.
The 1KG GDS file is mapped on hg38 (Lowy-Gallego et al. 2019).
This section can be skipped if you choose to use the pre-processed files.
Here is the output of sessionInfo()
in the environment
in which this document was compiled:
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 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.20.so; LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [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: UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Matrix_1.7-0 RAIDS_1.3.1 GENESIS_2.34.0 SNPRelate_1.38.0
## [5] gdsfmt_1.40.1 knitr_1.48 BiocStyle_2.32.1
##
## loaded via a namespace (and not attached):
## [1] jsonlite_1.8.8 shape_1.4.6.1
## [3] magrittr_2.0.3 jomo_2.7-6
## [5] GenomicFeatures_1.56.0 logistf_1.26.0
## [7] nloptr_2.1.1 rmarkdown_2.28
## [9] fs_1.6.4 BiocIO_1.14.0
## [11] zlibbioc_1.50.0 ragg_1.3.2
## [13] vctrs_0.6.5 memoise_2.0.1
## [15] minqa_1.2.8 Rsamtools_2.20.0
## [17] RCurl_1.98-1.16 quantsmooth_1.70.0
## [19] htmltools_0.5.8.1 S4Arrays_1.4.1
## [21] curl_5.2.2 broom_1.0.6
## [23] pROC_1.18.5 SparseArray_1.4.8
## [25] mitml_0.4-5 sass_0.4.9
## [27] bslib_0.8.0 htmlwidgets_1.6.4
## [29] desc_1.4.3 plyr_1.8.9
## [31] sandwich_3.1-0 zoo_1.8-12
## [33] cachem_1.1.0 GenomicAlignments_1.40.0
## [35] lifecycle_1.0.4 iterators_1.0.14
## [37] pkgconfig_2.0.3 R6_2.5.1
## [39] fastmap_1.2.0 GenomeInfoDbData_1.2.12
## [41] MatrixGenerics_1.16.0 digest_0.6.37
## [43] colorspace_2.1-1 GWASTools_1.50.0
## [45] AnnotationDbi_1.66.0 S4Vectors_0.42.1
## [47] textshaping_0.4.0 GenomicRanges_1.56.1
## [49] RSQLite_2.3.7 fansi_1.0.6
## [51] httr_1.4.7 abind_1.4-5
## [53] mgcv_1.9-1 compiler_4.4.1
## [55] bit64_4.0.5 backports_1.5.0
## [57] BiocParallel_1.38.0 DBI_1.2.3
## [59] highr_0.11 pan_1.9
## [61] MASS_7.3-61 quantreg_5.98
## [63] DelayedArray_0.30.1 rjson_0.2.22
## [65] DNAcopy_1.78.0 tools_4.4.1
## [67] lmtest_0.9-40 nnet_7.3-19
## [69] glue_1.7.0 restfulr_0.0.15
## [71] nlme_3.1-165 grid_4.4.1
## [73] generics_0.1.3 gtable_0.3.5
## [75] operator.tools_1.6.3 BSgenome_1.72.0
## [77] formula.tools_1.7.1 class_7.3-22
## [79] tidyr_1.3.1 ensembldb_2.28.1
## [81] data.table_1.16.0 utf8_1.2.4
## [83] XVector_0.44.0 BiocGenerics_0.50.0
## [85] GWASExactHW_1.2 stringr_1.5.1
## [87] foreach_1.5.2 pillar_1.9.0
## [89] splines_4.4.1 dplyr_1.1.4
## [91] lattice_0.22-6 survival_3.7-0
## [93] rtracklayer_1.64.0 bit_4.0.5
## [95] SparseM_1.84-2 BSgenome.Hsapiens.UCSC.hg38_1.4.5
## [97] tidyselect_1.2.1 SeqVarTools_1.42.0
## [99] Biostrings_2.72.1 bookdown_0.40
## [101] ProtGenerics_1.36.0 IRanges_2.38.1
## [103] SummarizedExperiment_1.34.0 stats4_4.4.1
## [105] xfun_0.47 Biobase_2.64.0
## [107] matrixStats_1.4.1 stringi_1.8.4
## [109] UCSC.utils_1.0.0 lazyeval_0.2.2
## [111] yaml_2.3.10 boot_1.3-30
## [113] evaluate_0.24.0 codetools_0.2-20
## [115] tibble_3.2.1 BiocManager_1.30.25
## [117] cli_3.6.3 rpart_4.1.23
## [119] systemfonts_1.1.0 munsell_0.5.1
## [121] jquerylib_0.1.4 Rcpp_1.0.13
## [123] GenomeInfoDb_1.40.1 png_0.1-8
## [125] XML_3.99-0.17 parallel_4.4.1
## [127] MatrixModels_0.5-3 ggplot2_3.5.1
## [129] pkgdown_2.1.0 blob_1.2.4
## [131] AnnotationFilter_1.28.0 bitops_1.0-8
## [133] lme4_1.1-35.5 glmnet_4.1-8
## [135] SeqArray_1.44.1 VariantAnnotation_1.50.0
## [137] scales_1.3.0 purrr_1.0.2
## [139] crayon_1.5.3 rlang_1.1.4
## [141] KEGGREST_1.44.1 mice_3.16.0