(phixonline)-->
procedure printcounts(sequence ss)
-- Given DNA sequence(s), report the sequence, length and base counts
for i=1 to length(ss) do
string dna = ss[i]
sequence acgt = repeat(0,6)
for j=1 to length(dna) do
acgt[find(dna[j],"ACGT")+1] += 1
end for
acgt[$] = sum(acgt)
string ncf = "Nucleotide counts for :"
printf(1,"%s%s\n",{ncf,join(split_by(dna,50),"\n"&repeat(' ',length(ncf)))})
printf(1,"Base counts: Other:%d, A:%d, C:%d, G:%d, T:%d, total:%d\n\n",acgt)
end for
end procedure
function deduplicate(sequence ss)
-- Remove any strings contained within a larger string from a set of strings
sequence filtered = {}
for i=1 to length(ss) do
string si = ss[i]
bool found = false
for j=1 to length(ss) do
if i!=j and match(si,ss[j]) then
found = true
exit
end if
end for
if not found then
filtered = append(filtered, si)
end if
end for
return filtered
end function
procedure shortest_common_superstring(sequence ss)
-- Returns shortest common superstring of a set of strings
ss = deduplicate(unique(ss))
sequence shortestsuper = {join(ss,"")}
integer shortest = length(shortestsuper[1])
for p=1 to factorial(length(ss)) do
sequence perm = permute(p,ss)
string sup = perm[1]
for i=2 to length(perm) do
string pi = perm[i]
for j=-min(length(pi),length(sup)) to 0 do
string overlap = sup[j..$]
if overlap = pi[1..length(overlap)] then
sup &= pi[length(overlap)+1..$]
pi = ""
exit
end if
end for
if length(pi) then ?9/0 end if -- (sanity chk)
end for
if length(sup) < shortest then
shortest = length(sup)
shortestsuper = {sup}
elsif length(sup) = shortest
and not find(sup,shortestsuper) then
shortestsuper = append(shortestsuper,sup)
end if
end for
printcounts(shortestsuper)
end procedure
constant tests = {
{"TA", "AAG", "TA", "GAA", "TA"},
{"CATTAGGG", "ATTAG", "GGG", "TA"},
{"AAGAUGGA", "GGAGCGCAUC", "AUCGCAAUAAGGA"},
{"ATGAAATGGATGTTCTGAGTTGGTCAGTCCCAATGTGCGGGGTTTCTTTTAGTACGTCGGGAGTGGTATTAT",
"GGTCGATTCTGAGGACAAAGGTCAAGATGGAGCGCATCGAACGCAATAAGGATCATTTGATGGGACGTTTCGTCGACAAAGT",
"CTATGTTCTTATGAAATGGATGTTCTGAGTTGGTCAGTCCCAATGTGCGGGGTTTCTTTTAGTACGTCGGGAGTGGTATTATA",
"TGCTTTCCAATTATGTAAGCGTTCCGAGACGGGGTGGTCGATTCTGAGGACAAAGGTCAAGATGGAGCGCATC",
"AACGCAATAAGGATCATTTGATGGGACGTTTCGTCGACAAAGTCTTGTTTCGAGAGTAACGGCTACCGTCTT",
"GCGCATCGAACGCAATAAGGATCATTTGATGGGACGTTTCGTCGACAAAGTCTTGTTTCGAGAGTAACGGCTACCGTC",
"CGTTTCGTCGACAAAGTCTTGTTTCGAGAGTAACGGCTACCGTCTTCGATTCTGCTTATAACACTATGTTCT",
"TGCTTTCCAATTATGTAAGCGTTCCGAGACGGGGTGGTCGATTCTGAGGACAAAGGTCAAGATGGAGCGCATC",
"CGTAAAAAATTACAACGTCCTTTGGCTATCTCTTAAACTCCTGCTAAATGCTCGTGC",
"GATGGAGCGCATCGAACGCAATAAGGATCATTTGATGGGACGTTTCGTCGACAAAGTCTTGTTTCGAGAGTAACGGCTACCGTCTTCGATT",
"TTTCCAATTATGTAAGCGTTCCGAGACGGGGTGGTCGATTCTGAGGACAAAGGTCAAGATGGAGCGCATC",
"CTATGTTCTTATGAAATGGATGTTCTGAGTTGGTCAGTCCCAATGTGCGGGGTTTCTTTTAGTACGTCGGGAGTGGTATTATA",
"TCTCTTAAACTCCTGCTAAATGCTCGTGCTTTCCAATTATGTAAGCGTTCCGAGACGGGGTGGTCGATTCTGAGGACAAAGGTCAAGA"}
}
papply(tests, shortest_common_superstring)