R/withingenednds.R
In im3sanger/dndscv: Poisson-based dN/dS models to quantify natural selection in somatic evolution
Defines functions
withingenednds
Documented in
withingenednds
#' withingenednds
#'
#' This function uses Poisson and Negative Binomial regression models at single-site level to study selection across different regions (coding and non-coding) within a gene.
#'
#' @author Inigo Martincorena (Wellcome Sanger Institute)
#' @details Martincorena I, et al. (2017) Universal patterns of selection in cancer and somatic tissues. Cell. 171(5):1029-1041.
#'
#' @param mutations Data frame with all the mutations detected in the study (5-column input table as for dndscv: sampleID, chr, pos, ref, mut).
#' @param gene Name of the gene of interest. This function is currently designed to work on a single gene, but combined analyses of multiple genes could be done using the sites output table generated by this function.
#' @param covtable Table with all sites of interest in the gene. This should be a data frame with one row per site and the following columns: chr, pos, dc (duplex depth). Additional columns will not be used.
#' @param dndsout dndscv output object for the dataset. This is mainly used for the MLEs of the substitution model. Running dndscv on all genes in the dataset is recommended unless the gene of interest is believed to have a different substitution model.
#' @param genomeFile Path to a reference fasta file for the genome assembly.
#' @param regionschr Optional data frame with user-defined regions of interest in the gene. This allows the user to define arbitrary regions within a gene (coding or non-coding) from which to calculate omega (selection or obs/exp) values (e.g. protein domains, splicing regulator regions, core promoters, etc). The table should contain the following columns: chr, start, end, wname (a unique name for the w parameter, e.g. wdomain1, wcorepromoter), impacts (e.g. Missense or Missense|Nonsense will restrict the w calculation with Missense or Missense and Nonsense mutations in the region, respectively), layered (1/0; using "0" removes other w parameters influencing the site, whereas using "1" models selection as relative to other w parameters active at these sites).
#' @param regionsaa Optional data frame with user-defined regions of interest in the gene, using aminoacid coordinates. The table should contain the following columns: gene, aa_start, aa_end, w feature name (e.g. wdomain1), impacts.
#' @param fixtheta Pre-calculated overdispersion (theta) parameter. This should be calculated using sitednds(., method="NB").
#' @param refdb Reference database (path to .rda file or a pre-loaded array object in the right format).
#' @param numcode NCBI genetic code number (default = 1; standard genetic code). To see the list of genetic codes supported use: ? seqinr::translate
#' @param normalisefromsyn Normalise the substitution rates based on the synonymous mutations in the gene. Using TRUE is recommended. Using FALSE uses the expected synonymous mutation rate of the gene from the dndscv negative binomial regression model (dndsout$genemuts).
#' @param syndrivers Vector of known synonymous driver sites defined by their aminoacid position, to be excluded from the background model (e.g. syndrivers = c("T125T","E224E","Q331Q") for TP53).
#' @param exon_flank_length Exon flank length in bp [default = 10]. Using a value higher than 0 will calculate a separate selection (w) coefficient for synonymous mutations in exon flanks.
#' @param intron_flank_length Intron flank length in bp [default = 10]. Intronic sites occurring within these flanks but not already classified as Essential_Splice will receive a separate w parameter.
#' @param sitefilename Optionally, provide a file name to save the table of all annotated sites in the gene. This table is also always contained in the output object.
#'
#' @return 'withingenednds' returns a list of objects:
#' @return - sites: Table with the annotation of all sites in the gene (from covtable), including all functional annotations in the "regions" input object as well as default annotations (Missense, Nonsense, Essential_Splice, Start_loss, Stop_loss, etc).
#' @return - par.pois: Poisson regression results (not recommended).
#' @return - par.nb: Negative binomial results fitting a new overdispersion parameter to the data (when fixtheta is not provided).
#' @return - par.nbfix: Negative binomial results using the input fixtheta value as recommended.
#' @return - model.pois: Poisson regression object.
#' @return - model.nb: Negative binomial regression object.
#' @return - model.nbfix: Negative binomial regression object.
#'
#' @export
withingenednds = function(mutations, gene, covtable, dndsout, genomeFile, regionschr = NULL, regionsaa = NULL, fixtheta = NULL, normalisefromsyn = TRUE, syndrivers = NULL, exon_flank_length = 10, intron_flank_length = 10, sitefilename = NULL, refdb = "hg19", numcode = 1) {
## 1. Environment
message("Running gene-level selection analyses...")
# [Input] Reference database
refdb_class = class(refdb)
if ("character" %in% refdb_class) {
if (refdb == "hg19") {
data("refcds_hg19", package="dndscv")
} else {
load(refdb)
}
} else if("array" %in% refdb_class) {
# use the user-supplied RefCDS object
RefCDS = refdb
} else {
stop("Expected refdb to be \"hg19\", a file path, or a RefCDS-formatted array object.")
}
# Annotating all possible coding changes in the positions provided in the input covtable
library("Rsamtools")
covtable$ref = as.vector(scanFa(genomeFile, GRanges(covtable$chr, IRanges(covtable$pos, covtable$pos))))
covtable$ref3 = as.vector(scanFa(genomeFile, GRanges(covtable$chr, IRanges(covtable$pos-1, covtable$pos+1))))
allsubs = data.frame(sampleID="allsubs", chr = rep(covtable$chr, each=3), pos = rep(covtable$pos, each=3), ref = rep(covtable$ref, each=3), mut = NA, ref3 = rep(covtable$ref3, each=3), mut3 = NA, ref_cod = NA, mut_cod = NA, ref3_cod = NA, mut3_cod = NA)
for (j in seq(1,nrow(allsubs),3)) {
allsubs$mut[j:(j+2)] = setdiff(c("A","C","G","T"),allsubs$ref[j])
}
allsubs$mut3 = paste(substr(allsubs$ref3,1,1), allsubs$mut, substr(allsubs$ref3,3,3), sep="")
aux = dndscv(allsubs, gene_list = gene, max_coding_muts_per_sample = Inf, max_muts_per_gene_per_sample = Inf, outp = 1)$annotmuts # Annotated mutations
# Adding duplex depth, functional impact annotation to all possible changes in the input table
aux$mstr = paste(aux$chr, aux$pos, aux$mut, sep=":") # mutation string ID
obssubs = dndsout$annotmuts[which(dndsout$annotmuts$ref %in% c("A","C","G","T") & dndsout$annotmuts$mut %in% c("A","C","G","T") & dndsout$annotmuts$gene == gene), ]
obssubs$mstr = paste(obssubs$chr, obssubs$pos, obssubs$mut, sep=":")
m$mstr = paste(m$chr, m$pos, m$mut, sep=":")
sites = cbind(data.frame(gene=gene, pid=aux$pid[1]), allsubs[,-1]) # Initialising the sites table
pos2dc = setNames(covtable$dc, covtable$pos)
sites$dc = pos2dc[as.character(sites$pos)]
sites$mstr = paste(sites$chr, sites$pos, sites$mut, sep=":")
sites$obs = table(m$mstr)[sites$mstr]; sites$obs[is.na(sites$obs)] = 0 # Annotating the number of observed mutations at each site
mstr2imp1 = setNames(aux$ntchange, aux$mstr)
mstr2imp2 = setNames(aux$aachange, aux$mstr)
mstr2imp3 = setNames(aux$impact, aux$mstr)
sites$ntchange = mstr2imp1[sites$mstr]
sites$aachange = mstr2imp2[sites$mstr]
sites$impact = mstr2imp3[sites$mstr]
# Adding strand and annotating the coding trinucleotide change
aux2 = unique(dndsout$annotmuts[,c("gene","strand")])
gene2strand = setNames(aux2[,2],aux2[,1])
sites$strand = gene2strand[sites$gene]
sites$ref_cod = sites$ref
sites$ref_cod[sites$strand==-1] = seqinr::comp(sites$ref[sites$strand==-1], forceToLower = F)
sites$mut_cod = sites$mut
sites$mut_cod[sites$strand==-1] = seqinr::comp(sites$mut[sites$strand==-1], forceToLower = F)
revcomp = function(seqvec) {
as.vector(sapply(seqvec, function(x) paste(seqinr::comp(rev(unlist(strsplit(x,split=""))), forceToLower=F), collapse=""))) # Reverse complement of a vector of sequence motifs
}
sites$ref3_cod = sites$ref3
sites$mut3_cod = sites$mut3
if (any(sites$strand==-1)) {
sites$ref3_cod[sites$strand==-1] = revcomp(sites$ref3[sites$strand==-1])
sites$mut3_cod[sites$strand==-1] = revcomp(sites$mut3[sites$strand==-1])
}
# Adding the trinucleotide substitution rates from the dndsout model
mle_submodel = setNames(dndsout$mle_submodel[,2], dndsout$mle_submodel[,1])
mle_submodel = c(mle_submodel, "TTT>TGT"=1) # Adding the missing rate (note that all rates in mle_submodel in dNdScv are relative to TTT>TGT)
mle_submodel = mle_submodel * mle_submodel["t"] # Absolute rates
sites$r = mle_submodel[paste(sites$ref3_cod, sites$mut3_cod, sep=">")]
sites$r = sites$r * sites$dc / mean(sites$dc) # Normalising the expected rates at a site by the duplex depth
# Normalising the global expected rates by the estimated mutation rate of the gene using one of two alternative models:
# 1. normalisefromsyn = TRUE. We normalise "r" using the synonymous mutations observed in the gene (excluding known syn driver sites)
# 2. normalisefromsyn = FALSE. We normalise "r" using the negative binomial model from dndscv.
if (normalisefromsyn == TRUE) {
mobs = sum(sites$obs[which(sites$impact=="Synonymous" & !(sites$aachange %in% syndrivers))])
mexp = sum(sites$r[which(sites$impact=="Synonymous" & !(sites$aachange %in% syndrivers))])
sites$rnorm = sites$r * mobs / mexp
} else {
message("Option not recommended: Normalising the mutation rate of the gene based on the negative regression model of dNdScv")
sites$rnorm = sites$r * dndsout$genemuts$exp_syn_cv[dndsout$genemuts$gene_name==gene] / dndsout$genemuts$exp_syn[dndsout$genemuts$gene_name==gene]
}
# Excluding sites with rate = 0 (e.g. due to lack of duplex coverage)
ratezero = which(sites$rnorm==0)
if (length(ratezero)>0) {
sites = sites[-ratezero, ]
message(sprintf("%0.0f sites have been excluded due to having a duplex depth and/or a predicted mutation rate of 0", length(ratezero)))
}
## 2. Creating an index regression matrix with all functional annotations (each row is a site and each column is a parameter in the selection model)
# Annotating Missense, Nonsense, Essential_Splice and Stop_loss mutations
mutmat = data.frame(t = rep(1,nrow(sites)), wmis = 0, wnon = 0, wspl = 0, wsyndriv = 0,
wintr = 0, woutcds = 0, wexfl = 0, winfl = 0, wnonlastex = 0, wstoploss = 0, wstartloss = 0)
mutmat$wmis[which(sites$impact=="Missense")] = 1
mutmat$wnon[which(sites$impact=="Nonsense")] = 1
mutmat$wspl[which(sites$impact=="Essential_Splice")] = 1
mutmat$wsyndriv[which(sites$impact=="Synonymous" & sites$aachange %in% syndrivers)] = 1
mutmat$wstoploss[which(sites$impact=="Stop_loss")] = 1
# Annotating Start_loss mutations (mutations in codon 1)
#start_loss = which(sites$impact != "Synonymous" & substr(sites$aachange,2,nchar(sites$aachange)-1) == "1") # There is no need to check for synonymous changes in ATG
start_loss = which(substr(sites$aachange,2,nchar(sites$aachange)-1) == "1") # There is no need to check for synonymous changes in ATG
mutmat$wstartloss[start_loss] = 1
sites$impact[start_loss] = "Start_loss"
# Annotating introns, exon flanks, and intron flanks
refcdsgene = RefCDS[[which(sapply(RefCDS, function(x) x$gene_name)==gene)]]
exons = refcdsgene$intervals_cds
esspl = refcdsgene$intervals_splice
gr_sites = GenomicRanges::GRanges(sites$gene, IRanges::IRanges(sites$pos,sites$pos))
gr_exons = GenomicRanges::GRanges(gene, IRanges::IRanges(exons[,1],exons[,2]))
gr_cds = GenomicRanges::GRanges(gene, IRanges::IRanges(min(exons),max(exons)))
gr_outcds = setdiff(gr_sites, gr_cds)
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_outcds, type="any", select="all"))
mutmat$woutcds[unique(ol[,1])] = 1 # Annotation of the intronic sites
# For genes with >1 exon, we also annotate essential splice sites, and intron and exon flanks.
if (length(esspl)>0) {
gr_esspl = GenomicRanges::GRanges(gene, IRanges::IRanges(esspl,esspl))
gr_introns = setdiff(setdiff(gr_cds, gr_exons), gr_esspl)
exfl = rbind(cbind(exons[2:nrow(exons),1],exons[2:nrow(exons),1]+exon_flank_length-1), cbind(exons[1:(nrow(exons)-1),2]-exon_flank_length+1,exons[1:(nrow(exons)-1),2]))
gr_exonflanks = GenomicRanges::GRanges(gene, IRanges::IRanges(exfl[,1],exfl[,2]))
gr_exonflanks = intersect(gr_exonflanks, gr_exons) # Ensuring that the exon flanks do not extend into introns
infl = rbind(cbind(exons[2:nrow(exons),1]-intron_flank_length,exons[2:nrow(exons),1]-1), cbind(exons[1:(nrow(exons)-1),2]+1,exons[1:(nrow(exons)-1),2]+intron_flank_length))
gr_intrflanks = GenomicRanges::GRanges(gene, IRanges::IRanges(infl[,1],infl[,2]))
gr_intrflanks = setdiff(gr_intrflanks, gr_exons) # Removing any overlaps with exonic sequences
gr_intrflanks = setdiff(gr_intrflanks, gr_esspl) # Removing any overlaps with essential splice site sequences
gr_intrflanks = intersect(gr_intrflanks, gr_cds) # Removing UTR sequences
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_introns, type="any", select="all"))
mutmat$wintr[unique(ol[,1])] = 1 # Annotation of the intronic sites
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_introns, type="any", select="all"))
mutmat$wintr[unique(ol[,1])] = 1 # Annotation of the intronic sites
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_exonflanks, type="any", select="all"))
mutmat$wexfl[seq(1,nrow(sites)) %in% ol[,1] & sites$impact=="Synonymous"] = 1 # Annotation of the exonic flank sites (only for synonymous mutations)
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_intrflanks, type="any", select="all"))
mutmat$winfl[unique(ol[,1])] = 1 # Annotation of the exonic flank sites
}
# Annotation of nonsense mutations in the last coding exon
if (nrow(exons)>1) {
if (refcdsgene$strand==1) {
lastexon = exons[nrow(exons),]
} else {
lastexon = exons[1,]
}
gr_lastexon = GenomicRanges::GRanges(gene, IRanges::IRanges(min(lastexon),max(lastexon)))
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_lastexon, type="any", select="all"))
mutmat$wnonlastex[seq(1,nrow(sites)) %in% ol[,1] & sites$impact=="Nonsense"] = 1
mutmat$wnon[seq(1,nrow(sites)) %in% ol[,1] & sites$impact=="Nonsense"] = 0 # Removing Nonsense mutations in the last exon from the wnon annotation
}
## 3. Other annotations from the user
# Regions defined by chr position
if (!is.null(regionschr)) {
wnames = unique(regionschr$wname)
badnames = intersect(wnames,unique(colnames(mutmat)))
if (length(badnames)>0) { stop(sprintf("The following w names are not allowed in the input regions as they match existing parameters: %s", paste(badnames,collapse = ","))) }
for (j in 1:length(wnames)) {
aux = regionschr[which(regionschr$wname==wnames[j]), ]
if (length(unique(aux$impacts))!=1) { stop("regionschr: different values found in the impacts column for a given feature, please correct your input object") }
imps = strsplit(unique(aux$impacts), split=",")[[1]]
if (any(!(imps %in% setdiff(unique(sites$impact),NA)))) { stop("regionschr: invalid impact terms used in the input object, please correct the impact column") }
gr_wname = GenomicRanges::GRanges(gene, IRanges::IRanges(aux$start,aux$end))
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_wname, type="any", select="all"))
mutmat[,wnames[j]] = 0 # Initialising this field in the data frame
if (length(imps)==0) {
indx = unique(ol[,1]) # If no impacts are indicated, we consider all sites independent of their impact
} else {
indx = intersect(unique(ol[,1]), which(sites$impact %in% imps))
}
mutmat[indx,wnames[j]] = 1
# Removing previous w parameters if layered=0
if (aux$layered[1]==0 | aux$layered[1]==FALSE) {
mutmat[indx, setdiff(colnames(mutmat),c("t",wnames[j]))] = 0 # Removing previous annotations at these sites
}
}
}
# Regions defined by aminoacid position
if (!is.null(regionsaa)) {
sites$aux = 1:nrow(sites)
aas = sites[which(!is.na(sites$aachange) & sites$aachange!="."),]
aas$aapos = as.numeric(substr(aas$aachange,2,nchar(aas$aachange)-1))
gr_aas = GenomicRanges::GRanges(gene, IRanges::IRanges(aas$aapos,aas$aapos))
wnames = unique(regionsaa$wname)
badnames = intersect(wnames,unique(colnames(mutmat)))
if (length(badnames)>0) { stop(sprintf("The following w names are not allowed in the input regions as they match existing parameters: %s", paste(badnames,collapse = ", "))) }
for (j in 1:length(wnames)) {
aux = regionsaa[which(regionsaa$wname==wnames[j]), ]
if (length(unique(aux$impacts))!=1) { stop("regionschr: different values found in the impacts column for a given feature, please correct your input object") }
imps = strsplit(unique(aux$impacts), split=",")[[1]]
if (any(!(imps %in% setdiff(unique(sites$impact),NA)))) { stop("regionschr: invalid impact terms used in the input object, please correct the impact column") }
gr_wname = GenomicRanges::GRanges(gene, IRanges::IRanges(aux$start,aux$end))
ol = as.data.frame(GenomicRanges::findOverlaps(gr_aas, gr_wname, type="any", select="all"))
mutmat[,wnames[j]] = 0 # Initialising this field in the data frame
if (length(imps)==0) {
indx = aas$aux[unique(ol[,1])] # If no impacts are indicated, we consider all sites independent of their impact
} else {
indx = aas$aux[intersect(unique(ol[,1]), which(aas$impact %in% imps))]
}
mutmat[indx,wnames[j]] = 1
# Removing previous w parameters if layered=0
if (aux$layered[1]==0 | aux$layered[1]==FALSE) {
mutmat[indx, setdiff(colnames(mutmat),c("t",wnames[j]))] = 0 # Removing previous annotations at these sites
}
}
sites = sites[, setdiff(colnames(sites),"aux")]
}
## 4. Poisson regression: fitting the selection model
model = glm(formula = sites$obs ~ offset(log(sites$rnorm)) + . -1, data=mutmat, family=poisson(link=log))
mle = exp(coefficients(model)) # Maximum-likelihood estimates for the rate params
pvals = coef(summary(model))[,4]
model.lrt = drop1(model, test= "Chisq")
pvals.lrt = setNames(model.lrt[[5]], row.names(model.lrt))
ci = exp(confint.default(model)) # Wald confidence intervals
par.pois = data.frame(name=gsub("\`","",rownames(ci)), mle=mle[rownames(ci)], cilow=ci[,1], cihigh=ci[,2], pval.wald=pvals[rownames(ci)], pval.lrt=pvals.lrt[rownames(ci)]) # MLEs and Wald CI95% for the selection parameters
# Adding obs/exp statistics to the regression outputs
nobs = apply(mutmat, 2, function(x) sum(sites$obs[x==1]))
nexp = apply(mutmat, 2, function(x) sum(sites$rnorm[x==1]))
par.pois$obs = nobs[par.pois$name]
par.pois$exp = nexp[par.pois$name]
## 5. Negative binomial regression
model.nbfix = model.nb = par.nbfix = par.nb = NULL
if (!is.null(fixtheta)) {
model.nbfix = glm(formula = sites$obs ~ offset(log(sites$rnorm)) + . -1, data=mutmat, family=MASS::negative.binomial(fixtheta))
mle = exp(coefficients(model.nbfix)) # Maximum-likelihood estimates for the rate params
pvals = coef(summary(model.nbfix))[,4]
model.lrt = drop1(model.nbfix, test= "Chisq")
pvals.lrt = setNames(model.lrt[[5]], row.names(model.lrt))
ci = exp(confint.default(model.nbfix)) # Wald confidence intervals
par.nbfix = data.frame(name=gsub("\`","",rownames(ci)), mle=mle[rownames(ci)], cilow=ci[,1], cihigh=ci[,2], pval.wald=pvals[rownames(ci)], pval.lrt=pvals.lrt[rownames(ci)]) # MLEs and Wald CI95% for the selection parameters
par.nbfix$obs = nobs[par.nbfix$name]
par.nbfix$exp = nexp[par.nbfix$name]
} else {
model.nb = MASS::glm.nb(formula = as.vector(sites$obs) ~ offset(log(sites$rnorm)) + . -1, data=mutmat)
mle = exp(coefficients(model.nb)) # Maximum-likelihood estimates for the rate params
pvals = coef(summary(model.nb))[,4]
model.lrt = drop1(model.nb, test= "Chisq")
pvals.lrt = setNames(model.lrt[[5]], row.names(model.lrt))
ci = exp(confint.default(model.nb)) # Wald confidence intervals
par.nb = data.frame(name=gsub("\`","",rownames(ci)), mle=mle[rownames(ci)], cilow=ci[,1], cihigh=ci[,2], pval.wald=pvals[rownames(ci)], pval.lrt=pvals.lrt[rownames(ci)]) # MLEs and Wald CI95% for the selection parameters
par.nb$obs = nobs[par.nb$name]
par.nb$exp = nexp[par.nb$name]
}
## 6. Outputs
# Annotated sites table
sites2 = cbind(sites, mutmat)
if (!is.null(sitefilename)) {
sites2 = apply(sites2,2,as.character)
write.table(sites2, file=sitefilename, col.names = T, row.names = F, quote = F, sep = "\t")
}
# Inform the user about the sites used as neutral reference by the model
neutral_sites = (rowSums(mutmat[,setdiff(colnames(mutmat),"t")])==0)
neutral_nsyn = sum(sites$obs[which(neutral_sites & sites$impact=="Synonymous")]) # Number of synonymous mutations used as background for dN/dS
neutral_nsynsites = length(which(neutral_sites & sites$impact=="Synonymous"))
neutral_othersites = length(which(neutral_sites & sites$impact!="Synonymous"))
nonneutral_nsyn = sum(sites$obs[which(neutral_sites==0 & sites$impact=="Synonymous")]) # Number of synonymous mutations excluded from the neutral backgroup by the "w" annotations
message(sprintf(" Sites used as neutral reference: %0.0f synonymous mutations across %0.0f synonymous sites.", neutral_nsyn, neutral_nsynsites))
if (neutral_othersites>0) {
message(sprintf(" %0.0f sites not classified as synonymous coding sites were used in the background model. Please ensure that this was intended.", neutral_othersites))
}
if (nonneutral_nsyn>0) {
message(sprintf(" %0.0f synonymous mutations were not used in the neutral background model. This may be due to excluding known synonymous driver mutations (when using syndrivers), excluding synonymous mutations in exon flanks (when using exon_flank_length>0) and/or due to annotations provided by the user. Please ensure that this is the desired behaviour.", nonneutral_nsyn))
}
# Output object
out = list(sites = sites2, dnds.pois = par.pois, dnds.nb = par.nb, dnds.nbfix = par.nbfix, model.pois = model, model.nb = model.nb, model.nbfix = model.nbfix)
}
im3sanger/dndscv documentation built on June 11, 2025, 8:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions withingenednds
Documented in withingenednds
#' withingenednds
#'
#' This function uses Poisson and Negative Binomial regression models at single-site level to study selection across different regions (coding and non-coding) within a gene.
#'
#' @author Inigo Martincorena (Wellcome Sanger Institute)
#' @details Martincorena I, et al. (2017) Universal patterns of selection in cancer and somatic tissues. Cell. 171(5):1029-1041.
#'
#' @param mutations Data frame with all the mutations detected in the study (5-column input table as for dndscv: sampleID, chr, pos, ref, mut).
#' @param gene Name of the gene of interest. This function is currently designed to work on a single gene, but combined analyses of multiple genes could be done using the sites output table generated by this function.
#' @param covtable Table with all sites of interest in the gene. This should be a data frame with one row per site and the following columns: chr, pos, dc (duplex depth). Additional columns will not be used.
#' @param dndsout dndscv output object for the dataset. This is mainly used for the MLEs of the substitution model. Running dndscv on all genes in the dataset is recommended unless the gene of interest is believed to have a different substitution model.
#' @param genomeFile Path to a reference fasta file for the genome assembly.
#' @param regionschr Optional data frame with user-defined regions of interest in the gene. This allows the user to define arbitrary regions within a gene (coding or non-coding) from which to calculate omega (selection or obs/exp) values (e.g. protein domains, splicing regulator regions, core promoters, etc). The table should contain the following columns: chr, start, end, wname (a unique name for the w parameter, e.g. wdomain1, wcorepromoter), impacts (e.g. Missense or Missense|Nonsense will restrict the w calculation with Missense or Missense and Nonsense mutations in the region, respectively), layered (1/0; using "0" removes other w parameters influencing the site, whereas using "1" models selection as relative to other w parameters active at these sites).
#' @param regionsaa Optional data frame with user-defined regions of interest in the gene, using aminoacid coordinates. The table should contain the following columns: gene, aa_start, aa_end, w feature name (e.g. wdomain1), impacts.
#' @param fixtheta Pre-calculated overdispersion (theta) parameter. This should be calculated using sitednds(., method="NB").
#' @param refdb Reference database (path to .rda file or a pre-loaded array object in the right format).
#' @param numcode NCBI genetic code number (default = 1; standard genetic code). To see the list of genetic codes supported use: ? seqinr::translate
#' @param normalisefromsyn Normalise the substitution rates based on the synonymous mutations in the gene. Using TRUE is recommended. Using FALSE uses the expected synonymous mutation rate of the gene from the dndscv negative binomial regression model (dndsout$genemuts).
#' @param syndrivers Vector of known synonymous driver sites defined by their aminoacid position, to be excluded from the background model (e.g. syndrivers = c("T125T","E224E","Q331Q") for TP53).
#' @param exon_flank_length Exon flank length in bp [default = 10]. Using a value higher than 0 will calculate a separate selection (w) coefficient for synonymous mutations in exon flanks.
#' @param intron_flank_length Intron flank length in bp [default = 10]. Intronic sites occurring within these flanks but not already classified as Essential_Splice will receive a separate w parameter.
#' @param sitefilename Optionally, provide a file name to save the table of all annotated sites in the gene. This table is also always contained in the output object.
#'
#' @return 'withingenednds' returns a list of objects:
#' @return - sites: Table with the annotation of all sites in the gene (from covtable), including all functional annotations in the "regions" input object as well as default annotations (Missense, Nonsense, Essential_Splice, Start_loss, Stop_loss, etc).
#' @return - par.pois: Poisson regression results (not recommended).
#' @return - par.nb: Negative binomial results fitting a new overdispersion parameter to the data (when fixtheta is not provided).
#' @return - par.nbfix: Negative binomial results using the input fixtheta value as recommended.
#' @return - model.pois: Poisson regression object.
#' @return - model.nb: Negative binomial regression object.
#' @return - model.nbfix: Negative binomial regression object.
#'
#' @export
withingenednds = function(mutations, gene, covtable, dndsout, genomeFile, regionschr = NULL, regionsaa = NULL, fixtheta = NULL, normalisefromsyn = TRUE, syndrivers = NULL, exon_flank_length = 10, intron_flank_length = 10, sitefilename = NULL, refdb = "hg19", numcode = 1) {
## 1. Environment
message("Running gene-level selection analyses...")
# [Input] Reference database
refdb_class = class(refdb)
if ("character" %in% refdb_class) {
if (refdb == "hg19") {
data("refcds_hg19", package="dndscv")
} else {
load(refdb)
}
} else if("array" %in% refdb_class) {
# use the user-supplied RefCDS object
RefCDS = refdb
} else {
stop("Expected refdb to be \"hg19\", a file path, or a RefCDS-formatted array object.")
}
# Annotating all possible coding changes in the positions provided in the input covtable
library("Rsamtools")
covtable$ref = as.vector(scanFa(genomeFile, GRanges(covtable$chr, IRanges(covtable$pos, covtable$pos))))
covtable$ref3 = as.vector(scanFa(genomeFile, GRanges(covtable$chr, IRanges(covtable$pos-1, covtable$pos+1))))
allsubs = data.frame(sampleID="allsubs", chr = rep(covtable$chr, each=3), pos = rep(covtable$pos, each=3), ref = rep(covtable$ref, each=3), mut = NA, ref3 = rep(covtable$ref3, each=3), mut3 = NA, ref_cod = NA, mut_cod = NA, ref3_cod = NA, mut3_cod = NA)
for (j in seq(1,nrow(allsubs),3)) {
allsubs$mut[j:(j+2)] = setdiff(c("A","C","G","T"),allsubs$ref[j])
}
allsubs$mut3 = paste(substr(allsubs$ref3,1,1), allsubs$mut, substr(allsubs$ref3,3,3), sep="")
aux = dndscv(allsubs, gene_list = gene, max_coding_muts_per_sample = Inf, max_muts_per_gene_per_sample = Inf, outp = 1)$annotmuts # Annotated mutations
# Adding duplex depth, functional impact annotation to all possible changes in the input table
aux$mstr = paste(aux$chr, aux$pos, aux$mut, sep=":") # mutation string ID
obssubs = dndsout$annotmuts[which(dndsout$annotmuts$ref %in% c("A","C","G","T") & dndsout$annotmuts$mut %in% c("A","C","G","T") & dndsout$annotmuts$gene == gene), ]
obssubs$mstr = paste(obssubs$chr, obssubs$pos, obssubs$mut, sep=":")
m$mstr = paste(m$chr, m$pos, m$mut, sep=":")
sites = cbind(data.frame(gene=gene, pid=aux$pid[1]), allsubs[,-1]) # Initialising the sites table
pos2dc = setNames(covtable$dc, covtable$pos)
sites$dc = pos2dc[as.character(sites$pos)]
sites$mstr = paste(sites$chr, sites$pos, sites$mut, sep=":")
sites$obs = table(m$mstr)[sites$mstr]; sites$obs[is.na(sites$obs)] = 0 # Annotating the number of observed mutations at each site
mstr2imp1 = setNames(aux$ntchange, aux$mstr)
mstr2imp2 = setNames(aux$aachange, aux$mstr)
mstr2imp3 = setNames(aux$impact, aux$mstr)
sites$ntchange = mstr2imp1[sites$mstr]
sites$aachange = mstr2imp2[sites$mstr]
sites$impact = mstr2imp3[sites$mstr]
# Adding strand and annotating the coding trinucleotide change
aux2 = unique(dndsout$annotmuts[,c("gene","strand")])
gene2strand = setNames(aux2[,2],aux2[,1])
sites$strand = gene2strand[sites$gene]
sites$ref_cod = sites$ref
sites$ref_cod[sites$strand==-1] = seqinr::comp(sites$ref[sites$strand==-1], forceToLower = F)
sites$mut_cod = sites$mut
sites$mut_cod[sites$strand==-1] = seqinr::comp(sites$mut[sites$strand==-1], forceToLower = F)
revcomp = function(seqvec) {
as.vector(sapply(seqvec, function(x) paste(seqinr::comp(rev(unlist(strsplit(x,split=""))), forceToLower=F), collapse=""))) # Reverse complement of a vector of sequence motifs
}
sites$ref3_cod = sites$ref3
sites$mut3_cod = sites$mut3
if (any(sites$strand==-1)) {
sites$ref3_cod[sites$strand==-1] = revcomp(sites$ref3[sites$strand==-1])
sites$mut3_cod[sites$strand==-1] = revcomp(sites$mut3[sites$strand==-1])
}
# Adding the trinucleotide substitution rates from the dndsout model
mle_submodel = setNames(dndsout$mle_submodel[,2], dndsout$mle_submodel[,1])
mle_submodel = c(mle_submodel, "TTT>TGT"=1) # Adding the missing rate (note that all rates in mle_submodel in dNdScv are relative to TTT>TGT)
mle_submodel = mle_submodel * mle_submodel["t"] # Absolute rates
sites$r = mle_submodel[paste(sites$ref3_cod, sites$mut3_cod, sep=">")]
sites$r = sites$r * sites$dc / mean(sites$dc) # Normalising the expected rates at a site by the duplex depth
# Normalising the global expected rates by the estimated mutation rate of the gene using one of two alternative models:
# 1. normalisefromsyn = TRUE. We normalise "r" using the synonymous mutations observed in the gene (excluding known syn driver sites)
# 2. normalisefromsyn = FALSE. We normalise "r" using the negative binomial model from dndscv.
if (normalisefromsyn == TRUE) {
mobs = sum(sites$obs[which(sites$impact=="Synonymous" & !(sites$aachange %in% syndrivers))])
mexp = sum(sites$r[which(sites$impact=="Synonymous" & !(sites$aachange %in% syndrivers))])
sites$rnorm = sites$r * mobs / mexp
} else {
message("Option not recommended: Normalising the mutation rate of the gene based on the negative regression model of dNdScv")
sites$rnorm = sites$r * dndsout$genemuts$exp_syn_cv[dndsout$genemuts$gene_name==gene] / dndsout$genemuts$exp_syn[dndsout$genemuts$gene_name==gene]
}
# Excluding sites with rate = 0 (e.g. due to lack of duplex coverage)
ratezero = which(sites$rnorm==0)
if (length(ratezero)>0) {
sites = sites[-ratezero, ]
message(sprintf("%0.0f sites have been excluded due to having a duplex depth and/or a predicted mutation rate of 0", length(ratezero)))
}
## 2. Creating an index regression matrix with all functional annotations (each row is a site and each column is a parameter in the selection model)
# Annotating Missense, Nonsense, Essential_Splice and Stop_loss mutations
mutmat = data.frame(t = rep(1,nrow(sites)), wmis = 0, wnon = 0, wspl = 0, wsyndriv = 0,
wintr = 0, woutcds = 0, wexfl = 0, winfl = 0, wnonlastex = 0, wstoploss = 0, wstartloss = 0)
mutmat$wmis[which(sites$impact=="Missense")] = 1
mutmat$wnon[which(sites$impact=="Nonsense")] = 1
mutmat$wspl[which(sites$impact=="Essential_Splice")] = 1
mutmat$wsyndriv[which(sites$impact=="Synonymous" & sites$aachange %in% syndrivers)] = 1
mutmat$wstoploss[which(sites$impact=="Stop_loss")] = 1
# Annotating Start_loss mutations (mutations in codon 1)
#start_loss = which(sites$impact != "Synonymous" & substr(sites$aachange,2,nchar(sites$aachange)-1) == "1") # There is no need to check for synonymous changes in ATG
start_loss = which(substr(sites$aachange,2,nchar(sites$aachange)-1) == "1") # There is no need to check for synonymous changes in ATG
mutmat$wstartloss[start_loss] = 1
sites$impact[start_loss] = "Start_loss"
# Annotating introns, exon flanks, and intron flanks
refcdsgene = RefCDS[[which(sapply(RefCDS, function(x) x$gene_name)==gene)]]
exons = refcdsgene$intervals_cds
esspl = refcdsgene$intervals_splice
gr_sites = GenomicRanges::GRanges(sites$gene, IRanges::IRanges(sites$pos,sites$pos))
gr_exons = GenomicRanges::GRanges(gene, IRanges::IRanges(exons[,1],exons[,2]))
gr_cds = GenomicRanges::GRanges(gene, IRanges::IRanges(min(exons),max(exons)))
gr_outcds = setdiff(gr_sites, gr_cds)
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_outcds, type="any", select="all"))
mutmat$woutcds[unique(ol[,1])] = 1 # Annotation of the intronic sites
# For genes with >1 exon, we also annotate essential splice sites, and intron and exon flanks.
if (length(esspl)>0) {
gr_esspl = GenomicRanges::GRanges(gene, IRanges::IRanges(esspl,esspl))
gr_introns = setdiff(setdiff(gr_cds, gr_exons), gr_esspl)
exfl = rbind(cbind(exons[2:nrow(exons),1],exons[2:nrow(exons),1]+exon_flank_length-1), cbind(exons[1:(nrow(exons)-1),2]-exon_flank_length+1,exons[1:(nrow(exons)-1),2]))
gr_exonflanks = GenomicRanges::GRanges(gene, IRanges::IRanges(exfl[,1],exfl[,2]))
gr_exonflanks = intersect(gr_exonflanks, gr_exons) # Ensuring that the exon flanks do not extend into introns
infl = rbind(cbind(exons[2:nrow(exons),1]-intron_flank_length,exons[2:nrow(exons),1]-1), cbind(exons[1:(nrow(exons)-1),2]+1,exons[1:(nrow(exons)-1),2]+intron_flank_length))
gr_intrflanks = GenomicRanges::GRanges(gene, IRanges::IRanges(infl[,1],infl[,2]))
gr_intrflanks = setdiff(gr_intrflanks, gr_exons) # Removing any overlaps with exonic sequences
gr_intrflanks = setdiff(gr_intrflanks, gr_esspl) # Removing any overlaps with essential splice site sequences
gr_intrflanks = intersect(gr_intrflanks, gr_cds) # Removing UTR sequences
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_introns, type="any", select="all"))
mutmat$wintr[unique(ol[,1])] = 1 # Annotation of the intronic sites
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_introns, type="any", select="all"))
mutmat$wintr[unique(ol[,1])] = 1 # Annotation of the intronic sites
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_exonflanks, type="any", select="all"))
mutmat$wexfl[seq(1,nrow(sites)) %in% ol[,1] & sites$impact=="Synonymous"] = 1 # Annotation of the exonic flank sites (only for synonymous mutations)
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_intrflanks, type="any", select="all"))
mutmat$winfl[unique(ol[,1])] = 1 # Annotation of the exonic flank sites
}
# Annotation of nonsense mutations in the last coding exon
if (nrow(exons)>1) {
if (refcdsgene$strand==1) {
lastexon = exons[nrow(exons),]
} else {
lastexon = exons[1,]
}
gr_lastexon = GenomicRanges::GRanges(gene, IRanges::IRanges(min(lastexon),max(lastexon)))
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_lastexon, type="any", select="all"))
mutmat$wnonlastex[seq(1,nrow(sites)) %in% ol[,1] & sites$impact=="Nonsense"] = 1
mutmat$wnon[seq(1,nrow(sites)) %in% ol[,1] & sites$impact=="Nonsense"] = 0 # Removing Nonsense mutations in the last exon from the wnon annotation
}
## 3. Other annotations from the user
# Regions defined by chr position
if (!is.null(regionschr)) {
wnames = unique(regionschr$wname)
badnames = intersect(wnames,unique(colnames(mutmat)))
if (length(badnames)>0) { stop(sprintf("The following w names are not allowed in the input regions as they match existing parameters: %s", paste(badnames,collapse = ","))) }
for (j in 1:length(wnames)) {
aux = regionschr[which(regionschr$wname==wnames[j]), ]
if (length(unique(aux$impacts))!=1) { stop("regionschr: different values found in the impacts column for a given feature, please correct your input object") }
imps = strsplit(unique(aux$impacts), split=",")[[1]]
if (any(!(imps %in% setdiff(unique(sites$impact),NA)))) { stop("regionschr: invalid impact terms used in the input object, please correct the impact column") }
gr_wname = GenomicRanges::GRanges(gene, IRanges::IRanges(aux$start,aux$end))
ol = as.data.frame(GenomicRanges::findOverlaps(gr_sites, gr_wname, type="any", select="all"))
mutmat[,wnames[j]] = 0 # Initialising this field in the data frame
if (length(imps)==0) {
indx = unique(ol[,1]) # If no impacts are indicated, we consider all sites independent of their impact
} else {
indx = intersect(unique(ol[,1]), which(sites$impact %in% imps))
}
mutmat[indx,wnames[j]] = 1
# Removing previous w parameters if layered=0
if (aux$layered[1]==0 | aux$layered[1]==FALSE) {
mutmat[indx, setdiff(colnames(mutmat),c("t",wnames[j]))] = 0 # Removing previous annotations at these sites
}
}
}
# Regions defined by aminoacid position
if (!is.null(regionsaa)) {
sites$aux = 1:nrow(sites)
aas = sites[which(!is.na(sites$aachange) & sites$aachange!="."),]
aas$aapos = as.numeric(substr(aas$aachange,2,nchar(aas$aachange)-1))
gr_aas = GenomicRanges::GRanges(gene, IRanges::IRanges(aas$aapos,aas$aapos))
wnames = unique(regionsaa$wname)
badnames = intersect(wnames,unique(colnames(mutmat)))
if (length(badnames)>0) { stop(sprintf("The following w names are not allowed in the input regions as they match existing parameters: %s", paste(badnames,collapse = ", "))) }
for (j in 1:length(wnames)) {
aux = regionsaa[which(regionsaa$wname==wnames[j]), ]
if (length(unique(aux$impacts))!=1) { stop("regionschr: different values found in the impacts column for a given feature, please correct your input object") }
imps = strsplit(unique(aux$impacts), split=",")[[1]]
if (any(!(imps %in% setdiff(unique(sites$impact),NA)))) { stop("regionschr: invalid impact terms used in the input object, please correct the impact column") }
gr_wname = GenomicRanges::GRanges(gene, IRanges::IRanges(aux$start,aux$end))
ol = as.data.frame(GenomicRanges::findOverlaps(gr_aas, gr_wname, type="any", select="all"))
mutmat[,wnames[j]] = 0 # Initialising this field in the data frame
if (length(imps)==0) {
indx = aas$aux[unique(ol[,1])] # If no impacts are indicated, we consider all sites independent of their impact
} else {
indx = aas$aux[intersect(unique(ol[,1]), which(aas$impact %in% imps))]
}
mutmat[indx,wnames[j]] = 1
# Removing previous w parameters if layered=0
if (aux$layered[1]==0 | aux$layered[1]==FALSE) {
mutmat[indx, setdiff(colnames(mutmat),c("t",wnames[j]))] = 0 # Removing previous annotations at these sites
}
}
sites = sites[, setdiff(colnames(sites),"aux")]
}
## 4. Poisson regression: fitting the selection model
model = glm(formula = sites$obs ~ offset(log(sites$rnorm)) + . -1, data=mutmat, family=poisson(link=log))
mle = exp(coefficients(model)) # Maximum-likelihood estimates for the rate params
pvals = coef(summary(model))[,4]
model.lrt = drop1(model, test= "Chisq")
pvals.lrt = setNames(model.lrt[[5]], row.names(model.lrt))
ci = exp(confint.default(model)) # Wald confidence intervals
par.pois = data.frame(name=gsub("\`","",rownames(ci)), mle=mle[rownames(ci)], cilow=ci[,1], cihigh=ci[,2], pval.wald=pvals[rownames(ci)], pval.lrt=pvals.lrt[rownames(ci)]) # MLEs and Wald CI95% for the selection parameters
# Adding obs/exp statistics to the regression outputs
nobs = apply(mutmat, 2, function(x) sum(sites$obs[x==1]))
nexp = apply(mutmat, 2, function(x) sum(sites$rnorm[x==1]))
par.pois$obs = nobs[par.pois$name]
par.pois$exp = nexp[par.pois$name]
## 5. Negative binomial regression
model.nbfix = model.nb = par.nbfix = par.nb = NULL
if (!is.null(fixtheta)) {
model.nbfix = glm(formula = sites$obs ~ offset(log(sites$rnorm)) + . -1, data=mutmat, family=MASS::negative.binomial(fixtheta))
mle = exp(coefficients(model.nbfix)) # Maximum-likelihood estimates for the rate params
pvals = coef(summary(model.nbfix))[,4]
model.lrt = drop1(model.nbfix, test= "Chisq")
pvals.lrt = setNames(model.lrt[[5]], row.names(model.lrt))
ci = exp(confint.default(model.nbfix)) # Wald confidence intervals
par.nbfix = data.frame(name=gsub("\`","",rownames(ci)), mle=mle[rownames(ci)], cilow=ci[,1], cihigh=ci[,2], pval.wald=pvals[rownames(ci)], pval.lrt=pvals.lrt[rownames(ci)]) # MLEs and Wald CI95% for the selection parameters
par.nbfix$obs = nobs[par.nbfix$name]
par.nbfix$exp = nexp[par.nbfix$name]
} else {
model.nb = MASS::glm.nb(formula = as.vector(sites$obs) ~ offset(log(sites$rnorm)) + . -1, data=mutmat)
mle = exp(coefficients(model.nb)) # Maximum-likelihood estimates for the rate params
pvals = coef(summary(model.nb))[,4]
model.lrt = drop1(model.nb, test= "Chisq")
pvals.lrt = setNames(model.lrt[[5]], row.names(model.lrt))
ci = exp(confint.default(model.nb)) # Wald confidence intervals
par.nb = data.frame(name=gsub("\`","",rownames(ci)), mle=mle[rownames(ci)], cilow=ci[,1], cihigh=ci[,2], pval.wald=pvals[rownames(ci)], pval.lrt=pvals.lrt[rownames(ci)]) # MLEs and Wald CI95% for the selection parameters
par.nb$obs = nobs[par.nb$name]
par.nb$exp = nexp[par.nb$name]
}
## 6. Outputs
# Annotated sites table
sites2 = cbind(sites, mutmat)
if (!is.null(sitefilename)) {
sites2 = apply(sites2,2,as.character)
write.table(sites2, file=sitefilename, col.names = T, row.names = F, quote = F, sep = "\t")
}
# Inform the user about the sites used as neutral reference by the model
neutral_sites = (rowSums(mutmat[,setdiff(colnames(mutmat),"t")])==0)
neutral_nsyn = sum(sites$obs[which(neutral_sites & sites$impact=="Synonymous")]) # Number of synonymous mutations used as background for dN/dS
neutral_nsynsites = length(which(neutral_sites & sites$impact=="Synonymous"))
neutral_othersites = length(which(neutral_sites & sites$impact!="Synonymous"))
nonneutral_nsyn = sum(sites$obs[which(neutral_sites==0 & sites$impact=="Synonymous")]) # Number of synonymous mutations excluded from the neutral backgroup by the "w" annotations
message(sprintf(" Sites used as neutral reference: %0.0f synonymous mutations across %0.0f synonymous sites.", neutral_nsyn, neutral_nsynsites))
if (neutral_othersites>0) {
message(sprintf(" %0.0f sites not classified as synonymous coding sites were used in the background model. Please ensure that this was intended.", neutral_othersites))
}
if (nonneutral_nsyn>0) {
message(sprintf(" %0.0f synonymous mutations were not used in the neutral background model. This may be due to excluding known synonymous driver mutations (when using syndrivers), excluding synonymous mutations in exon flanks (when using exon_flank_length>0) and/or due to annotations provided by the user. Please ensure that this is the desired behaviour.", nonneutral_nsyn))
}
# Output object
out = list(sites = sites2, dnds.pois = par.pois, dnds.nb = par.nb, dnds.nbfix = par.nbfix, model.pois = model, model.nb = model.nb, model.nbfix = model.nbfix)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.