Statistical inference of the extent at which Darwinian natural
selection had impacted on genetic data from multiple populations
commands a healthy portion of the phylogenetic literature (Jacques et al. 2023). Validation of these
codon-based models relies heavily on simulated data. A search of the
entries on the Bioconductor (Gentleman et al.
2004) platform, on 29 July 2024, with keywords
codon, mutation, selection,
simulate and simulation returned a total of 72
unique (out of the 2300 available Software) packages. None of the
retrieved entries was dedicated to codon data simulation for natural
selection analyses. Given the ever increasing diversity natural
selection inference models that exist, there is indeed need for
applicable packages on the platform.
Population genomic studies provided the mathematical foundation upon
which phylogenetics thrived (Wright 1931; Fisher
1922; Hardy 1908; Weinberg 1908; Darwin 1859). The thirst to
bridge the gap between these two genres of evolutionary biology continue
to drive the invention of more complex models of evolution (Aris-Brosou and Rodrigue 2012). Consequently,
there is need to develop codon sequence simulators to match the growth.
scoup is designed on the basis of the mutation-selection
(MutSel) framework (Halpern and Bruno
1998) as a contribution to this quest. Only a couple of existing
selection-focused simulators in the literature used the MutSel framework
(Spielman and Wilke 2015; Arenas and Posada
2014). This is most probably due to the perceived complexity of
the methodology. In scoup, the versatility of the Uhlenbeck and Ornstein (1930) algorithm as a
framework for evolutionary analyses (Bartoszek et
al. 2017) was used to circumvent the complexity.
In a bid to identify an appropriate quantifier that permits direct
comparison between the degree of selection signatures imposed during
simulation and that inferred, the ratio of the variance of the
non-synonymous to synonymous selection coefficients (vN/vS) was
discovered to be appropriate. The vN/vS statistic is consequently
posited as a quality selection discriminant metric. scoup
therefore represents an important contribution to the phylogenetic
modelling literature. Example code of how to successfully use the
package is presented below.
Use the following code to install scoup from the
Bioconductor platform.
Three primary evolutionary algorithms are available in
scoup. These include the frequency-dependent (Jones et al. 2017; Ayala and Campbell 1974),
the Ornstein-Uhlenbeck (OU) and the episodic algorithms. Example of
R (R Core Team 2024) code
where these functions are utilised are presented. The homogeneous (Muse and Gaut 1994; Goldman and Yang 1994) and
heterogeneous (site-wise and branch-wise) (Nielsen and Yang 1998; Kosakovsky Pond and Frost
2005) selection inference modelling contexts are explored. Data
quality is assessed by comparing the maximum likelihood inferred (ω) and the analytically calculated
(dN/dS) estimates to
the magnitude of the imposed selection pressure (measured by vN/vS).
Template code used to analyse the output data (to obtain ω) with PAML (Álvarez-Carretero, Kapli, and Yang 2023; Yang
2007) and FUBAR (Murrell et
al. 2013; Kosakovsky Pond, Frost, and Muse 2005) are presented in
the Appendix. The R code presented subsequently, require
that the user should have already installed the scoup
package.
## Loading required package: Matrix# Number of extant taxa
leaves <- 8
# Number of codon sites
sSize <- 15
# Number of data replications for each parameter combination
sims <- 1
# OU reversion parameter (Theta) value
eThta <- c(0.01)
# OU asymptotic variance value
eVary <- c(0.0001)
# OU landscape shift parameters
hbrunoStat <- hbInput(c(vNvS=1, nsynVar=0.01))
# Sequence alignment size information
seqStat <- seqDetails(c(nsite=sSize, ntaxa=leaves))
# Iterate over all listed OU variance values
for(g in seq(1,length(eVary))){
# Iterate over all listed OU reversion parameter values
for(h in seq(1,length(eThta))){
# Create appropriate simulation function ("ou") object
adaptStat <- ouInput(c(eVar=eVary[g],Theta=eThta[h]))
# Iterate over the specified number of replicates
for(i in seq(1,sims)){
# Execute simulation
simData <- alignsim(adaptStat, seqStat, hbrunoStat, NULL)
}
}
}
# Print simulated alignment
seqCOL(simData)## DNAStringSet object of length 8:
## width seq names
## [1] 45 CGACACTGCCAGTTTTTATCATTGGGTGTATTGTGCCGCACTGCA S001
## [2] 45 CGACGCTGCCAGTTTTTATCATTGGGTGTATTGTGCCGCACTGCA S002
## [3] 45 CGACACTGCCAGTTTTTATCATTGGGTGTATTGTACCGCACTGCA S003
## [4] 45 CGACACTGCCAGTTTTTATCATTGGGTGTATTGTACCGCACTGCA S004
## [5] 45 CAACACTGCCAGTTTTTATCATTGGGTGTATTATGCCGCACTGCA S005
## [6] 45 CAACACTGCCAGTTTTTATCATTGGGTGTATTATGCCGCACTGCA S006
## [7] 45 CAACACTGCCAGTTTTTATCATTGGGTATATTATGCCGTACTGCA S007
## [8] 45 CAACACTGCCAGTTTTTATCATTGGGTATATTATGCCGTACTGCA S008As expected, the correlation coefficient estimate was approximately 0.9974 when the means of the inferred (ω) and the calculated (dN/dS) selection effects were first compared. The correlation estimation included more iterations.
# Make package accessible in R session
library(scoup)
# Number of extant taxa
xtant <- 8
# Number of codon sites
siteSize <- 15
# Number of data replications for each parameter combination
simSize <- 1
# Variance of the non-synonymous selection coefficients
nsynVary <- c(0.001)
# Ratio of the variance of the non-synonymous to synonymous coeff.
vNvSvec <- c(1)
# Sequence alignment size information
seqStat <- seqDetails(c(nsite=siteSize, ntaxa=xtant))
# Iterate over all listed coefficient variance ratios
for(a in seq(1,length(vNvSvec))){
# Iterate over all listed non-synonymous coefficients variance
for(b in seq(1,length(nsynVary))){
# Create appropriate simulation function ("omega") object
adaptData <- wInput(list(vNvS=vNvSvec[a],nsynVar=nsynVary[b]))
# Iterate over the specified number of replicates
for(i in seq(1,simSize)){
# Execute simulation
simulateSeq <- alignsim(adaptData, seqStat, NA)
}
}
}
# Print simulated alignment
cseq(simulateSeq)## Sequence
## S001 TCGGACATGGATCCACTAAATCAGAGGACTTCTAAGTTAATTGTC
## S002 TCGGACATGGATCCACTAAATCAGAGGACTTCTAAGTTAATTGTC
## S003 TCGGACTTGGATCCACTAAATCAGAGGACTTCTAAGTTAATTGTC
## S004 TCGGACATGGATCCACTAAATCAGAGGACTTCTAAGTTAATTGTC
## S005 TCGGACATGGATCCACTAAATCAGAGGACTTCTAAGTTAATTGTC
## S006 TCGGACATGGATCCACTAAATCAGAGGACTTCTAAGTTAATTGTC
## S007 TCGGACATGGATCCACCAAATCAGAGGACTTCTAAGTTAGTTGCC
## S008 TCGGACATGGATCCACTAAATCAGAGGACTTCTAAGTTAGTTGTCSequences generated with the presented code (with more iterations included) produced strongly correlated selection inferences (correlation coefficient ≈ 0.9923) when the average dN/dS and the ω values were compared. This implementation is an example of how to execute the frequency-dependent evolutionary technique with the package.
# Make package accessible in R session
library(scoup)
# Number of codon sites
sitesize<- 15
# Variance of non-synonymous selection coefficients
nsynVary <- 0.01
# Number of extant taxa
## Commented value was used for results presented in article
taxasize <- 8
# Sequence alignment size information
seqsEntry <- seqDetails(c(nsite=sitesize, ntaxa=taxasize))
# Create the applicable ("ou") object for simulation function
## eVar= OU asymptotic variance, Theta=OU reversion parameter
adaptEntry <- ouInput(c(eVar=0.1,Theta=1))
# Ratio of the variance of the non-synonymous to synonymous coeff.
sratio <- c(1e-06)
# Iterate over all listed coefficient variance ratios
for(a0 in seq(1,length(sratio))){
# OU landscape shift parameters
mValues <- hbInput(c(vNvS=sratio[a0], nsynVar=nsynVary))
# Execute simulation
simSeq <- alignsim(adaptEntry, seqsEntry, mValues, NA)
}
# Print simulated codon sequence
cseq(simSeq)## Sequence
## S001 TACAGAATAATCTCACCCGTTCGGTTTATAGCCCGGCGTCTTCTC
## S002 TACGGAATAATCCTACCTGCCCGATTTAGAGACCGGCGACTTCTC
## S003 TACAGGATAATCTCACCCGCTCCGTCTGGAGATCGTCGACGCCTC
## S004 TACAGGCTAATCTCTCCCGCTGCATCTGGAAGCGGTCGACTCCTT
## S005 CATGAGATAATTCGATCATTTACTTTCAGAGGCCAGCCGCCCGCT
## S006 CACAAGATAATTAAATTTGTTACTGCCAGGCACCCGCCGCTCATT
## S007 TACAGGCTAATTCCATTTGTTGGCTTCAGGGACCAGCCACTCATT
## S008 TGCAGGACAATTCCATCAATTGGTTTCGGGGATCAGCCGCACACTThis is another example of how to call the OU shifting landscape evolutionary approach. The results obtained yielded a pairwise correlation coefficient estimate of approximately 0.9988 between the means of dN/dS and ω. The correlation coefficient estimates were approximately 0.8123 and 0.8305 when the averages were each compared to vN/vS, respectively.
# Make package accessible in R session
library(scoup)
# Number of internal nodes on the desired balanced tree
iNode <- 3
# Number of required codon sites
siteCount <- 15
# Variance of non-synonymous selection coefficients
nsnV <- 0.01
# Number of data replications for each parameter combination
nsim <- 1
# Ratio of the variance of the non-synonymous to synonymous coeff.
vNvSvec <- c(1e-06)
# Sequence alignment size information
seqsBwise <- seqDetails(c(nsite=siteCount, blength=0.10))
# Iterate over all listed coefficient variance ratios
for(h in seq(1,length(vNvSvec))){
# Iterate over the specified number of replicates
for(i in seq(1,nsim)){
# Create the parameter set applicable at each internal tree node
scInput <- rbind(vNvS=c(rep(0,iNode-1),vNvSvec[h]),
nsynVar=rep(nsnV,iNode))
# Create the applicable ("discrete") object for simulation function
adaptBranch <- discreteInput(list(p02xnodes=scInput))
# Execute simulation
genSeq <- alignsim(adaptBranch, seqsBwise, NULL)
}
}
# Print simulated sequence data
seqCOL(genSeq)## DNAStringSet object of length 8:
## width seq names
## [1] 45 AGCAAATCTAGAAACCGGCCAAGGGTAGCAGCAGTAGACAAATTA S001
## [2] 45 AGTAAATCGAGAAACCGGCCAAGGGTCACGGCAGTGGTCAGACTA S002
## [3] 45 AGCAAATCCAGAGATCTGCCTAGGATAACAGTAGTGGACAAATTG S003
## [4] 45 AGAAAATCTAGCAATCGGTCTAGGATAACAGTAGTGGACAGATTG S004
## [5] 45 AACAAATCTAGAAATCGGCCCAGAATAACGACAGTGGATAAATTG S005
## [6] 45 AACAAATCCATAAATGGGCCCCGGACAGTAGCAGTGCATAAGTTG S006
## [7] 45 AGCGAATCTAAGAATCAGCCCAGGATAACTGCAGCAGATATATTG S007
## [8] 45 AGCAAGTCTAAGAAGCGACCCAGGATAACTGCAGTAGATATATTG S008The correlation coefficient between the averages of the analytical dN/dS and the inferred ω estimates was approximately 0.9998, obtained from 50 independent iterations of the code for a set of vN/vS values. The correlation coefficient estimate was approximately 0.6349 for vN/vS vs ω and 0.6360 for vN/vS vs dN/dS.
Reference scoup code are presented to facilitate use of
the package. Although not explicitly presented, it is also possible to
generate data with signatures of directional selection by setting the
aaPlus element of the wInput entry of the
alignsim function accordingly. The capacity of the package
is expected to be extended in future versions.
A more appropriate citation will be provided for the package after it has been accepted for publication. In the meantime, to cite this package, use Sadiq, H. 2024. “scoup: Simulate Codon Sequences with Darwinian Selection Incorporated as an Ornstein-Uhlenbeck Process”. R Package. doi:10.18129/B9.bioc.scoup.
CODEML script executed in PAML to infer
single alignment-wide ω
estimates for data sets generated from 50 independent executions of the
sensitivity analyses code presented above. The same CODEML
script was used to analyse data (also 50 replicates) from the episodic
analyses code, with the model entry replaced with
2. The scoup simulated sequence data and tree
are seq.nex and seq.tre, respectively.
seqfile = seq.nex
treefile = seq.tre
outfile = seq.out
noisy = 0
verbose = 0
seqtype = 1
ndata = 1
icode = 0
cleandata = 0
model = 0
NSsites = 0
CodonFreq = 3
estFreq = 0
clock = 0
fix_omega = 0
omega = 1e-05FUBAR command executed with HyPhy through
the terminal in MacBook. Note that HyPhy was already
installed on the computer. The seq.nex input is the
scoup simulated codon sequence data that is saved in the
same NEXUS file with the tree data. The NEXUS file resides in the
working directory.
The output of sessionInfo() from the computer where this
file is generated is provided below.
## 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] scoup_1.1.4 Matrix_1.7-1 BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] crayon_1.5.3 httr_1.4.7 cli_3.6.3
## [4] knitr_1.49 rlang_1.1.4 xfun_0.50
## [7] UCSC.utils_1.3.0 generics_0.1.3 jsonlite_1.8.9
## [10] S4Vectors_0.45.2 buildtools_1.0.0 Biostrings_2.75.3
## [13] htmltools_0.5.8.1 maketools_1.3.1 sys_3.4.3
## [16] stats4_4.4.2 sass_0.4.9 rmarkdown_2.29
## [19] grid_4.4.2 evaluate_1.0.1 jquerylib_0.1.4
## [22] fastmap_1.2.0 GenomeInfoDb_1.43.2 IRanges_2.41.2
## [25] yaml_2.3.10 lifecycle_1.0.4 BiocManager_1.30.25
## [28] compiler_4.4.2 XVector_0.47.2 lattice_0.22-6
## [31] digest_0.6.37 R6_2.5.1 GenomeInfoDbData_1.2.13
## [34] bslib_0.8.0 tools_4.4.2 BiocGenerics_0.53.3
## [37] cachem_1.1.0