Pelagia.Rmd

---
title: "Pelagia Analysis"
author: "David Cerny"
date: "10/19/2017"
output: html_document
---

# Protocol summary

**Initial filtering step:** Exclude the 20% of loci with the shortest tree length from the 615 UCE pool and all of the loci in the top 20% of RF distance from the concatenated tree

**Track 1:** Run PartitionFinder on the sites remaining in the filtered UCE pool using *k*-means, then SortaDate the resulting partitions.

**Track 2:** Break up the loci into 150-bp chunks and run SortaDate on these. Rank the chunks by bipartition support, and keep adding the highest-support chunks until you reach the same total number of sites that is included in the alignment resulting from Track 3 (see below). Find out what value of bipartition support you end up at when this threshold is reached.

**Track 3:** Run SortaDate on intact loci. Concatenate those with bipartition support greater than 0.5, then run PartitionFinder on the resulting alignment (iteratively if needed) until you get 10 partitions or fewer. Run a PhyloBayes autocorrelation test on each of the selected loci.

**Subtrack 3.1:** Run the PhyloBayes model comparison (autocorrelation test) on each of the 138 loci whose gene trees have a bipartition support score of 0.5 or better.

**Subtrack 3.2:** See if the concatenated 138-locus analysis will start in BEAST.

**Subtrack 3.3:** In order to be able to perform a partitioned MCMCTree analysis, SortaDate the 21 partitions into which the 138 loci were grouped by PartitionFinder, and use the top 10.

# Initial filtering step

SortaDate was used to calculate the length of individual gene trees:

Step 1: convert the alignment files from Nexus to FASTA (Geneious -> Batch export).

Step 2: generate a script to perform a RAxML analysis with SH-like support values on each of the 615 loci:

```{r, eval = FALSE}
ucelist <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/ucelist.txt")
n <- length(unlist(ucelist))

# 1st part of the command: what to analyze

raxmlcommands <- vector(mode="character", length = n)
for(i in ucelist) {
  raxmlcommands[i] <- paste("ete3 build -w standard_raxml -n /Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/", i, sep = "")
}

# 2nd part of the command: names to give to the output files

withendings <- as.character(as.vector(as.matrix(ucelist)))
noendings <- vector()
for(i in withendings) {
  noendings <- append(noendings, sapply(strsplit(i, split='.', fixed=TRUE), function(x) (x[1])))
}
outputnames <- vector()
for(i in noendings) {
  outputnames <- append(outputnames, paste(" -o /Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/", i, sep=""))
}

length(outputnames)

# Put it together and print it to file:

commandlist <- paste(raxmlcommands, outputnames, "-tree &&", sep = "")
write(commandlist, "/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/615-gene-tree-analysis.sh")
```

Step 3: run the script as follows:

```
export PATH=~/anaconda_ete/bin:$PATH
chmod 755 /Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/615-gene-tree-analysis.sh
/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/615-gene-tree-analysis.sh
```

Step 4: move the resulting tree files from their nested subdirectories to the main directory:

```
find . -name '*.nw' -exec mv {} . \;
ls -lR ./*.nw | wc -l                     # Make sure you got all 615
```

Step 5: rename the trees to give them the `.tre` file ending that SortaDate will be looking for:

```
find . -name "*.fasta.final_tree.nw" -exec rename 's/.fasta.final_tree.nw$/.tre/' {} \;
```

Step 6: get_var_length

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_var_length.py /Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-loci --outg ceratoscopelus_warmingii,gadus_morhua
```

Two different outgroups were selected to increase the probability that at least one of them would be present in each alignment. Note that the script seems to be sensitive to the order in which the outgroups are passed to it: `--outg gadus_morhua,ceratoscopelus_warmingii` makes the program freeze and yields an empty output file, while `--outg ceratoscopelus_warmingii,gadus_morhua` works as intended, with values successfully estimated for all 615 loci (even those in which one of the outgroups was not present).

Step 7: get_bp_genetrees

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_bp_genetrees.py /Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete /Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-loci
```

Here, too, the congruence with the reference tree was successfully estimated for all 615 gene trees.

Step 8: combine_results

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/combine_results.py /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-loci /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-loci --outf /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-loci
```

Step 9: get Robinson-Foulds distances of the individual gene trees from the reference tree:

```
import os, uuid
from ete3 import Tree

t2 = Tree("/Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre", format = 9)
for file in os.listdir("/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/"):
    if file.endswith(".tre"):
        filename = "/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/" + file
        t1 = Tree(filename, format = 2)
        try:
            rf = t1.robinson_foulds(t2)
            print str(file), (rf[0])
        except:
            pass
```

```
python robinsonfoulds.py > /Users/David/Grive/Alfaro_Lab/Pelagia/rf-pelagia-loci
```

```{r, echo = FALSE}
sortadate_results <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-loci")
colnames(sortadate_results) <- c("Locus", "Root_to_tip_variance", "Tree_length", "Bipartition_support")
rf <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/rf-pelagia-loci")
colnames(rf) <- c("Locus", "Robinson_Foulds_distance")
par(mfrow=c(1,2), oma = c(0, 0, 2, 0))
hist(sortadate_results$Tree_length, breaks = 50, xlab = "Tree length (assumed proportional to rate)", main = "Distribution of gene tree lengths", col = "greenyellow")
hist(rf$Robinson_Foulds_distance, breaks = 50, xlab = "RF distance from the concatenated tree", main = "Distribution of RF distances", col = "lightcoral")
mtext("615 fully resolved gene trees", outer = TRUE, cex = 1.5)
```

Step 10: exclude the 20% of loci with the shortest tree length from the 615-UCE pool as well as all of the loci in the top 20% of RF distance from the concatenated tree:

```{r}
tree_quantile <- quantile(sortadate_results$Tree_length, probs = seq(0, 1, .2))
RF_quantile <- quantile(rf$Robinson_Foulds_distance, probs = seq(0, 1, .2))
# Merge the comb and RF files:
comb_collapsed <- merge(sortadate_results, rf, by = "Locus")
filtered_loci <- comb_collapsed[comb_collapsed$Tree_length > tree_quantile[2] & comb_collapsed$Robinson_Foulds_distance < RF_quantile[5],]
```
```{r, eval = FALSE}
write.table(filtered_loci, "/Users/David/Grive/Alfaro_Lab/Pelagia/filtered_loci.txt", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")
```

These selection criteria yielded a filtered dataset of 440 loci (71.5% of the original number), suggesting that the nearly invariant loci with the shortest tree length also performed the worst in terms of topology estimation.

Step 11: create a directory with just the selected loci:

Step 11.1: get a list of the names of the corresponding tree files:

```{r, eval = FALSE}
write.table(filtered_loci$Locus, "/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/filtered_loci_names.txt", row.names = FALSE, col.names = FALSE, quote = FALSE)
```

Step 11.2: change the endings so that the file lists alignments rather than trees:

```
gsed -i 's/tre/fasta/g' filtered_loci_names.txt
```

Step 11.3: copy the files listed to a new directory:

```
for file in `cat filtered_loci_names.txt`; do cp "$file" /Users/David/Grive/Alfaro_Lab/Pelagia/Filtered/ ; done
```

Step 12: subset the `comb-pelagia-loci` so that it contains only the 440 selected loci:

```{r, eval = FALSE}
comb_filtered <- sortadate_results[sortadate_results$Locus %in% filtered_loci$Locus,]
write.table(comb_filtered, "/Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-filtered-loci", row.names = FALSE, col.names = FALSE, quote = FALSE, sep = "\t")
```





# Track 1: PartitionFinder on sites

Step 1: change the alignment ending to make the files recognizable to FASconCAT:

```
find . -name "*.fasta" -exec rename 's/.fasta$/.fas/' {} \;
```

Step 2: concatenate the 440 loci using FASconCAT (make sure the Perl script is in the same directory as the FASTA files!):

```
perl FASconCAT_v1.1.pl -s
```

Step 3: convert the concatenated alignment from FASTA to PhyML to make it readable by PartitionFinder:

```
perl catfasta2phyml-master/catfasta2phyml.pl --verbose Pelagia/Filtered/FcC_smatrix.fas > pelagia-filtered-concat.phy
```

Step 4: run PartitionFinder on the following configuration file using the command below:

```{r comment='', echo = FALSE}
cat(readLines("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/partition_finder.cfg"), sep = '\n')
```

```
python PartitionFinder.py /home/analysis/partitionfinder-new/pelagia-filtered-loci/partition_finder.cfg --verbose --raxml --save-phylofiles
```

Step 5: extract the list of partitions from `best_scheme.txt` and read it into R:

```{r}
partitions <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/partitions.txt", sep = "|", header = TRUE, stringsAsFactors = FALSE, strip.white = TRUE)
```

Step 6: prepare a list of alignments and information files to download:

```{r}
alignments <- paste(gsub(" ", "", partitions[,4], fixed = TRUE), ".phy", sep = "")
infofiles <- paste("RAxML_info.", gsub(" ", "", partitions[,4], fixed = TRUE), "_GTR+G.txt", sep = "")
```
```{r, eval = FALSE}
write(cbind(alignments, infofiles), "/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/whattorsync.txt")
```

Step 7: download the relevant files as follows:

```
rsync --files-from=/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions/whattorsync.txt analysis@azathoth.eeb.ucla.edu:~/partitionfinder-new/pelagia-filtered-loci/analysis/phylofiles /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions
```

State 8: extract the branch length multipliers from all the info files:

```
for file in `cat whattorsync.txt`; do gsed -n -e 's/^.*Branch length scaler: //p' "$file"; done > multipliers.txt
```

State 9: add the rate info to your initial partition table:

```{r}
multipliers <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions/multipliers.txt")
partitions <- cbind(partitions, multipliers)
```

Step 10: plot partition rate against partition length:

```{r}
plot(partitions[,3], partitions[,6], xlab = "Partition length (bp)", ylab = "Relative rate", pch = 19, main = "All 37 partitions")
plot(partitions[,3], partitions[,6], xlim = c(0, 20000), xlab = "Partition length (bp)", ylab = "Relative rate", pch = 19, main = "(Longest partition excluded)")
```

```{r, echo = FALSE, eval = FALSE}
png("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/ratevslength.png", width = 1600, height = 2800, pointsize = 48)
par(mfrow=c(2,1))
plot(partitions[,3], partitions[,6], xlab = "Partition length (bp)", ylab = "Relative rate", pch = 19, main = "All 37 partitions")
plot(partitions[,3], partitions[,6], xlim = c(0, 20000), xlab = "Partition length (bp)", ylab = "Relative rate", pch = 19, main = "(Longest partition excluded)")
```

Step 11: convert the alignment files from Phylip to FASTA (Geneious -> Batch export) and simplify the file endings of the resulting files:

```
find . -name "*.phy.fasta" -exec rename 's/.phy.fasta$/.fasta/' {} \;
```

Step 12: run RAxML on each of the 37 partitions to find the maximum likelihood tree and its SH-like support values using the following script (`runthruraxml.R`):

```{r}
args <- commandArgs(trailingOnly = TRUE)

commandpart1 <- paste("ete3 build -w standard_raxml -n /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions/", args, sep = "")
commandpart2 <- paste(" -o /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions/", args, sep="")
wholecommand <- paste(commandpart1, commandpart2, "-tree", sep = "")
write(wholecommand, "/Users/David/Grive/Alfaro_Lab/Pelagia/37-pfinder-partition-analysis.sh")
```

```
cd /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions
Rscript runthruraxml.R *.fasta
export PATH=~/anaconda_ete/bin:$PATH
chmod 755 /Users/David/Grive/Alfaro_Lab/Pelagia/37-pfinder-partition-analysis.sh
/Users/David/Grive/Alfaro_Lab/Pelagia/37-pfinder-partition-analysis.sh
```

RAxML tree inference failed for the largest and nearly invariant 123-kbp partition.

Step 13: move and rename the tree files.

```
find . -name '*.nw' -exec mv {} . \;
find . -name "*.fasta.final_tree.nw" -exec rename 's/.fasta.final_tree.nw$/.tre/' {} \;
```

Step 14: SortaDate the 36 partitions for which RAxML trees are available.

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_var_length.py /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-kmeans --outg ceratoscopelus_warmingii,gadus_morhua
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_bp_genetrees.py /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/k-means/Partitions /Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-kmeans
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/combine_results.py /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-kmeans /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-kmeans --outf /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-kmeans
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_good_genes.py /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-kmeans --max 36 --order 3,1,2 --outf /Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-kmeans
```

Step 15: merge the SortaDate output with the list of partition sizes and branch rate multipliers:

```{r}
kmeans_parts <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-kmeans", header = TRUE, stringsAsFactors = FALSE)
colnames(kmeans_parts)
colnames(partitions)[4] <- "name"
colnames(partitions)[6] <- "multiplier"
kmeans_parts <- merge(kmeans_parts, partitions, by = "name")
```

Step 16: does tree length (estimated by SortaDate) correlate with branch rate multipliers (estimated by PartitionFinder)?

```{r}
library(ggplot2)

lm_eqn <- function(df){
    m <- lm(df[,3] ~ df[,9]);
    eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
         list(a = format(coef(m)[1], digits = 2), 
              b = format(coef(m)[2], digits = 2), 
             r2 = format(summary(m)$r.squared, digits = 3)))
    as.character(as.expression(eq));                 
}
# plot(kmeans_parts$treelength, kmeans_parts$multiplier, pch = 19)
# abline(lm(kmeans_parts$multiplier ~ kmeans_parts$treelength))

ggplot(kmeans_parts, aes(x=treelength, y = multiplier)) + geom_point(color = 'black') + xlab("Tree length (assumed proportional to rate)") +  ylab(expression("Branch rate multiplier")) + geom_smooth(method = 'lm', size = 0.5, color = 'black') + annotate("text", x = 30, y = 5, label = lm_eqn(kmeans_parts), parse = TRUE)
```

Step 17: bivariate plots showing the relationships among the three variables:

```{r, echo = FALSE}
par(oma = c(0, 0, 2, 0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(kmeans_parts$bipartition, kmeans_parts$treelength, xlab = "Bipartition support", ylab = "Tree length", main = "Tree length vs. bipartition support", col = "lightcoral")
abline(lm(kmeans_parts$treelength ~ kmeans_parts$bipartition))
plot(kmeans_parts$bipartition, kmeans_parts$root.to.tip_var, xlab = "Bipartition support", ylab = "Root-to-tip variance", main = "Clockiness vs. bipartition support", col = "greenyellow")
abline(lm(kmeans_parts$root.to.tip_var ~ kmeans_parts$bipartition))
plot(kmeans_parts$root.to.tip_var, kmeans_parts$treelength, xlab = "Root-to-tip variance", ylab = "Tree length", main = "Tree length vs. clockiness", col = "cornflowerblue")
abline(lm(kmeans_parts$treelength ~ kmeans_parts$root.to.tip_var))
mtext("36 k-means partitions", outer = TRUE, cex = 1.5)
```

Step 18: if you take the 10 partitions with the highest bipartition support, what will be the lowest bipartition score still included in the subset?

```{r}
kmeans_parts <- kmeans_parts[order(-kmeans_parts$bipartition),]
kmeans_parts$bipartition[10]
```

Step 19: if you keep adding partitions until the total length reaches that of the 138-locus alignment created in track 3 (see below), what will be the lowest bipartition score still included in the subset?

```{r}
totallength <- 0
for (i in 1:nrow(kmeans_parts)) {
  totallength <- totallength + kmeans_parts$No_sites[i]
  if (totallength >= 114462) {                         # see step 16 of track 3
    print(kmeans_parts$bipartition[i])
    break
  }
}
```





# Track 2: SortaDate on 150-bp chunks

Step 1: split the FASTA files of individual loci using a custom R script:

```
for file in `cat filtered_loci_names.txt`; do Rscript /Users/David/Grive/Alfaro_Lab/Pelagia/newandimprovedsplitter.R "$file" ; done
```

Step 2: add a file extension to the chunk files:

```
find . -type f -exec mv '{}' '{}'.fasta \;
```

### Fully resolved trees

Step 3: Generate a script to perform a RAxML analysis with SH-like support values on each of the 2397 chunks.

```{r, echo = FALSE, eval = FALSE}
chunklist <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/Chunks/chunklist.txt")
n <- length(unlist(chunklist))

# 1st part of the command: what to analyze

raxmlcommands <- vector(mode="character", length = n)
for(i in chunklist) {
  raxmlcommands[i] <- paste("ete3 build -w standard_raxml -n /Users/David/Grive/Alfaro_Lab/Pelagia/Chunks/", i, sep = "")
}

# 2nd part of the command: names to give to the output files

withendings <- as.character(as.vector(as.matrix(chunklist)))
noendings <- vector()
for(i in withendings) {
  noendings <- append(noendings, sapply(strsplit(i, split='.', fixed=TRUE), function(x) (x[1])))
}
outputnames <- vector()
for(i in noendings) {
  outputnames <- append(outputnames, paste(" -o /Users/David/Grive/Alfaro_Lab/Pelagia/Chunks/", i, sep=""))
}

length(outputnames)

# Put it together and print it to file:

commandlist <- paste(raxmlcommands, outputnames, "-tree", sep = "")
write(commandlist, "/Users/David/Grive/Alfaro_Lab/Pelagia/2397-chunk-analysis.sh")
```

The analysis was successful for 1180 chunks, corresponding to a success rate of 49.2%.

Step 4: move and rename the tree files.

```
find . -name '*.nw' -exec mv {} . \;
ls -lR ./*.nw | wc -l 
find . -name "*.fasta.final_tree.nw" -exec rename 's/.fasta.final_tree.nw$/.tre/' {} \;
```

Step 5: get_var_length

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_var_length.py /Users/David/Grive/Alfaro_Lab/Pelagia/Chunks --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-chunks --outg ceratoscopelus_warmingii,gadus_morhua
```

Step 6: get_bp_genetrees

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_bp_genetrees.py /Users/David/Grive/Alfaro_Lab/Pelagia/Chunks /Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-chunks
```

Step 7: combine_results

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/combine_results.py /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-chunks /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-chunks --outf /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-chunks
```

Step 8: get Robinson-Foulds distances from the ExaML tree

```
python robinsonfoulds.py > rf-pelagia-chunks
```

Step 9: run the `get_good_genes` script with default weights assigned to the three criteria (i.e., in order of descending priority: bipartition support, root-to-tip variance, tree length) to sort the chunks from best to worst:

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_good_genes.py /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-chunks --max 1180 --order 3,1,2 --outf /Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-chunks
```

Step 10: make bivariate plots to show the relationships between the three key variables:

```{r, echo = FALSE}
gg_pelagia_chunks <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-chunks", header = TRUE)
par(oma = c(0, 0, 2, 0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(gg_pelagia_chunks$bipartition, gg_pelagia_chunks$treelength, xlab = "Bipartition support", ylab = "Tree length", main = "Tree length vs. bipartition support", col = "lightcoral")
abline(lm(gg_pelagia_chunks$treelength ~ gg_pelagia_chunks$bipartition))
plot(gg_pelagia_chunks$bipartition, gg_pelagia_chunks$root.to.tip_var, xlab = "Bipartition support", ylab = "Root-to-tip variance", main = "Clockiness vs. bipartition support", col = "greenyellow")
abline(lm(gg_pelagia_chunks$root.to.tip_var ~ gg_pelagia_chunks$bipartition))
plot(gg_pelagia_chunks$root.to.tip_var, gg_pelagia_chunks$treelength, xlab = "Root-to-tip variance", ylab = "Tree length", main = "Tree length vs. clockiness", col = "cornflowerblue")
abline(lm(gg_pelagia_chunks$treelength ~ gg_pelagia_chunks$root.to.tip_var))
mtext("1180 chunks (fully resolved trees)", outer = TRUE, cex = 1.5)
```

### Collapsed trees

Step 11: collapse nodes with SH-like support values of less than 0.90 using the following Python script:

```
import os, uuid
from ete3 import Tree

for file in os.listdir("/Users/David/Grive/Alfaro_Lab/Pelagia/Chunks"):
    if file.endswith(".tre"):
        filename = "/Users/David/Grive/Alfaro_Lab/Pelagia/Chunks/" + file
        outname = "/Users/David/Grive/Alfaro_Lab/Pelagia/Min_SH-like_90_chunk_trees/" + file
        t = Tree(filename, format=2)

        print t.get_ascii(attributes=['support', 'name'])

        for node in t.get_descendants():
            if not node.is_leaf() and node.support <= 0.9:
                node.delete()

        print t.get_ascii(attributes=['support', 'name'])

        t.write(format=0, outfile=outname)
```

Step 12: Copy the locus alignments into the new folder:

```
cp Chunks/*.fasta Min_SH-like_90_chunk_trees
```

Step 13: Re-run SortaDate on the collapsed trees:

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_var_length.py /Users/David/Grive/Alfaro_Lab/Pelagia/Min_SH-like_90_chunk_trees --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-collapsed-chunks --outg ceratoscopelus_warmingii,gadus_morhua
```

Note: a new type of warning was printed to the screen when the script was run, presumably due to the presence of polytomies in the collapsed trees: "this really only works with nexus or newick". Consequently, the resulting file had to be cleaned up by deleting rows that contained no information at all (i.e., neither root-to-tip variance nor tree length) and those that contained NAs.

Step 14: clean up the resulting file:

```{r}
var_collapsed <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-collapsed-chunks",na.strings=c("", "NA"), sep = "\t")
var_collapsed <- na.omit(var_collapsed)
```
```{r, eval = FALSE}
write.table(var_collapsed, "/Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-collapsed-chunks", quote = F, col.names = F, row.names = F, sep = "\t")
```

The cleaned up file contains information about root-to-tip variance and tree length for only 558 chunks, suggesting that SortaDate was able to estimate these variables only in 47.3% of cases.

Step 15: get_bp_genetrees on chunks

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_bp_genetrees.py /Users/David/Grive/Alfaro_Lab/Pelagia/Min_SH-like_90_chunk_trees /Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-collapsed-chunks
```

The bipartition support was calculated successfully for all chunks. Therefore, the next step (combining the two files) was performed manually in R rather than using the provided script.

Step 16: manually combine the results with the `var-pelagia-collapsed-chunks` file:

```{r}
colnames(var_collapsed) <- c("Locus", "Variance", "Length")
bp_collapsed <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-collapsed-chunks")
colnames(bp_collapsed) <- c("Locus", "Support")
comb_collapsed <- merge(var_collapsed, bp_collapsed, by = "Locus")
```
```{r, eval = FALSE}
write.table(comb_collapsed, "/Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-collapsed-chunks", quote = F, col.names = F, row.names = F, sep = "\t")
```

Step 17: recalculate the RF distances from the reference tree:

```
python robinsonfoulds.py > rf-pelagia-collapsed-chunks
```

This task was successfully performed for 796 out of the 1180 chunks.

Step 18: get_good_genes

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_good_genes.py /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-collapsed-chunks --max 558 --order 3,1,2 --outf /Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-collapsed-chunks
```

Step 19: bivariate plots showing the relationships among the three variables:

```{r, echo = FALSE}
gg_pelagia_colchunks <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-collapsed-chunks", header = TRUE)
par(oma = c(0, 0, 2, 0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(gg_pelagia_colchunks$bipartition, gg_pelagia_colchunks$treelength, xlab = "Bipartition support", ylab = "Tree length", main = "Tree length vs. bipartition support", col = "lightcoral")
abline(lm(gg_pelagia_colchunks$treelength ~ gg_pelagia_colchunks$bipartition))
plot(gg_pelagia_colchunks$bipartition, gg_pelagia_colchunks$root.to.tip_var, xlab = "Bipartition support", ylab = "Root-to-tip variance", main = "Clockiness vs. bipartition support", col = "greenyellow")
abline(lm(gg_pelagia_colchunks$root.to.tip_var ~ gg_pelagia_colchunks$bipartition))
plot(gg_pelagia_colchunks$root.to.tip_var, gg_pelagia_colchunks$treelength, xlab = "Root-to-tip variance", ylab = "Tree length", main = "Tree length vs. clockiness", col = "cornflowerblue")
abline(lm(gg_pelagia_colchunks$treelength ~ gg_pelagia_colchunks$root.to.tip_var))
mtext("558 chunks (no nodes w/ SH-like < 0.9)", outer = TRUE, cex = 1.5)
```

Step 20: calculate chunk lengths using the following script:

```{r, eval = FALSE}
args <- commandArgs(trailingOnly = TRUE)
a <- read.table(paste("/Users/David/Grive/Alfaro_Lab/Pelagia/Chunks/", args, sep = ""), header = FALSE, stringsAsFactors = FALSE)
cat(c(paste(strsplit(args, ".", fixed = TRUE)[[1]][1], ".tre", sep = ""), nchar(a[2,])), "\n")
```

```
for file in *fasta; do Rscript chunklengthcounter.R "$file" ; done > /Users/David/Grive/Alfaro_Lab/Pelagia/chunklengths.txt
```

Step 21: merge this info with the list of bipartition support scores:

```{r}
comb_chunks <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-chunks", 
                          header = TRUE)
colnames(comb_chunks) <- c("Locus", "Root_to_tip_variance", "Tree_length", "Bipartition_support")
chunk_lengths <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/chunklengths.txt")
colnames(chunk_lengths) <- c("Locus", "Length")
allchunkinfo <- merge(comb_chunks, chunk_lengths, by = "Locus")
```

Step 22: sort chunks by bipartition support from highest to lowest; keep adding them until the total length reaches that of the 138-locus alignment created in track 3 (see below); then check the bipartition support value that you end up with:

```{r}
allchunkinfo <- allchunkinfo[order(-allchunkinfo$Bipartition_support),] 
totallength <- 0
for (i in 1:nrow(allchunkinfo)) {
  totallength <- totallength + allchunkinfo$Length[i]
  if (totallength >= 114462) {                         # see step 16 of track 3
    print(allchunkinfo$Bipartition_support[i])
    break
  }
}
```










# Track 3: SortaDate on loci

Step 1: get_good_genes

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_good_genes.py /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-filtered-loci --max 440 --order 3,1,2 --outf /Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-filtered
```

Step 2: make bivariate plots showing the relationships between each pair of variables:

```{r, echo = FALSE}
gg_pelagia <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-filtered", header = TRUE)
par(oma = c(0, 0, 2, 0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(gg_pelagia$bipartition, gg_pelagia$treelength, xlab = "Bipartition support", ylab = "Tree length", main = "Tree length vs. bipartition support", col = "lightcoral")
abline(lm(gg_pelagia$treelength ~ gg_pelagia$bipartition))
plot(gg_pelagia$bipartition, gg_pelagia$root.to.tip_var, xlab = "Bipartition support", ylab = "Root-to-tip variance", main = "Clockiness vs. bipartition support", col = "greenyellow")
abline(lm(gg_pelagia$root.to.tip_var ~ gg_pelagia$bipartition))
plot(gg_pelagia$root.to.tip_var, gg_pelagia$treelength, xlab = "Root-to-tip variance", ylab = "Tree length", main = "Tree length vs. clockiness", col = "cornflowerblue")
abline(lm(gg_pelagia$treelength ~ gg_pelagia$root.to.tip_var))
mtext("440 loci (fully resolved trees)", outer = TRUE, cex = 1.5)
```

### Collapsed trees

Step 3: collapse all nodes with SH-like support values below 0.90:

```
python collapsenodes.py
```

Step 4: Copy the locus alignments into the new folder:

```
cp 95_percent_complete/*.fasta Min_SH-like_90_trees
```

Step 5: Re-run SortaDate on the collapsed trees:

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_var_length.py /Users/David/Grive/Alfaro_Lab/Pelagia/Min_SH-like_90_trees --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-collapsed-loci --outg ceratoscopelus_warmingii,gadus_morhua
```

Step 6: clean up the resulting file.

```{r, echo = FALSE}
var_collapsed <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-collapsed-loci",na.strings=c("", "NA"), sep = "\t")
var_collapsed <- na.omit(var_collapsed)
```
```{r, echo = FALSE, eval = FALSE}
write.table(var_collapsed, "/Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-collapsed-loci", quote = F, col.names = F, row.names = F, sep = "\t")
```

This reduced the number of usable loci from 615 to 313 (success rate = 50.9%).

Step 7: get_bp_genetrees:

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_bp_genetrees.py /Users/David/Grive/Alfaro_Lab/Pelagia/Min_SH-like_90_trees /Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-collapsed-loci
```

The bipartition support scores were calculated successfully for all loci.

Step 8: manually combine the results with the `var-pelagia-collapsed-loci` file:

```{r}
colnames(var_collapsed) <- c("Locus", "Variance", "Length")
bp_collapsed <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-collapsed-loci")
colnames(bp_collapsed) <- c("Locus", "Support")
comb_collapsed <- merge(var_collapsed, bp_collapsed, by = "Locus")
```
```{r, eval = FALSE}
write.table(comb_collapsed, "/Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-collapsed-loci", quote = F, col.names = F, row.names = F, sep = "\t")
```

Step 9: Re-calculate the Robinson-Foulds distances from the ExaML tree:

```
python robinsonfoulds.py > /Users/David/Grive/Alfaro_Lab/Pelagia/rf-pelagia-collapsed-loci
```

Note that this step was successfully performed for 487 loci only. Moreover, not all of these overlap with the 313 loci for which root-to-tip variance and tree length are available.

Step 10: get_good_genes:

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_good_genes.py /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-collapsed-loci --max 313 --order 3,1,2 --outf /Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-collapsed
```

Step 11: show the relationships among the three variables:

```{r, echo = FALSE}
gg_pelagia_collapsed <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-collapsed", header = TRUE)
par(oma = c(0, 0, 2, 0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(gg_pelagia_collapsed$bipartition, gg_pelagia_collapsed$treelength, xlab = "Bipartition support", ylab = "Tree length", main = "Tree length vs. bipartition support", col = "lightcoral")
abline(lm(gg_pelagia_collapsed$treelength ~ gg_pelagia_collapsed$bipartition))
plot(gg_pelagia_collapsed$bipartition, gg_pelagia_collapsed$root.to.tip_var, xlab = "Bipartition support", ylab = "Root-to-tip variance", main = "Clockiness vs. bipartition support", col = "greenyellow")
abline(lm(gg_pelagia_collapsed$root.to.tip_var ~ gg_pelagia_collapsed$bipartition))
plot(gg_pelagia_collapsed$root.to.tip_var, gg_pelagia_collapsed$treelength, xlab = "Root-to-tip variance", ylab = "Tree length", main = "Tree length vs. clockiness", col = "cornflowerblue")
abline(lm(gg_pelagia_collapsed$treelength ~ gg_pelagia_collapsed$root.to.tip_var))
mtext("313 loci (no nodes w/ SH-like < 0.9)", outer = TRUE, cex = 1.5)
```

### Loci with the highest bipartition support

Step 12: create a list of all loci whose fully resolved trees have bipartition support values of 0.5 or higher:

```{r}
gg_pelagia <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-filtered", header = TRUE, stringsAsFactors = FALSE)
```
```{r, eval = FALSE}
write(gg_pelagia$name[gg_pelagia$bipartition >= 0.5], "/Users/David/Grive/Alfaro_Lab/Pelagia/95_percent_complete/highest-bipart-support-loci.txt", sep = "\n")
```

This criterion is satisfied by 138 (22.4%) loci.

Step 13: make sure the resulting file lists alignments rather than trees:

```
gsed -i 's/tre/fasta/g' highest-bipart-support-loci.txt
```

Step 14: copy the corresponding alignments into a new folder:

```
for file in `cat highest-bipart-support-loci.txt`; do cp "$file" /Users/David/Grive/Alfaro_Lab/Pelagia/138_loci/ ; done
```

Step 15: convert these from FASTA to Nexus (Geneious -> Batch export).

Step 16: concatenate these using a phyluce script:

```
phyluce_align_format_nexus_files_for_raxml \
    --alignments 138_loci \
    --output 138_loci_concat \
    --charsets
```

The total length of the concatenated file is 114,462 bp.

## Subtrack 3.1

Run the PhyloBayes model comparison ("autocorrelation test") on each of the 138 loci. To do this, the variance-covariance matrix for branch lengths needs to be computed using the `estbranches` program from the `multidistribute` package, which in turn requires a `baseml` estimate of the substitution model parameters. To perform this on the desired scale, the following script was used to generate `baseml` configuration files:

Step 1: convert the alignments from Nexus to Phylip (Geneious -> Batch export).

Step 2: in `baseml`, the guide tree must contain only those taxa that are also present in the alignment. Since not all alignments contain the full set of taxa present in the tree, the corresponding trees inferred using RAxML were copied into the same directory:

```
for file in `cat highest-bipart-support-trees.txt`; do cp "$file" /Users/David/Grive/Alfaro_Lab/Pelagia/138_loci/ ; done
```

Step 3: reformat the trees to strip them of branch length and node support information (which confuses `baseml`):

```
for file in *.tre; do gsed -i 's/[0-9]*\.[0-9]*//g' "$file" ; done
for file in *.tre; do gsed -i 's/:e-06//g' "$file" ; done
for file in *.tre; do gsed -i 's/:e-05//g' "$file" ; done
for file in *.tre; do gsed -i 's/:e-07//g' "$file" ; done
for file in *.tre; do gsed -i 's/://g' "$file" ; done
for file in *.tre; do gsed -i 's/)0)/))/g' "$file" ; done
for file in *.tre; do gsed -i 's/)1)/))/g' "$file" ; done
for file in *.tre; do gsed -i 's/)0,/),/g' "$file" ; done
for file in *.tre; do gsed -i 's/)1,/),/g' "$file" ; done
```

Step 4: run the following script to generate `baseml` configuration files for each locus:

```{r, eval = FALSE}
args <- commandArgs(trailingOnly = TRUE)
uce_name <- strsplit(args, ".", fixed = TRUE)[[1]][1]

seqfile <- paste("      seqfile = ", args, " * sequence data file name", sep = "")
outfile <- paste("      outfile = ", uce_name, " * main result file", sep = "")
treefile <- paste("     treefile = ", uce_name, ".tre * tree structure file name", sep = "")
rest <- paste("                                                                       ",
"        noisy = 3   * 0,1,2,3: how much rubbish on the screen",
"      verbose = 1   * 1: detailed output, 0: concise output",
"      runmode = 0   * 0: user tree; 1: semi-automatic; 2: automatic",
"                    * 3: StepwiseAddition; (4,5):PerturbationNNI",
"                                                                ",
"        model = 3   * 0:JC69, 1:K80, 2:F81, 3:F84, 4:HKY85",
"                    * 5:T92, 6:TN93, 7:REV, 8:UNREST, 9:REVu; 10:UNRESTu",
"        Mgene = 0   * 0:rates, 1:separate; 2:diff pi, 3:diff kapa, 4:all diff",
"                                                                             ",
"        ndata = 1   * number of data sets",
"        clock = 1   * 0:no clock, 1:clock; 2:local clock; 3:CombinedAnalysis",
"                                                                            ",
"    fix_alpha = 0   * 0: estimate alpha; 1: fix alpha at value below",
"        alpha = 0.1   * initial or fixed alpha, 0:infinity (constant rate)",
"       Malpha = 0   * 1: different alpha's for genes, 0: one alpha",
"        ncatG = 8   * # of categories in the dG, AdG, or nparK models of rates",
"                                                                              ",
"      fix_rho = 1   * 0: estimate rho; 1: fix rho at value below", 
"          rho = 0   * initial or fixed rho, 0:no correlation",
"        nparK = 0   * rate-class models. 1:rK, 2:rK&fK, 3:rK&MK(1/K), 4:rK&MK",
"                                                                             ",
"        nhomo = 1   * 0 & 1: homogeneous, 2: kappa for branches, 3: N1, 4: N2",
"        getSE = 0   * 0: don't want them, 1: want S.E.s of estimates",
"                                                                    ",
" RateAncestor = 0   * (0,1,2): rates (alpha>0) or ancestral states",
"                                                                  ",
"   Small_Diff = 1e-6",
"       method = 0   * 0: simultaneous; 1: one branch at a time", sep = "\n")
write.table(rbind(seqfile, outfile, rest), paste("/Users/David/Grive/Alfaro_Lab/Pelagia/138_loci/", uce_name, ".ctl", sep = ""), row.names = FALSE, col.names = FALSE, quote = FALSE)
```

The F84 model was selected for all loci to ensure compatibility with `estbranches`.

```
for file in *.phy; do Rscript generatebasemlctl.R "$file" ; done
```

Step 5: reformat the Phylip alignments so that they are readable to PAML:

```
for file in *.phy; do gsed -i 's/ /  /g' "$file" ; done
```

Step 6: upload the files to the cluster and start the analysis:

```
ls *.{phy,tre,ctl} > tosync.txt
rsync --files-from=/Users/David/Grive/Alfaro_Lab/Pelagia/138_loci/tosync.txt /Users/David/Grive/Alfaro_Lab/Pelagia/138_loci analysis@azathoth.eeb.ucla.edu:~/David/pelagiaPAML/138_loci
baseml *.ctl
```

**Note:** For UCE 1118, `baseml` kept returning the following error:

```
Seq/species #1 (acanthocybium_solandri_95) occurs more than once in the tree
```

This was not the case; the tree merely contained two similarly named taxa ("acanthocybium_solandri_95" and "acanthocybium_solandri_95b") that were also present in all the other trees. The sequences of these taxa were not identical. Moreover, even changing the order of the two taxa in the tree, renaming one or both of them, or even deleting them did not solve the problem. Therefore, this locus was omitted from the analysis. The same approach was taken when the same error was encountered in analyzing UCEs 1199, 1239, and 1243.

## Subtrack 3.2

Step 1: rerun `phyluce` on the 138 loci to produce a concatenated alignment in the Nexus format:

```
phyluce_align_format_nexus_files_for_raxml \
    --alignments 138_loci \
    --output 138_loci_concat_nexus \
    --nexus \
    --charsets
```

Step 2: calculate the means of exponential node priors such that the 95th percentile equals the Hedman-based upper bound:

```{r, echo = FALSE}
aipichthys <- (143.01 - 98.0)/qexp(0.95)
paste("Aipichthys", aipichthys, sep = " = ")
anenchelum <- (65.16 - 41.2)/qexp(0.95)
paste("Anenchelum", anenchelum, sep = " = ")
argestichthys <- (75.14 - 54.17)/qexp(0.95)
paste("Argestichthys", argestichthys, sep = " = ")
ariomma <- (56.37 - 5.94)/qexp(0.95)
paste("Ariomma", ariomma, sep = " = ")
bannikovichthys <- (69.47 - 11.9)/qexp(0.95)
paste("Bannikovichthys", bannikovichthys, sep = " = ")
cretatriacanthus <- (98.08 - 69.71)/qexp(0.95)
paste("Cretatriacanthus", cretatriacanthus, sep = " = ")
eocoelopoma <- (67.61 - 54.17)/qexp(0.95)
paste("Eocoelopoma", eocoelopoma, sep = " = ")
gasterorhamphosus <- (98.08 - 69.71)/qexp(0.95)
paste("Gasterorhamphosus", gasterorhamphosus, sep = " = ")
hemithyristes <- (45.46 - 15.97)/qexp(0.95)
paste("Hemithyristes", hemithyristes, sep = " = ")
landanichthys <- (86.02 - 62.29)/qexp(0.95)
paste("Landanichthys", landanichthys, sep = " = ")
palaeocybium <- (61.68 - 50.5)/qexp(0.95)
paste("Palaeocybium", palaeocybium, sep = " = ")
pinichthys <- (71.9 - 32.02)/qexp(0.95)
paste("Pinichthys", pinichthys, sep = " = ")
scomber_colias <- (48.12 - 5.94)/qexp(0.95)
paste("Scomber_colias", scomber_colias, sep = " = ")
scomber_sadii <- (59.06 - 33.9)/qexp(0.95)
paste("Scomber_sadii", scomber_sadii, sep = " = ")
scomberomorus <- (53.29 - 15.1)/qexp(0.95)
paste("Scomberomorus", scomberomorus, sep = " = ")
stichocentrus <- (128.83 - 98)/qexp(0.95)
paste("Stichocentrus", stichocentrus, sep = " = ")
thrysitoides <- (56.17 - 33.9)/qexp(0.95)
paste("Thrysitoides", thrysitoides, sep = " = ")
thunnus <- (59.67 - 38.62)/qexp(0.95)
paste("Thunnus", thunnus, sep = " = ")
urosphenopsis <- (74.88 - 54.17)/qexp(0.95)
paste("Urosphenopsis", urosphenopsis, sep = " = ")
mene <- (84.67 - 55.2)/qexp(0.95)
paste("Mene", mene, sep = " = ")
```

Step 3: create a BEAST 1.8 XML file with the following settins:

- substitution models: a separate GTR+$\Gamma$ model with estimated base frequencies and 4 discrete rate categories for each of the 138 loci/partitions
- node priors: see above, plus a truncated exponential root prior from 98 to 143 Ma with the mean at the midpoint (22.5)
- `ucld.mean` prior: lognormal, mean = 0.005 in real space, log(st. dev.) = 1, initial value = 0.005
- `ucld.stdev` prior: truncated exponential from 0 to 1, mean = 0.3, initial value = 0.1
- `allMus` prior: gamma, shape = 0.001, scale = 1000, initial value = 1
- `ucld.mean` operator: tuning = 0.9, weight = 6.0
- `ucld.stdev` operator: tuning = 0.9, weight = 6.0
- MCMC settings: 800 million generations, log every 1000 generations

## Subtrack 3.3

Step 1: run PartitionFinder on the following configuration file using the command below:

```{r comment='', echo = FALSE}
cat(readLines("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/partition_finder.cfg"), sep = '\n')
```

```
python PartitionFinder.py /home/analysis/partitionfinder-new/pelagia-138-loci/partition_finder.cfg --verbose --raxml --save-phylofiles
```

The algorithm finished after 118 steps with 21 final partitions. The respective charsets were extracted from the `best_scheme` file and plugged as data blocks into a new round of partitioning. However, the second analysis was unable to improve on the results of the first, and did not cluster the 21 partitions any further.

Step 2: prepare a list of alignments and information files to download:

```{r}
new_partitions <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/partitions.txt", sep = "|", header = TRUE, stringsAsFactors = FALSE, strip.white = TRUE)
new_alignments <- paste(gsub(" ", "", new_partitions[,4], fixed = TRUE), ".phy", sep = "")
infofiles <- paste("RAxML_info.", gsub(" ", "", new_partitions[,4], fixed = TRUE), "_GTR+G.txt", sep = "")
```
```{r, eval = FALSE}
write(cbind(new_alignments, infofiles), "/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/whattorsync.txt")
```

Step 3: convert the partition alignments from Phylip to FASTA (Geneious -> Batch export) and infer partition trees using RAxML:

```{r}
args <- commandArgs(trailingOnly = TRUE)

commandpart1 <- paste("ete3 build -w standard_raxml -n /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/Partitions/", args, sep = "")
commandpart2 <- paste(" -o /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/Partitions/", args, sep="")
wholecommand <- paste(commandpart1, commandpart2, "-tree", sep = "")
```
```{r, eval = FALSE}
write(wholecommand, "/Users/David/Grive/Alfaro_Lab/Pelagia/21-pfinder-partition-analysis.sh")
```

Step 4: extract the 21 trees from their subdirectories and SortaDate them:

```
find . -name '*.nw' -exec mv {} . \;
find . -name "*.fasta.final_tree.nw" -exec rename 's/.fasta.final_tree.nw$/.tre/' {} \;
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_var_length.py /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/Partitions --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-21-partitions --outg ceratoscopelus_warmingii,gadus_morhua
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_bp_genetrees.py /Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/Partitions /Users/David/Grive/Alfaro_Lab/Pelagia/ExaML_pelagia-75p.tre --flend .tre --outf /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-21-partitions
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/combine_results.py /Users/David/Grive/Alfaro_Lab/Pelagia/var-pelagia-21-partitions /Users/David/Grive/Alfaro_Lab/Pelagia/bp-pelagia-21-partitions --outf /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-21-partitions
```

Step 5: does tree length (estimated by SortaDate) correlate with branch rate multipliers (estimated by PartitionFinder)?

```{r, echo = FALSE}
library(ggplot2)

new_multipliers <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/PartitionFinder/138_loci/multipliers.txt")
partitions <- cbind(new_partitions, new_multipliers)
colnames(partitions)[4] <- "name"
colnames(partitions)[6] <- "multiplier"
sortadated_partitions <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-21-partitions", header = TRUE, stringsAsFactors = FALSE)
colnames(sortadated_partitions) <- c("name", "root.to.tip_var", "treelength", "bipartition")
sortadated_partitions <- merge(sortadated_partitions, partitions, by = "name")

ggplot(sortadated_partitions, aes(x=treelength, y = multiplier)) + geom_point(color = 'black') + xlab("Tree length (assumed proportional to rate)") +  ylab(expression("Branch rate multiplier")) + geom_smooth(method = 'lm', size = 0.5, color = 'black') + annotate("text", x = 10, y = 1.75, label = lm_eqn(sortadated_partitions), parse = TRUE)
```

Step 6: bivariate plots showing the relationships among the three variables:

```{r, echo = FALSE}
par(oma = c(0, 0, 2, 0))
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(sortadated_partitions$bipartition, sortadated_partitions$treelength, xlab = "Bipartition support", ylab = "Tree length", main = "Tree length vs. bipartition support", col = "lightcoral")
abline(lm(sortadated_partitions$treelength ~ sortadated_partitions$bipartition))
plot(sortadated_partitions$bipartition, sortadated_partitions$root.to.tip_var, xlab = "Bipartition support", ylab = "Root-to-tip variance", main = "Clockiness vs. bipartition support", col = "greenyellow")
abline(lm(sortadated_partitions$root.to.tip_var ~ sortadated_partitions$bipartition))
plot(sortadated_partitions$root.to.tip_var, sortadated_partitions$treelength, xlab = "Root-to-tip variance", ylab = "Tree length", main = "Tree length vs. clockiness", col = "cornflowerblue")
abline(lm(sortadated_partitions$treelength ~ sortadated_partitions$root.to.tip_var))
mtext("21 pfinder partitions of the 138-locus alignment", outer = TRUE, cex = 1.5)
```

Step 7: select the 10 best partitions

```
python /Users/David/Grive/Alfaro_Lab/SortaDate/src/get_good_genes.py /Users/David/Grive/Alfaro_Lab/Pelagia/comb-pelagia-21-partitions --max 10 --order 3,1,2 --outf /Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-21-partitions
```

Step 8: how many basepairs does this selection account for?

```{r, echo = FALSE}
best_10_part <- read.table("/Users/David/Grive/Alfaro_Lab/Pelagia/gg-pelagia-21-partitions", header = TRUE)
sum(sortadated_partitions[paste(sortadated_partitions$name, ".phy.tre", sep = "") %in% best_10_part$name, 7])
```





# MCMCTree concatenated analysis (615 loci, unpartitioned)

Step 1: delete node numbers from the ExaML tree, as the tree file used by MCMCTree should contain no information except for topology and calibrations:

```
gsed -i 's/)[0-9]*:/):/g' ExaML_pelagia-75p-examl.tre
```

Step 2: delete the branch lengths.

```
gsed -i 's/:[^:]*,/,/g' ExaML_pelagia-75p-examl.tre
gsed -i 's/:[^:]*)/)/g' ExaML_pelagia-75p-examl.tre
```

Step 3: make the phylip alignment file compatible with PAML by ensuring that there are at least two spaces between the name of a taxon and the corresponding sequence:

```
gsed -i 's/ /  /g' raxml_95_alignments.phylip
```

Step 4: run a baseml analysis to find appropriate priors to be placed on the `rgene_gamma` and `sigma2_gamma` hyperparameters. The analysis was run under the strict clock model conditioned on the root age, which was set to 113.02 Ma (the mean of the exponential distribution whose offset equals 98 Ma and whose 95th percentile equals 143 Ma):

```{r comment='', echo = FALSE}
cat(readLines("/Users/David/Grive/Alfaro_Lab/Pelagia/pelagia_concatenated_baseml.ctl"), sep = '\n')
```

The baseml analysis finished up in 97:33:54, with the final substitution rate estimate of 0.145074 substitution per 100 million years. Following Alfaro et al., $\alpha$ (the gamma shape parameter) was set to 2 and $\beta$ (the corresponding scale parameter) was chosen so that the mean ($\frac{\alpha}{\beta}$) would equal to the rate when rescaled to a time unit of 10 million years:

```{r}
beta_pelagia <- 2 / (0.145074 / 10)
beta_pelagia
```

This corresponds to the following `rgene_gamma` prior (i.e., the prior on the mean of the lognormal distribution from which branch are drawn):

```{r, echo = FALSE}
curve(dgamma(x, shape = 2, rate = 137.86), from = 0, to = 0.5, main = "Rgene_gamma (mean substitution rate) prior", xlab = "Substitutions per 10 million years", ylab = "Probability density")
```

Following Alfaro et al., the `sigma2_prior`, which specifies the variance of the logarithm of the rate of evolution, was parameterized as the gamma distribution with $\alpha$ = 2 and $\beta$ = 5. This yielded the following rate distribution:

```{r, echo = FALSE}
curve(dlnorm(x, meanlog = log(2/137.86), sdlog = sqrt(2/5)), from = 0, to = 0.5, main = "Relaxed clock rate distribution", xlab = "Substitutions per 10 million years", ylab = "Probability density")
```

Step 5: start 15 MCMCTree analyses under the settings described in the configuration file below:

```{r comment='', echo = FALSE, message = FALSE, warning = FALSE}
cat(readLines("/Users/David/Grive/Alfaro_Lab/Pelagia/pelagia_concatenated_mcmctree.ctl"), sep = '\n')
```

Chain      | Run time
-----------|---------
pelagia 1  | 13:31:44
pelagia 2  | 13:14:21
pelagia 3  | 13:12:01
pelagia 4  | 13:27:41
pelagia 5  | 13:16:03
pelagia 6  | 13:24:40
pelagia 7  | 13:19:08
pelagia 8  | 13:15:22
pelagia 9  | 13:16:06
pelagia 10 | 13:15:47
pelagia 11 | 13:27:14
pelagia 12 | 13:21:35
pelagia 13 | 13:24:51
pelagia 14 | 13:26:08
pelagia 15 | 13:15:59