[med-svn] [r-cran-seqinr] 12/15: New upstream version 3.4-5
Andreas Tille
tille at debian.org
Sat Sep 23 06:55:23 UTC 2017
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-seqinr.
commit c7de922e1b94dc138d99eb4d47f3ceabe0beda8c
Author: Andreas Tille <tille at debian.org>
Date: Sat Sep 23 08:36:04 2017 +0200
New upstream version 3.4-5
---
DESCRIPTION | 21 ++--
MD5 | 86 ++++++++-------
NAMESPACE | 4 +-
NEWS | 1 +
R/dotchart.uco.R | 75 +++++++++++++
R/kaks.R | 12 ++-
R/read.alignment.R | 17 +--
R/uco.R | 78 --------------
data/AnoukResult.RData | Bin 319 -> 319 bytes
data/ECH.RData | Bin 38356 -> 38432 bytes
data/EXP.RData | Bin 1066 -> 1066 bytes
data/JLO.RData | Bin 39576 -> 39576 bytes
data/SEQINR.UTIL.RData | Bin 1775 -> 1775 bytes
data/aacost.RData | Bin 638 -> 638 bytes
data/aaindex.RData | Bin 81845 -> 81830 bytes
data/caitab.RData | Bin 941 -> 941 bytes
data/chargaff.RData | Bin 351 -> 351 bytes
data/clustal.RData | Bin 422 -> 422 bytes
data/datalist | 1 +
data/dinucl.RData | Bin 65740 -> 65740 bytes
data/ec999.RData | Bin 317767 -> 317766 bytes
data/fasta.RData | Bin 1028 -> 1029 bytes
data/gcO2.rda | Bin 6076 -> 6098 bytes
data/gcT.rda | Bin 17828 -> 17860 bytes
data/gs500liz.RData | Bin 176 -> 176 bytes
data/identifiler.RData | Bin 554 -> 554 bytes
data/kaksTorture.RData | Bin 0 -> 2220 bytes
data/m16j.RData | Bin 418892 -> 418884 bytes
data/mase.RData | Bin 543 -> 543 bytes
data/msf.RData | Bin 838 -> 838 bytes
data/pK.RData | Bin 322 -> 322 bytes
data/phylip.RData | Bin 276 -> 276 bytes
data/prochlo.RData | Bin 412468 -> 412468 bytes
data/revaligntest.RData | Bin 435 -> 435 bytes
data/sysdata.rda | Bin 1775 -> 1775 bytes
data/toyaa.RData | Bin 164 -> 164 bytes
data/toycodon.RData | Bin 195 -> 195 bytes
data/waterabs.RData | Bin 568 -> 568 bytes
inst/sequences/kaks-torture.fasta | 145 +++++++++++++++++++++++++
man/dotchart.uco.Rd | 10 +-
man/gcO2.Rd | 2 +-
man/kaks-torture.Rd | 36 +++++++
man/kaks.Rd | 46 +++++---
man/oriloc.Rd | 4 +-
src/getzlibsock.c | 3 +-
src/kaks.c | 222 ++++++++++++++++++++------------------
src/packagename_init.c | 43 ++++++++
47 files changed, 536 insertions(+), 270 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index 08633e5..c06e9c1 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,7 +1,7 @@
Encoding: UTF-8
Package: seqinr
-Version: 3.3-3
-Date: 2016-10-12
+Version: 3.4-5
+Date: 2017-08-01
Title: Biological Sequences Retrieval and Analysis
Authors at R: c(person("Delphine", "Charif", role = "aut"),
person("Olivier", "Clerc", role = "ctb"),
@@ -13,20 +13,21 @@ Authors at R: c(person("Delphine", "Charif", role = "aut"),
person("Guy", "Perrière", role = "ctb"))
Author: Delphine Charif [aut], Olivier Clerc [ctb], Carolin Frank [ctb], Jean R. Lobry [aut, cph], Anamaria Necşulea [ctb], Leonor Palmeira [ctb], Simon Penel [cre], Guy Perrière [ctb]
Maintainer: Simon Penel <simon.penel at univ-lyon1.fr>
-BugReports:
- http://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/seqinr-forum
+BugReports: http://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/seqinr-forum
Imports: ade4,segmented
Depends: R (>= 2.10.0)
Description: Exploratory data analysis and data visualization
- for biological sequence (DNA and protein) data. Includes also
- utilities for sequence data management under the ACNUC system.
+ for biological sequence (DNA and protein) data. Seqinr includes
+ utilities for sequence data management under the ACNUC system
+ described in Gouy, M. et al. (1984) Nucleic Acids Res.
+ 12:121-127 <doi:10.1093/nar/12.1Part1.121>.
License: GPL (>= 2)
SystemRequirements: zlib headers and library.
URL: http://seqinr.r-forge.r-project.org/
NeedsCompilation: yes
Repository: CRAN
Repository/R-Forge/Project: seqinr
-Repository/R-Forge/Revision: 2020
-Repository/R-Forge/DateTimeStamp: 2016-10-12 17:32:34
-Date/Publication: 2016-10-13 20:37:13
-Packaged: 2016-10-12 17:45:23 UTC; rforge
+Repository/R-Forge/Revision: 2078
+Repository/R-Forge/DateTimeStamp: 2017-08-01 09:25:54
+Date/Publication: 2017-08-01 13:18:36 UTC
+Packaged: 2017-08-01 10:05:09 UTC; rforge
diff --git a/MD5 b/MD5
index 0d7f920..a7a8ad8 100644
--- a/MD5
+++ b/MD5
@@ -1,5 +1,6 @@
-572ca39a8b8fae2be5f07aafed8c509d *DESCRIPTION
-f4eeeb197363233719e576484e515136 *NAMESPACE
+8779c4725ac84bff043854cea8005e2a *DESCRIPTION
+6376fe6e77b9c49dbb5d6d26434834ed *NAMESPACE
+d2ff218f55f61084bf070ee2b5341928 *NEWS
dd981a8bc19bea5e91c256f75ea8b5c9 *R/AAstat.R
7175d6708d717dca65176df83a41f9bd *R/ClassSeq.R
1e91dc3520dc9fd00e6e78da768c4b2e *R/GC.R
@@ -31,6 +32,7 @@ e354a7785d3f515d4311ddf75b54fc71 *R/countsubseqs.R
1394355049f1bb135f481a50508994d5 *R/dia.bactgensize.R
d47075fe840da6f000af83ece511af2c *R/dist.alignment.R
3dc6f329eea0d8db614bfb6ded99c056 *R/dotPlot.R
+35bf83b40232e876bb2e619749056162 *R/dotchart.uco.R
e7c8d78234ac7a398766943c2b179d91 *R/draw.oriloc.R
c5d3818fb836288a75c41ce6dd7c19e0 *R/draw.rearranged.oriloc.R
4f7b1e44bac0700d4d313876c1f19639 *R/draw.recstat.R
@@ -55,7 +57,7 @@ b582f5a822b09825f9f4dabc2e220cc4 *R/getType.R
e2af7c88e0e7c42cbb6e6208aa6364d8 *R/gfrag.R
1b9650fd96f5cbce29b1b79015ad99c5 *R/ghelp.R
ade42b1d2a6034de880e6676cf1f5da5 *R/isenum.R
-22d9950d657c07d38f245a908f615150 *R/kaks.R
+8834450f10c4d06cc047966749dda721 *R/kaks.R
106b5a8d04534e8441081bdadc23959f *R/knowndbs.R
3ea277d7f032f1d31004cda14e46fbab *R/lseqinr.R
6e4f1e677e5d870e8c5b3bf57fa25b91 *R/modifylist.R
@@ -74,7 +76,7 @@ ae47cf06593fc89cb8a105dfb69b6e93 *R/pmw.R
cb9f8a17e0bdafa1fe7b8a4584eb8ab0 *R/query.r
126d5ade572e95511b2b66b7189807d9 *R/quitacnuc.R
5c829cc5536f04c051d686cf1b553859 *R/read.abif.R
-568dca4f1a96d5619854a9e30ee60f58 *R/read.alignment.R
+33fa1c3e2f5f20e0f3c75792269e1e66 *R/read.alignment.R
bbed196d15ef9fb2b8d4da9d294ab5ed *R/read.fasta.R
8a46bbfddc96af2244eb21d7fa7157b6 *R/readBins.R
7fa1d45bb89c1229cd494c007ce94d5f *R/readPanels.R
@@ -99,42 +101,43 @@ b06cb8115a0f2f568cbc8946908a52ee *R/stutterabif.R
ab86b9a3f2568f102986b6fb93a545bc *R/test.li.recstat.R
a3350838cdb114e256e2441990c0dfa1 *R/translate.R
db37533bde4cb23943363ec1649c395b *R/trimSpace.R
-30f245251279c3e71f13f4e96ef30079 *R/uco.R
+8e180cdee361f0de8b4854c0428eb166 *R/uco.R
3c550b2b0e6b7f012c5884ec97a5f2b3 *R/util.R
5c6a726318b034002ec8478005d0ec54 *R/where.is.this.acc.R
4e4e196dab3926727e78e6af44616ebd *R/words.R
a56eb0ff10c2eb36431885fa78e0fdb8 *R/words.pos.R
204830bfeeef59878373d10b947e6fc2 *R/write.fasta.R
6daeda4870b681a9def5e87b0c2a692c *R/zscore.R
-fb5ea3868c5553ee8200b6ccaa8be14e *data/AnoukResult.RData
-89c880c7792e1fd0f247a675b00eac7e *data/ECH.RData
-9cc885d5878894f15e4180c920ad27bd *data/EXP.RData
-0b616d65d5b57d1386cf6fc4ebf9ee15 *data/JLO.RData
-a5062809357a2fe2539c70140d930245 *data/SEQINR.UTIL.RData
-1c88524ef8cd3758aa248f064deae5ce *data/aacost.RData
-30ce52aad5219c69e83336b93f7a162f *data/aaindex.RData
-f68ee04f58aadf44274a164adbe00821 *data/caitab.RData
-3af9fcdfaaf9d05c69a752710834a5c5 *data/chargaff.RData
-a83b249863be4cd4318d8e88922f7186 *data/clustal.RData
-806f1a3796b1e80cb87cb5b538b2586f *data/datalist
-6d345d4fefb0514cca31ace5fbc88730 *data/dinucl.RData
-39a1ebbe4ecb0e5acd1be77c70b2b700 *data/ec999.RData
-2813531b2d6fbd35e92467f338135660 *data/fasta.RData
-a68302621e3d01c9b556ad0181a9c70f *data/gcO2.rda
-76dd231427b0f6094387b1fd881f6cc7 *data/gcT.rda
-46a8298a82862a0a4b314f2769fa45c8 *data/gs500liz.RData
-3289605c76c50413d7ec66a0e7eae1bc *data/identifiler.RData
-50245edb26f8168405f878587e41f8c2 *data/m16j.RData
-ac915f82c50f06f8b16dc57cee87fd18 *data/mase.RData
-3d22106f5b1db1c0b0932367624d6124 *data/msf.RData
-e9b35c92761a0bba2f6ed13c3809cf90 *data/pK.RData
-b5aff4e0546952e316c8adeb40e7fffe *data/phylip.RData
-1828730dd85cef8c8c39ad05af6d6a01 *data/prochlo.RData
-d42f8673a6bbfde5ed82e67309303f3c *data/revaligntest.RData
-a5062809357a2fe2539c70140d930245 *data/sysdata.rda
-ed4eef2ae0c20040b249777ae83f36f7 *data/toyaa.RData
-7e5270892f1ac2b68a267f525a0f188b *data/toycodon.RData
-f8df825ad5220bae76d3a9378ecf99b9 *data/waterabs.RData
+033d80a8adda07cac603ce6b611d608f *data/AnoukResult.RData
+2b0a3f963e0280158c600e0a835c68f0 *data/ECH.RData
+1865b10477052ef2fb80a197be5cef8b *data/EXP.RData
+c9817ee41aba54c1210233bd7c545e9e *data/JLO.RData
+62e519367c3f3cf6be1d57af71bb1d31 *data/SEQINR.UTIL.RData
+aa077468c90bc3bc98c26b3fd261eda5 *data/aacost.RData
+2257b0c567b40276e9e210a5c8af8b3c *data/aaindex.RData
+8acad0cc4d855222639ae95acbcd6123 *data/caitab.RData
+060289fc3bc5f4625f08d3b957510f95 *data/chargaff.RData
+6cce89d032b04f46f559fd5bbb27af3d *data/clustal.RData
+831178c5bec48693c19f47db9f513546 *data/datalist
+a73c91d0f67b8597c69ab04a8c8b59df *data/dinucl.RData
+da07577d8b3ea0dc2cb2c8e0f368b74b *data/ec999.RData
+eba93070a5d7fb24a834da2218927e55 *data/fasta.RData
+a4a6a8e82355a655d37391ef30b42203 *data/gcO2.rda
+d8e3c74857a632fb70f87d494b808067 *data/gcT.rda
+47ba374c4260e17e6976b055a605ba37 *data/gs500liz.RData
+9266be8ec69dce388b8f44140f893db7 *data/identifiler.RData
+74030fda1c4203ba18b3969e2b106f55 *data/kaksTorture.RData
+3f255b9c54dbbf89586ce62076237ff8 *data/m16j.RData
+ea792211e964f9b87891dc5bf10a2687 *data/mase.RData
+ad734da5e5cb63a40d4b91559a0528cd *data/msf.RData
+8143331b64f077f8b75754ca2103eb94 *data/pK.RData
+140a09b00e84b2d8c3e07c63aa1d340e *data/phylip.RData
+cdbd5b6634836c95a8a103bb9e2fdc33 *data/prochlo.RData
+139282991e6521aa3106295a165829dc *data/revaligntest.RData
+62e519367c3f3cf6be1d57af71bb1d31 *data/sysdata.rda
+25253e4bac590ba3844a7967db019c32 *data/toyaa.RData
+74d31c280517b21cb57f291ea490eea4 *data/toycodon.RData
+bba3aacde26c6235c5725427af6cc482 *data/waterabs.RData
df797f2cdf697b330d68ac8194db7a53 *inst/CITATION
79179e6254bc1f46eb07303d85aa084e *inst/abif/1_0000206138_C01_005.fsa
45ad0f3fea5bf2143fc74a810a519945 *inst/abif/1_FAC321_0000205983_B02_004.fsa
@@ -171,6 +174,7 @@ cf31afbaea651a30e8f2b530483830bb *inst/sequences/hannah.txt
e03266d6ed9efa1dad9755a88e024565 *inst/sequences/humanMito.fasta
c56d8223a10995c70a6f965e46a06b44 *inst/sequences/input.dat
f82ee7968ac65ab690bdc1cfd10f055d *inst/sequences/input.out
+dc806adb55559a8486727292e56d58cf *inst/sequences/kaks-torture.fasta
1b714fae98378e5142d7ec658127ef36 *inst/sequences/legacy.fasta
d84b82c9cbc37133cb4d1f4e3c651878 *inst/sequences/louse.fasta
49138a950cfee1b2ebce2f95e6c4569c *inst/sequences/louse.names
@@ -229,7 +233,7 @@ d8c994f0a0bc546ae8456ad895772671 *man/countsubseqs.Rd
9463b95f58a2e0eae9c9d09f692fc607 *man/dinucl.Rd
e116cf59ab497b198e0cf53313493b0c *man/dist.alignment.Rd
fb1715453475d9bc554fa3803e9d1ed9 *man/dotPlot.Rd
-635d745b4d9d002b0136ed326ce302e4 *man/dotchart.uco.Rd
+83a53b12a5ef9a3717b8aeb2d9abc621 *man/dotchart.uco.Rd
7d235f8acb82c7c2cec55a7ee31b32eb *man/draw.oriloc.Rd
47ceb185e4d32c72952960cce8137126 *man/draw.rearranged.oriloc.Rd
cb50f090cefbd8f07526064d3760525e *man/draw.recstat.Rd
@@ -247,7 +251,7 @@ d2b9d3566289986cbded9e3ec42b349b *man/fasta.Rd
ebd1dd082d26b71029c3d5c478c34b47 *man/gb2fasta.Rd
dc73697f22c67c0762ceae6cc04f4f1f *man/gbk2g2.Rd
3013a472ea7f9eb0e83e58eccd8915b2 *man/gbk2g2.euk.Rd
-23ce44a155027f1adf77018b5000fb84 *man/gcO2.Rd
+ca71f60ee685df319cd6d10ac4050436 *man/gcO2.Rd
4683c2678cd9bac7e8da85772869fa18 *man/gcT.Rd
43e1ba4a463958d3fc30dd834cf84f51 *man/get.db.growth.Rd
97a1bc320a34ee6919ae9d5a537d8ed2 *man/get.ncbi.Rd
@@ -267,7 +271,8 @@ dedec3250e5f737066c28cbc819f9a75 *man/ghelp.Rd
79f3dc6846a328f749391c2333b840c5 *man/gs500liz.Rd
9e88c23480ad91dc30ad2e5f4a452932 *man/identifiler.Rd
193b78a2d66647ef32a26d3a6a9cb570 *man/isenum.Rd
-d1c3109d8cb096356e2a03d1e82c4693 *man/kaks.Rd
+d8dd4aad6b30a96ccb2a120382866b93 *man/kaks-torture.Rd
+a0cf3f443ed4063c266998eb9b209e99 *man/kaks.Rd
99cc2e82b216d68651eae627f741ed68 *man/knowndbs.Rd
a4394470a6facecb904714bfd79bbbec *man/lseqinr.Rd
e65af94d33f2668d6be2742b85657558 *man/m16j.Rd
@@ -276,7 +281,7 @@ e65af94d33f2668d6be2742b85657558 *man/m16j.Rd
4f32da53eef0d1cadb1547bf25d463cf *man/move.Rd
d28add4f994e4b9b07aa0d2322ec555f *man/msf.Rd
82c40f0f71a304c5778f2e5a271f18e5 *man/n2s.Rd
-35b13ee74557aa68fd8c6536b2b58f60 *man/oriloc.Rd
+b11b56d61b05e9b547540f283c6096f5 *man/oriloc.Rd
333457a72f8a13ed0558924f0adfa49d *man/pK.Rd
56afbbbf5ea360c7f53d43eba732aee8 *man/parser.socket.Rd
723e1c897ddc1cd527e4b68e72b61d16 *man/peakabif.Rd
@@ -337,7 +342,8 @@ cb407d8c24b063ecc44bb3bb49ab6216 *src/Makevars.win
436a2de8b172239b4ccf21ad72641759 *src/alignment.c
72d7cf9e818c492b578f390f794f10f7 *src/alignment.h
320b43b2d16f79ed0b12b471bad1c13c *src/fastacc.c
-ffbe5a8a0eb4f0468a358e69160f8c28 *src/getzlibsock.c
-f6eb92933aec8516b6fa19a6db75747c *src/kaks.c
+f001eff6e29fdf85842e50945fc1c0e6 *src/getzlibsock.c
+7d056d33ca14fe8ba3c7e805f37c2d00 *src/kaks.c
+fcee75c26a5363b9f427e319ab9e9b71 *src/packagename_init.c
3604f3611d34f9864f85ed904739cbe3 *src/util.c
833983b499aae68695d2852d4373abe4 *src/zsockr.c
diff --git a/NAMESPACE b/NAMESPACE
index 4fb4f9f..704a814 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,6 @@
# Export all names
-useDynLib(seqinr, .registration = TRUE)
-exportPattern(".")
+useDynLib(seqinr, .registration = TRUE,.fixes = "C_")
+exportPattern("^[^C_]")
importFrom("ade4", "dudi.coa")
importFrom("grDevices", "col2rgb", "grey", "rainbow", "rgb")
importFrom("graphics", "abline", "axis", "box", "dotchart", "hist",
diff --git a/NEWS b/NEWS
new file mode 100644
index 0000000..6846908
--- /dev/null
+++ b/NEWS
@@ -0,0 +1 @@
+Release notes are available at: http://seqinr.r-forge.r-project.org/src/appendix/releasenotes.pdf
\ No newline at end of file
diff --git a/R/dotchart.uco.R b/R/dotchart.uco.R
new file mode 100644
index 0000000..dbce359
--- /dev/null
+++ b/R/dotchart.uco.R
@@ -0,0 +1,75 @@
+dotchart.uco <- function(x, numcode = 1, aa3 = TRUE, pt.cex = 0.7,
+ alphabet = s2c("tcag"), pch = 21, gpch = 20, bg = par("bg"), cex = 0.7,
+ color = "black", gcolor = "black", lcolor = grey(0.9), xlim, ...)
+{
+ if( is.null(names(x)) ) names(x) <- words( alphabet = alphabet )
+ bcknames <- names(x)
+ x <- as.numeric(x)
+ names(x) <- bcknames
+#
+# General sorting
+#
+ x <- sort(x)
+ labels <- names(x)
+ stringlabel = paste(labels, sep = "", collapse = "")
+ groups <- as.factor(translate(s2c(stringlabel), numcode = numcode))
+ gdata <- sapply(split(x, groups), sum)
+#
+# Now, sorting by aa order
+#
+ gordered <- rank(gdata)
+ xidx <- numeric(64)
+
+ for( i in seq_len(64) )
+ {
+ xidx[i] <- -0.01*i + gordered[groups[i]]
+ }
+
+ x <- x[order(xidx)]
+ labels <- names(x)
+ stringlabel = paste(labels, sep = "", collapse = "")
+ aa <- translate(s2c(stringlabel), numcode = numcode)
+ groups <- factor(aa, levels = unique(aa))
+ gdata <- sapply(split(x, groups), sum)
+
+ if( missing(xlim) ) xlim <- c(0, max(gdata))
+ if( aa3 ) levels(groups) <- aaa(levels(groups))
+
+ dotchart(x = x, labels = labels, groups = groups, gdata = gdata,
+ pt.cex = pt.cex, pch = pch, gpch = gpch, bg = bg, color = color,
+ gcolor = gcolor, lcolor = lcolor, cex = cex, xlim, ...)
+#
+# Return invisibly for further plots
+#
+ result <- list(0)
+ result$x <- x
+ result$labels <- labels
+ result$groups <- groups
+ result$gdata <- gdata
+
+ ypg <- numeric( length(levels(groups)) )
+ i <- 1
+ for( aa in levels(groups) )
+ {
+ ypg[i] <- length(which(groups == aa)) + 2
+ i <- i + 1
+ }
+ ypg <- rev(cumsum(rev(ypg))) - 1
+ names(ypg) <- levels(groups)
+ result$ypg <- ypg
+
+ ypi <- numeric( length(x) )
+ for( i in seq_len(length(x)) )
+ {
+ ypi[i] <- ypg[groups[i]]
+ }
+ antirank <- function(x)
+ {
+ return( seq(length(x),1,by=-1 ))
+ }
+ ypi <- ypi - unlist(sapply(split(x, groups),antirank))
+ names(ypi) <- labels
+ result$ypi <- ypi
+
+ return( invisible(result) )
+}
diff --git a/R/kaks.R b/R/kaks.R
index 6038f7e..82bdbf0 100644
--- a/R/kaks.R
+++ b/R/kaks.R
@@ -1,4 +1,4 @@
-kaks <- function(x, verbose=FALSE, debug = FALSE, forceUpperCase = TRUE){
+kaks <- function(x, verbose = FALSE, debug = FALSE, forceUpperCase = TRUE, rmgap = TRUE){
#
# Check argument class:
#
@@ -41,9 +41,17 @@ kaks <- function(x, verbose=FALSE, debug = FALSE, forceUpperCase = TRUE){
}
}
#
+ # Choose option for gap removal
+ #
+ if(rmgap) {
+ gaprm = 0 # positions with at least one gap are removed
+ } else {
+ gaprm = 1 # only all gap positions are removed
+ }
+ #
# Call internal C function:
#
- l <- .Call("kaks", x$seq, x$nb, debug, PACKAGE = "seqinr")
+ l <- .Call("kaks", x$seq, x$nb, debug, gaprm, PACKAGE = "seqinr")
if(debug){
cat("<--- Result l storage is --->\n")
print(str(l))
diff --git a/R/read.alignment.R b/R/read.alignment.R
index b7c6b67..f092443 100644
--- a/R/read.alignment.R
+++ b/R/read.alignment.R
@@ -9,21 +9,22 @@ read.alignment <- function(file, format, forceToLower = TRUE)
file <- path.expand(file)
if(file.access(file, mode = 4) != 0) stop(paste("File", file, "is not readable"))
- ali <- switch( format,
- fasta = .Call("read_fasta_align", file, PACKAGE = "seqinr"),
- FASTA = .Call("read_fasta_align", file, PACKAGE = "seqinr"),
+ fasta2ali <- function(file){
+ tmp <- read.fasta(file, as.string = TRUE)
+ list(length(tmp), getName(tmp), unlist(getSequence(tmp, as.string = TRUE)))
+ }
+ ali <- switch( tolower(format),
+ fasta = fasta2ali(file),
mase = .Call("read_mase", file, PACKAGE = "seqinr"),
- MASE = .Call("read_mase", file, PACKAGE = "seqinr"),
phylip = .Call("read_phylip_align", file, PACKAGE = "seqinr"),
- PHYLIP = .Call("read_phylip_align", file, PACKAGE = "seqinr"),
msf = .Call("read_msf_align", file, PACKAGE = "seqinr"),
- MSF = .Call("read_msf_align", file, PACKAGE = "seqinr"),
- CLUSTAL = .Call("read_clustal_align", file, PACKAGE = "seqinr"),
clustal = .Call("read_clustal_align", file, PACKAGE = "seqinr"),
- stop("Wrong format name: Format available are fasta,mase,phylip,msf,clustal")
+ stop("Wrong format name: Format available are fasta, mase, phylip, msf, clustal")
)
ali <- lapply(ali, as.character)
+ #cleaning \r char
+ ali <- lapply(ali, function (x ){gsub ('\r','',x)})
if(forceToLower) ali[[3]] <- lapply(ali[[3]], tolower)
if(format == "mase"){
ali <- list(nb = as.numeric(ali[[1]]), nam = ali[[2]], seq = ali[[3]], com = ali[[4]])
diff --git a/R/uco.R b/R/uco.R
index 4eca45e..02a8cfe 100755
--- a/R/uco.R
+++ b/R/uco.R
@@ -41,81 +41,3 @@ as.data.frame = FALSE, NA.rscu = NA)
}
}
-
-dotchart.uco <- function(x, numcode = 1, aa3 = TRUE, cex = 0.7,
- alphabet = s2c("tcag"), pch = 21, gpch = 20, bg = par("bg"),
- color = par("fg"), gcolor = par("fg"), lcolor = "gray", xlim, ...)
-{
- if( is.null(names(x)) ) names(x) <- words( alphabet = alphabet )
-#
-# General sorting
-#
- x <- sort(x)
- labels <- names(x)
- stringlabel = paste(labels, sep="", collapse="")
- groups <- as.factor(translate(s2c(stringlabel), numcode = numcode))
- gdata <- sapply(split(x, groups), sum)
-#
-# Now, sorting by aa order
-#
- gordered <- rank(gdata)
- xidx <- numeric(64)
-
- for( i in seq_len(64) )
- {
- xidx[i] <- -0.01*i + gordered[groups[i]]
- }
-
- x <- x[order(xidx)]
- labels <- names(x)
- stringlabel = paste(labels, sep="", collapse="")
- aa <- translate(s2c(stringlabel), numcode = numcode)
- groups <- factor(aa, levels = unique(aa))
- gdata <- sapply(split(x, groups), sum)
-
- if( missing(xlim) ) xlim <- c(0, max(gdata))
- if( aa3 )
- {
- levels(groups) <- aaa(levels(groups))
- }
- dotchart(x = x, labels = labels, groups = groups, gdata = gdata,
- cex = cex, pch = pch, gpch = gpch, bg = bg, color = color,
- gcolor = gcolor, lcolor = lcolor, xlim, ...)
-#
-# Return invisibly for further plots
-#
- result <- list(0)
- result$x <- x
- result$labels <- labels
- result$groups <- groups
- result$gdata <- gdata
-
- ypg <- numeric( length(levels(groups)) )
- i <- 1
- for( aa in levels(groups) )
- {
- ypg[i] <- length(which(groups == aa)) + 2
- i <- i + 1
- }
- ypg <- rev(cumsum(rev(ypg))) - 1
- names(ypg) <- levels(groups)
- result$ypg <- ypg
-
- ypi <- numeric( length(x) )
- for( i in seq_len(length(x)) )
- {
- ypi[i] <- ypg[groups[i]]
- }
- antirank <- function(x)
- {
- return( seq(length(x),1,by=-1 ))
- }
- ypi <- ypi - unlist(sapply(split(x, groups),antirank))
- names(ypi) <- labels
- result$ypi <- ypi
-
- return( invisible(result) )
-}
-
-
-
diff --git a/data/AnoukResult.RData b/data/AnoukResult.RData
index 7842380..dbe4e1b 100644
Binary files a/data/AnoukResult.RData and b/data/AnoukResult.RData differ
diff --git a/data/ECH.RData b/data/ECH.RData
index b0cbd74..c9036eb 100644
Binary files a/data/ECH.RData and b/data/ECH.RData differ
diff --git a/data/EXP.RData b/data/EXP.RData
index 859d313..9cfb52a 100644
Binary files a/data/EXP.RData and b/data/EXP.RData differ
diff --git a/data/JLO.RData b/data/JLO.RData
index 452be8e..9e739e5 100644
Binary files a/data/JLO.RData and b/data/JLO.RData differ
diff --git a/data/SEQINR.UTIL.RData b/data/SEQINR.UTIL.RData
index df81124..3af358e 100644
Binary files a/data/SEQINR.UTIL.RData and b/data/SEQINR.UTIL.RData differ
diff --git a/data/aacost.RData b/data/aacost.RData
index 3e14ea2..d4dd180 100644
Binary files a/data/aacost.RData and b/data/aacost.RData differ
diff --git a/data/aaindex.RData b/data/aaindex.RData
index 4239ca6..89f6b10 100644
Binary files a/data/aaindex.RData and b/data/aaindex.RData differ
diff --git a/data/caitab.RData b/data/caitab.RData
index 101f9e8..cd84e97 100644
Binary files a/data/caitab.RData and b/data/caitab.RData differ
diff --git a/data/chargaff.RData b/data/chargaff.RData
index de8a70c..07c4605 100644
Binary files a/data/chargaff.RData and b/data/chargaff.RData differ
diff --git a/data/clustal.RData b/data/clustal.RData
index 8d65c4f..2a7d11d 100644
Binary files a/data/clustal.RData and b/data/clustal.RData differ
diff --git a/data/datalist b/data/datalist
index 6d53430..d71c85f 100644
--- a/data/datalist
+++ b/data/datalist
@@ -13,6 +13,7 @@ ec999
fasta
gs500liz
identifiler
+kaksTorture
m16j
mase
msf
diff --git a/data/dinucl.RData b/data/dinucl.RData
index 7d4b56c..4fef760 100644
Binary files a/data/dinucl.RData and b/data/dinucl.RData differ
diff --git a/data/ec999.RData b/data/ec999.RData
index f6fbb8e..abf8bec 100644
Binary files a/data/ec999.RData and b/data/ec999.RData differ
diff --git a/data/fasta.RData b/data/fasta.RData
index 0624d1c..78d8033 100644
Binary files a/data/fasta.RData and b/data/fasta.RData differ
diff --git a/data/gcO2.rda b/data/gcO2.rda
index 4d5a7c3..7d59626 100644
Binary files a/data/gcO2.rda and b/data/gcO2.rda differ
diff --git a/data/gcT.rda b/data/gcT.rda
index a8255ba..f8927d7 100644
Binary files a/data/gcT.rda and b/data/gcT.rda differ
diff --git a/data/gs500liz.RData b/data/gs500liz.RData
index 9b49439..d6ac880 100644
Binary files a/data/gs500liz.RData and b/data/gs500liz.RData differ
diff --git a/data/identifiler.RData b/data/identifiler.RData
index 054fd9c..5d20d0c 100644
Binary files a/data/identifiler.RData and b/data/identifiler.RData differ
diff --git a/data/kaksTorture.RData b/data/kaksTorture.RData
new file mode 100644
index 0000000..33f85f0
Binary files /dev/null and b/data/kaksTorture.RData differ
diff --git a/data/m16j.RData b/data/m16j.RData
index 6de552d..eb62692 100644
Binary files a/data/m16j.RData and b/data/m16j.RData differ
diff --git a/data/mase.RData b/data/mase.RData
index 63291fd..a6af9ad 100644
Binary files a/data/mase.RData and b/data/mase.RData differ
diff --git a/data/msf.RData b/data/msf.RData
index 4b2771e..feefbbc 100644
Binary files a/data/msf.RData and b/data/msf.RData differ
diff --git a/data/pK.RData b/data/pK.RData
index b2841b7..b60eac6 100644
Binary files a/data/pK.RData and b/data/pK.RData differ
diff --git a/data/phylip.RData b/data/phylip.RData
index 5fcca0d..fc64bd7 100644
Binary files a/data/phylip.RData and b/data/phylip.RData differ
diff --git a/data/prochlo.RData b/data/prochlo.RData
index 2b8d9c6..447acf3 100644
Binary files a/data/prochlo.RData and b/data/prochlo.RData differ
diff --git a/data/revaligntest.RData b/data/revaligntest.RData
index 24850e2..e3ad060 100644
Binary files a/data/revaligntest.RData and b/data/revaligntest.RData differ
diff --git a/data/sysdata.rda b/data/sysdata.rda
index df81124..3af358e 100644
Binary files a/data/sysdata.rda and b/data/sysdata.rda differ
diff --git a/data/toyaa.RData b/data/toyaa.RData
index b808c37..2ee2f6b 100644
Binary files a/data/toyaa.RData and b/data/toyaa.RData differ
diff --git a/data/toycodon.RData b/data/toycodon.RData
index 90cc976..b364136 100644
Binary files a/data/toycodon.RData and b/data/toycodon.RData differ
diff --git a/data/waterabs.RData b/data/waterabs.RData
index 06f2f97..34c1c87 100644
Binary files a/data/waterabs.RData and b/data/waterabs.RData differ
diff --git a/inst/sequences/kaks-torture.fasta b/inst/sequences/kaks-torture.fasta
new file mode 100644
index 0000000..0383754
--- /dev/null
+++ b/inst/sequences/kaks-torture.fasta
@@ -0,0 +1,145 @@
+;
+; This is a test file for the kaks() function from seqinr package.
+; This is an alignment of the 64 possible codons so that kaks
+; values are computed for possible pairs of one-codon sequences.
+;
+; Only finite values should be returned. Expected values are given
+; in data(kaks-torture).
+;
+; This file was generated with the following R script :
+;
+; x <- words() # The 64 codons
+; aa1 <- sapply(x, function(x) translate(s2c(x)))
+; aa3 <- aaa(aa1)
+; names <- paste(x, aa1, aa3, sep = "-")
+; sequences <- as.list(toupper(x))
+; write.fasta(sequences, names, file.out = "kaks-torture.fasta")
+;
+>aaa-K-Lys
+AAA
+>aac-N-Asn
+AAC
+>aag-K-Lys
+AAG
+>aat-N-Asn
+AAT
+>aca-T-Thr
+ACA
+>acc-T-Thr
+ACC
+>acg-T-Thr
+ACG
+>act-T-Thr
+ACT
+>aga-R-Arg
+AGA
+>agc-S-Ser
+AGC
+>agg-R-Arg
+AGG
+>agt-S-Ser
+AGT
+>ata-I-Ile
+ATA
+>atc-I-Ile
+ATC
+>atg-M-Met
+ATG
+>att-I-Ile
+ATT
+>caa-Q-Gln
+CAA
+>cac-H-His
+CAC
+>cag-Q-Gln
+CAG
+>cat-H-His
+CAT
+>cca-P-Pro
+CCA
+>ccc-P-Pro
+CCC
+>ccg-P-Pro
+CCG
+>cct-P-Pro
+CCT
+>cga-R-Arg
+CGA
+>cgc-R-Arg
+CGC
+>cgg-R-Arg
+CGG
+>cgt-R-Arg
+CGT
+>cta-L-Leu
+CTA
+>ctc-L-Leu
+CTC
+>ctg-L-Leu
+CTG
+>ctt-L-Leu
+CTT
+>gaa-E-Glu
+GAA
+>gac-D-Asp
+GAC
+>gag-E-Glu
+GAG
+>gat-D-Asp
+GAT
+>gca-A-Ala
+GCA
+>gcc-A-Ala
+GCC
+>gcg-A-Ala
+GCG
+>gct-A-Ala
+GCT
+>gga-G-Gly
+GGA
+>ggc-G-Gly
+GGC
+>ggg-G-Gly
+GGG
+>ggt-G-Gly
+GGT
+>gta-V-Val
+GTA
+>gtc-V-Val
+GTC
+>gtg-V-Val
+GTG
+>gtt-V-Val
+GTT
+>taa-*-Stp
+TAA
+>tac-Y-Tyr
+TAC
+>tag-*-Stp
+TAG
+>tat-Y-Tyr
+TAT
+>tca-S-Ser
+TCA
+>tcc-S-Ser
+TCC
+>tcg-S-Ser
+TCG
+>tct-S-Ser
+TCT
+>tga-*-Stp
+TGA
+>tgc-C-Cys
+TGC
+>tgg-W-Trp
+TGG
+>tgt-C-Cys
+TGT
+>tta-L-Leu
+TTA
+>ttc-F-Phe
+TTC
+>ttg-L-Leu
+TTG
+>ttt-F-Phe
+TTT
diff --git a/man/dotchart.uco.Rd b/man/dotchart.uco.Rd
index 64c7ad8..cdf9ef0 100644
--- a/man/dotchart.uco.Rd
+++ b/man/dotchart.uco.Rd
@@ -5,20 +5,22 @@
Draw a Cleveland dot plot for codon usage tables
}
\usage{
-dotchart.uco(x, numcode = 1, aa3 = TRUE, cex = 0.7, alphabet = s2c("tcag"),
- pch = 21, gpch = 20, bg = par("bg"), color = par("fg"), gcolor = par("fg"),
-lcolor = "gray", xlim, ...)
+dotchart.uco(x, numcode = 1, aa3 = TRUE, pt.cex = 0.7, alphabet =
+ s2c("tcag"), pch = 21, gpch = 20, bg = par("bg"), cex
+ = 0.7, color = "black", gcolor = "black", lcolor =
+ grey(0.9), xlim, ...)
}
\arguments{
\item{x}{table of codon usage as computed by \code{uco}. }
\item{numcode}{the number of the code to be used by \code{translate}.}
\item{aa3}{logical. If TRUE use the three-letter code for amino-
acids. If FALSE use the one-letter code for amino-acids. }
- \item{cex}{the character size to be used. }
+ \item{pt.cex}{the character size to be used for points. }
\item{alphabet}{character for codons labels}
\item{pch}{the plotting character or symbol to be used.}
\item{gpch}{the plotting character or symbol to be used for group values. }
\item{bg}{the background color to be used. }
+ \item{cex}{the character expansion size passed to \code{\link{dotchart}}. }
\item{color}{the color(s) to be used for points an labels. }
\item{gcolor}{the single color to be used for group labels and values.}
\item{lcolor}{the color(s) to be used for the horizontal lines.}
diff --git a/man/gcO2.Rd b/man/gcO2.Rd
index 9bb073f..fd727bc 100644
--- a/man/gcO2.Rd
+++ b/man/gcO2.Rd
@@ -11,7 +11,7 @@
\source{
Naya, H., Romero, H., Zavala, A., Alvarez, B. and Musto, H. (2002) Aerobiosis increases the Genomic Guanine Plus Cytosine Content (GC%) in prokaryotes. \emph{Journal of Molecular Evolution}, \bold{55}:260-264.
-Data imported into seqinr by J.R. Lobry on 09-OCT-2016. Original source location given in the article was \code{http://oeg.fcien.edu.uy/GCprok/} but is no more active. Data were copied at \url{http://pbil.univ-lyon1.fr/R/donnees/gcO2.txt} (\emph{cf.} section 2.1 in Lobry, J.R (2004) Life history traits and genome structure: aerobiosis and G+C content in bacteria. \emph{Lecture Notes in Computer Sciences}, \bold{3039}:679-686). Import was from this last ressource.
+Data imported into seqinr by J.R. Lobry on 09-OCT-2016. Original source location given in the article was \code{http://oeg.fcien.edu.uy/GCprok/} but is no more active. Data were copied at \url{http://pbil.univ-lyon1.fr/R/donnees/gcO2.txt} (\emph{cf.} section 2.1 in Lobry, J.R (2004) Life history traits and genome structure: aerobiosis and G+C content in bacteria. \emph{Lecture Notes in Computer Sciences}, \bold{3039}:679-686). Import was from this last ressource. There are 130 aerobic ge [...]
}
\references{
\code{citation("seqinr")}
diff --git a/man/kaks-torture.Rd b/man/kaks-torture.Rd
new file mode 100644
index 0000000..67b8014
--- /dev/null
+++ b/man/kaks-torture.Rd
@@ -0,0 +1,36 @@
+\name{kaksTorture}
+\alias{kaksTorture}
+\docType{data}
+\title{Expected numeric results for Ka and Ks in extreme cases}
+\description{
+This data set is what should be obtained when runing \code{kaks()}
+on the test file kaks-torture.fasta in the sequences directory of the
+seqinR package.
+}
+\usage{data(kaksTorture)}
+\format{
+ A list with 4 components of class dist.
+ \describe{
+ \item{ka}{Ka}
+ \item{ks}{Ks}
+ \item{vka}{variance for Ka}
+ \item{vks}{variance for Ks}
+ }
+}
+\source{
+See comments in kaks-torture.fasta for R code used to produce it.
+}
+\references{
+\code{citation("seqinr")}
+}
+\examples{
+data(kaksTorture)
+kaks.torture <- read.alignment(file = system.file("sequences/kaks-torture.fasta",
+ package = "seqinr"), format = "fasta")
+#
+# Failed on windows :
+#
+# stopifnot(identical(kaksTorture, kaks(kaks.torture)))
+# stopifnot(identical(kaksTorture, kaks(kaks.torture, rmgap = FALSE)))
+}
+\keyword{datasets}
diff --git a/man/kaks.Rd b/man/kaks.Rd
index 47cccaa..5522d80 100644
--- a/man/kaks.Rd
+++ b/man/kaks.Rd
@@ -4,15 +4,16 @@
\description{ Ks and Ka are, respectively, the number of substitutions per synonymous site and per non-synonymous site between two protein-coding genes. They are also denoted as ds and dn in the literature. The ratio of nonsynonymous (Ka) to synonymous (Ks) nucleotide substitution rates is an indicator of selective pressures on genes. A ratio significantly greater than 1 indicates positive selective pressure. A ratio around 1 indicates either neutral evolution at the protein level or an [...]
}
\usage{
-kaks(x,verbose = FALSE, debug = FALSE, forceUpperCase = TRUE)
+kaks(x, verbose = FALSE, debug = FALSE, forceUpperCase = TRUE, rmgap = TRUE)
}
\arguments{
\item{x}{ An object of class \code{alignment}, obtained for instance by importing into R the data from an alignment file with the \code{\link{read.alignment}} function. This is typically a set of coding sequences aligned at the protein level, see \code{\link{reverse.align}}.}
- \item{verbose}{ If TRUE add to the results the value of L0, L2, L4 (respectively the number of non-synonymous sites, of 2-fold synonymous sites, of 4-fold synonymous sites), A0, A2, A4 (respectively the number of transitional changes at non-synonymous, 2-fold, and 4-fold synonymous sites ) and B0, B2, B4 (respectively the number of transversional changes at non-synonymous, 2-fold, and 4-fold synonymous sites).}
+ \item{verbose}{ If TRUE add to the results the value of L0, L2, L4 (respectively the frequency of non-synonymous sites, of 2-fold synonymous sites, of 4-fold synonymous sites), A0, A2, A4 (respectively the number of transitional changes at non-synonymous, 2-fold, and 4-fold synonymous sites ) and B0, B2, B4 (respectively the number of transversional changes at non-synonymous, 2-fold, and 4-fold synonymous sites).}
\item{debug}{ If TRUE turns debug mode on.}
\item{forceUpperCase}{ If TRUE, the default value, all character in sequences are forced to the upper case
if at least one 'a', 'c', 'g', or 't' is found in the sequences.
Turning it to FALSE if the sequences are already in upper case will save time.}
+ \item{rmgap}{ If TRUE all positions with at least one gap are removed. If FALSE only positions with nothing else than gaps are removed.}
}
\value{
\item{ ks }{ matrix of Ks values }
@@ -21,16 +22,19 @@ kaks(x,verbose = FALSE, debug = FALSE, forceUpperCase = TRUE)
\item{ vka }{ variance matrix of Ka }
}
\references{
-Li, W.-H. (1993) Unbiased estimation of the rates of synonymous and nonsynonymous substitution.
-\emph{J. Mol. Evol.}, \bold{36}:96-99.\cr
+Li, W.-H., Wu, C.-I., Luo, C.-C. (1985) A new method for estimating synonymous and nonsynonymous rates of nucleotide substitution considering the relative likelihood of nucleotide and codon changes. \emph{Mol. Biol. Evol}, \bold{2}:150-174\cr
+
+Li, W.-H. (1993) Unbiased estimation of the rates of synonymous and nonsynonymous substitution. \emph{J. Mol. Evol.}, \bold{36}:96-99.\cr
+
+Pamilo, P., Bianchi, N.O. (1993) Evolution of the \emph{Zfx} and \emph{Zfy} genes: Rates and interdependence between genes. \emph{Mol. Biol. Evol}, \bold{10}:271-281\cr
Hurst, L.D. (2002) The Ka/Ks ratio: diagnosing the form of sequence evolution.
\emph{Trends Genet.}, \bold{18}:486-486.\cr
The C programm implementing this method was provided by Manolo Gouy. More info is
needed here to trace back the original C source so as to credit correct source.
-The original FORTRAN-77 code by Chung-I Wu modified by Ken Wolfe is available
-here \url{http://wolfe.gen.tcd.ie/lab/pub/li93/}.\cr
+The original FORTRAN-77 code by Chung-I Wu modified by Ken Wolfe was available
+here \url{http://wolfe.gen.tcd.ie/lab/pub/li93/} but this is no more true as 2017-07-01.\cr
For a more recent discussion about the estimation of Ka and Ks see:\cr
@@ -50,20 +54,21 @@ As from seqinR 2.0-3, when there is at least one non ACGT base in a codon, this
Gap-codons (\code{---}) are not used for computations.
-When the alignment does not contain enough information (\emph{i.e.} close to saturation), the Ka and Ks values are forced to 10.
+When the alignment does not contain enough information (\emph{i.e.} close to saturation), the Ka and Ks values are forced to 10 (more exactly to 9.999999).
Negative values indicate that Ka and Ks can not be computed.
-According to Li (1993) J. Mol. Evol. 36(1):96-99,
+According to Li (1993) and Pamilo and Bianchi (1993),
the rate of synonymous substitutions Ks is computed as:
Ks = (L2.A2 + L4.A4) / (L2 + L4) + B4
and the rate of non-synonymous substitutions Ka is computed as:
Ka = A0 + (L0.B0 + L2.B2) / (L0 + L2)
- }
+}
\author{D. Charif, J.R. Lobry}
-\seealso{\code{\link{read.alignment}} to import alignments from files, \code{\link{reverse.align}} to align CDS at the aa level.}
+\seealso{\code{\link{read.alignment}} to import alignments from files, \code{\link{reverse.align}} to align CDS at the aa level,
+\code{\link{kaksTorture}} for test on one-codon CDS.}
\examples{
#
# Simple Toy example:
@@ -77,11 +82,11 @@ Ka = A0 + (L0.B0 + L2.B2) / (L0 + L2)
data(AnoukResult)
Anouk <- read.alignment(file = system.file("sequences/Anouk.fasta", package = "seqinr"),
format = "fasta")
-## if( ! all.equal(kaks(Anouk), AnoukResult) ) {
-## warning("Poor numeric results with respect to AnoukResult standard")
-## } else {
-## print("Results are consistent with AnoukResult standard")
-## }
+ if( ! all.equal(kaks(Anouk), AnoukResult) ) {
+ warning("Poor numeric results with respect to AnoukResult standard")
+ } else {
+ print("Results are consistent with AnoukResult standard")
+ }
#
# As from seqinR 2.0-3 the following alignment with out-of-frame gaps
# should return a zero Ka value.
@@ -98,4 +103,15 @@ Ka = A0 + (L0.B0 + L2.B2) / (L0 + L2)
Darren <- read.alignment(file = system.file("sequences/DarrenObbard.fasta", package = "seqinr"),
format = "fasta")
stopifnot( all.equal(kaks(Darren)$ka[1], 0) )
+#
+# As from seqinR 3.4-0, non-finite values should never be returned for
+# Ka and Ks even for small sequences. The following test checks that this
+# is true for an alignement of the 64 codons, so that we compute Ka and
+# Ks for all possible pairs of codons.
+#
+wrd <- as.alignment(nb = 64, nam = words(), seq = words())
+res <- kaks(wrd)
+if(any(!is.finite(res$ka))) stop("Non finite value returned for Ka")
+if(any(!is.finite(res$ks))) stop("Non finite value returned for Ks")
+
}
diff --git a/man/oriloc.Rd b/man/oriloc.Rd
index 9be6030..8112480 100644
--- a/man/oriloc.Rd
+++ b/man/oriloc.Rd
@@ -76,12 +76,12 @@ The original paper for oriloc:\cr
Frank, A.C., Lobry, J.R. (2000) Oriloc: prediction of replication
boundaries in unannotated bacterial chromosomes. \emph{Bioinformatics},
\bold{16}:566-567.\cr
-\url{http://bioinformatics.oupjournals.org/cgi/reprint/16/6/560}\cr\cr
+\url{https://doi.org/10.1093/bioinformatics/16.6.560}\cr\cr
A simple informal introduction to DNA-walks:\cr
Lobry, J.R. (1999) Genomic landscapes. \emph{Microbiology Today},
\bold{26}:164-165.\cr
-\url{http://www.microbiologysociety.org/download.cfm/docid/C66F92D0-6292-4CB0-83E344560770384C}\cr\cr
+\url{http://seqinr.r-forge.r-project.org/MicrTod_1999_26_164.pdf}\cr\cr
An early and somewhat historical application of DNA-walks:\cr
Lobry, J.R. (1996) A simple vectorial representation of DNA sequences
diff --git a/src/getzlibsock.c b/src/getzlibsock.c
index 2acc09b..92e8362 100644
--- a/src/getzlibsock.c
+++ b/src/getzlibsock.c
@@ -36,10 +36,9 @@ SEXP getzlibsock(SEXP sock, SEXP nmax, SEXP debug)
int numsoc;
- int i,j, n, nn, nnn, ok, warn, nread, c;
+ int i, n, nn, nnn, nread;
int itest,itestd;
- char *test;
char *res;
int flagend =0;
diff --git a/src/kaks.c b/src/kaks.c
index 61c1cdc..def988b 100644
--- a/src/kaks.c
+++ b/src/kaks.c
@@ -8,13 +8,13 @@ void prefastlwl(double **, double **, double **, double **, double **, double **
int fastlwl(char **, int, int, double **, double **, double **, double **, double **, double **, double **, double **,
double **, double **, double **, double **, double **,double **, double **,double **,double **, double **,double **,double **, double **,double **);
-SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks)
+SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks, SEXP gaprm)
{
char **seqIn; /* local working copy of sequences */
char **seq; /* pointer to original sequences from R object */
double *tl0[64], *tl1[64], *tl2[64], *tti0[64], *tti1[64], *tti2[64], *ttv0[64], *ttv1[64], *ttv2[64];
int i, j, totseqs, lgseq, n;
- int debugon;
+ int debugon, option;
double *rl[21];
double **ka, **ks, **vka, **vks;
double **l0, **l2,**l4;
@@ -75,6 +75,7 @@ SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks)
debugon = INTEGER_VALUE(debugkaks);
totseqs = INTEGER_VALUE(nbseq);
+ option = INTEGER_VALUE(gaprm);
if(debugon) Rprintf("C> mode degug is on at C level with %d sequences\n", totseqs);
@@ -196,7 +197,7 @@ SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks)
PROTECT(ra0 = NEW_NUMERIC(totseqs*totseqs));
PROTECT(ra2 = NEW_NUMERIC(totseqs*totseqs));
PROTECT(ra4 = NEW_NUMERIC(totseqs*totseqs));
- PROTECT(rb0 = NEW_NUMERIC(totseqs*totseqs));
+ PROTECT(rb0 = NEW_NUMERIC(totseqs*totseqs));
PROTECT(rb2 = NEW_NUMERIC(totseqs*totseqs));
PROTECT(rb4 = NEW_NUMERIC(totseqs*totseqs));
@@ -258,7 +259,7 @@ SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks)
/* */
/******************************************************************************/
- reresh(seqIn,totseqs,0);
+ reresh(seqIn, totseqs, option); /* seqIn est modifié par reresh */
for(i = 0 ; i < totseqs ; i++){
if(debugon) Rprintf("reresh-->%s<--\n", seqIn[i]);
@@ -380,38 +381,36 @@ SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks)
SET_ELEMENT(res, 10, rb0);
SET_ELEMENT(res, 11, rb2);
SET_ELEMENT(res, 12, rb4);
+
+ if(debugon) Rprintf("C> %s", "End of C level....................\n");
+
UNPROTECT(14);
- if(debugon) Rprintf("C> %s", "End of C level....................\n");
return(res);
}
int num(char *cod)
{
- int n1, n2, n3;
-
- n1 = n2 = n3 = 0;
- if (cod[0] == 'C')
- n1 = 1;
- if (cod[1] == 'C')
- n2 = 1;
- if (cod[2] == 'C')
- n3 = 1;
- if (cod[0] == 'G')
- n1 = 2;
- if (cod[1] == 'G')
- n2 = 2;
- if (cod[2] == 'G')
- n3 = 2;
- if (cod[0] == 'T')
- n1 = 3;
- if (cod[1] == 'T')
- n2 = 3;
- if (cod[2] == 'T')
- n3 = 3;
-
- return 16 * n1 + 4 * n2 + n3;
+ int n1, n2, n3;
+//MG
+ static const char bases[] = "ACGT";
+ if(strchr(bases, cod[0]) == NULL ||
+ strchr(bases, cod[1]) == NULL ||
+ strchr(bases, cod[2]) == NULL) return 64;
+//MG
+ n1 = n2 = n3 = 0;
+ if (cod[0] == 'C') n1 = 1;
+ if (cod[1] == 'C') n2 = 1;
+ if (cod[2] == 'C') n3 = 1;
+ if (cod[0] == 'G') n1 = 2;
+ if (cod[1] == 'G') n2 = 2;
+ if (cod[2] == 'G') n3 = 2;
+ if (cod[0] == 'T') n1 = 3;
+ if (cod[1] == 'T') n2 = 3;
+ if (cod[2] == 'T') n3 = 3;
+
+ return 16 * n1 + 4 * n2 + n3;
}
@@ -425,7 +424,7 @@ int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
aaa[3], bb[3], flgseq, va[3], vb[3];
char cod1[3], cod2[3];
int i, j, ii, num1, num2, sat, sat1, sat2;
- sat = sat1 = sat2 = 2;
+ sat = sat1 = sat2 = 2;
/*
Internal check at C level: this should be no more be necessary, I'll keep it just in case. JRL - 26-APR-2009
*/
@@ -434,19 +433,21 @@ int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
REprintf("Fatal error: the number of nucleotide after gap removal is not a multiple of 3.\nPlease report this bug on the seqinr diffusion list.\n");
return(0); /* Should be R's NA but an int is returned by fastlwl */
}
+
for (i = 0; i < nbseq - 1; i++) {
for (j = i + 1; j < nbseq; j++) {
l[0] = l[1] = l[2] = 0;
ti[0] = ti[1] = ti[2] = tv[0] = tv[1] = tv[2] = 0;
for (ii = 0; ii < lgseq / 3; ii++) {
- cod1[0] = *(seq[i] + 3 * ii);
- cod1[1] = *(seq[i] + 3 * ii + 1);
+ cod1[0] = *(seq[i] + 3 * ii);
+ cod1[1] = *(seq[i] + 3 * ii + 1);
cod1[2] = *(seq[i] + 3 * ii + 2);
cod2[0] = *(seq[j] + 3 * ii);
cod2[1] = *(seq[j] + 3 * ii + 1);
cod2[2] = *(seq[j] + 3 * ii + 2);
num1 = num(cod1);
num2 = num(cod2);
+ if(num1 == 64 || num2 == 64) continue;//MG ignore - or N-containing codons
l[0] += *(tl0[num1] + num2);
l[1] += *(tl1[num1] + num2);
l[2] += *(tl2[num1] + num2);
@@ -457,7 +458,7 @@ int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
tv[1] += *(ttv1[num1] + num2);
tv[2] += *(ttv2[num1] + num2);
}
- l0[i][j]=l[0];
+ l0[i][j]=l[0];
l2[i][j]=l[1];
l4[i][j]=l[2];
for (ii = 0; ii < 3; ii++) {
@@ -466,13 +467,15 @@ int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
aaa[ii] = 1 / (1 - 2 * p[ii] - q[ii]);
bb[ii] = 1 / (1 - 2 * q[ii]);
cc[ii] = (aaa[ii] + bb[ii]) / 2;
- if (bb[ii] <= 0) {
- b[ii] = 10;
+ /* adding the isfinite condition - JLO JUL 2017 */
+ if (bb[ii] <= 0 || !isfinite(bb[ii])) {
+ b[ii] = 10.0;
} else {
b[ii] = 0.5 * (double) log(bb[ii]);
}
- if ((aaa[ii] <= 0) || (bb[ii] <= 0)) {
- a[ii] = 10;
+ /* adding the isfinite condition - JLO JUL 2017 */
+ if ((aaa[ii] <= 0) || (bb[ii] <= 0) || !isfinite(aaa[ii]) || !isfinite(bb[ii])) {
+ a[ii] = 10.0;
} else {
a[ii] = 0.5 * (double) log(aaa[ii]) - 0.25 * log(bb[ii]);
}
@@ -480,24 +483,24 @@ int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
vb[ii] = bb[ii] * bb[ii] * q[ii] * (1 - q[ii]) / l[ii];
}
- if ((a[1] != 10) && (a[2] != 10) && (b[2] != 10)){
+ if ((a[1] < 10) && (a[2] < 10) && (b[2] < 10)){
ks[i][j] = (l[1] * a[1] + l[2] * a[2]) / (l[2] + l[1]) + b[2];
vks[i][j] = (l[1] * l[1] * va[1] + l[2] * l[2] * va[2]) / ((l[1] + l[2]) * (l[1]+l[2])) + vb[2] - bb[2] * q[2] * (2 * aaa[2] * p[2] - cc[2] * (1 - q[2]))/(l[1]+l[2]);
} else {
sat1 = 1;
vks[i][j]=ks[i][j] = 9.999999;
}
- if ((a[0] != 10) && (b[0] != 10) && (b[1] != 10)){
+ if ((a[0] < 10) && (b[0] < 10) && (b[1] < 10)){
ka[i][j] = a[0] + (l[0] * b[0] + l[1] * b[1]) / (l[0] + l[1]);
vka[i][j] = (l[0] * l[0] * vb[0] + l[1] * l[1] * vb[1]) / ((l[1] + l[0]) * (l[1]+l[0])) + va[0] - bb[0] * q[0] * (2 * aaa[0] * p[0] - cc[0] * (1 - q[0]))/(l[1]+l[0]);
} else {
vka[i][j]=ka[i][j] = 9.999999;
sat2 = 1;
}
- a0[i][j]=a[0];
+ a0[i][j]=a[0];
a2[i][j]=a[1];
- a4[i][j]=a[2];
- b0[i][j]=b[0];
+ a4[i][j]=a[2];
+ b0[i][j]=b[0];
b2[i][j]=b[1];
b4[i][j]=b[2];
@@ -511,10 +514,10 @@ int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
B0, B2, B4: # of transversional changes at non-synonymous, 2-fold, and 4-fold synonymous sites
- Ces quantit�s sont les suivantes dans la fonction fastlwl():
- L0, L2, l4 correspondent � l[0], l[1], l[2]
- A0, A2, A4 correspondent � a[0], a[1], a[2]
- B0, B2, B4 correspondent � b[0], b[1], b[2]
+ Ces quantités sont les suivantes dans la fonction fastlwl():
+ L0, L2, l4 correspondent à l[0], l[1], l[2]
+ A0, A2, A4 correspondent à a[0], a[1], a[2]
+ B0, B2, B4 correspondent à b[0], b[1], b[2]
*/
@@ -761,44 +764,48 @@ code_mt: des TI syno en site 2-fold qui ont ete comptees normalement
void titv2(char *cod1, char *cod2, double *ti, double *tv, double* l, int *aa, double **rl, int* pos)
{
- char codint1[3], codint2[3];
- int i, j, n = 0;
- double l1, l2, p1, p2;
+ char codint1[4], codint2[4];
+ int i, j, n, aa1, aa2, aaint1, aaint2;
+ double l1, l2, p1, p2;
void titv1(char *, char *, double, double *, double *,double*);
- memcpy(codint1, cod1, 3);
- memcpy(codint2, cod1, 3);
+ memcpy(codint1, cod1, 3);
+ memcpy(codint2, cod1, 3); /* codint_2_ <-- cod_1_ : no problem */
for (i = 0; i < 2; i++) {
- if (cod1[i] != cod2[i])
+ if (cod1[i] != cod2[i]){
codint1[i] = cod2[i];
- if (cod1[i] != cod2[i])
break;
+ }
}
for (j = i + 1; j <= 2; j++) {
- if (cod1[j] != cod2[j])
+ if (cod1[j] != cod2[j]){
codint2[j] = cod2[j];
- if (cod1[j] != cod2[j])
break;
+ }
}
-
- l1 = *(rl[aa[num(cod1)]] + aa[num(codint1)]) * *(rl[aa[num(codint1)]] + aa[num(cod2)]);
- l2 = *(rl[aa[num(cod1)]] + aa[num(codint2)]) * *(rl[aa[num(codint2)]] + aa[num(cod2)]);
- p1 = l1 / (l1 + l2);
- p2 = 1 - p1;
+ aa1=aa[num(cod1)]; aa2=aa[num(cod2)];
+ aaint1=aa[num(codint1)]; aaint2=aa[num(codint2)];
+
+ l1 = *(rl[aa1] + aaint1) * *(rl[aaint1] + aa2);
+ l2 = *(rl[aa1] + aaint2) * *(rl[aaint2] + aa2);
+ p1 = (l1+l2)? l1 / (l1 + l2) : 0.;
+ p2 = (l1+l2)? 1.-p1 : 0.;
for (i=0;i<3;i++) if (pos[i]==0) n=i+1;
l[catsite(cod1[0], cod1[1] ,cod1[2], n)]+=0.333333;
l[catsite(cod2[0], cod2[1] ,cod2[2], n)]+=0.333333;
l[catsite(codint1[0], codint1[1] ,codint1[2], n)]+=0.333333*p1;
l[catsite(codint2[0], codint2[1] ,codint2[2], n)]+=0.333333*p2;
+
+
+
titv1(cod1, codint1, p1, ti, tv,l);
titv1(cod2, codint1, p1, ti, tv,l);
titv1(cod1, codint2, p2, ti, tv,l);
titv1(cod2, codint2, p2, ti, tv,l);
-
}
void titv3(char *cod1, char *cod2, double *ti, double *tv, double* l, int *aa, double **rl)
@@ -1068,58 +1075,61 @@ return;
void reresh(char** seq, int nbseq, int option){
-/* Si option = 0, toutes les positions avec au moins un gap sont eliminees */
-
+/* Si option = 0, toutes les positions avec au moins un gap sont éliminées.
+ Sinon, seules les positions avec uniquement des gaps sont éliminées */
int lgseq, l, drapeau, i, j, k;
char **seqref;
- seqref = (char **) R_alloc(nbseq, sizeof(char *));
-
- lgseq = strlen(seq[1]);
+/* Allocation dynamique du tableau seqref de l'alignement */
- for(i = 0 ; i < nbseq ; i++){
- seqref[i] = (char*) R_alloc(lgseq + 1, sizeof(char));
- }
+ seqref = (char **) R_alloc(nbseq, sizeof(char *));
+ lgseq = strlen(seq[1]);
+ for(i = 0 ; i < nbseq ; i++){
+ seqref[i] = (char*) R_alloc(lgseq + 1, sizeof(char));
+ }
+
+ l = -1; /* position de la colonne courante dans seqref */
+ if (option == 0){
+ for(i = 0 ; i < lgseq ; i++){
+ drapeau = 0; /* 0 si pas de gap */
+ for(j = 0 ; j < nbseq; j++){
+ if (*(seq[j] + i) == '-') drapeau = 1;
+ }
+ if (drapeau == 0){ /* on recopie la colonne i de seq dans la colonne l de seqref */
+ l++;
+ for(k = 0 ; k < nbseq ; k++) *(seqref[k] + l) = *(seq[k] + i);
+ }
+ }
+ }
+ else{
+ for(i = 0 ; i < lgseq ; i++){
+ drapeau = 0; /* 1 au premier non gap */
+ for(j = 0 ; j < nbseq ; j++){
+ if (*(seq[j] + i) != '-') {
+ drapeau = 1;
+ break;
+ }
+ }
+ if (drapeau == 1){ /* on recopie la colonne i de seq dans la colonne l de seqref */
+ l++;
+ for(k = 0 ; k < nbseq ; k++) *(seqref[k] + l) = *(seq[k] + i);
+ }
+ }
+ }
+/* Ajout de caractères nuls en fin d'alignement dans seqref */
+ for(i = 0 ; i < nbseq ; i++){
+ for(j = l + 1 ; j < lgseq ; j++) {
+ *(seqref[i] + j) = '\0';
+ }
+ }
- l=-1;
- if (option==0){
- for(i=0;i<lgseq;i++){
- drapeau=0;
- for(j=0;j<nbseq;j++){
- if (*(seq[j]+i)=='-') drapeau=1;
- }
- if (drapeau==0){
- l++;
- for(k=0;k<nbseq;k++) *(seqref[k]+l)=*(seq[k]+i);
- }
- }
- }
- else{
- for(i=0;i<lgseq;i++){
- drapeau=0;
- for(j=0;j<nbseq;j++){
- if (*(seq[j]+i)!='-') {
- drapeau=1;
- break;
- }
- }
- if (drapeau==1){
- l++;
- for(k=0;k<nbseq;k++) *(seqref[k]+l)=*(seq[k]+i);
- }
- }
- }
- for(i=0;i<nbseq;i++){
- for (j=l+1;j<lgseq;j++) {
- *(seqref[i]+j)='\0';
- }
- }
- for (i=0;i<nbseq;i++) {
- for (j=0;j<lgseq;j++){
- *(seq[i]+j)=*(seqref[i]+j);
- }
- }
+/* Recopie de seqref dans seq */
+ for(i = 0 ; i < nbseq ; i++) {
+ for(j = 0 ; j < lgseq ; j++){
+ *(seq[i] + j) = *(seqref[i] + j);
+ }
+ }
}
diff --git a/src/packagename_init.c b/src/packagename_init.c
new file mode 100644
index 0000000..33e1f38
--- /dev/null
+++ b/src/packagename_init.c
@@ -0,0 +1,43 @@
+#include <R.h>
+#include <Rinternals.h>
+#include <stdlib.h> // for NULL
+#include <R_ext/Rdynload.h>
+
+/* FIXME:
+ Check these declarations against the C/Fortran source code.
+*/
+
+/* .Call calls */
+extern SEXP distance(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP fastacc(SEXP, SEXP, SEXP, SEXP, SEXP);
+extern SEXP getzlibsock(SEXP, SEXP, SEXP);
+extern SEXP is_a_protein_seq(SEXP);
+extern SEXP kaks(SEXP, SEXP, SEXP);
+extern SEXP read_clustal_align(SEXP);
+extern SEXP read_fasta_align(SEXP);
+extern SEXP read_mase(SEXP);
+extern SEXP read_msf_align(SEXP);
+extern SEXP read_phylip_align(SEXP);
+extern SEXP s2c(SEXP);
+
+static const R_CallMethodDef CallEntries[] = {
+ {"distance", (DL_FUNC) &distance, 5},
+ {"fastacc", (DL_FUNC) &fastacc, 5},
+ {"getzlibsock", (DL_FUNC) &getzlibsock, 3},
+ {"is_a_protein_seq", (DL_FUNC) &is_a_protein_seq, 1},
+ {"kaks", (DL_FUNC) &kaks, 4},
+ {"read_clustal_align", (DL_FUNC) &read_clustal_align, 1},
+ {"read_fasta_align", (DL_FUNC) &read_fasta_align, 1},
+ {"read_mase", (DL_FUNC) &read_mase, 1},
+ {"read_msf_align", (DL_FUNC) &read_msf_align, 1},
+ {"read_phylip_align", (DL_FUNC) &read_phylip_align, 1},
+ {"s2c", (DL_FUNC) &s2c, 1},
+ {NULL, NULL, 0}
+};
+
+void R_init_seqinr(DllInfo *dll)
+{
+ R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+ R_useDynamicSymbols(dll, FALSE);
+}
+
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-seqinr.git
More information about the debian-med-commit
mailing list