#Run after LepCosall.R or similar tssv scoring script
setwd("/media/bob/backup data/Fastq")
library(dplyr)
library(seqinr)
loci <- c("pim4", "pim8", "pim10", "pim11", "pim12", "pim13", "pim15", "pim16", "pim17", "pim20")
for (i in loci){
  #for each locus read in sequence alignment
j <- paste("aln", i, ".fas.phy.fas", sep = "")
k <- read.alignment(j, "fasta")
l <- data.frame(k$nam, seq=unlist(k$seq))
#replace "0" and "1" with "a" and "c" for POFAD
m <- gsub('0', 'a', l$seq)
n <- gsub('1', 'c', m)
o <-data.frame(l$k.nam, n)  
#output edited fasta
write.fasta(as.list(o$n), as.list(o$l.k.nam), paste(i, "out.fas", sep = ""), open = "w", nbchar = 1000, as.string = FALSE )
#count "taxa"
ntaxa = nrow(o)
#count characters
nchar = nchar(o$n[1])
#make some nexus formating lines
header1 <- as.character("#NEXUS")
header2 <- as.character("Begin data;")
header3 <- as.character(paste(    "DIMENSIONS NTAX=", ntaxa, " NCHAR=", nchar, ";", sep = "")  )
header4 <- as.character("     FORMAT MISSING=? DATATYPE=DNA GAP=-;")
header5 <- as.character("MATRIX")
footer <- as.character("; end;")
#write nexus headers, then the name/sequences, then end characters for nexus 
  write.table(header1, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE)
  write.table(header2, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE, append = TRUE)
  write.table(header3, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE, append = TRUE)
  write.table(header4, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE, append = TRUE)
  write.table(header5, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE, append = TRUE)
  write.table(o, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE, append = TRUE)
  write.table(footer, file=(paste(i, ".nex", sep = "")), row.names = FALSE, quote = FALSE, col.names = FALSE, append = TRUE)
}
#run raxml.sh then use Tree2Fasta to separate loci within each contig group.