Re: [R-sig-phylo] ape package in R

2020-09-19 Thread Emmanuel Paradis 
Hi,

Actually, Ranjita can use the number of gene copies if it is reasonable to 
assume that these numbers evolve among samples, maybe a stepwise model with 
transitions: 

n -> n + 1
n -> n - 1

In that case, a simple approach is to calculate the Euclidean distances among 
samples and do an NJ or ME tree. Another approach is to use phangorn and a 
parsimony or maximum likelihood (ML) method depending on the model you want to 
assume.

The functions needed are in different packages (as usual you are advised to 
check the help pages and the options):

utils::read.csv
stats::dist
ape::nj
ape::fastme.bal
phangorn::phyDat
phangorn::parsimony

And as suggested by Salvador, if you have the sequences of the gene copies for 
each sample (which would be a very nice study), you can get the evolutionary 
distances with:

ape::dist.dna

or fit an evolutionary model by ML with phangorn::pml.

HTH.

Best,

Emmanuel

- Le 20 Sep 20, à 3:14, Salvador Espada Hinojosa salvador.esp...@gmail.com 
a écrit :

> If you want to use ape for making a phylogenetic tree, I think you need the
> actual sequences of the genes, and not only the copy numbers
> 
> El sáb., 19 sept. 2020 a las 21:40, Ranjita Thapa ()
> escribió:
> 
>> Hi,
>>
>> I want to construct a phylogenetic tree with ape package. I have a csv
>> file with gene ids and copy number for each gene across different samples.
>> I want to construct phylogenetic tree. I am wondering how could I convert
>> the csv file to input file for ape. Could you please suggest me how can I
>> convert csv file to input file that is acceptable by ape?
>>
>> Thanks
>> Ranjita___
>> R-sig-phylo mailing list - R-sig-phylo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>> Searchable archive at
>> http://www.mail-archive.com/r-sig-phylo@r-project.org/
>>
> 
> 
> --
> Salvador Espada Hinojosa
> Monika Bright's Lab
> Department of Functional and Evolutionary Ecology
> University of Vienna
> Althanstraße, 14
> 1090 Vienna (Austria)
> Mobile: (+43) 660 5952724
> skype: salva_e
> http://salvae.net 
> 
>   [[alternative HTML version deleted]]
> 
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] ape package in R

2020-09-19 Thread Salvador Espada Hinojosa 
If you want to use ape for making a phylogenetic tree, I think you need the
actual sequences of the genes, and not only the copy numbers

El sáb., 19 sept. 2020 a las 21:40, Ranjita Thapa ()
escribió:

> Hi,
>
> I want to construct a phylogenetic tree with ape package. I have a csv
> file with gene ids and copy number for each gene across different samples.
> I want to construct phylogenetic tree. I am wondering how could I convert
> the csv file to input file for ape. Could you please suggest me how can I
> convert csv file to input file that is acceptable by ape?
>
> Thanks
> Ranjita___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at
> http://www.mail-archive.com/r-sig-phylo@r-project.org/
>


-- 
Salvador Espada Hinojosa
Monika Bright's Lab
Department of Functional and Evolutionary Ecology
University of Vienna
Althanstraße, 14
1090 Vienna (Austria)
Mobile: (+43) 660 5952724
skype: salva_e
http://salvae.net 

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] ape package in R

2020-09-19 Thread Ranjita Thapa 
Hi,

I want to construct a phylogenetic tree with ape package. I have a csv file 
with gene ids and copy number for each gene across different samples. I want to 
construct phylogenetic tree. I am wondering how could I convert the csv file to 
input file for ape. Could you please suggest me how can I convert csv file to 
input file that is acceptable by ape? 

Thanks
RanjitaGENE,SAMN05064989,SAMN05064990,SAMN05064991,SAMN05064992,SAMN05064993,SAMN05064994,SAMN05064995,SAMN05064996,SAMN05064997,SAMN05064998,SAMN05064999,SAMN05065000,SAMN05065001
MD00G1000200,1.999678,2.644484,1.843331,1.868703,2.034242,1.752502,1.859093,2.039209,2.441564,2.711497,2.184576,1.473055,1.876445
MD00G1000300,1.415342,1.48209,1.549128,1.706571,2.046146,1.158263,1.646863,1.173129,1.59711,1.738894,1.419168,1.371778,1.568272
MD00G1000400,1.398499,1.797129,1.802143,1.166775,1.390045,1.236685,1.314259,1.58712,1.303329,1.782833,1.305629,1.461014,1.712468
MD00G1000500,2.41541,4.003027,2.883078,2.829228,2.876289,2.686079,2.253276,2.52543,2.45406,3.053208,3.668503,3.119604,2.609322
MD00G1000600,1.357889,1.514805,1.754632,1.905162,1.853981,1.248514,1.889641,1.832417,1.537696,1.39775,1.808635,1.389049,1.623491
MD00G1000700,1.629476,2.228273,2.048752,1.74234,1.974563,1.961737,2.069386,1.617159,1.82523,1.496132,2.092336,2.065827,1.946314
MD00G1000800,1.896418,2.22733,1.985247,1.983065,1.83869,1.9923,1.68222,1.750361,1.779146,2.209676,1.973086,1.981803,1.630689
MD00G1000900,3.134369,2.734057,2.550392,3.401735,3.493074,2.421759,2.86019,3.21655,2.768818,2.675039,3.344572,2.324965,2.332387
MD00G1001000,1.933826,1.876517,2.021223,1.739334,1.900607,1.864082,1.791559,1.847408,1.798716,2.155571,1.955044,2.377852,1.829788
MD00G1001100,0.945261,1.069836,1.105532,1.39809,1.147328,1.434332,1.195645,0.999553,0.836034,1.329002,1.742932,0.969243,1.431592
MD00G1001200,2.231896,2.332669,3.101326,2.176223,3.526095,2.718395,2.382203,2.327535,2.215922,1.528227,2.62628,2.409097,3.385232
MD00G1001300,4.502762,3.921681,4.116981,3.10278,3.640254,3.509639,3.559179,3.580477,3.159958,3.558723,4.421366,4.787455,4.684205
MD00G1001400,2.236153,2.537318,2.032665,2.0862,2.509157,2.348267,2.240651,2.156838,1.998623,2.466158,2.506128,2.508376,2.115701
MD00G1001500,2.238276,2.20046,2.405819,2.438955,2.578238,2.362307,2.68222,2.044836,2.145149,2.553174,2.299561,2.901258,2.219329
MD00G1001700,12.567245,11.491921,12.917022,9.731756,10.496296,10.999282,10.04299,11.116823,11.920293,10.890662,10.974454,9.573012,10.618516
MD00G1001800,2.255836,2.143793,2.275005,2.017084,2.705015,2.045186,2.490839,2.174137,2.108093,2.173543,2.182403,1.969187,2.330892
MD00G1001900,2.611557,3.243386,2.359358,2.605919,2.039701,2.751724,2.535086,2.946174,2.491261,2.874563,2.854939,2.685266,3.34816
MD00G1002000,2.41494,2.463868,2.470049,2.613261,2.576815,2.296066,2.277826,2.382715,2.329795,2.857351,2.16088,2.585115,2.3057
MD00G1002100,1.856245,1.766623,2.192718,1.9758,1.372465,1.346592,1.294279,2.434389,2.401637,1.943977,2.171607,2.156644,1.889808
MD00G1002200,2.415919,2.234855,2.118146,2.487616,1.265053,1.655353,1.578375,1.438708,1.571695,1.564428,2.190602,2.621737,2.232404
MD00G1002500,1.794938,2.574585,2.19974,2.300694,2.203994,2.157548,2.21563,2.38745,1.942388,2.597048,2.321994,2.174149,1.843977
MD00G1002600,2.511032,2.147372,2.178689,2.033836,2.076008,1.957138,3.11344,2.450917,3.184336,2.385136,3.445098,2.460429,2.427365
MD00G1002700,2.466562,2.200628,2.402447,1.768002,2.576956,2.776423,2.992975,2.044115,3.389897,2.034527,3.470035,4.297166,2.281435
MD00G1002800,3.994084,3.585871,4.357335,4.134188,3.843827,4.029738,5.410013,4.366464,5.32427,5.007784,5.071291,4.958608,4.356433
MD00G1003400,2.028538,2.072182,2.080781,2.25311,2.311936,2.206668,2.107647,1.86438,1.809568,2.150551,2.07075,2.415351,2.620719
MD00G1003500,2.226612,1.528872,1.259334,1.752868,1.8164,1.521194,1.386067,1.855044,2.288861,1.823343,1.259023,1.369945,1.69697
MD00G1003600,1.073772,1.715179,2.027966,1.970396,2.462572,1.751237,2.013319,1.483089,1.99988,1.443287,2.069001,1.939983,1.855052
MD00G1003800,2.390463,2.873044,2.403015,2.404447,2.018727,2.537998,2.189413,2.712841,2.135105,2.285238,2.239648,2.636597,2.004905
MD00G1003900,1.703446,1.970978,1.510863,1.712658,1.421343,2.311724,1.868775,2.210896,1.416742,1.720507,1.567725,1.078489,1.795938
MD00G1004000,1.703446,1.970978,1.510863,1.712658,1.421343,2.311724,1.868775,2.210896,1.416742,1.720507,1.567725,1.078489,1.795938
MD00G1004100,2.077799,1.755617,1.755502,1.791851,1.969101,2.207733,1.756074,2.178037,1.699324,1.789315,1.634722,1.565095,2.105419
MD00G1004200,2.548076,2.641978,2.604369,2.44473,1.836607,1.962447,2.378719,2.61517,2.478992,2.56008,2.582974,1.634768,1.673753
MD00G1004300,1.612036,1.967232,2.218305,1.745579,1.786646,1.702544,2.079315,1.724835,2.201695,1.814631,2.220065,2.921447,2.336502
MD00G1004400,2.164274,3.437788,2.550175,2.197593,2.42354,2.774165,2.647565,2.981612,3.011943,2.766472,2.750117,2.558302,2.101304

Re: [R-sig-phylo] ape-package: unknown states in ace analyses?

2018-06-25 Thread Théo Léger 
Dear Emmanuel and Liam,


thank you for your answers. I tried ace with the modifications suggested by 
Emmanuel and make.simmap from phytools suggested by Liam, both worked well for 
me. I think it would be worth to add these modifications to the ace function, 
or at least a note in the manual stating that missing info should be coded "NA" 
(which wasn't obvious to me). Surprisingly I got an error message when running 
an ace with "NA" instead of "?" in my character file:

Error in matexpo(Q * EL[i]) : NA/NaN/Inf in foreign function call (arg 1)
In addition: Warning message:
In ace(char3_NA, tree_cramb, model = "ER", type = "discrete") :
model fit suspicious: gradients apparently non-finite

However, I did not pay more attention to it as first two solutions work fine.

Cheers,
Th�o






De : Emmanuel Paradis 
Envoy� : vendredi, 22 juin 2018 15:45:05
� : Th�o L�ger; r-sig-phylo@r-project.org
Objet : Re: [R-sig-phylo] ape-package: unknown states in ace analyses?

Two more comments:

1/ The first insertion in the code of ace() below is not needed if the
user replaces "?" (or any coding for unknown/uncertain state) with NA --
which seems reasonable to do.

2/ This approach of taking uncertainty into account in likelihood
ancestral state reconstruction is also used for DNA sequences and
implemented in phangorn where uncertainty is (implicitly) taken from the
IUPAC ambiguity code. Because there are here four states, there is more
than one possibility of ambiguity, but the same logic applies: if the
tip "may be" in state x, then set column x to one in the matrix of
likelihoods:

1 0 0 0 # if the base observed is A
0 0 1 0 # if the base observed is G
1 0 1 0 # if the base observed is R
...
1 1 1 1 # if the base observed is N

This can be applied to any character with more than two states.

Best,

Emmanuel

Le 22/06/2018 � 11:41, Emmanuel Paradis a �crit :
> Hi Th�o,
>
> It is possible to modify the code of ace() to take uncertain character
> states into account. If you edit the code with fix(ace), after these
> four lines:
>
> if (method != "ML")
> stop("only ML estimation is possible for discrete characters.")
> if (any(phy$edge.length <= 0))
> stop("some branches have length zero or negative")
>
> you insert:
>
>  if (any(x == "?")) x[x == "?"] <- NA
>
> Then around 30 lines below, after this line:
>
> liks[cbind(TIPS, x)] <- 1
>
> you insert:
>
> if (anyNA(x)) liks[which(is.na(x)), ] <- 1
>
> Save and close. I tested this with a simple example:
>
> tr <- compute.brlen(stree(8, "b"), 1)
> x <- rep(0:1, each = 4)
> plot(tr)
> tiplabels(x)
>
> So there are two clades in state 0 and state 1, respectively. With no
> uncertainty, the node likelihoods are well unbalanced except at the root:
>
> R> ace(x, tr, "d")$lik.anc
> 0   1
> [1,] 0.5 0.5
> [2,] 0.953204905 0.046795095
> [3,] 0.995642789 0.004357211
> [4,] 0.995642789 0.004357211
> [5,] 0.046795095 0.953204905
> [6,] 0.004357211 0.995642789
> [7,] 0.004357211 0.995642789
>
> Let's say the first value is unknown:
>
> x[1] <- "?"
>
> Now the node likelihoods are slightly "in favour" of state 1 because of
> the possibility of this state within the "state 0" clade:
>
> R> ace(x, tr, "d")$lik.anc
> 0   1
> [1,] 0.479307333 0.520692667
> [2,] 0.904647485 0.095352515
> [3,] 0.943101903 0.056898097
> [4,] 0.990270355 0.009729645
> [5,] 0.053105869 0.946894131
> [6,] 0.005881435 0.994118565
> [7,] 0.005881435 0.994118565
>
> For comparison, here's what is returned by the current version in ape:
>
> R> ace(x, tr, "d")$lik.anc
> ?   0  1
> [1,] 0.093607676 0.404434328 0.50195800
> [2,] 0.096246928 0.787370632 0.11638244
> [3,] 0.201389158 0.750401358 0.04820948
> [4,] 0.011648214 0.974955230 0.01339656
> [5,] 0.017176009 0.050038984 0.93278501
> [6,] 0.003422534 0.006275985 0.99030148
> [7,] 0.003422534 0.006275985 0.99030148
>
> You may try this modification with your data and give feed-back: if it
> useful, I'll fix that in ape.
>
> Cheers,
>
> Emmanuel
>
> Le 21/06/2018 � 23:32, Th�o L�ger a �crit :
>> Hello,
>>
>>
>> I am working on the phylogeny of a subfamily of moths and I use ace
>> from the ape-package to reconstruct the ancestral state of a bunch of
>> morphological characters.
>>
>> I encountered a problem with the few unknown states I have on my
>> matrix (coded "?", either because

Re: [R-sig-phylo] ape-package: unknown states in ace analyses?

2018-06-22 Thread Jacob Berv 
This is so cool. Thanks Liam!
J

> On Jun 21, 2018, at 8:22 PM, Liam J. Revell  wrote:
> 
> Dear Theo.
> 
> This can be done with the function rerootingMethod as described here 
> http://blog.phytools.org/2013/04/estimating-ancestral-states-when.html, 
> although it is *very* important to note that rerootingMethod is only valid if 
> the fitted transition model is symmetric (e.g., "ER", "SYM", or any custom 
> model in which i->j == j->i for all i & j). It can also be done with 
> make.simmap for both symmetric & non-symmetric transition models. If you 
> sample enough stochastic maps & then run summary on the set of maps the 
> posterior probabilities will converge on the marginal ancestral states under 
> the defaults for make.simmap, so long as a flat prior is used on internal 
> nodes. This is also on my blog here: 
> http://blog.phytools.org/2013/03/estimating-ancestral-states-when-tips.html.
> 
> All the best, Liam
> 
> Liam J. Revell, Associate Professor of Biology
> University of Massachusetts Boston
> & Profesor Asociado, Programa de Biología
> Universidad del Rosario
> web: http://faculty.umb.edu/liam.revell/
> 
> On 6/21/2018 5:32 PM, Théo Léger wrote:
>> Hello,
>> I am working on the phylogeny of a subfamily of moths and I use ace from the 
>> ape-package to reconstruct the ancestral state of a bunch of morphological 
>> characters.
>> I encountered a problem with the few unknown states I have on my matrix 
>> (coded "?", either because material for examination was missing or the state 
>> could not fit in any of the categories): they are treated as other 
>> characters. Is there a way to ignore them? Is there a way to estimate the 
>> state for the species with missing information? I know it is possible to 
>> estimate the state at a tip with missing information in functions like 
>> AncTresh (using bayesian reconstruction) from the phytools package, but I am 
>> not sure if ML reconstruction can do it.
>> Thanks in advance to those who can enlighten me!
>> Cheers,
>> Th�o L�ger
>>  [[alternative HTML version deleted]]
>> ___
>> R-sig-phylo mailing list - R-sig-phylo@r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
>> 
> 
> ___
> R-sig-phylo mailing list - R-sig-phylo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
> Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] ape-package: unknown states in ace analyses?

2018-06-22 Thread Emmanuel Paradis 

Two more comments:

1/ The first insertion in the code of ace() below is not needed if the 
user replaces "?" (or any coding for unknown/uncertain state) with NA -- 
which seems reasonable to do.


2/ This approach of taking uncertainty into account in likelihood 
ancestral state reconstruction is also used for DNA sequences and 
implemented in phangorn where uncertainty is (implicitly) taken from the 
IUPAC ambiguity code. Because there are here four states, there is more 
than one possibility of ambiguity, but the same logic applies: if the 
tip "may be" in state x, then set column x to one in the matrix of 
likelihoods:


1 0 0 0 # if the base observed is A
0 0 1 0 # if the base observed is G
1 0 1 0 # if the base observed is R
...
1 1 1 1 # if the base observed is N

This can be applied to any character with more than two states.

Best,

Emmanuel

Le 22/06/2018 à 11:41, Emmanuel Paradis a écrit :

Hi Théo,

It is possible to modify the code of ace() to take uncertain character 
states into account. If you edit the code with fix(ace), after these 
four lines:


    if (method != "ML")
    stop("only ML estimation is possible for discrete characters.")
    if (any(phy$edge.length <= 0))
    stop("some branches have length zero or negative")

you insert:

     if (any(x == "?")) x[x == "?"] <- NA

Then around 30 lines below, after this line:

    liks[cbind(TIPS, x)] <- 1

you insert:

    if (anyNA(x)) liks[which(is.na(x)), ] <- 1

Save and close. I tested this with a simple example:

tr <- compute.brlen(stree(8, "b"), 1)
x <- rep(0:1, each = 4)
plot(tr)
tiplabels(x)

So there are two clades in state 0 and state 1, respectively. With no 
uncertainty, the node likelihoods are well unbalanced except at the root:


R> ace(x, tr, "d")$lik.anc
    0   1
[1,] 0.5 0.5
[2,] 0.953204905 0.046795095
[3,] 0.995642789 0.004357211
[4,] 0.995642789 0.004357211
[5,] 0.046795095 0.953204905
[6,] 0.004357211 0.995642789
[7,] 0.004357211 0.995642789

Let's say the first value is unknown:

x[1] <- "?"

Now the node likelihoods are slightly "in favour" of state 1 because of 
the possibility of this state within the "state 0" clade:


R> ace(x, tr, "d")$lik.anc
    0   1
[1,] 0.479307333 0.520692667
[2,] 0.904647485 0.095352515
[3,] 0.943101903 0.056898097
[4,] 0.990270355 0.009729645
[5,] 0.053105869 0.946894131
[6,] 0.005881435 0.994118565
[7,] 0.005881435 0.994118565

For comparison, here's what is returned by the current version in ape:

R> ace(x, tr, "d")$lik.anc
    ?   0  1
[1,] 0.093607676 0.404434328 0.50195800
[2,] 0.096246928 0.787370632 0.11638244
[3,] 0.201389158 0.750401358 0.04820948
[4,] 0.011648214 0.974955230 0.01339656
[5,] 0.017176009 0.050038984 0.93278501
[6,] 0.003422534 0.006275985 0.99030148
[7,] 0.003422534 0.006275985 0.99030148

You may try this modification with your data and give feed-back: if it 
useful, I'll fix that in ape.


Cheers,

Emmanuel

Le 21/06/2018 à 23:32, Théo Léger a écrit :

Hello,


I am working on the phylogeny of a subfamily of moths and I use ace 
from the ape-package to reconstruct the ancestral state of a bunch of 
morphological characters.


I encountered a problem with the few unknown states I have on my 
matrix (coded "?", either because material for examination was missing 
or the state could not fit in any of the categories): they are treated 
as other characters. Is there a way to ignore them? Is there a way to 
estimate the state for the species with missing information? I know it 
is possible to estimate the state at a tip with missing information in 
functions like AncTresh (using bayesian reconstruction) from the 
phytools package, but I am not sure if ML reconstruction can do it.


Thanks in advance to those who can enlighten me!


Cheers,

Théo Léger


___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at 
http://www.mail-archive.com/r-sig-phylo@r-project.org/








___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] ape-package: unknown states in ace analyses?

2018-06-22 Thread Emmanuel Paradis 

Hi Théo,

It is possible to modify the code of ace() to take uncertain character 
states into account. If you edit the code with fix(ace), after these 
four lines:


   if (method != "ML")
   stop("only ML estimation is possible for discrete characters.")
   if (any(phy$edge.length <= 0))
   stop("some branches have length zero or negative")

you insert:

if (any(x == "?")) x[x == "?"] <- NA

Then around 30 lines below, after this line:

   liks[cbind(TIPS, x)] <- 1

you insert:

   if (anyNA(x)) liks[which(is.na(x)), ] <- 1

Save and close. I tested this with a simple example:

tr <- compute.brlen(stree(8, "b"), 1)
x <- rep(0:1, each = 4)
plot(tr)
tiplabels(x)

So there are two clades in state 0 and state 1, respectively. With no 
uncertainty, the node likelihoods are well unbalanced except at the root:


R> ace(x, tr, "d")$lik.anc
   0   1
[1,] 0.5 0.5
[2,] 0.953204905 0.046795095
[3,] 0.995642789 0.004357211
[4,] 0.995642789 0.004357211
[5,] 0.046795095 0.953204905
[6,] 0.004357211 0.995642789
[7,] 0.004357211 0.995642789

Let's say the first value is unknown:

x[1] <- "?"

Now the node likelihoods are slightly "in favour" of state 1 because of 
the possibility of this state within the "state 0" clade:


R> ace(x, tr, "d")$lik.anc
   0   1
[1,] 0.479307333 0.520692667
[2,] 0.904647485 0.095352515
[3,] 0.943101903 0.056898097
[4,] 0.990270355 0.009729645
[5,] 0.053105869 0.946894131
[6,] 0.005881435 0.994118565
[7,] 0.005881435 0.994118565

For comparison, here's what is returned by the current version in ape:

R> ace(x, tr, "d")$lik.anc
   ?   0  1
[1,] 0.093607676 0.404434328 0.50195800
[2,] 0.096246928 0.787370632 0.11638244
[3,] 0.201389158 0.750401358 0.04820948
[4,] 0.011648214 0.974955230 0.01339656
[5,] 0.017176009 0.050038984 0.93278501
[6,] 0.003422534 0.006275985 0.99030148
[7,] 0.003422534 0.006275985 0.99030148

You may try this modification with your data and give feed-back: if it 
useful, I'll fix that in ape.


Cheers,

Emmanuel

Le 21/06/2018 à 23:32, Théo Léger a écrit :

Hello,


I am working on the phylogeny of a subfamily of moths and I use ace from the 
ape-package to reconstruct the ancestral state of a bunch of morphological 
characters.

I encountered a problem with the few unknown states I have on my matrix (coded 
"?", either because material for examination was missing or the state could not 
fit in any of the categories): they are treated as other characters. Is there a way to 
ignore them? Is there a way to estimate the state for the species with missing 
information? I know it is possible to estimate the state at a tip with missing 
information in functions like AncTresh (using bayesian reconstruction) from the phytools 
package, but I am not sure if ML reconstruction can do it.

Thanks in advance to those who can enlighten me!


Cheers,

Théo Léger


___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] ape-package: unknown states in ace analyses?

2018-06-21 Thread Théo Léger 
Hello,


I am working on the phylogeny of a subfamily of moths and I use ace from the 
ape-package to reconstruct the ancestral state of a bunch of morphological 
characters.

I encountered a problem with the few unknown states I have on my matrix (coded 
"?", either because material for examination was missing or the state could not 
fit in any of the categories): they are treated as other characters. Is there a 
way to ignore them? Is there a way to estimate the state for the species with 
missing information? I know it is possible to estimate the state at a tip with 
missing information in functions like AncTresh (using bayesian reconstruction) 
from the phytools package, but I am not sure if ML reconstruction can do it.

Thanks in advance to those who can enlighten me!


Cheers,

Th�o L�ger

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] ape package: improvement suggestion in color handling for mosaic in phydataplot?

2017-05-23 Thread Elizabeth Purdom 
Dear all,

I am trying to use the phydataplot function in the ape package with the style 
mosaic. I am having difficulty getting the colors of the rectangles to be the 
color I desire in an automatic way. I have come to the conclusion that I can’t 
do and it would have to be an enhancement to the package, though I would be 
glad to know I am wrong. 

Here is the example in phydataplot for how to set the colors

set.seed(1378923)
## use type = "mosaic" on a 30x5 matrix:
tr <- rtree(n <- 30)
p <- 5
x <- matrix(sample(3, size = n*p, replace = TRUE), n, p)
dimnames(x) <- list(paste0("t", 1:n), LETTERS[1:p])
plot(tr, x.lim = 35, align.tip = TRUE, adj = 1)
phydataplot(x, tr, "m", 2, width = 2, border = "white", lwd = 3, legend = 
"side")
f <- function(n) c("yellow", "blue", "red")
phydataplot(x, tr, "m", 18, width = 2, border = "white", lwd = 3,
legend = "side", funcol = f)

Notice we set the colors with the function  c("yellow", "blue", "red”). 
However, we see when we plot them, that this assignment resulted in 2=yellow, 
3=blue, 1=red. If instead I want 1=yellow, 2=blue, 3=red, I can use this plot I 
just made to see how they should be ordered to get my desired color assignment, 
and rearrange the colors: 

plot(tr, x.lim = 35, align.tip = TRUE, adj = 1)
phydataplot(x, tr, "m", 2, width = 2, border = "white", lwd = 3, legend = 
"side")
f <- function(n) c("blue", "red", "yellow")
phydataplot(x, tr, "m", 18, width = 2, border = "white", lwd = 3,
legend = "side", funcol = f)

However, I see no way to do this without first plotting them, manually looking 
at the arbitrary order created, and then fixing it. Yet, I would like to do 
this in an automatic fashion — I’m trying to build a function on top of this 
function — so I need to know ahead of time how phydataplot will assign the 
colors (and so if my tree changes it will still give the desired colors). 
Looking at the internal function, the relevant stream of code is the following 
(i.e. cutting away the parts not needed for determining the colors) and shows 
that the assignment of colors to data is entirely dependent on the order of the 
data after the function internally matches it to the tree:

lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
x <- .matchDataPhylo(x, phy)
n <- length(phy$tip.label)
one2n <- seq_len(n)
lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
y1 <- lastPP$yy[one2n]
o <- order(y1)
x <- if (style == "image") x[o, o]
else if (is.vector(x)) x[o]
else x[o, ]
nux <- length(ux <- unique.default(x))
x <- match(x, ux)
co <- funcol(nux)
rect(xl, yb, xr, yt, col = co[x], xpd = TRUE, ...)

This means that the colors provide by the user are assigned to the values of 
unique(x) *after* x is reordered internally by .matchDataPhylo and then again 
reordered x[o,], where o is the order of the tips of the tree in the plot, I 
believe. This is not something returned to the user, unfortunately, by 
plot.phylo, but is assigned to an hidden environment variable :
assign("last_plot.phylo", c(L, list(edge = xe, xx = xx, yy = yy)), 
envir = .PlotPhyloEnv)

These seems to make it next to impossible to get the right colors without just 
trying it out interactively. If anyone had any other alternatives I’d love to 
hear them.

I’ve tried making a little function:
getColFun<-function(x,phy,namedColors){
x <- ape:::.matchDataPhylo(x, phy)
n <- length(phy$tip.label)
one2n <- seq_len(n)
lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
y1 <- lastPP$yy[one2n]
o <- order(y1)
ux<-unique.default(x[o])
m<-match(as.character(ux),names(namedColors))
function(n){namedColors[m]}
}
But it’s a bit awkward because it relies on .PlotPhyloEnv which is not always 
available if you haven’t loaded the package, as well as the internal function 
.matchDataPhylo. This is problematic for putting it into a package. Is there a 
function I’m not aware of that would be helpful?

In terms of enhancement, it would be nice if instead phydataplot accepted a 
named vector of colors, where the vector of colors is required to have names 
that match the values of the data in x, and then the phydataplot internally 
would figure out how to match them together. Less desirable would be if the 
output of plot.phylo returned the information in ‘L’ that is saved in the 
.PlotPhyloEnv, then I could at least recreate a more stable version of my 
little function above.

(it would also be nice if the legend of the colors was in some reasonable 
sorted order, i.e. 1,2,3, and not in the arbitrary order)

Thanks,
Elizabeth Purdom






[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org

Re: [R-sig-phylo] ape package gives 'Error in FI[i]:LA[i] : NA/NaN argument'

2015-09-12 Thread Uwe Menzel 
Emmanuel Paradis  writes:

> 
> Hi Matt,
> 
> This is because this accession number does not point to the sequence. 
> For this particular one, you could use:
> 
> seq1 <- read.GenBank("U00096.3")
> 
> Best,
> 
> Emmanuel
> 
> Le 27/02/2014 01:35, Matt Curcio a écrit :
> > Greetings all,,
> > I received this error while using R version 2.15.1 and ape 3.0.11.
> > Error in FI[i]:LA[i] : NA/NaN argument
> >
> > I have tried several approaches:
> >
> > ref = "NC_000913.3"
> > seq1 = read.GenBank(ref)
> > Error in FI[i]:LA[i] : NA/NaN argument
> >
> > ref = read.csv(file="test2.csv", head=T)
> >
> > ref = str(ref)
> >
> > # Where test2.csv is in two rows
> >
> > 'data.frame':   1 obs. of  1 variable:
> > $ seq: Factor w/ 1 level "NC+AF8-000913": 1
> >
> > seq2 = read.GenBank(ref)
> > Error in FI[i]:LA[i] : NA/NaN argument
> >
> >
> > I was wondering if someone could lend a hand based on their own experience.
> >
> >
> >
> 
> 

Hi,

I do have similar problems. 
I try downloading 1814 sequences from K. flaccidum.

I just try with two of them:

accessions = c("DF238762", "DF238763")
flaccidum_scaffolds <- read.GenBank(accessions)

# Error in FI[i]:LA[i] : NA/NaN argument

If I try out these accesions directly in the Browser at
http://www.ncbi.nlm.nih.gov/genbank/ ,
it works very well.

So, how do I infer the "correct" accessions to be used in read.GenBank?
Or, how do I come from "NC_000913.3" to "U00096.3" in the example above?
That doesn't seem straightforward 

Thanks, Uwe









___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] Ape package.

2014-04-22 Thread Frederico Faleiro 
You need put the path of the directory with this simbol / and not this
\, so the correct path is F:*/*tetrapods_tree.txt.

Cheers!


On Sun, Apr 20, 2014 at 5:17 PM, Alexey Fomin lesh...@rambler.ru wrote:

  Hi!

 Whether someone work with ape package?
 I have long tree (see attached, it is in newick format, it is ultra
 metric, it is from http://www.onezoom.org/user/tetrapods_tree.js ).
 I need dates of nodes.
 Someone advised me to do:

 1) library(ape)
 2) tree-read.tree(?yourtree.nex?)
 3) br.times-branching.times(tree)

 But I can not not do the second step.

 I tried (after library(ape))
 tree-read.tree(F:\tetrapods_tree.txt)
 But it is gives:
 Error: unexpected input in tree-read.tree(F:\

 What can I do?

 Alexey.

 ___
 R-sig-phylo mailing list - R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
 Searchable archive at
 http://www.mail-archive.com/r-sig-phylo@r-project.org/


[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] Ape package.

2014-04-20 Thread Klaus Schliep 
The tree file is not proper formatted, it is missing a ; in the end.

Regards,
Klaus


On Sun, Apr 20, 2014 at 10:17 PM, Alexey Fomin lesh...@rambler.ru wrote:

  Hi!

 Whether someone work with ape package?
 I have long tree (see attached, it is in newick format, it is ultra
 metric, it is from http://www.onezoom.org/user/tetrapods_tree.js ).
 I need dates of nodes.
 Someone advised me to do:

 1) library(ape)
 2) tree-read.tree(?yourtree.nex?)
 3) br.times-branching.times(tree)

 But I can not not do the second step.

 I tried (after library(ape))
 tree-read.tree(F:\tetrapods_tree.txt)
 But it is gives:
 Error: unexpected input in tree-read.tree(F:\

 What can I do?

 Alexey.

 ___
 R-sig-phylo mailing list - R-sig-phylo@r-project.org
 https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
 Searchable archive at
 http://www.mail-archive.com/r-sig-phylo@r-project.org/




-- 
Klaus Schliep

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] ape package.

2014-04-11 Thread Alexey Fomin 
Whether anybody help?
I can non start the ape package.
 when I do

 library(ape)

it give:

Error in library(ape) : there is no package called 'ape' 

(I use R-3.0.3-win)

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] ape package gives 'Error in FI[i]:LA[i] : NA/NaN argument'

2014-02-27 Thread Emmanuel Paradis 

Hi Matt,

This is because this accession number does not point to the sequence. 
For this particular one, you could use:


seq1 - read.GenBank(U00096.3)

Best,

Emmanuel

Le 27/02/2014 01:35, Matt Curcio a écrit :

Greetings all,,
I received this error while using R version 2.15.1 and ape 3.0.11.
Error in FI[i]:LA[i] : NA/NaN argument

I have tried several approaches:

ref = NC_000913.3
seq1 = read.GenBank(ref)
Error in FI[i]:LA[i] : NA/NaN argument

ref = read.csv(file=test2.csv, head=T)

ref = str(ref)

# Where test2.csv is in two rows

'data.frame':   1 obs. of  1 variable:
$ seq: Factor w/ 1 level NC+AF8-000913: 1

seq2 = read.GenBank(ref)
Error in FI[i]:LA[i] : NA/NaN argument


I was wondering if someone could lend a hand based on their own experience.





___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

[R-sig-phylo] ape package gives 'Error in FI[i]:LA[i] : NA/NaN argument'

2014-02-26 Thread Matt Curcio 
Greetings all,,
I received this error while using R version 2.15.1 and ape 3.0.11.
Error in FI[i]:LA[i] : NA/NaN argument

I have tried several approaches:

ref = NC_000913.3
seq1 = read.GenBank(ref)
Error in FI[i]:LA[i] : NA/NaN argument

ref = read.csv(file=test2.csv, head=T)

ref = str(ref)

# Where test2.csv is in two rows

'data.frame':   1 obs. of  1 variable:
$ seq: Factor w/ 1 level NC+AF8-000913: 1

seq2 = read.GenBank(ref)
Error in FI[i]:LA[i] : NA/NaN argument


I was wondering if someone could lend a hand based on their own experience.



-- 
Thank You,

Matt Curcio
Bioinformatics Dept.
Worcester Polytechnic Institute
Worcester, MA, USA
M: 401-316-5358
E: matt.curcio...@gmail.com

[[alternative HTML version deleted]]

___
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/

Re: [R-sig-phylo] ape-package/ possibility to get the transformed residuals (phylogentic corrected residuals)?

2010-10-12 Thread Simon Blomberg 
Actually, the function residuals.gls calculates 3 different types of 
residuals: response or raw residuals, pearson or standardized 
residuals, and normalized residuals. For diagnostic purposes, you 
probably want the normalized residuals. The default is the response  
residuals (although the ?residuals.gls says that the pearson residuals 
are the default. This seems to be a bug in the help file). See 
?residuals.gls


Also, it's better programming practice to use the extractor function 
(residuals.gls in this case), rather than using result$residuals 
directly, since the contents of slots can change on the whim of the 
developer, but the extractor function should, in general, return what 
you ask for.


Cheers,

Simon.

On 13/10/10 03:40, Liam J. Revell wrote:

Hi Jan:

Yes.  gls() returns an object of class gls which contains a vector 
containing the residuals as one of its components.

So for instance you might compute:
 
result-gls(y~x,data=test.data,correlation=corPagel(phy=tree,value=1,fixed=FALSE)) 


in which case:
 result$residuals
contains the model residuals.
You can also get the residuals using the function resid(), i.e.:
 resid(result)

Residual standard error is:
 result$sigma
so residual sum of squares should just be:
 resid.SS-n*result$sigma^2
for n species.

Hope this is helpful. - Liam

Liam J. Revell
NESCent, Duke University
web: http://anolis.oeb.harvard.edu/~liam/
NEW email: lrev...@nescent.org



Werner, Jan wrote:

Dear R-sig-phylo members,
for some analyses I used the gls and gnls functions (package nlme) 
and controlled for phylogenetic effects by the phylogenetic 
correlation structure corPagel provided by the ape-package. My 
question is, is there a possibility to get the transformed residuals 
(phylogentic corrected

residuals) or the phylogenetic corrected residual sum of squares?

Many thanks for your efforts.

Best regards

Jan Werner

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo


--
Simon Blomberg, BSc (Hons), PhD, MAppStat, AStat
Lecturer and Consultant Statistician
School of Biological Sciences
The University of Queensland
St. Lucia Queensland 4072
Australia
T: +61 7 3365 2506
email: S.Blomberg1_at_uq.edu.au
http://www.uq.edu.au/~uqsblomb/

Policies:
1.  I will NOT analyse your data for you.
2.  Your deadline is your problem

Statistics is the grammar of science - Karl Pearson.

___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo