-
Notifications
You must be signed in to change notification settings - Fork 0
/
metazoan_hypothesis_test_Luis.Rev
103 lines (80 loc) · 4.63 KB
/
metazoan_hypothesis_test_Luis.Rev
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
###########################################################################
#
# RevBayes Script: Testing specific hypothesis for metazoans
#
# To run this analysis, call `rb metazoan_hypothesis_test.Rev --args xxx_run_1.trees xxx_run_1.txt`
#
# author: Sebastian Hoehna
#
###########################################################################
## we extract the arguments and store them
INPUT_FILENAME = args[1]
OUTPUT_FILENAME = args[2]
## read in the tree trace
#trace = readTreeTrace(INPUT_FILENAME, treetype="non-clock", outgroup=clade("SCPO","SACE","MAOR","CRNE","ALMA","SPPU"))
#trace = readTreeTrace(INPUT_FILENAME, treetype="non-clock", nexus=TRUE)
trace <- readTrees(INPUT_FILENAME, treetype="non-clock")
## Define the clades of interest
Porifera = clade("Oscarella", "Amphimedon", "Sycon", "Tethya")
Ctenophora = clade("Pleurobrachia", "Mnemiopsis")
Placozoa = clade("Hoilungia", "Trichoplax")
Cnidaria = clade("Acropora", "Hydra", "Nematostella")
Xenacoelomorpha = clade("Nemertoderma", "Symsagittifera", "Meara", "Xenoturbella", "Pseudoaphanostoma")
Hemichordata = clade("Saccoglossus", "Ptychodera")
Echinodermata = clade("Acanthaster", "Australostichopus", "Ophionereis", "Strongylocentrotus")
Chordata = clade("Branchiostoma", "Gallus", "Xenopus", "Ciona", "Homo", "Danio")
Lophotrochozoa = clade("Schistosoma", "Helobdella", "Lottia", "Capitella")
Ecdysozoa = clade("Pristionchus", "Drosophila", "Ixodes", "Caenorhabditis", "Daphnia")
## compute the posterior probability for each of these clades
pp_Porifera = trace.cladeProbability(Porifera)
pp_Ctenophora = trace.cladeProbability(Ctenophora)
pp_Cnidaria = trace.cladeProbability(Cnidaria)
pp_Xenacoelomorpha = trace.cladeProbability(Xenacoelomorpha)
pp_Hemichordata = trace.cladeProbability(Hemichordata)
pp_Echinodermata = trace.cladeProbability(Echinodermata)
pp_Chordata = trace.cladeProbability(Chordata)
pp_Lophotrochozoa = trace.cladeProbability(Lophotrochozoa)
pp_Ecdysozoa = trace.cladeProbability(Ecdysozoa)
## Test the monophyly of metazoans without Ctenophora and Porifera
metazoans_ingroup = clade(Ecdysozoa, Lophotrochozoa, Chordata, Echinodermata, Hemichordata, Xenacoelomorpha, Cnidaria, Placozoa)
pp_metazoans_ingroup = trace.cladeProbability(metazoans_ingroup)
## Test the monophyly of metazoans with Ctenophora and Porifera
metazoans = clade(Ecdysozoa, Lophotrochozoa, Chordata, Echinodermata, Hemichordata, Xenacoelomorpha, Cnidaria, Placozoa, Porifera, Ctenophora)
pp_metazoans = trace.cladeProbability(metazoans)
## Test the Ctenophora-Sister hypothesis
Ctenophora_Sister = clade(Ecdysozoa, Lophotrochozoa, Chordata, Echinodermata, Hemichordata, Xenacoelomorpha, Cnidaria, Placozoa, Porifera)
pp_Ctenophora_Sister = trace.cladeProbability(Ctenophora_Sister)
## Test the Porifera-Sister hypothesis
Porifera_Sister = clade(Ecdysozoa, Lophotrochozoa, Chordata, Echinodermata, Hemichordata, Xenacoelomorpha, Cnidaria, Placozoa, Ctenophora)
pp_Porifera_Sister = trace.cladeProbability(Porifera_Sister)
## Test the Deuterostome hypothesis (or paraphyly, depending on null hypothesis, see below)
Deuterostome = clade(Hemichordata, Echinodermata, Chordata)
pp_Deuterostome = trace.cladeProbability(Deuterostome)
## Test the Xenambulacraria hypothesis
Nephrozoa = clade(Ecdysozoa, Lophotrochozoa, Chordata, Echinodermata, Hemichordata)
Xenambulacraria = clade(Xenacoelomorpha, Hemichordata, Echinodermata)
pp_Nephrozoa = trace.cladeProbability(Nephrozoa)
pp_Xenambulacraria = trace.cladeProbability(Xenambulacraria)
## Test the Coelenterata hypothesis
Coelenterata = clade(Cnidaria, Ctenophora)
pp_Coelenterata = trace.cladeProbability(Coelenterata)
## Now write the results into a file
write(file=OUTPUT_FILENAME,trace.size(),"\n",append=FALSE)
write(file=OUTPUT_FILENAME,pp_Porifera,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Ctenophora,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Cnidaria,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Xenacoelomorpha,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Hemichordata,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Echinodermata,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Lophotrochozoa,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Ecdysozoa,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_metazoans,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_metazoans_ingroup,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Ctenophora_Sister,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Porifera_Sister,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Deuterostome,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Nephrozoa,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Xenambulacraria,"\n",append=TRUE)
write(file=OUTPUT_FILENAME,pp_Coelenterata,"\n",append=TRUE)
## quit
q()