[Seqinr-commits] r1767 - in pkg: R man src
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Feb 17 15:29:26 CET 2012
Author: simonpenel
Date: 2012-02-17 15:29:26 +0100 (Fri, 17 Feb 2012)
New Revision: 1767
Modified:
pkg/R/kaks.R
pkg/man/kaks.Rd
pkg/src/kaks.c
Log:
Modification of kaks: option verbose display values of L0,L2,L4,A0,A2,A4,B0,B2,B4
Modified: pkg/R/kaks.R
===================================================================
--- pkg/R/kaks.R 2012-02-15 09:30:35 UTC (rev 1766)
+++ pkg/R/kaks.R 2012-02-17 14:29:26 UTC (rev 1767)
@@ -1,4 +1,4 @@
-kaks <- function(x, debug = FALSE, forceUpperCase = TRUE){
+kaks <- function(x, verbose=FALSE, debug = FALSE, forceUpperCase = TRUE){
#
# Check argument class:
#
@@ -47,6 +47,7 @@
if(debug){
cat("<--- Result l storage is --->\n")
print(str(l))
+ print (l)
cat("<--- Result l storage is --->\n")
}
#
@@ -61,8 +62,20 @@
tmp <- matrix( k, x$nb, x$nb, byrow = TRUE, dimnames = list(x$nam, x$nam))
as.dist(t(tmp))
}
- result <- lapply(l[seq_len(4)], mkresult)
- names(result) <- c("ka", "ks", "vka", "vks")
- return(result)
+ #result <- lapply(l[seq_len(4)], mkresult)
+ if (verbose)
+ {
+ result <- lapply(l[seq_len(13)], mkresult)
+ check <- result[[5]] + result[[6]] + result[[7]]
+ result[[14]] <-check
+ names(result) <- c("ka", "ks", "vka", "vks","l0","l2","l4","a0","a2","a4", "b0","b2","b4","checksuml")
+ return(result)
+ }
+ else
+ {
+ result <- lapply(l[seq_len(4)], mkresult)
+ names(result) <- c("ka", "ks", "vka", "vks")
+ return(result)
+ }
}
Modified: pkg/man/kaks.Rd
===================================================================
--- pkg/man/kaks.Rd 2012-02-15 09:30:35 UTC (rev 1766)
+++ pkg/man/kaks.Rd 2012-02-17 14:29:26 UTC (rev 1767)
@@ -4,14 +4,15 @@
\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 averaging of sites under positive and negative selective pressures. A ratio less than 1 indicates pressures to conserve protein sequence (\emph{i.e.} purifying selection). This function estimates the Ka and Ks values for a set of aligned sequences using the method published by Li (1993) and gives the associated variance matrix.
}
\usage{
-kaks(x, debug = FALSE, forceUpperCase = TRUE)
+kaks(x,verbose = FALSE, debug = FALSE, forceUpperCase = 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{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{debug}{ If TRUE turns debug mode on.}
}
\value{
\item{ ks }{ matrix of Ks values }
@@ -53,6 +54,13 @@
Negative values indicate that Ka and Ks can not be computed.
+According to Li (1993) J. Mol. Evol. 36(1):96-99,
+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.}
@@ -67,11 +75,11 @@
#
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.
Modified: pkg/src/kaks.c
===================================================================
--- pkg/src/kaks.c 2012-02-15 09:30:35 UTC (rev 1766)
+++ pkg/src/kaks.c 2012-02-17 14:29:26 UTC (rev 1767)
@@ -6,7 +6,7 @@
void reresh(char **, int, int);
void prefastlwl(double **, double **, double **, double **, 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 **, double **,double **,double **, double **,double **);
SEXP kaks(SEXP sequences, SEXP nbseq, SEXP debugkaks)
{
@@ -17,8 +17,17 @@
int debugon;
double *rl[21];
double **ka, **ks, **vka, **vks;
+ double **l0, **l2,**l4;
+ double **a0, **a2,**a4;
+ double **b0, **b2,**b4;
+
+ double *xl0,*xl2,*xl4;
+ double *xa0,*xa2,*xa4;
+ double *xb0,*xb2,*xb4;
+
double *xka, *xks, *xvka, *xvks;
+
double mat[19][19] = {{.382, .382, .343, .382, .382, .382, .382, .128, .040, .128, .040, .128, .040, .128, .040, .128, .343, .128, .040 },
{ .382, .382, .128, .343, .343, .343, .343, .128, .040, .128, .040, .128, .040, .128, .040, .128, .128, .040, .040 },
{ .343, .128, .343, .382, .382, .382, .343, .128, .040, .128, .128, .343, .128, .343, .128, .343, .343, .128, .040 },
@@ -44,6 +53,24 @@
SEXP rvka;
SEXP rvks;
SEXP res;
+
+ /* -- addition fevrier 2012 --*/
+ SEXP rl0;
+ SEXP rl2;
+ SEXP rl4;
+
+ SEXP ra0;
+ SEXP ra2;
+ SEXP ra4;
+
+ SEXP rb0;
+ SEXP rb2;
+ SEXP rb4;
+
+
+
+ /* --------------------------- */
+
/* SEXP lsequtil; The effective number of sites used, not used yet */
debugon = INTEGER_VALUE(debugkaks);
@@ -102,11 +129,36 @@
vka = (double **) R_alloc(totseqs, sizeof(double *));
vks = (double **) R_alloc(totseqs, sizeof(double *));
+ l0 = (double **) R_alloc(totseqs, sizeof(double *));
+ l2 = (double **) R_alloc(totseqs, sizeof(double *));
+ l4 = (double **) R_alloc(totseqs, sizeof(double *));
+
+ a0 = (double **) R_alloc(totseqs, sizeof(double *));
+ a2 = (double **) R_alloc(totseqs, sizeof(double *));
+ a4 = (double **) R_alloc(totseqs, sizeof(double *));
+
+
+ b0 = (double **) R_alloc(totseqs, sizeof(double *));
+ b2 = (double **) R_alloc(totseqs, sizeof(double *));
+ b4 = (double **) R_alloc(totseqs, sizeof(double *));
+
+
for (i = 0; i < totseqs; i++) {
ka[i] = (double *) R_alloc(totseqs, sizeof(double));
vka[i] = (double *) R_alloc(totseqs, sizeof(double));
ks[i] = (double *) R_alloc(totseqs, sizeof(double));
vks[i] = (double *) R_alloc(totseqs, sizeof(double));
+ l0[i] = (double *) R_alloc(totseqs, sizeof(double));
+ l2[i] = (double *) R_alloc(totseqs, sizeof(double));
+ l4[i] = (double *) R_alloc(totseqs, sizeof(double));
+ a0[i] = (double *) R_alloc(totseqs, sizeof(double));
+ a2[i] = (double *) R_alloc(totseqs, sizeof(double));
+ a4[i] = (double *) R_alloc(totseqs, sizeof(double));
+ b0[i] = (double *) R_alloc(totseqs, sizeof(double));
+ b2[i] = (double *) R_alloc(totseqs, sizeof(double));
+ b4[i] = (double *) R_alloc(totseqs, sizeof(double));
+
+
}
/******************************************************************************/
@@ -130,14 +182,33 @@
/******************************************************************************/
- PROTECT(res = NEW_LIST(5));
+ PROTECT(res = NEW_LIST(14));
PROTECT(rka = NEW_NUMERIC(totseqs*totseqs));
PROTECT(rks = NEW_NUMERIC(totseqs*totseqs));
PROTECT(rvka = NEW_NUMERIC(totseqs*totseqs));
PROTECT(rvks = NEW_NUMERIC(totseqs*totseqs));
-
-
- for (i = 2; i < 21; i++) {
+
+ /* -- addition fevrier 2012 --*/
+
+ PROTECT(rl0 = NEW_NUMERIC(totseqs*totseqs));
+ PROTECT(rl2 = NEW_NUMERIC(totseqs*totseqs));
+ PROTECT(rl4 = NEW_NUMERIC(totseqs*totseqs));
+ 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(rb2 = NEW_NUMERIC(totseqs*totseqs));
+ PROTECT(rb4 = NEW_NUMERIC(totseqs*totseqs));
+
+
+
+
+/* PROTECT(rl024 = NEW_NUMERIC(3));
+ PROTECT(ra024 = NEW_NUMERIC(3));
+ PROTECT(rb024 = NEW_NUMERIC(3));*/
+ /* --------------------------- */
+
+ for (i = 2; i < 21; i++) {
for (j = 1; j < i; j++) {
*(rl[i] + j) = mat[j-1][i-2] ;
}
@@ -198,6 +269,18 @@
ks[i][j] = -1;
vka[i][j] = -1;
vks[i][j] = -1;
+ l0[i][j] = 0;
+ l2[i][j] = 0;
+ l4[i][j] = 0;
+ a0[i][j] = 0;
+ a2[i][j] = 0;
+ a4[i][j] = 0;
+ b0[i][j] = 0;
+ b2[i][j] = 0;
+ b4[i][j] = 0;
+
+
+
}
}
@@ -231,11 +314,14 @@
}
*/
- lgseq = strlen(seqIn[0]);
- fastlwl(seqIn, totseqs, lgseq, ka, ks, tti0, tti1, tti2, ttv0, ttv1, ttv2, tl0, tl1, tl2, vka, vks);
+ lgseq = strlen(seqIn[0]);
+ /* l024 = NUMERIC_POINTER(rl024);
+ a024 = NUMERIC_POINTER(ra024);
+ b024 = NUMERIC_POINTER(rb024);*/
+ fastlwl(seqIn, totseqs, lgseq, ka, ks, tti0, tti1, tti2, ttv0, ttv1, ttv2, tl0, tl1, tl2, vka, vks,l0,l2,l4,a0,a2,a4,b0,b2,b4);
for(i = 0 ; i < totseqs ; i++){
- if(debugon) Rprintf("fastlwl-->%s<--\n", seqIn[i]);
+ if(debugon) Rprintf(" -->%s<--\n", seqIn[i]);
}
/******************************************************************************/
@@ -250,25 +336,52 @@
xks = NUMERIC_POINTER(rks);
xvka = NUMERIC_POINTER(rvka);
xvks = NUMERIC_POINTER(rvks);
-
+ xl0 = NUMERIC_POINTER(rl0);
+ xl2 = NUMERIC_POINTER(rl2);
+ xl4 = NUMERIC_POINTER(rl4);
+ xa0 = NUMERIC_POINTER(ra0);
+ xa2 = NUMERIC_POINTER(ra2);
+ xa4 = NUMERIC_POINTER(ra4);
+ xb0 = NUMERIC_POINTER(rb0);
+ xb2 = NUMERIC_POINTER(rb2);
+ xb4 = NUMERIC_POINTER(rb4);
+
+
+
for(i = 0 ; i < totseqs ; i++){
for(j = 0 ; j < totseqs ; j++){
xka[n] = ka[i][j];
xks[n] = ks[i][j];
xvka[n] = vka[i][j];
xvks[n] = vks[i][j];
- if(debugon) Rprintf("C> i = %d, j = %d, n = %d, ka = %lf, ks = %lf, vka = %lf, vks = %lf\n", i, j, n, ka[i][j], ks[i][j], vka[i][j], vks[i][j]);
+ xl0[n] = l0[i][j];
+ xl2[n] = l2[i][j];
+ xl4[n] = l4[i][j];
+ xa0[n] = a0[i][j];
+ xa2[n] = a2[i][j];
+ xa4[n] = a4[i][j];
+ xb0[n] = b0[i][j];
+ xb2[n] = b2[i][j];
+ xb4[n] = b4[i][j];
+ if(debugon) Rprintf("C> i = %d, j = %d, n = %d, ka = %lf, ks = %lf, vka = %lf, vks = %lf, l0 = %lf, l2 = %lf, l4 = %lf, a0 = %lf, a2 = %lf, a4 = %lf, b0 = %lf, b2 = %lf, b4 = %lf\n", i, j, n, ka[i][j], ks[i][j], vka[i][j], vks[i][j],l0[i][j],l2[i][j],l4[i][j],a0[i][j],a2[i][j],a4[i][j],b0[i][j],b2[i][j],b4[i][j] );
n++;
}
- }
-
+ }
SET_ELEMENT(res, 0, rka);
SET_ELEMENT(res, 1, rks);
SET_ELEMENT(res, 2, rvka);
SET_ELEMENT(res, 3, rvks);
+ SET_ELEMENT(res, 4, rl0);
+ SET_ELEMENT(res, 5, rl2);
+ SET_ELEMENT(res, 6, rl4);
+ SET_ELEMENT(res, 7, ra0);
+ SET_ELEMENT(res, 8, ra2);
+ SET_ELEMENT(res, 9, ra4);
+ SET_ELEMENT(res, 10, rb0);
+ SET_ELEMENT(res, 11, rb2);
+ SET_ELEMENT(res, 12, rb4);
+ UNPROTECT(14);
- UNPROTECT(5);
-
if(debugon) Rprintf("C> %s", "End of C level....................\n");
return(res);
}
@@ -305,14 +418,13 @@
int fastlwl(char **seq, int nbseq, int lgseq, double **ka, double **ks,
double **tti0, double **tti1, double **tti2, double **ttv0,
double **ttv1, double **ttv2, double **tl0, double **tl1,
- double **tl2, double **vka, double **vks)
+ double **tl2, double **vka, double **vks, double **l0, double **l2,double **l4, double **a0, double **a2,double **a4, double **b0, double **b2,double **b4)
{
const double trois = 3.0;
double l[3], a[3], b[3], p[3], q[3], ti[3], tv[3], cc[3],
aaa[3], bb[3], flgseq, va[3], vb[3],es1,es2;
char cod1[3], cod2[3];
int i, j, ii, num1, num2, sat, sat1, sat2;
-
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
@@ -345,6 +457,9 @@
tv[1] += *(ttv1[num1] + num2);
tv[2] += *(ttv2[num1] + num2);
}
+ l0[i][j]=l[0];
+ l2[i][j]=l[1];
+ l4[i][j]=l[2];
for (ii = 0; ii < 3; ii++) {
p[ii] = ti[ii] / l[ii];
q[ii] = tv[ii] / l[ii];
@@ -381,8 +496,43 @@
vka[i][j]=ka[i][j] = 9.999999;
sat2 = 1;
}
+ a0[i][j]=a[0];
+ a2[i][j]=a[1];
+ a4[i][j]=a[2];
+ b0[i][j]=b[0];
+ b2[i][j]=b[1];
+ b4[i][j]=b[2];
+
}
}
+
+ /* -- addition fevrier 2012 --*/
+ /* L0, L2, L4: # of non-synonymous sites, of 2-fold synonymous sites, of 4-fold synonymous sites
+
+ A0, A2, A4: # of transitional changes at non-synonymous, 2-fold, and 4-fold synonymous sites
+
+ 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]
+ */
+
+ /* l024[0]=l[0];
+ l024[1]=l[1];
+ l024[2]=l[2];
+
+ a024[0]=a[0];
+ a024[1]=a[1];
+ a024[2]=a[2];
+
+ b024[0]=b[0];
+ b024[1]=b[1];
+ b024[2]=b[2];*/
+
+
+
if (sat1 == 1)
sat = 1;
if (sat2 == 1)
@@ -446,7 +596,7 @@
char transf(char nt1, char nt2)
{
if (nt1 == nt2) {
- printf("Same nt, patate.\n");
+ Rprintf("Same nt, patate.\n");
return 'S';
}
if ((nt1 == 'A') && (nt2 == 'C'))
More information about the Seqinr-commits
mailing list