I must say that I did understand completely your problem.

However, if you have 15 assemblies (meaning 15 draft genomes) I would suggest you to annotate them and then run a program for pan/core genomes analysis. You could use PROKKA, for annotation it is quite fast, or RAST.

After that you can use a program for calculating the orthologues gene families in your genomes, and consequently the pan-genome. There are other discussion on this topic: A: Pan- and core- genome in bacteria

I personally use get_homologues

If you than want to align your genomes to a reference you can use something like MAUVE, but there are tons of tools for this

Hope this helps.

Hey.

Not sure if I understand what you are trying to accomplish. Are you trying to build a reference genome containing all the "pan-genomic" variation, so you can use that as a baseline for describing genomic variation?

If so, you can achieve this by:

1. Running get_homologues on the dataset as suggested by dago. (I advise using both -M and -G to intersect later, and also consider using -t 2 or 3 to remove singletons, some of which may be artifacial)
2. Using the downstream tool "compare_clusters.pl" included in get_homologues to produce a pan-intersection syntax: compare_clusters.pl -d (dir to MCL),(dir to COG) -m -o (output). You can use the flag -n if you want nucleotides instead of amino acids.
3. Read in each individual .fna/.faa file from 2), align them, and produce a consensus from the alignment. This can be performed using an R-script with the packages "DECIPHER" and "BioStrings", which have functions for aligning nucleotides/amino acids and for producing a consensus sequence from alignments.

R-script could be a loop that looks something like this (probably better ways of doing it). The example is for .fna (using -n in compare_clusters.pl). Script should (could) work if you have the packages installed and set the working directory to the output made by compare_clusters.pl

library("DECIPHER")
library("BioStrings")

library("seqinr") #For write.fasta

files <- list.files(pattern = ".fna")

Consensus.Seq.Vector <- NULL
for(a in 1:length(files)) {
  Seq <- readDNAStringSet(files[a])
  Aligned.seq <- AlignSeqs(Seq)
  Consensus.Seq.Vector[a] <- paste(consensusString(Aligned.seq, ambiguityMap="?", threshold=0.7))
}
Consensus.Seq.String <- paste(Consensus.Seq.Vector, sep="", collapse="")
write.fasta(sequences = Consensus.Seq.String, file.out="ConsensusGenome.fasta", open="a")

Hope this helps.