# MSA2dist Vignette

Kristian K Ullrich

#### 2022-04-26

#### Abstract

MSA2dist calculates pairwise distances between all sequences of a DNAStringSet or a AAStringSet using a custom score matrix and conducts codon based analysis. It uses scoring matrices to be used in these pairwise distance calcualtions which can be adapted to any scoring for DNA or AA characters. E.g. by using literal distances MSA2dist calcualtes pairwise IUPAC distances. DNAStringSet alignments can be analysed as codon alignments to look for synonymous and nonsynonymous substitutions in a parallelised fashion.

# Contents

# 1 Introduction

Calculating pairwise distances of either DNA or AA sequences is a common task for evolutionary biologist. The distance calculations are either based on specific nucleotide, codon or amino acid models or on a scoring matrix.

Note: Sequences need to be pre-aligned into so called multiple sequence alignments (MSA), which can be done with a multitude of existing software. Just to mention for example mafft, muscle or the R package `msa`.

The R package `ape` for example offers the `ape::dist.dna()` function, which has implemented a collection of different evolutionary models. `MSA2dist` extends the possibility to directly calculate pairwise nucloetide distances of an `Biostrings::DNAStringSet` object or pairwise amino acid distances of an `Biostrings::AAStringSet` object. The scoring matrix based calculations are implemented in `c++` with `RcppThread` to parallelise pairwise combinations.

It is a non-trivial part to resolve haploid (1n) sequences from a diploid (2n) individual (aka phasing) to further use the haploid sequences for distance calculations. To cope with this situation, `MSA2dist` uses a literal distance (Chang et al. 2017) which can be directly applied on `IUPAC` nucleotide ambiguity encoded sequences with the `dnastring2dist()` function. `IUPAC` sequences can be for example obtained directly from mapped `BAM` files and the angsd `-doFasta 4` option (Korneliussen, Albrechtsen, and Nielsen 2014).

The Grantham’s score (Grantham 1974) attempts to predict the distance between two amino acids, in an evolutionary sense considering the amino acid composition, polarity and molecular volume. `MSA2dist` offers with the `aastring2dist()` function the possibility to obtain pairwise distances of all sequences in an `Biostrings::AAStringSet` (needs to be pre-aligned). The resulting distance matrix can be used to calculate neighbor-joining trees via the R package `ape`.

Calculating synonymous (Ks) and nonsynonymous (Ka) substitutions from coding sequences and its ratio Ka/Ks can be used as an indicator of selective pressure acting on a protein. The `dnastring2kaks()` function can be applied on pre-aligned `Biostrings::DNAStringSet()` objects to calculate these values either according to (Li 1993) via the R package `seqinr` or according to the model of (Nei and Gojobori 1986).

Further, all codons can be evaluated among the coding sequence alignment and be plotted to for example protein domains with substitutions or indels with the `codonmat2xy()` function.

# 2 Installation

To install this package, start R (version “4.1”) and enter:

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")

    BiocManager::install("MSA2dist")

# 3 Load MSA2dist

    # load MSA2dist
    library(MSA2dist)
    # load example data
    data(hiv, package="MSA2dist")
    data(AAMatrix, package="MSA2dist")
    data(woodmouse, package="ape")

# 4 Sequence Format conversion

To be able to use distance calculation functions from other R packages, like `ape` or `seqinr`, it is necessary to have dedicated sequence format conversion functions. Here, some examples are shown, how to convert from and to a `Biostrings::DNAStringSet` object.

`?Biostrings::DNAStringSet()` >>> `?seqinr::as.alignment()`

    ## define two cds sequences
    cds1 <- Biostrings::DNAString("ATGCAACATTGC")
    cds2 <- Biostrings::DNAString("ATG---CATTGC")
    cds1.cds2.aln <- c(Biostrings::DNAStringSet(cds1),
        Biostrings::DNAStringSet(cds2))
    ## define names
    names(cds1.cds2.aln) <- c("seq1", "seq2")
    ## convert into alignment
    cds1.cds2.aln |> dnastring2aln()

    ## $nb
    ## [1] 2
    ##
    ## $nam
    ## [1] "seq1" "seq2"
    ##
    ## $seq
    ## [1] "atgcaacattgc" "atg---cattgc"
    ##
    ## $com
    ## [1] NA
    ##
    ## attr(,"class")
    ## [1] "alignment"

`?seqinr::as.alignment()` >>> `?Biostrings::DNAStringSet()`

    ## convert back into DNAStringSet
    cds1.cds2.aln |> dnastring2aln() |> aln2dnastring()

    ## DNAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]    12 ATGCAACATTGC                                      seq1
    ## [2]    12 ATG---CATTGC                                      seq2

`?Biostrings::DNAStringSet()` >>> `?ape::DNAbin()`

    ## convert into alignment
    cds1.cds2.aln |> dnastring2dnabin()

    ## 2 DNA sequences in binary format stored in a matrix.
    ##
    ## All sequences of same length: 12
    ##
    ## Labels:
    ## seq1
    ## seq2
    ##
    ## Base composition:
    ##     a     c     g     t
    ## 0.286 0.238 0.190 0.286
    ## (Total: 24 bases)

`?ape::DNAbin()` >>> `?Biostrings::DNAStringSet()`

    ## convert back into DNAStringSet
    cds1.cds2.aln |> dnastring2dnabin() |> dnabin2dnastring()

    ## DNAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]    12 ATGCAACATTGC                                      seq1
    ## [2]    12 ATG---CATTGC                                      seq2

    ## use woodmouse data
    woodmouse |> dnabin2dnastring()

    ## DNAStringSet object of length 15:
    ##      width seq                                              names
    ##  [1]   965 NTTCGAAAAACACACCCACTACT...GCCCAATTACTCAGACCCTATA No305
    ##  [2]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCTGTN No304
    ##  [3]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCTATA No306
    ##  [4]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAATACNNNN No0906S
    ##  [5]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCNNNN No0908S
    ##  ...   ... ...
    ## [11]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCNNNN No1007S
    ## [12]   965 NNNNNNNNNNNNNNNNNNNNNNN...GCCCAATTACTCAAACCCNNNN No1114S
    ## [13]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCNNNN No1202S
    ## [14]   965 ATTCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCNNNN No1206S
    ## [15]   965 NNNCGAAAAACACACCCACTACT...GCCCAATTACTCAAACCCNNNN No1208S

`?Biostrings::AAStringSet()` >>> `?seqinr::as.alignment()`

    ## translate cds into aa
    aa1.aa2.aln <- cds1.cds2.aln |> cds2aa()
    ## convert into alignment
    aa1.aa2.aln |> aastring2aln()

    ## $nb
    ## [1] 2
    ##
    ## $nam
    ## [1] "seq1" "seq2"
    ##
    ## $seq
    ## [1] "mqhc" "mxhc"
    ##
    ## $com
    ## [1] NA
    ##
    ## attr(,"class")
    ## [1] "alignment"

`?seqinr::as.alignment()` >>> `?Biostrings::AAStringSet()`

    ## convert back into AAStringSet
    aa1.aa2.aln |> aastring2aln() |> aln2aastring()

    ## AAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]     4 MQHC                                              seq1
    ## [2]     4 MXHC                                              seq2

`?Biostrings::AAStringSet()` >>> `?ape::as.AAbin()`

    ## convert into AAbin
    aa1.aa2.aln |> aastring2aabin()

    ## 2 amino acid sequences in a list
    ##
    ## All sequences of the same length: 4

`?ape::as.AAbin()` >>> `?Biostrings::AAStringSet()`

    ## convert back into AAStringSet
    aa1.aa2.aln |> aastring2aabin() |> aabin2aastring()

    ## AAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]     4 MQHC                                              seq1
    ## [2]     4 MXHC                                              seq2

# 5 Frame aware `Biostrings::DNAStringSet` translation (`cds2aa()`)

To be able to translate a coding sequence into amino acids, sometimes the sequences do not start at the first frame. The `cds2aa` function can take an alternative codon start site into account (`frame = 1` or `frame = 2` or `frame = 3`). However, sometimes it is also necessary that the resulting coding sequence length is a multiple of three. This can be forced by using the `shorten = TRUE` option.

Simple translation:

    ## define two cds sequences
    cds1 <- Biostrings::DNAString("ATGCAACATTGC")
    cds2 <- Biostrings::DNAString("ATG---CATTGC")
    cds1.cds2.aln <- c(Biostrings::DNAStringSet(cds1),
        Biostrings::DNAStringSet(cds2))
    ## define names
    names(cds1.cds2.aln) <- c("seq1", "seq2")
    ## translate cds into aa
    cds1.cds2.aln |> cds2aa()

    ## AAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]     4 MQHC                                              seq1
    ## [2]     4 MXHC                                              seq2

    aa1.aa2.aln <- cds1.cds2.aln |> cds2aa()

Translation keeping multiple of three sequence length:

    ## translate cds into aa using frame = 2
    ## result is empty due to not multiple of three
    cds1.cds2.aln |> cds2aa(frame=2)

    ## AAStringSet object of length 0

    ## translate cds into aa using frame = 2 and shorten = TRUE
    cds1.cds2.aln |> cds2aa(frame=2, shorten=TRUE)

    ## AAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]     3 CNI                                               seq1
    ## [2]     3 XXI                                               seq2

    ## translate cds into aa using frame = 3 and shorten = TRUE
    cds1.cds2.aln |> cds2aa(frame=3, shorten=TRUE)

    ## AAStringSet object of length 2:
    ##     width seq                                               names
    ## [1]     3 ATL                                               seq1
    ## [2]     3 XXL                                               seq2

    ## use woodmouse data
    woodmouse |> dnabin2dnastring() |> cds2aa(shorten=TRUE)

    ## AAStringSet object of length 15:
    ##      width seq                                              names
    ##  [1]   321 XRKTHPLLKXISHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTL No305
    ##  [2]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTL No304
    ##  [3]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTL No306
    ##  [4]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQIX No0906S
    ##  [5]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTX No0908S
    ##  ...   ... ...
    ## [11]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTX No1007S
    ## [12]   321 XXXXXXXXXXXXXXXIDLPAPSN...LLPFLHTSKQRSLIFRPITQTX No1114S
    ## [13]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTX No1202S
    ## [14]   321 IRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTX No1206S
    ## [15]   321 XRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLIFRPITQTX No1208S

Translation using alternative genetic code:

As you can see from the above example, the initial amino acids `I` will change into `M` due to the mitochondrial translation code and also some `*` stop codons will change into a `W` amino acid.

    ## alternative genetic code
    ## use woodmouse data
    woodmouse |> dnabin2dnastring() |> cds2aa(shorten=TRUE,
        genetic.code=Biostrings::getGeneticCode("2"))

    ## AAStringSet object of length 15:
    ##      width seq                                              names
    ##  [1]   321 XRKTHPLLKXISHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTL No305
    ##  [2]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTL No304
    ##  [3]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTL No306
    ##  [4]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQMX No0906S
    ##  [5]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTX No0908S
    ##  ...   ... ...
    ## [11]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTX No1007S
    ## [12]   321 XXXXXXXXXXXXXXXIDLPAPSN...LLPFLHTSKQRSLMFRPITQTX No1114S
    ## [13]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTX No1202S
    ## [14]   321 MRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTX No1206S
    ## [15]   321 XRKTHPLLKIINHSFIDLPAPSN...LLPFLHTSKQRSLMFRPITQTX No1208S

# 6 Pairwise sequence comparison

## 6.1 Calculate pairwise AA distances (`aastring2dist()`)

### 6.1.1 Grantham’s distance

    ## calculate pairwise AA distances based on Grantham's distance
    aa.dist <- hiv |> cds2aa() |> aastring2dist(score=granthamMatrix())

    ##
    Computing: [========================================] 100% (done)

    ## obtain distances
    head(aa.dist$distSTRING)

    ##           U68496   U68497    U68498    U68499    U68500   U68501   U68502
    ## U68496  0.000000  4.43956 11.516484  9.879121 13.098901 17.05495 19.71429
    ## U68497  4.439560  0.00000 12.681319 10.406593 13.626374 18.21978 19.35165
    ## U68498 11.516484 12.68132  0.000000  8.769231  8.571429 13.25275 16.23077
    ## U68499  9.879121 10.40659  8.769231  0.000000  4.780220 15.58242 15.21978
    ## U68500 13.098901 13.62637  8.571429  4.780220  0.000000 15.38462 14.35165
    ## U68501 17.054945 18.21978 13.252747 15.582418 15.384615  0.00000 16.43956
    ##          U68503    U68504   U68505   U68506   U68507   U68508
    ## U68496 17.82418 13.857143 14.24176 14.92308 18.01099 19.01099
    ## U68497 19.80220 13.890110 14.27473 14.82418 18.54945 18.90110
    ## U68498 16.15385  9.516484 10.27473 11.04396 12.15385 15.56044
    ## U68499 16.58242 10.483516 12.47253 12.01099 14.24176 13.52747
    ## U68500 14.95604 10.571429 12.56044 12.06593 13.84615 10.89011
    ## U68501 17.50549 10.758242 13.08791 11.78022 14.91209 13.86813

    ## obtain pairwise sites used
    head(aa.dist$sitesUsed)

    ##        U68496 U68497 U68498 U68499 U68500 U68501 U68502 U68503 U68504 U68505
    ## U68496     91     91     91     91     91     91     91     91     91     91
    ## U68497     91     91     91     91     91     91     91     91     91     91
    ## U68498     91     91     91     91     91     91     91     91     91     91
    ## U68499     91     91     91     91     91     91     91     91     91     91
    ## U68500     91     91     91     91     91     91     91     91     91     91
    ## U68501     91     91     91     91     91     91     91     91     91     91
    ##        U68506 U68507 U68508
    ## U68496     91     91     91
    ## U68497     91     91     91
    ## U68498     91     91     91
    ## U68499     91     91     91
    ## U68500     91     91     91
    ## U68501     91     91     91

    ## create and plot bionj tree
    aa.dist.bionj <- ape::bionj(as.dist(aa.dist$distSTRING))
    plot(aa.dist.bionj)

To use a different score matrix, here as an example the `AAMatrix` from the R package `alakazam` is used:

    ## use AAMatrix data
    head(AAMatrix)

    ##   A B C D E F G H I J K L M N P Q R S T V W X Y Z * - .
    ## A 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0
    ## B 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 0 1 1 1 1 0 1 1 1 0 0
    ## C 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0
    ## D 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0
    ## E 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 0 0
    ## F 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0

    aa.dist.AAMatrix <- hiv |> cds2aa() |> aastring2dist(score=AAMatrix)

    ##
    Computing: [========================================] 100% (done)

    head(aa.dist.AAMatrix$distSTRING)

    ##            U68496     U68497    U68498     U68499     U68500    U68501
    ## U68496 0.00000000 0.06593407 0.1758242 0.17582418 0.23076923 0.2417582
    ## U68497 0.06593407 0.00000000 0.2197802 0.19780220 0.25274725 0.2857143
    ## U68498 0.17582418 0.21978022 0.0000000 0.13186813 0.16483516 0.1868132
    ## U68499 0.17582418 0.19780220 0.1318681 0.00000000 0.08791209 0.2307692
    ## U68500 0.23076923 0.25274725 0.1648352 0.08791209 0.00000000 0.2637363
    ## U68501 0.24175824 0.28571429 0.1868132 0.23076923 0.26373626 0.0000000
    ##           U68502    U68503    U68504    U68505    U68506    U68507    U68508
    ## U68496 0.2857143 0.2527473 0.2417582 0.2307692 0.2527473 0.2747253 0.3186813
    ## U68497 0.2967033 0.2747253 0.2527473 0.2417582 0.2637363 0.3076923 0.3186813
    ## U68498 0.2087912 0.2417582 0.1758242 0.1758242 0.2087912 0.1868132 0.2637363
    ## U68499 0.2197802 0.2417582 0.1758242 0.1758242 0.2087912 0.2087912 0.2197802
    ## U68500 0.2527473 0.2527473 0.2087912 0.2087912 0.2307692 0.2417582 0.2197802
    ## U68501 0.2527473 0.2637363 0.2087912 0.2307692 0.2197802 0.2417582 0.2747253

## 6.2 Calculate pairwise DNA distances (`dnastring2dist()`)

### 6.2.1 `ape::dist.dna` models

    ## use hiv data
    dna.dist <- hiv |> dnastring2dist(model="K80")

    ##
    Computing: [========================================] 100% (done)

    ## obtain distances
    head(dna.dist$distSTRING)

    ##            U68496     U68497     U68498     U68499     U68500    U68501
    ## U68496 0.00000000 0.03381189 0.07731910 0.08135801 0.11058044 0.1022214
    ## U68497 0.03381189 0.00000000 0.09396372 0.08960527 0.11916567 0.1194893
    ## U68498 0.07731910 0.09396372 0.00000000 0.05333488 0.07342797 0.0773191
    ## U68499 0.08135801 0.08960527 0.05333488 0.00000000 0.04143570 0.1022214
    ## U68500 0.11058044 0.11916567 0.07342797 0.04143570 0.00000000 0.1237391
    ## U68501 0.10222140 0.11948926 0.07731910 0.10222140 0.12373909 0.0000000
    ##           U68502    U68503     U68504     U68505    U68506     U68507    U68508
    ## U68496 0.1503488 0.1282034 0.11948926 0.11991759 0.1157063 0.13779195 0.1684152
    ## U68497 0.1594427 0.1280155 0.13725343 0.12865927 0.1381014 0.15620300 0.1682063
    ## U68498 0.1105804 0.1067010 0.09823725 0.09004207 0.1110678 0.08987206 0.1366493
    ## U68499 0.1150839 0.1157063 0.09808005 0.09412740 0.1108796 0.10254762 0.1232278
    ## U68500 0.1457404 0.1286593 0.12357099 0.11948926 0.1323187 0.12415474 0.1234296
    ## U68501 0.1150839 0.1241547 0.10670096 0.10711515 0.1108796 0.11570627 0.1368233

    ## obtain pairwise sites used
    head(dna.dist$sitesUsed)

    ##        U68496 U68497 U68498 U68499 U68500 U68501 U68502 U68503 U68504 U68505
    ## U68496    273    273    273    273    273    273    273    273    273    273
    ## U68497    273    273    273    273    273    273    273    273    273    273
    ## U68498    273    273    273    273    273    273    273    273    273    273
    ## U68499    273    273    273    273    273    273    273    273    273    273
    ## U68500    273    273    273    273    273    273    273    273    273    273
    ## U68501    273    273    273    273    273    273    273    273    273    273
    ##        U68506 U68507 U68508
    ## U68496    273    273    273
    ## U68497    273    273    273
    ## U68498    273    273    273
    ## U68499    273    273    273
    ## U68500    273    273    273
    ## U68501    273    273    273

    ## create and plot bionj tree
    dna.dist.bionj <- ape::bionj(as.dist(dna.dist$distSTRING))

It is also possible to compare the amino acid and nucleotide based trees:

    ## creation of the association matrix:
    association <- cbind(aa.dist.bionj$tip.label, aa.dist.bionj$tip.label)
    ## cophyloplot
    ape::cophyloplot(aa.dist.bionj,
        dna.dist.bionj,
        assoc=association,
        length.line=4,
        space=28,
        gap=3,
        rotate=FALSE)

### 6.2.2 `IUPAC` distance

    ## use hiv data
    hiv.dist.iupac <- head(hiv |> dnastring2dist(model="IUPAC"))

    ##
    Computing: [========================================] 100% (done)

    head(hiv.dist.iupac$distSTRING)

    ##            U68496     U68497     U68498     U68499     U68500     U68501
    ## U68496 0.00000000 0.03296703 0.07326007 0.07692308 0.10256410 0.09523810
    ## U68497 0.03296703 0.00000000 0.08791209 0.08424908 0.10989011 0.10989011
    ## U68498 0.07326007 0.08791209 0.00000000 0.05128205 0.06959707 0.07326007
    ## U68499 0.07692308 0.08424908 0.05128205 0.00000000 0.04029304 0.09523810
    ## U68500 0.10256410 0.10989011 0.06959707 0.04029304 0.00000000 0.11355311
    ## U68501 0.09523810 0.10989011 0.07326007 0.09523810 0.11355311 0.00000000
    ##           U68502    U68503     U68504     U68505    U68506     U68507    U68508
    ## U68496 0.1355311 0.1172161 0.10989011 0.10989011 0.1062271 0.12454212 0.1501832
    ## U68497 0.1428571 0.1172161 0.12454212 0.11721612 0.1245421 0.13919414 0.1501832
    ## U68498 0.1025641 0.0989011 0.09157509 0.08424908 0.1025641 0.08424908 0.1245421
    ## U68499 0.1062271 0.1062271 0.09157509 0.08791209 0.1025641 0.09523810 0.1135531
    ## U68500 0.1318681 0.1172161 0.11355311 0.10989011 0.1208791 0.11355311 0.1135531
    ## U68501 0.1062271 0.1135531 0.09890110 0.09890110 0.1025641 0.10622711 0.1245421

    ## run multi-threaded
    system.time(hiv |> dnastring2dist(model="IUPAC", threads=1))

    ##
    Computing: [========================================] 100% (done)

    ##    user  system elapsed
    ##       0       0       0

    system.time(hiv |> dnastring2dist(model="IUPAC", threads=2))

    ##
    Computing: [========================================] 100% (done)

    ##    user  system elapsed
    ##    0.02    0.00    0.02

Woodmouse data example:

    ## use woodmouse data
    woodmouse.dist <- woodmouse |> dnabin2dnastring() |> dnastring2dist()

    ##
    Computing: [========================================] 100% (done)

    head(woodmouse.dist$distSTRING)

    ##              No305       No304       No306     No0906S    No0908S    No0909S
    ## No305   0.00000000 0.016684046 0.013541667 0.018789144 0.01670146 0.01670146
    ## No304   0.01668405 0.000000000 0.005208333 0.013555787 0.01147028 0.01564129
    ## No306   0.01354167 0.005208333 0.000000000 0.009384776 0.00729927 0.01147028
    ## No0906S 0.01878914 0.013555787 0.009384776 0.000000000 0.01248699 0.01664932
    ## No0908S 0.01670146 0.011470282 0.007299270 0.012486993 0.00000000 0.01456816
    ## No0909S 0.01670146 0.015641293 0.011470282 0.016649324 0.01456816 0.00000000
    ##             No0910S     No0912S     No0913S     No1103S     No1007S    No1114S
    ## No305   0.017745303 0.014613779 0.018789144 0.012526096 0.016701461 0.01531729
    ## No304   0.012513034 0.013555787 0.005213764 0.011470282 0.015641293 0.01642935
    ## No306   0.008342023 0.009384776 0.005213764 0.007299270 0.011470282 0.01533406
    ## No0906S 0.009365245 0.014568158 0.012486993 0.012486993 0.016649324 0.02076503
    ## No0908S 0.011446410 0.012486993 0.012486993 0.010405827 0.014568158 0.02076503
    ## No0909S 0.015608741 0.010405827 0.016649324 0.008324662 0.002081165 0.02076503
    ##             No1202S     No1206S     No1208S
    ## No305   0.016701461 0.016701461 0.018828452
    ## No304   0.011470282 0.012513034 0.017782427
    ## No306   0.007299270 0.008342023 0.013598326
    ## No0906S 0.008324662 0.011446410 0.018789144
    ## No0908S 0.010405827 0.009365245 0.016701461
    ## No0909S 0.014568158 0.015608741 0.002087683

# 7 Coding sequences

## 7.1 Calculating synonymous and nonsynonymous substitutions (`dnastring2kaks()`)

    ## use hiv data
    ## model Li
    head(hiv |> dnastring2kaks(model="Li"))

    ## Joining, by = c("Comp1", "Comp2")
    ## Joining, by = c("Comp1", "Comp2")
    ## Joining, by = c("Comp1", "Comp2")

    ##    Comp1  Comp2         ka         ks          vka          vks
    ## 1 U68496 U68497 0.03026357 0.03170319 0.0003051202 0.0004007730
    ## 2 U68496 U68498 0.09777332 0.01761416 0.0009314173 0.0005091970
    ## 3 U68496 U68499 0.10295875 0.01767311 0.0009595527 0.0005787387
    ## 4 U68496 U68500 0.13461355 0.04639690 0.0013731885 0.0020497481
    ## 5 U68496 U68501 0.12607831 0.02844294 0.0013277282 0.0006310195
    ## 6 U68496 U68502 0.17441037 0.10926532 0.0017871934 0.0040766687

    ## model NG86
    head(hiv |> dnastring2kaks(model="NG86", threads=1))

    ##          Comp1 Comp2   seq1   seq2 Codons Compared Ambigiuous Indels Ns
    ## result.1     1     2 U68496 U68497     91       91          0      0  0
    ## result.2     1     3 U68496 U68498     91       91          0      0  0
    ## result.3     1     4 U68496 U68499     91       91          0      0  0
    ## result.4     1     5 U68496 U68500     91       91          0      0  0
    ## result.5     1     6 U68496 U68501     91       91          0      0  0
    ## result.6     1     7 U68496 U68502     91       91          0      0  0
    ##                        Sd               Sn                S                N
    ## result.1                3                6               57              216
    ## result.2              1.5             18.5             57.5            215.5
    ## result.3              1.5             19.5             56.5            216.5
    ## result.4              2.5             25.5 56.1666666666667 216.833333333333
    ## result.5              2.5             23.5 57.3333333333333 215.666666666667
    ## result.6 5.83333333333333 31.1666666666667             57.5            215.5
    ##                          ps                 pn             pn/ps
    ## result.1 0.0526315789473684 0.0277777777777778 0.527777777777778
    ## result.2 0.0260869565217391   0.08584686774942   3.2907965970611
    ## result.3 0.0265486725663717 0.0900692840646651  3.39260969976905
    ## result.4 0.0445103857566766  0.117601844734819  2.64212144504228
    ## result.5 0.0436046511627907  0.108964451313756  2.49891808346213
    ## result.6  0.101449275362319  0.144624903325599  1.42558833278091
    ##                          ds                 dn             dn/ds
    ## result.1 0.0545695157118212 0.0283052459871352 0.518700699793888
    ## result.2  0.026551445288187 0.0911703470876381  3.43372445823884
    ## result.3 0.0270299523623977 0.0959537646568416  3.54990505977829
    ## result.4 0.0458858673563109  0.127915513677956  2.78768869474935
    ## result.5 0.0449236061858017  0.117741219733924  2.62092093067847
    ## result.6  0.108999740568461  0.160668709585171  1.47402836692312

## 7.2 Codon comparison

As an example for the codon comparison data from the Human Immunodeficiency Virus Type 1 is used (Ganeshan et al. 1997), (Yang et al. 2000).

The window plots are constructed with the R package `ggplot2`.

### 7.2.1 Create codon matrix (`dnastring2codonmat()`)

    ## define two cds sequences
    cds1 <- Biostrings::DNAString("ATGCAACATTGC")
    cds2 <- Biostrings::DNAString("ATG---CATTGC")
    cds1.cds2.aln <- c(Biostrings::DNAStringSet(cds1),
        Biostrings::DNAStringSet(cds2))
    ## convert into alignment
    cds1.cds2.aln |> dnastring2codonmat()

    ##      [,1]  [,2]
    ## [1,] "ATG" "ATG"
    ## [2,] "CAA" "---"
    ## [3,] "CAT" "CAT"
    ## [4,] "TGC" "TGC"

Like the `cds2aa()` function, also the `dnastring2codonmat()` function is implemented to handle different frames.

    ## use frame 2 and shorten to circumvent multiple of three error
    cds1 <- Biostrings::DNAString("-ATGCAACATTGC-")
    cds2 <- Biostrings::DNAString("-ATG---CATTGC-")
    cds1.cds2.aln <- c(Biostrings::DNAStringSet(cds1),
        Biostrings::DNAStringSet(cds2))
    cds1.cds2.aln |> dnastring2codonmat(frame=2, shorten=TRUE)

    ##      [,1]  [,2]
    ## [1,] "ATG" "ATG"
    ## [2,] "CAA" "---"
    ## [3,] "CAT" "CAT"
    ## [4,] "TGC" "TGC"

### 7.2.2 Calculate average behavior of each codon (`codonmat2xy()`)

    ## use hiv data
    ## calculate average behavior
    hiv.xy <- hiv |> dnastring2codonmat() |> codonmat2xy()

    ## Joining, by = "Codon"
    ## Joining, by = "Codon"
    ## Joining, by = "Codon"

### 7.2.3 Plot average behavior of each codon

    print(hiv.xy |> dplyr::select(Codon,SynMean,NonSynMean,IndelMean) |>
        tidyr::gather(variable, values, -Codon) |>
        ggplot2::ggplot(aes(x=Codon, y=values)) +
        ggplot2::geom_line(aes(colour=factor(variable))) +
        ggplot2::geom_point(aes(colour=factor(variable))) +
        ggplot2::ggtitle("HIV-1 sample 136 patient 1 from
            Sweden envelope glycoprotein (env) gene"))

# 8 References

Chang, Peter L, Emily Kopania, Sara Keeble, Brice AJ Sarver, Erica Larson, Annie Orth, Khalid Belkhir, et al. 2017. “Whole Exome Sequencing of Wild-Derived Inbred Strains of Mice Improves Power to Link Phenotype and Genotype.” Mammalian Genome 28 (9): 416–25.

Ganeshan, Shanthi, Ruth E Dickover, BT Korber, Yvonne J Bryson, and Steven M Wolinsky. 1997. “Human Immunodeficiency Virus Type 1 Genetic Evolution in Children with Different Rates of Development of Disease.” Journal of Virology 71 (1): 663–77.

Grantham, Richard. 1974. “Amino Acid Difference Formula to Help Explain Protein Evolution.” Science 185 (4154): 862–64.

Harr, Bettina, Emre Karakoc, Rafik Neme, Meike Teschke, Christine Pfeifle, Željka Pezer, Hiba Babiker, et al. 2016. “Genomic Resources for Wild Populations of the House Mouse, Mus Musculus and Its Close Relative Mus Spretus.” Scientific Data 3 (1): 1–14.

Korneliussen, Thorfinn Sand, Anders Albrechtsen, and Rasmus Nielsen. 2014. “ANGSD: Analysis of Next Generation Sequencing Data.” BMC Bioinformatics 15 (1): 1–13.

Li, Wen-Hsiung. 1993. “Unbiased Estimation of the Rates of Synonymous and Nonsynonymous Substitution.” Journal of Molecular Evolution 36 (1): 96–99.

Nei, Masatoshi, and Takashi Gojobori. 1986. “Simple Methods for Estimating the Numbers of Synonymous and Nonsynonymous Nucleotide Substitutions.” Molecular Biology and Evolution 3 (5): 418–26.

Yang, Ziheng, Rasmus Nielsen, Nick Goldman, and Anne-Mette Krabbe Pedersen. 2000. “Codon-Substitution Models for Heterogeneous Selection Pressure at Amino Acid Sites.” Genetics 155 (1): 431–49.

# 9 Session Info

    sessionInfo()

    ## R version 4.2.0 RC (2022-04-19 r82224 ucrt)
    ## Platform: x86_64-w64-mingw32/x64 (64-bit)
    ## Running under: Windows Server x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] ggplot2_3.3.5       tidyr_1.2.0         dplyr_1.0.8
    ##  [4] ape_5.6-2           Biostrings_2.64.0   GenomeInfoDb_1.32.0
    ##  [7] XVector_0.36.0      IRanges_2.30.0      S4Vectors_0.34.0
    ## [10] BiocGenerics_0.42.0 MSA2dist_1.0.0      BiocStyle_2.24.0
    ##
    ## loaded via a namespace (and not attached):
    ##  [1] Rcpp_1.0.8.3           lattice_0.20-45        assertthat_0.2.1
    ##  [4] digest_0.6.29          foreach_1.5.2          utf8_1.2.2
    ##  [7] R6_2.5.1               evaluate_0.15          highr_0.9
    ## [10] pillar_1.7.0           zlibbioc_1.42.0        rlang_1.0.2
    ## [13] jquerylib_0.1.4        magick_2.7.3           rmarkdown_2.14
    ## [16] labeling_0.4.2         stringr_1.4.0          RCurl_1.98-1.6
    ## [19] munsell_0.5.0          compiler_4.2.0         xfun_0.30
    ## [22] pkgconfig_2.0.3        htmltools_0.5.2        tidyselect_1.1.2
    ## [25] tibble_3.1.6           GenomeInfoDbData_1.2.8 bookdown_0.26
    ## [28] codetools_0.2-18       fansi_1.0.3            withr_2.5.0
    ## [31] crayon_1.5.1           MASS_7.3-57            bitops_1.0-7
    ## [34] grid_4.2.0             nlme_3.1-157           jsonlite_1.8.0
    ## [37] gtable_0.3.0           lifecycle_1.0.1        DBI_1.1.2
    ## [40] magrittr_2.0.3         scales_1.2.0           cli_3.3.0
    ## [43] stringi_1.7.6          farver_2.1.0           doParallel_1.0.17
    ## [46] seqinr_4.2-8           bslib_0.3.1            ellipsis_0.3.2
    ## [49] generics_0.1.2         vctrs_0.4.1            iterators_1.0.14
    ## [52] tools_4.2.0            ade4_1.7-19            glue_1.6.2
    ## [55] purrr_0.3.4            parallel_4.2.0         fastmap_1.1.0
    ## [58] yaml_2.3.5             colorspace_2.0-3       BiocManager_1.30.17
    ## [61] GenomicRanges_1.48.0   knitr_1.38             sass_0.4.1