Main Steps
This is an overview of the genetic ancestry inference from cancer-derived molecular data:
The main steps are:
Step 1. Set-up and provide population reference files
Step 2 Sample the reference data for donor genotypes to be used for synthesis and optimize ancestry inference parameters
Step 3 Infer ancestry for the subjects of the external study
Step 4 Present and interpret the results
These steps are described in detail in the following.
Step 1. Set-up and provide population reference files
1.1 Create a directory structure
First, a specific directory structure should be created. The structure must correspond to this:
#############################################################################
## Working directory structure
#############################################################################
workingDirectory/
data/
refGDS
profileGDS
This following code creates a temporary working directory structure where the example will be run.
#############################################################################
## Create a temporary working directory structure
#############################################################################
pathWorkingDirectory <- file.path(tempdir(), "workingDirectory")
pathWorkingDirectoryData <- file.path(pathWorkingDirectory, "data")
if (!dir.exists(pathWorkingDirectory)) {
dir.create(pathWorkingDirectory)
dir.create(pathWorkingDirectoryData)
dir.create(file.path(pathWorkingDirectoryData, "refGDS"))
dir.create(file.path(pathWorkingDirectoryData, "profileGDS"))
}
1.2 Download the population reference files
The population reference files should be downloaded in the data/refGDS sub-directory. This following code downloads the complete pre-processed files for 1000 Genomes (1KG), in hg38. The size of the 1KG GDS file is 15GB.
#############################################################################
## How to download the pre-processed files for 1000 Genomes (1KG) (15 GB)
#############################################################################
cd workingDirectory
cd data/refGDS
wget https://labshare.cshl.edu/shares/krasnitzlab/aicsPaper/matGeno1000g.gds
wget https://labshare.cshl.edu/shares/krasnitzlab/aicsPaper/matAnnot1000g.gds
cd -
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.
In this running example, the demonstration files are copied in the required data/refGDS directory.
#############################################################################
## Load RAIDS package
#############################################################################
library(RAIDS)
#############################################################################
## The population reference GDS file and SNV Annotation GDS file
## need to be located in the same sub-directory.
## Note that the population reference GDS file used for this example is a
## simplified version and CANNOT be used for any real analysis
#############################################################################
## Path to the demo 1KG GDS file is located in this package
dataDir <- system.file("extdata", package="RAIDS")
pathReference <- file.path(dataDir, "tests")
fileGDS <- file.path(pathReference, "ex1_good_small_1KG.gds")
fileAnnotGDS <- file.path(pathReference, "ex1_good_small_1KG_Annot.gds")
file.copy(fileGDS, file.path(pathWorkingDirectoryData, "refGDS"))
## [1] TRUE
file.copy(fileAnnotGDS, file.path(pathWorkingDirectoryData, "refGDS"))
## [1] TRUE
Step 2 Ancestry inference with RAIDS
2.1 Set-up required directories
All required directories need to be created. In addition, the path to the reference files are kept in variables that will be used later.
#############################################################################
## The file path to the population reference GDS file
## is required (refGenotype will be used as input later)
## The file path to the population reference SNV Annotation GDS file
## is also required (refAnnotation will be used as input later)
#############################################################################
pathReference <- file.path(pathWorkingDirectoryData, "refGDS")
refGenotype <- file.path(pathReference, "ex1_good_small_1KG.gds")
refAnnotation <- file.path(pathReference, "ex1_good_small_1KG_Annot.gds")
#############################################################################
## The output profileGDS directory, inside workingDirectory/data, must be
## created (pathProfileGDS will be used as input later)
#############################################################################
pathProfileGDS <- file.path(pathWorkingDirectoryData, "profileGDS")
if (!dir.exists(pathProfileGDS)) {
dir.create(pathProfileGDS)
}
2.2 Sample reference donor profiles from the reference data
With the 1KG reference, we recommend sampling 30 donor profiles per population. For reproducibility, be sure to use the same random-number generator seed.
In the following code, only 2 profiles per population are sampled from the demo population GDS file:
#############################################################################
## Fix 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 would use 30 profiles for each subcontinental-level.
## The 1KG files in this example only have 6 profiles for each
## subcontinental-level (for demo purpose only).
#############################################################################
dataRef <- select1KGPopForSynthetic(fileReferenceGDS=refGenotype,
nbProfiles=2L)
The output object is going to be used later.
2.3 Perform the ancestry inference
Ancestry inference can be done in one function call. Within a single function call, data synthesis is performed, the synthetic data are used to optimize the inference parameters and, with these, the ancestry of the input profile donor is inferred.
According to the type of input data (RNA or DNA), a specific function should be called. The inferAncestry() function is used for DNA profiles while the inferAncestryGeneAware() function is RNA specific.
The inferAncestry() function requires a specific profile input format. The format is set by the genoSource parameter.
One of those formats is in a VCF format (genoSource=c(âVCFâ)). This format follows the VCF standard with at least those genotype fields: GT, AD and DP. The SNVs must be germline variants and should include the genotype of the wild-type homozygous at the selected positions in the reference. The VCF file must be gzipped.
A generic SNP file can replace the VCF file (genoSource=c(âgenericâ)). The format is coma separated and the mandatory columns are:
- Chromosome: The name of the chromosome
- Position: The position on the chromosome
- Ref: The reference nucleotide
- Alt: The aternative nucleotide
- Count: The total count
- File1R: The count for the reference nucleotide
- File1A: The count for the alternative nucleotide
Beware that the starting position in the population reference GDS file is zero (like BED files). The generic SNP file should also start at position zero.
In this example, the profile is from DNA source and requires the use of the inferAncestry() function.
###########################################################################
## 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 is required
## chr23 is chrX, chr24 is chrY and chrM is 25
#######################################################################
genome <- BSgenome.Hsapiens.UCSC.hg38::Hsapiens
chrInfo <- GenomeInfoDb::seqlengths(genome)[1:25]
#######################################################################
## The demo SNP VCF file of the DNA profile donor
#######################################################################
fileDonorVCF <- file.path(dataDir, "example", "snpPileup", "ex1.vcf.gz")
#######################################################################
## The ancestry inference call
#######################################################################
resOut <- inferAncestry(profileFile=fileDonorVCF,
pathProfileGDS=pathProfileGDS,
fileReferenceGDS=refGenotype,
fileReferenceAnnotGDS=refAnnotation,
chrInfo=chrInfo,
syntheticRefDF=dataRef,
genoSource=c("VCF"))
}
A profile GDS file is created in the pathProfileGDS directory while all the ancestry and optimal parameters information are integrated in the output object.
At last, all temporary files created in this example should be deleted.
#######################################################################
## Remove temporary files created for this demo
#######################################################################
unlink(pathWorkingDirectory, recursive=TRUE, force=TRUE)
Step 3. Examine the value of the inference call
The inferred ancestry and the optimal parameters are present in the list object generated by the inferAncestry() and inferAncestryGeneAware() functions.
###########################################################################
## The output is a list object with multiple entries
###########################################################################
class(resOut)
## [1] "list"
names(resOut)
## [1] "pcaSample" "paraSample" "KNNSample" "KNNSynthetic" "Ancestry"
3.1 Inspect the inference and the optimal parameters
For the global ancestry inference using PCA followed by nearest neighbor classification these parameters are D (the number of the top principal directions retained) and k (the number of nearest neighbors).
The information is stored in the Ancestry entry as a data.frame object. It is a contains those columns:
- sample.id: The unique identifier of the sample
- D: The optimal PCA dimension value used to infer the ancestry
- k: The optimal number of neighbors value used to infer the ancestry
- SuperPop: The inferred ancestry
###########################################################################
## The ancestry information is stored in the 'Ancestry' entry
###########################################################################
print(resOut$Ancestry)
## sample.id D K SuperPop
## 171 ex1 14 4 AFR
3.2 Visualize the RAIDS performance for the synthetic data
The createAUROCGraph() function enable the visualization of RAIDS performance for the synthetic data, as a function of D and k.
###########################################################################
## Create a graph showing the perfomance for the synthetic data
## The output is a ggplot object
###########################################################################
createAUROCGraph(dfAUROC=resOut$paraSample$dfAUROC, title="Example ex1")
In this specific example, the performances are lower than expected with a real profile and a complete reference population file.