# load up libs
require("ape") #tested using 3.2
require("spider") #tested using 1.3-0
require("phyloch") #tested using 1.5-3

## open genbank data spreadsheet 
tab <- read.table("genbank_scrape.csv", header=TRUE, sep = ",", stringsAsFactors=FALSE, na.strings="NA", quote="\"")
# subset just the 'yes' rows, i.e. change the yes/no answer in the "Use" column to choose that species to be downloaded 
tab2 <- tab[tab$Use=="yes",]
# get gene names from row headers
tags <- names(tab2)[7:11]

## downloads all the sequences and writes them into a fasta file
jj <- list()
gb <- list()
sp <- list()
ac <- list()
for(i in 1:length(tags)){
  jj[[i]] <- tab2[,tags[i]][!is.na(tab2[,tags[i]])]
  gb[[i]] <- read.GB(jj[[i]])
  sp[[i]] <- sapply(tags[i], function(x) tab2[,"Species"][which(tab2[,tags[i]] %in% jj[[i]])])
  ac[[i]] <- attr(gb[[i]], "accession_num")
  names(gb[[i]]) <- paste(sp[[i]], "_", ac[[i]], "_", tags[i], sep="") 
  nam <- paste(tags[i], ".fas", sep="")
  write.dna(gb[[i]], file=nam, format="fasta", colw=10000) 
}#

## concatenate the sequences AFTER alignment and trimming
# read in the aligned sequences
rag <- read.dna(file="rag1_ALIGNED.fas", format="fasta")
sh3 <- read.dna(file="SH3PX3_ALIGNED.fas", format="fasta")
sreb <- read.dna(file="sreb2_ALIGNED.fas", format="fasta")
zic <- read.dna(file="zic1_ALIGNED.fas", format="fasta")
plag <- read.dna(file="plagl2_ALIGNED.fas", format="fasta")

# edit the names so all are identical for concatenation: strip away gene name and genbank number 
dimnames(rag)[[1]] <- sapply(strsplit(labels(rag), "_"), function(x) paste(x[1], x[2], sep="_"))
dimnames(sh3)[[1]] <- sapply(strsplit(labels(sh3), "_"), function(x) paste(x[1], x[2], sep="_"))
dimnames(sreb)[[1]] <- sapply(strsplit(labels(sreb), "_"), function(x) paste(x[1], x[2], sep="_"))
dimnames(zic)[[1]] <- sapply(strsplit(labels(zic), "_"), function(x) paste(x[1], x[2], sep="_"))
dimnames(plag)[[1]] <- sapply(strsplit(labels(plag), "_"), function(x) paste(x[1], x[2], sep="_"))

# concatenate the genes into supermatrix - full alignment (5 genes)
all <- c.genes(rag, sh3, sreb, zic, plag, match=FALSE)
# sort the names into alphabetical order
ord <- all[sort(dimnames(all)[[1]]),]
# replace "-" with "n"
ord <- as.DNAbin(gsub("-", "n", ord))

## write concatenated files
# fasta format
write.dna(ord, file="full_alignment.fas", format="fasta", colw=10000)
# phylip format
write.dna(ord, file="full_alignment.phy", format="sequential", colw=10000)
# nexus format (note that there is a bug in ape 3.2 that does not print the "=" between "DATATYPE" and "DNA" so need to change this by hand)
write.nexus.data(ord, file="full_alignment.nex", interleaved=FALSE)