Skip to content

Commit ad99ede

Browse files
hpageshpages
committed
igblastr 1.1.8: Add extract_codons() and translate_codons()
Also rename argument 'max_codons' of function translate_fwr4() to 'max.codons'.
1 parent 5f5e6df commit ad99ede

7 files changed

Lines changed: 224 additions & 38 deletions

File tree

‎DESCRIPTION‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ Description: The igblastr package provides functions to conveniently install
77
biocViews: Immunology, Immunogenetics, ImmunoOncology, CellBiology
88
URL: https://bioconductor.org/packages/igblastr
99
BugReports: https://github.com/HyrienLab/igblastr/issues
10-
Version: 1.1.7
10+
Version: 1.1.8
1111
License: Artistic-2.0
1212
Encoding: UTF-8
1313
Authors@R: c(
@@ -31,6 +31,7 @@ VignetteBuilder: knitr
3131
Collate:
3232
utils.R
3333
long_to_wide_airr.R
34+
translate_codons.R
3435
file-utils.R
3536
loci-utils.R
3637
db-utils.R

‎NAMESPACE‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -71,6 +71,10 @@ export(
7171
###
7272

7373
export(
74+
## translate_codons.R:
75+
extract_codons,
76+
translate_codons,
77+
7478
## get_igblast_root.R:
7579
get_igblast_root, set_igblast_root,
7680

‎R/auxdata-utils.R‎

Lines changed: 16 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -138,16 +138,6 @@ load_igblast_auxiliary_data <- function(...)
138138
setNames(auxdata_col[match(J_names, auxdata_allele_name)], J_names)
139139
}
140140

141-
.translate_codons <- function(dna, start)
142-
{
143-
stopifnot(is(dna, "DNAStringSet"), is.integer(start),
144-
length(dna) == length(start))
145-
end <- -1L - (width(dna) - start + 1L) %% 3L
146-
codons <- narrow(dna, start=start, end=end)
147-
aa <- translate(codons, no.init.codon=TRUE, if.fuzzy.codon="solve")
148-
as.character(aa)
149-
}
150-
151141
### Translate the coding frame.
152142
### Only needs access to the "coding_frame_start" column in 'auxdata'.
153143
### Returns the amino acid sequences in a named character vector that
@@ -161,9 +151,10 @@ translate_J_alleles <- function(J_alleles, auxdata)
161151
ans <- rep.int(NA_character_, length(J_alleles))
162152
selection_idx <- which(!is.na(coding_frame_start))
163153
if (length(selection_idx) != 0L) {
164-
alleles <- J_alleles[selection_idx]
165-
start <- coding_frame_start[selection_idx] + 1L # 0-based to 1-based
166-
ans[selection_idx] <- .translate_codons(alleles, start=start)
154+
dna <- J_alleles[selection_idx]
155+
offset <- coding_frame_start[selection_idx]
156+
aa <- translate_codons(dna, offset=offset)
157+
ans[selection_idx] <- as.character(aa)
167158
}
168159
setNames(ans, names(J_alleles))
169160
}
@@ -186,23 +177,24 @@ J_allele_has_stop_codon <- function(J_alleles, auxdata)
186177
### is parallel to 'J_alleles' and has the allele names on it.
187178
### The returned vector will contain an NA for any allele that is not
188179
### annotated in 'auxdata' or for which 'auxdata$cdr3_end' has an NA.
189-
translate_fwr4 <- function(J_alleles, auxdata, max_codons=NA)
180+
translate_fwr4 <- function(J_alleles, auxdata, max.codons=NA)
190181
{
191-
if (!isSingleNumberOrNA(max_codons))
192-
stop(wmsg("'max_codons' must be a single number or NA"))
193-
if (!is.integer(max_codons))
194-
max_codons <- as.integer(max_codons)
182+
if (!isSingleNumberOrNA(max.codons))
183+
stop(wmsg("'max.codons' must be a single number or NA"))
184+
if (!is.integer(max.codons))
185+
max.codons <- as.integer(max.codons)
195186

196-
cdr3_end <- .query_auxdata(auxdata, J_alleles, "cdr3_end")
187+
cdr3_end <- .query_auxdata(auxdata, J_alleles, "cdr3_end") # 0-based
197188
ans <- rep.int(NA_character_, length(J_alleles))
198189
selection_idx <- which(!is.na(cdr3_end))
199190
if (length(selection_idx) != 0L) {
200-
alleles <- J_alleles[selection_idx]
201-
start <- cdr3_end[selection_idx] + 2L # 1-based FWR4 start
202-
ans[selection_idx] <- .translate_codons(alleles, start=start)
191+
dna <- J_alleles[selection_idx]
192+
offset <- cdr3_end[selection_idx] + 1L # 0-based FWR4 start
193+
aa <- translate_codons(dna, offset=offset)
194+
ans[selection_idx] <- as.character(aa)
203195
}
204-
if (!is.na(max_codons))
205-
ans <- substr(ans, 1L, max_codons)
196+
if (!is.na(max.codons))
197+
ans <- substr(ans, 1L, max.codons)
206198
setNames(ans, names(J_alleles))
207199
}
208200

‎R/translate_codons.R‎

Lines changed: 70 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,70 @@
1+
### =========================================================================
2+
### extract_codons() and translate_codons()
3+
### -------------------------------------------------------------------------
4+
###
5+
6+
### Normalize a single offset.
7+
.normarg_offset1 <- function(offset, dna_len)
8+
{
9+
if (!isSingleNumber(offset))
10+
stop(wmsg("'offset' must be a single number"))
11+
if (!is.integer(offset))
12+
offset <- as.integer(offset)
13+
if (offset < 0L || offset > dna_len)
14+
stop(wmsg("'offset' must be >= 0 and <= 'length(dna)'"))
15+
offset
16+
}
17+
18+
### Normalize a vector of offsets.
19+
.normarg_offset2 <- function(offset, widths)
20+
{
21+
if (!is.numeric(offset))
22+
stop(wmsg("'offset' must be numeric"))
23+
if (!is.integer(offset))
24+
offset <- as.integer(offset)
25+
noffset <- length(offset)
26+
nseq <- length(widths)
27+
if (noffset == 1L) {
28+
if (is.na(offset) || offset < 0L)
29+
stop(wmsg("'offset' cannot be NA or negative"))
30+
offset <- rep.int(offset, nseq)
31+
if (any(offset > widths))
32+
stop(wmsg("'offset' must be less than the length ",
33+
"of the shortest sequence in 'dna'"))
34+
} else {
35+
if (noffset != nseq)
36+
stop(wmsg("'offset' must be a single value or ",
37+
"a vector with one value per sequence in 'dna'"))
38+
if (anyNA(offset) || any(offset < 0L))
39+
stop(wmsg("'offset' cannot contain NAs or negative values"))
40+
if (any(offset > widths))
41+
stop(wmsg("each value in 'offset' must be less than the ",
42+
"length of the corresponding sequence in 'dna'"))
43+
}
44+
offset
45+
}
46+
47+
extract_codons <- function(dna, offset=0)
48+
{
49+
if (is(dna, "DNAStringSet")) {
50+
dna_nchar <- nchar(dna) # same as width(dna)
51+
offset <- .normarg_offset2(offset, dna_nchar)
52+
} else if (is(dna, "DNAString")) {
53+
dna_nchar <- nchar(dna) # same as length(dna)
54+
offset <- .normarg_offset1(offset, dna_nchar)
55+
} else {
56+
stop(wmsg("'dna' must be DNAStringSet or DNAString object"))
57+
}
58+
width <- dna_nchar - offset
59+
width <- width - width %% 3L
60+
subseq(dna, start=offset+1L, width=width)
61+
}
62+
63+
translate_codons <- function(dna, offset=0, with.init.codon=FALSE)
64+
{
65+
if (!isTRUEorFALSE(with.init.codon))
66+
stop(wmsg("'with.init.codon' must be TRUE or FALSE"))
67+
codons <- extract_codons(dna, offset)
68+
translate(codons, no.init.codon=!with.init.codon, if.fuzzy.codon="solve")
69+
}
70+

‎man/auxdata-utils.Rd‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,7 @@ load_auxdata(organism, which=c("live", "original"))
2626

2727
translate_J_alleles(J_alleles, auxdata)
2828
J_allele_has_stop_codon(J_alleles, auxdata)
29-
translate_fwr4(J_alleles, auxdata, max_codons=NA)
29+
translate_fwr4(J_alleles, auxdata, max.codons=NA)
3030

3131
compute_auxdata(J_alleles)
3232
}
@@ -57,9 +57,9 @@ compute_auxdata(J_alleles)
5757
A data.frame as returned by \code{load_auxdata()} or
5858
\code{compute_auxdata()}.
5959
}
60-
\item{max_codons}{
60+
\item{max.codons}{
6161
The maximum number of FWR4 codons to translate. By default (i.e. when
62-
\code{max_codons} is \code{NA}) all the FWR4 codons are translated.
62+
\code{max.codons} is \code{NA}) all the FWR4 codons are translated.
6363
}
6464
}
6565

@@ -134,8 +134,8 @@ compute_auxdata(J_alleles)
134134
\item \link[Biostrings]{DNAStringSet} and \link[Biostrings]{AAStringSet}
135135
objects in the \pkg{Biostrings} package.
136136
137-
\item The \code{\link[Biostrings]{translate}} function in the
138-
\pkg{Biostrings} package.
137+
\item The \code{\link{translate_codons}} function on which
138+
\code{translate_J_alleles()} and \code{translate_fwr4()} are based.
139139
140140
\item The \code{\link{igblastn}} function to run the \code{igblastn}
141141
\emph{standalone executable} included in IgBLAST from R. This
@@ -178,13 +178,13 @@ stopifnot(setequal(names(J_alleles), rabbit_auxdata$allele_name))
178178
179179
## Note that this might change with future IMGT releases.
180180
181-
## Let's put the allele sequences in 'J_alleles' in the same order than
181+
## Let's put the allele sequences in 'J_alleles' in the same order as
182182
## in 'rabbit_auxdata':
183183
J_alleles <- J_alleles[rabbit_auxdata$allele_name]
184184
stopifnot(identical(names(J_alleles), rabbit_auxdata$allele_name))
185185

186186
## The 'coding_frame_start' column in 'rabbit_auxdata' contains integer
187-
## values that are >= 0 and < 3. They indicate how many nucleotides
187+
## values that are >= 0 and <= 2. They indicate how many nucleotides
188188
## precede the first codon on each allele sequence. In other words,
189189
## this is the number of nucleotides that we need to trim on the 5'
190190
## end of the germline J allele sequence before we start translating
@@ -193,7 +193,7 @@ stopifnot(identical(names(J_alleles), rabbit_auxdata$allele_name))
193193
J_aa <- translate_J_alleles(J_alleles, rabbit_auxdata)
194194
J_aa
195195

196-
## No sequence in 'J_aa' should contain the letter "*" that is used
196+
## No sequence in 'J_aa' should contain the letter "*" which is used
197197
## by translate_J_alleles() to represent a stop codon. However, one
198198
## J allele in IMGT-202531-1.Oryctolagus_cuniculus.IGH+IGK+IGL seems
199199
## to disobey:
@@ -213,7 +213,7 @@ J_aa[has_stop_codon] # indeed!
213213

214214
## Let's use translate_fwr4() to extract and translate the first 4
215215
## codons of the FWR4 region:
216-
fwr4_head <- translate_fwr4(J_alleles, rabbit_auxdata, max_codons=4)
216+
fwr4_head <- translate_fwr4(J_alleles, rabbit_auxdata, max.codons=4)
217217

218218
## We expect to see the "WGXG" and "FGXG" motifs here, and most of the
219219
## time we do:

‎man/translate_codons.Rd‎

Lines changed: 119 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,119 @@
1+
\name{translate_codons}
2+
3+
\alias{translate_codons}
4+
\alias{extract_codons}
5+
6+
\title{Extract and translate codons from a set of DNA sequences}
7+
8+
\description{
9+
The \code{translate_codons()} function extracts and translates codons
10+
from a set of DNA sequences.
11+
}
12+
13+
\usage{
14+
translate_codons(dna, offset=0, with.init.codon=FALSE)
15+
16+
## Used internally by translate_codons():
17+
extract_codons(dna, offset=0)
18+
}
19+
20+
\arguments{
21+
\item{dna}{
22+
A \link[Biostrings]{DNAStringSet} (or \link[Biostrings]{DNAString})
23+
object containing the codons to translate.
24+
}
25+
\item{offset}{
26+
The number of nucleotides that precede the first codon to translate.
27+
This must be supplied as a numeric vector with one value per
28+
sequence in \code{dna}, or as a single value. If the latter,
29+
then the same offset is used for all sequences.
30+
}
31+
\item{with.init.codon}{
32+
Is the first codon to translate in each DNA sequence the initiation
33+
codon? By default, \code{with.init.codon} is set to \code{FALSE}, in
34+
which case \code{translate_codons()} assumes that the first codon
35+
to translate in each DNA sequence is \emph{not} the initiation codon.
36+
37+
See the \code{no.init.codon} argument in
38+
\code{?\link[Biostrings]{translate}} in the \pkg{Biostrings} package
39+
for more information.
40+
}
41+
}
42+
43+
\value{
44+
\code{translate_codons()} returns an \link[Biostrings]{AAStringSet}
45+
object with one amino acid sequence per input sequence.
46+
47+
\code{extract_codons()} returns a \link[Biostrings]{DNAStringSet}
48+
object with one sequence per input sequence. The output sequences
49+
are obtained by trimming the original sequences as follow:
50+
\itemize{
51+
\item On their 5' end, sequences are trimmed by the amount of
52+
nucleotides specified in \code{offset}.
53+
\item On their 3' end, sequences are trimmed by the smallest
54+
amount of nucleotides that makes the length of the trimmed
55+
sequence a multiple of 3. Note that this will always be 0, 1,
56+
or 2 nucleotides.
57+
}
58+
}
59+
60+
\seealso{
61+
\itemize{
62+
\item \link[Biostrings]{DNAStringSet} and \link[Biostrings]{AAStringSet}
63+
objects in the \pkg{Biostrings} package.
64+
65+
\item The \code{\link[Biostrings]{translate}} function in the
66+
\pkg{Biostrings} package on which \code{translate_codons()}
67+
is based.
68+
69+
\item \code{\link{list_germline_dbs}} to list all the \emph{cached
70+
germline dbs}, that is, all the germline databases currently
71+
installed in \pkg{igblastr}'s persistent cache.
72+
}
73+
}
74+
75+
\examples{
76+
## ---------------------------------------------------------------------
77+
## translate_codons()
78+
## ---------------------------------------------------------------------
79+
80+
## Load germline V gene allele sequences for human:
81+
list_germline_dbs()
82+
db_name <- "_AIRR.human.IGH+IGK+IGL.202410"
83+
V_alleles <- load_germline_db(db_name, region_types="V")
84+
V_alleles # DNAStringSet object
85+
86+
## Translate them:
87+
V_aa <- translate_codons(V_alleles)
88+
V_aa # AAStringSet object
89+
90+
## Some human germline V gene allele sequences have a stop codon:
91+
has_stop_codon <- grepl("*", as.character(V_aa), fixed=TRUE)
92+
V_aa[has_stop_codon]
93+
94+
## Handling of initiation codons (see '?translate' in the Biostrings
95+
## package for how initiation codons are handled):
96+
dna2 <- DNAStringSet(c("TTGTCCTTTATA", "GAATCATTTATC", "CTGTCGTTTATT"))
97+
translate_codons(dna2)
98+
translate_codons(dna2, with.init.codon=TRUE)
99+
100+
## ---------------------------------------------------------------------
101+
## extract_codons()
102+
## ---------------------------------------------------------------------
103+
dna <- DNAStringSet(c("CCCAAAGGGTTT",
104+
"CCAAAGGGTTT",
105+
"CAAAGGGTTT",
106+
"AAAGGGTTT"))
107+
extract_codons(dna)
108+
extract_codons(dna, offset=1)
109+
extract_codons(dna, offset=2)
110+
extract_codons(dna, offset=3)
111+
extract_codons(dna, offset=4)
112+
113+
extract_codons(dna, offset=3:0)
114+
extract_codons(dna, offset=4:1)
115+
extract_codons(dna, offset=5:2)
116+
extract_codons(dna, offset=6:3)
117+
}
118+
119+
\keyword{utilities}

‎tests/testthat/test-auxdata-utils.R‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -128,7 +128,7 @@ test_that("translate_fwr4()", {
128128
expect_true(is.character(fwr4_aa))
129129
expect_identical(names(fwr4_aa), names(J_alleles))
130130
expect_false(anyNA(fwr4_aa))
131-
fwr4_head <- translate_fwr4(J_alleles, auxdata, max_codons=4L)
131+
fwr4_head <- translate_fwr4(J_alleles, auxdata, max.codons=4L)
132132
expect_true(is.character(fwr4_head))
133133
expect_identical(names(fwr4_head), names(J_alleles))
134134
expect_false(anyNA(fwr4_head))
@@ -138,7 +138,7 @@ test_that("translate_fwr4()", {
138138
db_name <- "IMGT-202531-1.Homo_sapiens.IGH+IGK+IGL"
139139
J_alleles <- load_germline_db(db_name, region_types="J")
140140
allele_is_known <- names(J_alleles) %in% auxdata$allele_name
141-
fwr4_head <- translate_fwr4(J_alleles, auxdata, max_codons=4L)
141+
fwr4_head <- translate_fwr4(J_alleles, auxdata, max.codons=4L)
142142
expect_identical(unname(is.na(fwr4_head)), !allele_is_known)
143143
expect_true(all(grepl("^[WF]G.G$", fwr4_head[allele_is_known])))
144144

@@ -154,7 +154,7 @@ test_that("translate_fwr4()", {
154154
expect_true(sum(!allele_is_known) <= 3L)
155155
## Get rid of the "unknown" J alleles.
156156
known_J_alleles <- J_alleles[allele_is_known]
157-
fwr4_head <- translate_fwr4(known_J_alleles, auxdata, max_codons=4L)
157+
fwr4_head <- translate_fwr4(known_J_alleles, auxdata, max.codons=4L)
158158
expect_false(anyNA(fwr4_head))
159159
## 3 "known" mouse J alleles in IMGT release 202531-1 don't have
160160
## the expected motif at the beginning of their FWR4 region.
@@ -169,7 +169,7 @@ test_that("translate_fwr4()", {
169169

170170
db_name <- "IMGT-202531-1.Rattus_norvegicus.IGH+IGK+IGL"
171171
J_alleles <- load_germline_db(db_name, region_types="J")
172-
fwr4_head <- translate_fwr4(J_alleles, auxdata, max_codons=4L)
172+
fwr4_head <- translate_fwr4(J_alleles, auxdata, max.codons=4L)
173173
## translate_fwr4() uses 'auxdata$cdr3_end' to get the position of
174174
## the first FWR4 codons, but this column has an NA for IGKJ3*01.
175175
## This could change in the future.

0 commit comments

Comments
 (0)