%%%%%%% 
% oligo.m
% Author: Lance Harden, Davidson College
% July 2006
% 
% Description: Oligo tells the user where to make cuts in order to
% optimize the melting temperatures of each individual oligo.  Optimization
% is defined as the maximin of the oligo melting temperatures.  This
% program allows the user to synthesize shorter oligos and then anneal them
% together.

%% CGCGGCCGCTTCTAGATTCGGGTGTCAACAATTGACCAAAATATCGATTTACAGCGTAATGCGCTTTCTAGTGCAAATTGTGACCGCATTTTACTAGTTGCGGCCGCC

seq = input('Enter your top strand sequence: ', 's');
seq = lower(seq);
total = length(seq);

tic

% 4 Oligos

temps4 = [];
oligos4 = [];
tempcount4 = 0;
runningcount4 = 0;
toadd4 = 0;

if total/3 <= 80
    for iA4 = 20:80
        BC4 = total - iA4;
        if BC4 < 40
            break
        end
        for iB4 = 20:(min(BC4 - 20,80))
            C4 = total - iA4 - iB4;
            if C4 < 20 | C4 > 80
                break
            end
            oligoA4 = seq(1:iA4);
            oligoB4 = seq((iA4 + 1):(iA4 + iB4));
            oligoC4 = seq((iA4 + iB4 + 1):total);
            oligos4(iB4 - 19 + runningcount4,:) = [length(oligoA4) length(oligoB4) length(oligoC4)];
            A4_gc_content = ( length(find(oligoA4 == 'g')) + length(find(oligoA4 == 'c')) )/ length(oligoA4);
            B4_gc_content = ( length(find(oligoB4 == 'g')) + length(find(oligoB4 == 'c')) )/ length(oligoB4);
            C4_gc_content = ( length(find(oligoC4 == 'g')) + length(find(oligoC4 == 'c')) )/ length(oligoC4);
            meltA4 = 81.5 + (16.6*log10(.1)) + (41*A4_gc_content) - (675/length(oligoA4));
            meltB4 = 81.5 + (16.6*log10(.1)) + (41*B4_gc_content) - (675/length(oligoB4));
            meltC4 = 81.5 + (16.6*log10(.1)) + (41*C4_gc_content) - (675/length(oligoC4));
            temps4(iB4 - 19 + runningcount4,:) = [meltA4 meltB4 meltC4];
        end
        toadd4 = total - 59 - tempcount4;
        runningcount4 = runningcount4 + toadd4;
        tempcount4 = tempcount4 + 1;
    end
    oligos4;
    temps4;
    minmelt4 = min(temps4,[],2);
    [maximin4,index4] = max(minmelt4);
    bestcuts4 = oligos4(index4,:)
    besttemps4 = temps4(index4,:)
end

% 5 Oligos

temps5 = [];
oligos5 = [];
Acount5 = 0;
runningcount5 = 0;
toadd5 = 0;

if total/4 <= 80
    for iA5 = 20:80
        Bcount5 = 0;
        BCD5 = total - iA5;
        if BCD5 < 60
            break
        end
        for iB5 = 20:(min(BCD5 - 40,80))
            CD5 = total - iA5 - iB5;
            if CD5 < 40
                break
            end
            for iC5 = 20:(min(CD5 - 20,80))
                D5 = total - iA5 - iB5 - iC5;
                if D5 < 20 | D5 > 80
                    break
                end
                oligoA5 = seq(1:iA5);
                oligoB5 = seq((iA5 + 1):(iA5 + iB5));
                oligoC5 = seq((iA5 + iB5 + 1):(iA5 + iB5 + iC5));
                oligoD5 = seq((iA5 + iB5 + iC5 + 1):total);
                oligos5(iC5 - 19 + runningcount5,:) = [length(oligoA5) length(oligoB5) length(oligoC5) length(oligoD5)];
                A5_gc_content = ( length(find(oligoA5 == 'g')) + length(find(oligoA5 == 'c')) )/ length(oligoA5);
                B5_gc_content = ( length(find(oligoB5 == 'g')) + length(find(oligoB5 == 'c')) )/ length(oligoB5);
                C5_gc_content = ( length(find(oligoC5 == 'g')) + length(find(oligoC5 == 'c')) )/ length(oligoC5);
                D5_gc_content = ( length(find(oligoD5 == 'g')) + length(find(oligoD5 == 'c')) )/ length(oligoD5);
                meltA5 = 81.5 + (16.6*log10(.1)) + (41*A5_gc_content) - (675/length(oligoA5));
                meltB5 = 81.5 + (16.6*log10(.1)) + (41*B5_gc_content) - (675/length(oligoB5));
                meltC5 = 81.5 + (16.6*log10(.1)) + (41*C5_gc_content) - (675/length(oligoC5));
                meltD5 = 81.5 + (16.6*log10(.1)) + (41*D5_gc_content) - (675/length(oligoD5));
                temps5(iC5 - 19 + runningcount5,:) = [meltA5 meltB5 meltC5 meltD5];
            end
            toadd5 = total - 79 - Acount5 - Bcount5;
            runningcount5 = runningcount5 + toadd5;
            Bcount5 = Bcount5 + 1;
        end
        Acount5 = Acount5 + 1;
    end
    oligos5;
    temps5;
    minmelt5 = min(temps5,[],2);
    [maximin5,index5] = max(minmelt5);
    bestcuts5 = oligos5(index5,:)
    besttemps5 = temps5(index5,:)
end

% 6 Oligos

temps6 = [];
oligos6 = [];
Acount6 = 0;
runningcount6 = 0;
toadd6 = 0;

if total/5 <= 80
    for iA6 = 20:80
        Bcount6 = 0;
        BCDE6 = total - iA6;
        if BCDE6 < 80
            break
        end
        for iB6 = 20:(min(BCDE6 - 60,80)
            Ccount6 = 0;
            CDE6 = total - iA6 - iB6;
            if CDE6 < 60
                break
            end
            for iC6 = 20:(min(CDE6 - 40,80)
                DE6 = total - iA6 - iB6 - iC6;
                if DE6 < 40
                    break
                end
                for iD6 = 20:(min(DE6 - 20,80)
                    E6 = total - iA6 - iB6 - iC6 - iD6;
                    if E6 < 20 | E6 > 80
                        break
                    end
                    oligoA6 = seq(1:iA6);
                    oligoB6 = seq((iA6 + 1):(iA6 + iB6));
                    oligoC6 = seq((iA6 + iB6 + 1):(iA6 + iB6 + iC6));
                    oligoD6 = seq((iA6 + iB6 + iC6 + 1):(iA6 + iB6 + iC6 + iD6));
                    oligoE6 = seq((iA6 + iB6 + iC6 + iD6 + 1):total);
                    oligos6(iD6 - 19 + runningcount6,:) = [length(oligoA6) length(oligoB6) length(oligoC6) length(oligoD6) length(oligoE6)];
                    A6_gc_content = ( length(find(oligoA6 == 'g')) + length(find(oligoA6 == 'c')) )/ length(oligoA6);
                    B6_gc_content = ( length(find(oligoB6 == 'g')) + length(find(oligoB6 == 'c')) )/ length(oligoB6);
                    C6_gc_content = ( length(find(oligoC6 == 'g')) + length(find(oligoC6 == 'c')) )/ length(oligoC6);
                    D6_gc_content = ( length(find(oligoD6 == 'g')) + length(find(oligoD6 == 'c')) )/ length(oligoD6);
                    E6_gc_content = ( length(find(oligoE6 == 'g')) + length(find(oligoE6 == 'c')) )/ length(oligoE6);
                    meltA6 = 81.5 + (16.6*log10(.1)) + (41*A6_gc_content) - (675/length(oligoA6));
                    meltB6 = 81.5 + (16.6*log10(.1)) + (41*B6_gc_content) - (675/length(oligoB6));
                    meltC6 = 81.5 + (16.6*log10(.1)) + (41*C6_gc_content) - (675/length(oligoC6));
                    meltD6 = 81.5 + (16.6*log10(.1)) + (41*D6_gc_content) - (675/length(oligoD6));
                    meltE6 = 81.5 + (16.6*log10(.1)) + (41*E6_gc_content) - (675/length(oligoE6));
                    temps6(iD6 - 19 + runningcount6,:) = [meltA6 meltB6 meltC6 meltD6 meltE6];
                end
                toadd6 = total - 99 - Acount6 - Bcount6 - Ccount6;
                runningcount6 = runningcount6 + toadd6;
                Ccount6 = Ccount6 + 1;
            end
            Bcount6 = Bcount6 + 1;
        end
        Acount6 = Acount6 + 1;
    end
    oligos6;
    temps6;
    minmelt6 = min(temps6,[],2);
    [maximin6,index6] = max(minmelt6);
    bestcuts6 = oligos6(index6,:)
    besttemps6 = temps6(index6,:)
end

% 7 Oligos

temps7 = [];
oligos7 = [];
Acount7 = 0;
runningcount7 = 0;
toadd7 = 0;

if total/6 <= 80
    for iA7 = 20:80
        Bcount7 = 0;
        BCDEF7 = total - iA7;
        if BCDEF7 < 100
            break
        end
        for iB7 = 20:(min(BCDEF7 - 80,80))
            Ccount7 = 0;
            CDEF7 = total - iA7 - iB7;
            if CDEF7 < 80
                break
            end
            for iC7 = 20:(min(CDEF7 - 60,80))
                Dcount7 = 0;
                DEF7 = total - iA7 - iB7 - iC7;
                if DEF7 < 60
                    break
                end
                for iD7 = 20:(min(DEF7 - 40,80))
                    EF7 = total - iA7 - iB7 - iC7 - iD7;
                    if EF7 < 40
                        break
                    end
                    for iE7 = 20:(min(EF7 - 20,80))
                        F7 = total - iA7 - iB7 - iC7 - iD7 - iE7;
                        if F7 < 20 | F7 > 80
                            break
                        end
                        oligoA7 = seq(1:iA7);
                        oligoB7 = seq((iA7 + 1):(iA7 + iB7));
                        oligoC7 = seq((iA7 + iB7 + 1):(iA7 + iB7 + iC7));
                        oligoD7 = seq((iA7 + iB7 + iC7 + 1):(iA7 + iB7 + iC7 + iD7));
                        oligoE7 = seq((iA7 + iB7 + iC7 + iD7 + 1):(iA7 + iB7 + iC7 + iD7 + iE7));
                        oligoF7 = seq((iA7 + iB7 + iC7 + iD7 + iE7 + 1):total);
                        oligos7(iE7 - 19 + runningcount7,:) = [length(oligoA7) length(oligoB7) length(oligoC7) length(oligoD7) length(oligoE7) length(oligoF7)];
                        A7_gc_content = ( length(find(oligoA7 == 'g')) + length(find(oligoA7 == 'c')) )/ length(oligoA7);
                        B7_gc_content = ( length(find(oligoB7 == 'g')) + length(find(oligoB7 == 'c')) )/ length(oligoB7);
                        C7_gc_content = ( length(find(oligoC7 == 'g')) + length(find(oligoC7 == 'c')) )/ length(oligoC7);
                        D7_gc_content = ( length(find(oligoD7 == 'g')) + length(find(oligoD7 == 'c')) )/ length(oligoD7);
                        E7_gc_content = ( length(find(oligoE7 == 'g')) + length(find(oligoE7 == 'c')) )/ length(oligoE7);
                        F7_gc_content = ( length(find(oligoF7 == 'g')) + length(find(oligoF7 == 'c')) )/ length(oligoF7);
                        meltA7 = 81.5 + (16.6*log10(.1)) + (41*A7_gc_content) - (675/length(oligoA7));
                        meltB7 = 81.5 + (16.6*log10(.1)) + (41*B7_gc_content) - (675/length(oligoB7));
                        meltC7 = 81.5 + (16.6*log10(.1)) + (41*C7_gc_content) - (675/length(oligoC7));
                        meltD7 = 81.5 + (16.6*log10(.1)) + (41*D7_gc_content) - (675/length(oligoD7));
                        meltE7 = 81.5 + (16.6*log10(.1)) + (41*E7_gc_content) - (675/length(oligoE7));
                        meltF7 = 81.5 + (16.6*log10(.1)) + (41*F7_gc_content) - (675/length(oligoF7));
                        temps7(iE7 - 19 + runningcount7,:) = [meltA7 meltB7 meltC7 meltD7 meltE7 meltF7];
                    end
                    toadd7 = total - 119 - Acount7 - Bcount7 - Ccount7 - Dcount7;
                    runningcount7 = runningcount7 + toadd7;
                    Dcount7 = Dcount7 + 1;
                end
                Ccount7 = Ccount7 + 1;
            end
            Bcount7 = Bcount7 + 1;
        end
        Acount7 = Acount7 + 1;
    end
    oligos7;
    temps7;
    minmelt7 = min(temps7,[],2);
    [maximin7,index7] = max(minmelt7);
    bestcuts7 = oligos7(index7,:)
    besttemps7 = temps7(index7,:)
end

% 8 Oligos

temps8 = [];
oligos8 = [];
Acount8 = 0;
runningcount8 = 0;
toadd8 = 0;

if total/7 <= 80
    for iA8 = 20:80
        Bcount8 = 0;
        BCDEFG8 = total - iA8;
        if BCDEFG8 < 120
            break
        end
        for iB8 = 20:(min(BCDEFG8 - 100,80))
            Ccount8 = 0;
            CDEFG8 = total - iA8 - iB8;
            if CDEFG8 < 100
                break
            end
            for iC8 = 20:(min(CDEFG8 - 80,80))
                Dcount8 = 0;
                DEFG8 = total - iA8 - iB8 - iC8;
                if DEFG8 < 80
                    break
                end
                for iD8 = 20:(min(DEFG8 - 60,80))
                    Ecount8 = 0;
                    EFG8 = total - iA8 - iB8 - iC8 - iD8;
                    if EFG8 < 60
                        break
                    end
                    for iE8 = 20:(min(EFG8 - 40,80))
                        FG8 = total - iA8 - iB8 - iC8 - iD8 - iE8;
                        if FG8 < 40
                            break
                        end
                        for iF8 = 20:(min(FG8 - 20,80))
                            G8 = total - iA8 - iB8 - iC8 - iD8 - iE8 - iF8;
                            if G8 < 20 | G8 > 80
                                break
                            end
                            oligoA8 = seq(1:iA8);
                            oligoB8 = seq((iA8 + 1):(iA8 + iB8));
                            oligoC8 = seq((iA8 + iB8 + 1):(iA8 + iB8 + iC8));
                            oligoD8 = seq((iA8 + iB8 + iC8 + 1):(iA8 + iB8 + iC8 + iD8));
                            oligoE8 = seq((iA8 + iB8 + iC8 + iD8 + 1):(iA8 + iB8 + iC8 + iD8 + iE8));
                            oligoF8 = seq((iA8 + iB8 + iC8 + iD8 + iE8 + 1):(iA8 + iB8 + iC8 + iD8 + iE8 + iF8));
                            oligoG8 = seq((iA8 + iB8 + iC8 + iD8 + iE8 + iF8 + 1):total);
                            oligos8(iF8 - 19 + runningcount8,:) = [length(oligoA8) length(oligoB8) length(oligoC8) length(oligoD8) length(oligoE8) length(oligoF8) length(oligoG8)];
                            A8_gc_content = ( length(find(oligoA8 == 'g')) + length(find(oligoA8 == 'c')) )/ length(oligoA8);
                            B8_gc_content = ( length(find(oligoB8 == 'g')) + length(find(oligoB8 == 'c')) )/ length(oligoB8);
                            C8_gc_content = ( length(find(oligoC8 == 'g')) + length(find(oligoC8 == 'c')) )/ length(oligoC8);
                            D8_gc_content = ( length(find(oligoD8 == 'g')) + length(find(oligoD8 == 'c')) )/ length(oligoD8);
                            E8_gc_content = ( length(find(oligoE8 == 'g')) + length(find(oligoE8 == 'c')) )/ length(oligoE8);
                            F8_gc_content = ( length(find(oligoF8 == 'g')) + length(find(oligoF8 == 'c')) )/ length(oligoF8);
                            G8_gc_content = ( length(find(oligoG8 == 'g')) + length(find(oligoG8 == 'c')) )/ length(oligoG8);
                            meltA8 = 81.5 + (16.6*log10(.1)) + (41*A8_gc_content) - (675/length(oligoA8));
                            meltB8 = 81.5 + (16.6*log10(.1)) + (41*B8_gc_content) - (675/length(oligoB8));
                            meltC8 = 81.5 + (16.6*log10(.1)) + (41*C8_gc_content) - (675/length(oligoC8));
                            meltD8 = 81.5 + (16.6*log10(.1)) + (41*D8_gc_content) - (675/length(oligoD8));
                            meltE8 = 81.5 + (16.6*log10(.1)) + (41*E8_gc_content) - (675/length(oligoE8));
                            meltF8 = 81.5 + (16.6*log10(.1)) + (41*F8_gc_content) - (675/length(oligoF8));
                            meltG8 = 81.5 + (16.6*log10(.1)) + (41*G8_gc_content) - (675/length(oligoG8));
                            temps8(iF8 - 19 + runningcount8,:) = [meltA8 meltB8 meltC8 meltD8 meltE8 meltF8 meltG8];
                        end
                        toadd8 = total - 139 - Acount8 - Bcount8 - Ccount8 - Dcount8 - Ecount8;
                        runningcount8 = runningcount8 + toadd8;
                        Ecount8 = Ecount8 + 1;
                    end
                    Dcount8 = Dcount8 + 1;
                end
                Ccount8 = Ccount8 + 1;
            end
            Bcount8 = Bcount8 + 1;
        end
        Acount8 = Acount8 + 1;
    end
    oligos8;
    temps8;
    minmelt8 = min(temps8,[],2);
    [maximin8,index8] = max(minmelt8);
    bestcuts8 = oligos8(index8,:)
    besttemps8 = temps8(index8,:)
end

time = toc