Monday, September 9, 2013

To see if the geometric mean could be used to predict body size, I performed a phylogenetic regression.

Here is the following script used to generate the data for the paper:
(start with body_size.csv).

library (caper)
#make a comparative data frame
data<-comparative.data(phy=tree,data=collapsed.babbler,names.col=matched.tip )

#make two different models; one using maximum likelihood and one using browning motion
mod1<-pgls(log(mean.mass)~log(mean.gm),data, lambda = "ML")
mod2<-pgls(log(mean.mass)~log(mean.gm),data)

#compare the two models using log likelihood
logLik(mod1)
logLik(mod2)

#results of ML model

#results of brownian motion model


#when comparing the two models, compare the LL--since very close, use simpler model (Brownian motion)

Here is the ggplot2 code to generate the scatterplot:
library(ggplot2)
qplot(log(mean.gm),log(mean.mass), data=new,xlab = "Geometric Mean (log)", ylab = "Body Mass (log)",  size = I(7))+theme_bw()+theme(text = element_text(size=30))+ geom_abline(intercept=-1.82487, slope=2.28979)

Monday, September 2, 2013

babblersPCA.csv is the data to be used for analysis for publication

babblers.csv was the file used without any imputation or editing yet 

Monday, August 12, 2013

Going back to the species distribution project, I am now considering not choosing my BIOCLIM layers based on correlation coefficients but rather running a PCA of random samples of points from each layer.

For the Alcippe peracensis extent, the first four principal components will be used:
PC1-->BIO4: Temperature seasonality
PC2-->altitude
PC3-->BIO12: Annual precipitation
PC4-->BIO19:Precipitation of the coldest quarter

I am currently at La Selva biological research station collecting bats for a transcriptomic analysis. The purpose of this post is to record what we did for specimen preparation for future reference:


Lab set up for dissections of tissues:

Checklist:
  • isofluorane (anesthetic)
  • cotton balls
  • RNAlater (lots of bottles!)
  • RNAase AWAY
  • cryovials (lots!)
  • plastic droppers (lots!)
  • boxes
  • disposable vial rack
  • sharpies (small and large)
  • ethanol (to clean)
  • bleach (to clean)
  • plastic petri dishes for dissection
  • dissection kit (I especially used sharp and dull forceps, scissors, small scissors, curved forceps)
  • dissection scope
  • cutting board
  • light source for scope
  • camera and lens
  • tape
  • organized binder with data sheets
  • calipers/ruler and pesolas for measurements
  • sponges
  • bottles for sharps
  • plastic squeezy bottles for RNA later
  • Cooler (for ice)
  • Little bucket for ice
  • Ice
  • Cooler (for dry ice if you want to collect and preserve the whole brain)
  • Dry Ice (for brain)
  • Wing punches
  • scalpels
  • liquid nitrogen and dry shipper
  • tiny plastic baggies to put cryovials in
  • string to hang bat bags on
  • RNA protect
  • DNA away
  • ALL protect
  • aluminum foil
Next time, I will get a voice recorder to record observations during dissection (which will also be useful in the field) and transcribe the notes later on. 
  1. Make sheet before beginning with labeled tissues in respective order
  2. Bleach entire bench and wash down to prevent bleach from contaminating the tissue
  3. Wash everything with RNAase AWAY
  4. Final solution on top of everything with RNAlater
  5. Label cryovials
  6. Remove bat from bat bag
  7. Euthanize with anesthetic
  8. Immediately take measurements and record on data sheet before rigor mortis sets in
  9. Place cotton ball in mouth to prevent rigor mortis
  10. Separate head from body and have two people working on the dissection simultaneously to remove the different tissues
  11. Place tissues in respective labeled tubes
  12. Top off the tubes with RNAlater
  13. Place tubes on ice as soon as tissue is placed in RNAlater
  14. If harvesting the brain, gently scoop the whole brain into a pocket of powdery dry ice and cover with dry ice until the brain hardens. Be careful to not damage the olfactory bulb and keep in mind the strong optic nerve. 
  15. Remove brain from dry ice. We had been placing it in RNA later but not sure if this is the best way to preserve the actual brain. 
  16. Place all cryovials in tiny plastic baggies and label. Place the entire plastic bag in dry shipper. 
  17. Place the remaining specimen in the foil to be frozen at -20 and later prepared in formadehyde
  18. Wash and clean bench in between each specimen. 


The most organized way we found to do this was to perform the actions in the following order
To prepare the final specimen:

  • Plastic sandwich baggies
  • Bottles of formaldehyde 
  • WriteInTheRain paper (will not degrade in formaldehyde)
  • Sharpies
  • Gloves
  • Masks
LANDSAT is a huge fail Downloading 50m resolution PALSAR data for the Vietnam land use project suggested by Dong, et al (2012). I downloaded PALSAR 50m Orthorectified Mosaic Product - Indochina from the following page :http://www.eorc.jaxa.jp/ALOS/en/kc_mosaic/kc_50_indochina.htm

However, I am only using tiles A03, A04, A06, A07, and A09 because I only need the data for Vietnam.

Sunday, June 16, 2013

Okay, so I have a sample of 500 dated trees from my posterior. How to use this method for all of them? 
I haven't come across anyone actually doing this in practice. So let me work out the code right now.

This is probably wrong but I am going to try it:
#declare variables
i<-1
count<-0
final.data.frame<-data.frame()
other.data.frame<-data.frame()

#for every tree in the list
while (i <= length(mini.list)){
  #for that tree run the functions;#10000for testings sake
result<-evol.rate.mcmc(mini.list[[i]],x,ngen=10000,control=list(sd2=2.0,sdlnr=3.0))

  min_split<-minSplit(mini.list[[i]],result$mcmc[21:101,c("node","bp")]) 
  mcmc<-posterior.evolrate(mini.list[[i]],min_split,result$mcmc[21:101,],result$tips[201:1001])
  this.mcmc.data<-as.data.frame(mcmc)
  
#if the nodes produce the same edge, combine the data
if((mini.list[[1]]$edge[this.mcmc.data$node,]==mini.list[[i]]$edge[this.mcmc.data$node,])==TRUE)
  {
    final.data.frame<-rbind(final.data.frame, this.mcmc.data)
  }
#if not put in a different data frame
  else{
    other.data.frame<-rbind(other.data.frame, this.mcmc.data)
    count<-rbind(count, i)
  }
  
  i<-i+1
}
cat("count:", count, "\n")
To plot a nifty plot with the rates:

ave.rates(tree.trim,min_split,extract.clade(tree.trim, node=min_split$node)$tip.label,colMeans(mcmc)["sig1"], colMeans(mcmc)["sig2"],min_split) 

Saturday, June 15, 2013

WOOT--figured out how to date a random sample of my posterior of trees. Super straightforward:

#read in tree file to make multiphylo R object--ignoroing the first 25% as burn-in
tree.list<-read.tree(file.choose(),keep.multi=TRUE, skip = 3374)
#take random sample of trees
sample.list<-sample(tree.list,500)

#date the trees using lapply() and parameters established from Moyle, et al (2012)
newmulti<-lapply(unclass(sample.list),chronopl,lambda=0.5,age.min=c(8.4, 13, 27.1), age.max=c(11.4,17, 37.3), node = c(336, 297, 293))
#rename the object's class
class(newmulti)<-"multiPhylo"

As per Revell's comment on this post: http://comments.gmane.org/gmane.comp.lang.r.phylo/2773

Friday, June 14, 2013

Finally have some meaningful(?) output. This is going to be lengthy but it will be really helpful in organizing my thought. Because I had been including morphological data for very "distantly-related babblers", all of the nodes for the rate shifts occurred, of course on the root node. I removed some more taxa to only include measurements for birds within the babbler family. Therefore, I reran the PCA. Still similar trends:

Count the number of species in my analysis and the number of specimens per species.
> babblers.species.count<-ddply(.data=new.babblers, .(Genus, Species), species.count=length(Wing.Length), summarize)
> babblers.species.count
             Genus         Species species.count
1       Actinodura        egertoni             6
2       Actinodura    morrisoniana             3
3       Actinodura      nipalensis             3
4       Actinodura         waldeni             5
5          Alcippe    brunneicauda             4
6          Alcippe          grotei            11
7          Alcippe      morrisonia             5
8          Alcippe      peracensis             2
9          Alcippe    poioicephala            11
10           Babax     lanceolatus             4
11           Cutia      nipalensis             7
12         Dumetia      hyperythra             6
13  Gampsorhynchus         rufulus             8
14        Garrulax         affinis            11
15        Garrulax     albogularis             5
16        Garrulax         canorus            14
17        Garrulax       chinensis            40
18        Garrulax      cineraceus             4
19        Garrulax erythrocephalus             1
20        Garrulax        formosus             4
21        Garrulax     leucolophus            41
22        Garrulax        lineatus             6
23        Garrulax           maesi             3
24        Garrulax       merulinus             1
25        Garrulax         milleti             1
26        Garrulax          milnei             4
27        Garrulax        mitratus             6
28        Garrulax       monileger             8
29        Garrulax        nuchalis             3
30        Garrulax       ocellatus             7
31        Garrulax       palliatus             4
32        Garrulax poecilorhynchus             4
33        Garrulax          sannio            12
34        Garrulax       squamatus             3
35        Garrulax        striatus             6
36        Garrulax     subunicolor             6
37        Garrulax      variegatum             3
38        Garrulax        virgatus             2
39    Heterophasia     auricularis             3
40    Heterophasia      capistrata             6
41    Heterophasia        gracilis             4
42    Heterophasia    maelanoleuca             9
43    Heterophasia       picaoides             7
44    Heterophasia       pulchella             3
45     Illadopsis       albipectus             4
46     Illadopsis         cleaveri             4
47     Illadopsis       fulvescens             8
48     Illadopsis           puveli             2
49     Illadopsis      pyrrhoptera             2
50     Illadopsis        rufescens             3
51     Illadopsis       rufipennis             6
52         Kenopia         striata             2
53       Leiothrix     argentauris            15
54       Leiothrix           lutea             4
55       Liocichla       phoenicea             1
56       Liocichla          steeri             3
57        Macronus         gularis            25
58        Macronus        ptilosus             3
59        Macronus     striaticeps            12
60    Malacocincla         abbotti             4
61    Malacocincla     cinereiceps             2
62    Malacocincla     malaccensis             6
63    Malacocincla        sepiaria             3
64    Malacopteron          affine             3
65    Malacopteron      albogulare             4
66    Malacopteron        cinereum             1
67    Malacopteron     magnirostre             7
68    Malacopteron          magnum             5
69    Malacopteron    palawanenese             3
70           Minla      ignotincta             6
71           Minla        strigula             7
72       Napothera    brevicaudata             4
73       Napothera          crassa             2
74       Napothera     crispifrons             7
75       Napothera     epilepidota            10
76      Pellorneum      albiventre             9
77      Pellorneum    capristratum            15
78      Pellorneum      pyrrogenys             6
79      Pellorneum        ruficeps             5
80      Pellorneum        tickelli             1
81     Phyllanthus      atripennis             6
82    Pomatorhinus    erythrogenys            18
83    Pomatorhinus    ferruginosus             4
84    Pomatorhinus      horsfieldi             6
85    Pomatorhinus      hypoleucos             4
86    Pomatorhinus        montanus             4
87    Pomatorhinus      ruficollis             9
88    Pomatorhinus        swinhoei             2
89     Ptilocichla         falcata             2
90     Ptilocichla   leucogrammica             3
91     Ptilocichla     mindanensis             3
92       Ptyrticus        turdinus             6
93    Rhopocichla         atriceps             4
94         Rimator    malacoptilus             3
95    Schoeniparus         brunnea            10
96    Schoeniparus     castaneceps             1
97    Schoeniparus         cinerea             3
98    Schoeniparus           dubia             6
99    Schoeniparus     rufogularis             9
100    Spelaeornis    chocolatinus             5
101    Spelaeornis   troglodytoide             3
102   Sphenocichla           humei             3
103      Stachyris         ambigua             4
104      Stachyris        chrysaea            11
105      Stachyris    erythroptera             8
106      Stachyris      grammiceps             2
107      Stachyris        leucotis             3
108      Stachyris        maculata             3
109      Stachyris    melanothorax             6
110      Stachyris       nigriceps             7
111      Stachyris     nigricollis             4
112      Stachyris           oglei             3
113      Stachyris    poliocephala             4
114      Stachyris        pyrrhops             3
115      Stachyris        ruficeps            10
116      Stachyris       rufifrons             3
117      Stachyris       striolata             5
118      Stachyris       thoracica             4
119        Timalia         pileata             8
120    Trichastoma         bicolor             4
121    Trichastoma       celebense            11
122    Trichastoma       rostratum             5
123      Turdoides         bicolor             3
124      Turdoides         gularis             4
125      Turdoides       jardineii            10
126      Turdoides        plebejus            16
127      Turdoides      reinwardii             3
128  Xiphirhynchus   superciliaris             3

I had to match up the names of my species with the names on the tips (huge pain) and I ended up doing it by hand--I could have probably scripted it but I think it would have taken me just as long. 

I collapsed the PCA values into means for each species but I believe I already posted the script for that. So now I have a .csv that contains all the mean PCA that corresponds to each species with the species name that matches the tip. YAY.

Now I finally have evol.rate.mcmc() running and I ran the output for the first four principle components:
#read in the PCA vector:
temp<-read.csv(file.choose())
x<-as.vector(temp$Tarsus) #vector of tarsus measurement traits

names(x)<-as.vector(temp$Species) #names the measurements
#i have to date the trees
undated.tree<-read.tree(file.choose())
#the tree was based on two calibration points (see Moyle, et al. 2012) and the root

tree<-chronopl(undated.tree, lambda=.5, age.min=c(8.4, 13, 27.1), age.max=c(11.4,17, 37.3), node = c(336, 297, 293)) 

#now time for the juicy part--using default parameters to start with as suggested by Revell, et al. 2011

result<-evol.rate.mcmc(tree,x,ngen=100000,control=list(sd2=2.0,sdlnr=3.0))
min_split<-minSplit(tree,result$mcmc[21:101,c("node","bp")]) 
mcmc<-posterior.evolrate(tree,min_split,result$mcmc[21:101,],result$tips[201:1001])
mcmc.data.1<-as.data.frame(mcmc)
#nicely summarize the mcmc for Tarsus
output.summary.1<-ddply(.data=mcmc.data.1, .(node), node.count=length(bp), mean.bp=mean(bp),mean.sig1=mean(sig1), mean.sig2=mean(sig2), summarize)

#let's visualize the nodes! first trim the tree to visualize only nodes that were used in analysis
trim.tree<-drop.tip(tree,tree$tip.label[-(match(names(x),tree$tip.label))]) 

Nodes 3 and 4 are probably root nodes and have high probability. Let's look at 186, 185, and 188.

#now extract the clades we want to look at

this.extract<-extract.clade(tree.trim, tree.trim$edge[186])
this.extract<-extract.clade(tree.trim, tree.trim$edge[58])
#perhaps the most interesting because it contains nearly the whole subfamily
this.extract<-extract.clade(tree.trim, tree.trim$edge[4])
Now time for the bill length:
#read in the PCA vector:
temp<-read.csv(file.choose())
x<-as.vector(temp$Bill.Length) #vector of bill length measurement traits

#you have to reset the names
names(x)<-as.vector(temp$Species) #names the measurements

#use the same dated tree as before
result<-evol.rate.mcmc(tree,x,ngen=100000,control=list(sd2=2.0,sdlnr=3.0))
min_split<-minSplit(tree,result$mcmc[21:101,c("node","bp")]) 
mcmc<-posterior.evolrate(tree,min_split,result$mcmc[21:101,],result$tips[201:1001])
mcmc.data.2<-as.data.frame(mcmc)
#nicely summarize the mcmc for Bill Length
output.summary.2<-ddply(.data=mcmc.data.2, .(node), node.count=length(bp), mean.bp=mean(bp),mean.sig1=mean(sig1), mean.sig2=mean(sig2), summarize)

Hmm...very strange. 250 and 251 are only tips of two sister taxa.
> this.extract<-extract.clade(tree.trim, tree.trim$edge[122])
> plot(this.extract)

Now for the Hallux:
x<-as.vector(temp$Hallux) #vector of hallux measurement traits

#you have to reset the names
names(x)<-as.vector(temp$Species) #names the measurements

#use the same dated tree as before
result<-evol.rate.mcmc(tree,x,ngen=100000,control=list(sd2=2.0,sdlnr=3.0))
min_split<-minSplit(tree,result$mcmc[21:101,c("node","bp")]) 
mcmc<-posterior.evolrate(tree,min_split,result$mcmc[21:101,],result$tips[201:1001])
mcmc.data.2<-as.data.frame(mcmc)
#nicely summarize the mcmc for hallux
output.summary.2<-ddply(.data=mcmc.data.3, .(node), node.count=length(bp), mean.bp=mean(bp),mean.sig1=mean(sig1), mean.sig2=mean(sig2), summarize)

Node 9 has high probability (low frequency):
> this.extract<-extract.clade(tree.trim, tree.trim$edge[9])
> plot(this.extract)

> this.extract<-extract.clade(tree.trim, tree.trim$edge[29])
> plot(this.extract)

> this.extract<-extract.clade(tree.trim, tree.trim$edge[235])
> plot(this.extract)

Enough of that, let's finally look at Tail Length and see if there is anything interesting:
x<-as.vector(temp$Tail.Length) #vector of tail measurement traits

#you have to reset the names
names(x)<-as.vector(temp$Species) #names the measurements

#use the same dated tree as before
result<-evol.rate.mcmc(tree,x,ngen=100000,control=list(sd2=2.0,sdlnr=3.0))
min_split<-minSplit(tree,result$mcmc[21:101,c("node","bp")]) 
mcmc<-posterior.evolrate(tree,min_split,result$mcmc[21:101,],result$tips[201:1001])
mcmc.data.2<-as.data.frame(mcmc)
#nicely summarize the mcmc for tail length
output.summary.4<-ddply(.data=mcmc.data.4, .(node), node.count=length(bp), mean.bp=mean(bp),mean.sig1=mean(sig1), mean.sig2=mean(sig2), summarize)



> this.extract<-extract.clade(tree.trim, tree.trim$edge[8])
> plot(this.extract)

Okay I am kind of petering out on uploading all of these. Gives me something to compare to though.

Important questions:
Next step is to start making my slides for my presentation and date random samples of trees and get feedback from adviser to make sure I am on the right track. 

Thursday, June 6, 2013

Okay checked all my data and things were looking fishy and realized there were still some typos (missing a decimal point, etc.). Everything was carefully reexamined and I think (fingers crossed) these are official PCA results.

There are 145 unique species of babblers to be used in my analysis.

> library(plyr)
> babblers.species.count<-ddply(.data=new.babblers, .(Genus, Species), species.count=length(Wing.Length), summarize)
> babblers.species.count

                   Genus         Species species.count
1             Actinodura        egertoni             6
2             Actinodura    morrisoniana             3
3             Actinodura      nipalensis             3
4             Actinodura         waldeni             5
5                Alcippe    brunneicauda             4
6                Alcippe          grotei            11
7                Alcippe      morrisonia             5
8                Alcippe      peracensis             2
9                Alcippe    poioicephala            11
10               Alcippe     rufogularis             5
11                 Babax      lanceolata             4
12             Chrysomma         sinense             6
13  Chrysomma (Moupinia)      poecilotis             4
14                 Cutia      nipalensis             7
15              Dumentia      hyperythra             6
16        Gampsorhynchus         rufulus             8
17              Garrulax         affinis            11
18              Garrulax     albogularis             5
19              Garrulax         canorus            14
20              Garrulax       chinensis            40
21              Garrulax      cineraceus             4
22              Garrulax erythrocephalus             1
23              Garrulax        formosus             4
24              Garrulax     leucolophus            41
25              Garrulax        lineatus             6
26              Garrulax           maesi             3
27              Garrulax       merulinus             1
28              Garrulax         milleti             1
29              Garrulax          milnei             4
30              Garrulax        mitratus             6
31              Garrulax       monileger             8
32              Garrulax        nuchalis             3
33              Garrulax       ocellatus             7
34              Garrulax       palliatus             4
35              Garrulax      pectoralis             1
36              Garrulax  perspicillatus            14
37              Garrulax poecilorhynchus             4
38              Garrulax          sannio            12
39              Garrulax         sinesis             9
40              Garrulax       squamatus             3
41              Garrulax        striatus             6
42              Garrulax     subunicolor             6
43              Garrulax      variegatum             3
44              Garrulax         vassali             7
45              Garrulax        virgatus             2
46         Heterophaesia       annectens             4
47         Heterophaesia     auricularis             3
48         Heterophaesia      capistrata             6
49         Heterophaesia      desgodinsi             1
50         Heterophaesia        gracilis             4
51         Heterophaesia    maelanoleuca             9
52         Heterophaesia       picaoides             7
53         Heterophaesia       pulchella             3
54          Heterophasia      desgodinsi             7
55           Illadopsis       albipectus             4
56           Illadopsis         cleaveri             4
57           Illadopsis       fulvescens             8
58           Illadopsis           puveli             2
59           Illadopsis      pyrrhoptera             2
60           Illadopsis        rufescens             3
61           Illadopsis       rufipennis             6
62               Kenopia         striata             2
63             Leiothrix     argentauris            15
64             Leiothrix           lutea             4
65            Leonardina           woodi             1
66             Liocichla       phoenicea             1
67             Liocichla         steerii             3
68              Macronus         gularis            25
69              Macronus        ptilosus             3
70              Macronus     striaticeps            12
71          Malacocincla         abbotti             4
72          Malacocincla     cinereiceps             2
73          Malacocincla     malaccensis             6
74          Malacocincla        sepiaria             3
75          Malacopteron          affine             3
76          Malacopteron      albogulare             4
77          Malacopteron        cinereum             1
78          Malacopteron     magnirostre             7
79          Malacopteron          magnum             5
80          Malacopteron    palawanenese             3
81                 Malia           grata             6
82         Micromacronus        sordidus             2
83                 Minla      ignotincta             6
84                 Minla        strigula             7
85             Napothera    brevicaudata             4
86             Napothera          crassa             2
87             Napothera     crispifrons             7
88             Napothera     epilepidota            10
89            Pellorneum      albiventre             9
90            Pellorneum    capristratum            15
91            Pellorneum      pyrrogenys             6
92            Pellorneum        ruficeps             5
93            Pellorneum        tickelli             1
94           Phyllanthus      atripennis             6
95              Pnoepyga         pusilla             1
96          Pomatorhinus    erythrogenys            18
97          Pomatorhinus    ferruginosus             4
98          Pomatorhinus      horsfieldi             6
99          Pomatorhinus      hypoleucos             4
100         Pomatorhinus        montanus             4
101         Pomatorhinus      ruficollis             9
102         Pomatorhinus        swinhoei             2
103          Ptilocichla         falcata             2
104          Ptilocichla   leucogrammica             3
105          Ptilocichla     mindanensis             3
106            Ptyrticus        turdinus             6
107         Rhopocichla         atriceps             4
108              Rimator    malacoptilus             3
109         Schoeniparus         brunnea            10
110         Schoeniparus     castaneceps             1
111         Schoeniparus         cinerea             3
112         Schoeniparus           dubia             6
113         Schoeniparus     rufogularis             4
114          Spelaeornis    chocolatinus             5
115          Spelaeornis        formosus             2
116          Spelaeornis   troglodytoide             3
117         Sphenocichla           humei             3
118            Stachyris         ambigua             4
119            Stachyris        chrysaea            11
120            Stachyris    erythroptera             8
121            Stachyris      grammiceps             2
122            Stachyris        leucotis             3
123            Stachyris        maculata             3
124            Stachyris    melanothorax             6
125            Stachyris       nigriceps             7
126            Stachyris     nigricollis             4
127            Stachyris           oglei             3
128            Stachyris    poliocephala             4
129            Stachyris        pyrrhops             3
130            Stachyris        ruficeps            10
131            Stachyris       rufifrons             3
132            Stachyris       striolata             5
133            Stachyris       thoracica             4
134              Timalia         pileata             8
135          Trichastoma         bicolor             4
136          Trichastoma       celebense            11
137          Trichastoma       rostratum             5
138          Trichastoma        tickelli            10
139          Trichastona        tickelli             2
140            Turdoides         bicolor             3
141            Turdoides         gularis             4
142            Turdoides       jardineii            10
143            Turdoides        plebejus            16
144            Turdoides      reinwardii             3
145        Xiphirhynchus   superciliaris             3