monkey2.0.tar

monkey2.0/0000755000163100010400000000000010457514562011044 5ustar  alanallmonkey2.0/all.tree0000644000163100010400000000014510432332204012455 0ustar  alanallt 7
r 6	n 5	l 0
n 5	n 4	l 2
n 4	l 1	l 3
b
4	0.09787	0.13253	
5	0.12665	0.23934	
6	0.154675	0.154675	
monkey2.0/bg.freq0000644000163100010400000000003010432332204012264 0ustar  alanallA 0.3
C 0.2
G 0.2
T 0.3
monkey2.0/eg.list0000644000163100010400000000025710432332204012320 0ustar  alanallexample/YBR018C_upstream.aln.fasta
example/YBR019C_upstream.aln.fasta
example/YBR020W_upstream.aln.fasta
example/YML051W_upstream.aln.fasta
example/YOR120W_upstream.aln.fasta
monkey2.0/example/0000755000163100010400000000000010457514562012477 5ustar  alanallmonkey2.0/example/YBR018C_upstream.aln.fasta0000644000163100010400000000635010432332204017164 0ustar  alanall>YBR018C_Sbay
TTTG-T-CCGCTCAGCTTTCTT---AAAG-AAAAATATACTTGATATA-CTTAGTTAGTTCATT-CCACAACACATATACAAT-TATATAACGAACCTTATGTCTACTAATTCTCCTTTTAATGAATATAACTAACTACAATAACTTTTAAATGGAG--TCTT-CT-GGCAGTAAAGCACCAAAATCTCCTCGAAA--TTTTACGGTTTTTTTATGTACAAACGAAAAGAAATTGATAGAGCAGAGCTTTTAGATGAACGGGAACTGTAGTAAAAAGAATA-AACATGTGCCGAGGA----ATTTAACGCT-CTCTTGAAA------------AC-GCTACG--CTA-ACTAAAAAAATTCCTTAAA-CTCA----TCTTCTAGAGATT-TCCATTGAA------------CTGTAGATGAGAAAGT-GTACACAATTC---CCCTGAACCATGTGCGTACCTTTTT-AGATTACTTAGCATCCC---TTGGAAGGTTGCTCTTCGGAGCACTGTTGGACGAAGGCTGATCGGATCA-TTTTCCTGTTATTTTCCGTCTGCACCAGTACAATAGGAAAGGGGGCCCTTGACCGCTCGGACAGCAGTCGACCGTGATCCGACATATTGACCATGCAGTGTCATGAAAAA--A----AAGAG-GCAATTTAGTATATCTAATT---G-----TTAATCACTGCAGATATATAAGCAG----ATCGATCAGAGGCAGTTATGAGTATCGGAGTGTGTGGAGTATCAGCATGACAAGAAAGATATAACCA-CCAACTAAACTTTTCTTTTCAATA
>YBR018C_Scer
TTTG-C-CAGCTTACTATCCT-TCT--TGAAAATAT--GCACTCTATATCTT--TTAGTTCTTAATTGCAACACATAGATTTGCTGTATAACGAATTTTATGCTATTT-TTTA--AATT---TGGAGTTCAGTGA---TAAA-AGTGTCACAGCGAA--TTTC-CT-CACATGTAGGGAC-----------CG-AA--TTGTTTACAAGTTCTCTGTACCACCATGGAGACATCAAAAATTGAAAATCTATGGAAAGATATGGACGGTAGCAACAAGAATATAGCACGAGCCGCGGAGTTCATTT--CGTTACTTTTGATATCACTCACAACTATTGCGAAGCGCTTCAGTGAAAAAAT--CATAAGGA-AA----AGTTGTAAATATTATTGGTAGTA-TTCGTTTGGTAAAGTAGA-----GGGGG-TA-ATTTTTC---CCCTTTATTTTGTTCATACATTCTT-AAATTGCTTTGCCTCTCCTTTTGGAAAGCTATACTTCGGAGCACTGTTGAGCGAAGGCTCATTAGATATATTTT-CTGTCATTTTCC--TTAACCCA-AAAATAAGGGAAAGGGTCCAAAAAGCGCTCGGACAACTGTTGACCGTGATCCGAAGGACTGGCTATACAGTGTTCACAAAAT--AGCCAA-GCTGAAAATAATGTGTAGCTATGTTCAGTTAGTTTGGCTAG--CAAAGATATAAA-AGCAG-GTCGG----AAATATTTATGGGCAT--TATTATGCAGAGCATCAACATGATA--AAA-A-A-AA-A---CAGTTGAAT-ATTC-CCTCAAAA
>YBR018C_Smik
TTTCTCACA-CT-GGTATCATTACCAATG-AAGAAT-GACACTCTATACCTTGATTAGTTTAT--TTGCGACACATATATGAT-TATATAACGAATTTTATGTTTAGT-TTTC----TT---TCA--TTCAATGA---TGAG-AGTATTTT--GGAATGTTTCTATGCACCTTT-GATGTAGCATTACATTCTTAA--TTGGGCAAACATTTCATTTCTTATCAAAGCGATACTAGAAATTCAAAATTTATAAAATAAGGAGAATTTTAGTAACAAGACCGATGC--GAGCAAC--AGAT-GGTT---GCTATTTTTGATATTTCTCACGAATATTGTAGAGTATTCGGGCGAAGACAT--CGTACAGGCAAGAATAATTGTAAGTGTTGGCCACAACAGTTTTCCTAGTGAAGTAAAGGAGAAGTTCGCA-TTTTT---TCCCCTTGA-------CACACATTCAT-AGATTACTTTACCTATCCTCT-GGA---CCATACTTCGGACCACTGTTAAGCGAAGATCTATCGGATAA-TGTTCCTGTCATTTCCCACTTGTCCCA-CGCACAAGGAAAAGGGTCCTA-AGGCGCTCGGACAACAGTTGACCGTGATCCGTCAGACTGGCTGTGCAGTGTCAATAAAGA----CCGATGCTCGTAATAGAGTGTGTCTACTT---G-----TTAGTTAATTTAGATATAAAAACAGCACCGT-GG----AGATACTGATGGGTAT--TACTAGTTATAGTATTAACTTAACA-GAAA-A-A-TA-AAATCAATTGAAC-TAT---TTCGATA
>YBR018C_Spar
TTTG-C-CAGCTTAGTATCCT-ACA--A-AAACTAT--ACAATCTATATCTTGATCAGTTTTTGATCGCAACACATATATTAGCTGTATAACTTATTTCATATTATGT-TATT----TT---TGCTTTTTAGTGA---TAAT-AGTGTCATAAGGAA--TTTC-CT-CACATTTCGGTGCGTGATTACACTCGAAAGATTGTTTGCAGGTTTTTTGTACCACGAAAGAGATATTAAAAATTCAAAGGCTATAGAAGGATGCGAACGGTAGTAACAACAATA--GCGCGAACAGCGGGGGTCCCTT--CGTTACCTTTGATATT---------GATTGCAAGGTGCTCCGGTGAAAACAT--CATAAAGG-AA----AATTGTAAAAGTTATCCACAGTA-TTCTTTTGGTGAAGGAGA-----TGGAACTA-ATTTTTCTTCCCCTTGATTAGGTACATACATCCCTCAGTTTACTTTACCTCTCCTCTTGGAAAGCTACACTTCGGAGCACTGTTGAGCGAAAGCTCATCAGGTAT-TTTTCCTGTCATTTCCCATTTGTCCCA-CAAATAAAGGAAAGGGTCCAAGAGCCGCTCGGACAACAGTCGACCGTGATCCGACAGACTGGCTGTGCAGTGTCAATAAAAAGAATCCGA-ATTGACAATAGAGTGTACCTATTT---G-----TTGATTAA--TAAATATATAAA-AGCAG-GTCGGT---AGATAATTATGGAAAT--TACTATGTAGGGCATCGACTTAAT--GAAA-A-A-AA-TAGTCTGCTGAAC-TTTC-TTTCAATA
monkey2.0/example/YBR019C_upstream.aln.fasta0000644000163100010400000000604010432332204017161 0ustar  alanall>YBR019C_Sbay
TATATGTTGTTGTT-GTTTTTT--TGG--GATAAA--------CAAATGTTGAATACTTGTTATTGATGTTTTGATAATTGGAGACAAAGTAAAAC----CGGAAATAAC-AAATCTAAGAAGTACAAACGGAATAAAGAAAAACTATACAATTGAGGCGATCTTACTTTCATATTTATATGTTTGCAAAAGGAACAGAAAACCATTACTTT--CTGCATGAGCCACTACAGAAGTAAAACGG--AAAA--CGAGTAG-GCT-TGACCGCAACCCCATCATTGTTCTATCCAGTGTGCTGTACACGGAACAA-AGTATATTGGGGAGGGGCAAGGAATTGTCATTTCACGTTCA-CCATAAC--TGCCACTACAGGCATTTGTTC-TTGCATTTGCATTGTTCGGGGTATGGCGGGGCGAGACGCTTCAGTA-ATCAGGGACGTGATCCGCACAGACGACTTCGCACGGAGGACAGTTTTCCGCCGGAGTACTGTCATCCGGTCGGCAGCTTCTGATCCGTATTTCAACATTGCGAAGTAATGGACAATAAGACAAAATAGGGAAAGCCCTACAGAA-AATCCCCTTCTTTTGT-TCGAGGTAGCACGTTTATTGGCATCTCTCTA-TAATATATAAGAAAGATCAGAAATATCTATATATGTGATAGTGA---CAGATAGCATGATGTGACTACACTTTCTTTT-ATAAT-TT-TGAAGTTAAAGAAAGAAAACTAAAGAAATCAAAAAAATAAACGTAT
>YBR019C_Scer
TATA-GT--TTTTTCT--CCTTG----ACGTTAAA-GTATAGAGGTATATTAACAATTTTTTGTTGATACTTTTATGAC----ATTTGAATAAGAAGTAATACAAACTGA-AAATGTTGAAAGTAT-----TAGTTAA-------------AGTG--GTTATGCAGCTTTTCCATTTATATATCTGTTAATAGATCA-AAAATCATCGCTTC--GCTGATTAATTACCCCAGA---AATAAGGCTAAAAAACTAATCGCATTATCATC------CTATGGTTG------TTAATTTGATTCGTTAATTTGAA-GGTTTG-TGGGGCCAGGTT----ACTGCCAATTTTTCCTCT-TCATAACCATAAAAG--CTAGTATT----G-TAGAATCTTTATTGTTCGGAGCAGTGCGGCGCGAGGCACATCTGCGTTTCAGGAACGCGACCGGTGAAGACGAGGACGCACGGAGGAGAGTCTTCCGTCGGAGGGCTGTCGCCCGCTCGGCGGCTTCTAATCCGTACTTCAATATAGC-----AATGAGCAGTTAAGCGTATTACTGAAAGTTCCAAAGAGAAG-GT---TTTTTTAGGCTAAGATAATGGGGCTCTTTACAT-T-TCCA-CAACATATAAGTAAGATTAGATATGGATATGTATATGGTGGTAATGCCATGTAATATGATT--ATTAAAC-TTCTTTGCGTCCATCCAAAAAAAAAGT-AAGAATTTTTGAAAA-TTCAATATA---------A
>YBR019C_Smik
--TA-GT--TTTTTTT---CTTGATGAACGATGTATGTATAG----ATACTGACTATTTTTTTTTGGTATTTCTATAACTCCAAAGAAAAC---TAGTACTAAATA-TAATTAGTCTCAAAGGTTCGAATGGAATGAA-------------T-TGA-TCTCTGCGACTTTTGCATTTATATGTTTGCTAATAGTCCT-AAA-CCATCATTTC--GGTGATGAATTATCACAGA---GA-AAGG-TGAGAACCTAAT-ATTTTGCTACCAAAATTCTATGCTAG------CCAGTTTTACGCATCGATTTGAA-GGTATG-TGGGGCCGGGGG----ACTGTTTATTTTGCGTTGGGCATCACCGTAGAAA--TTAGTACCTCTTTATGCCAGTATTATAGTTCGGACCACGGTGGTGCGAGCCGTACCTATG-TTCAGGGACGCGACCAACGTAGAGGACGCCGCGCGGAGGAGAATGTTCCGTCGGGGAATTGTCGTCCGTTCGGCGGCTTCTCATCCGTATTTCAATATC--------TTGAACAGCTAAACATAT-AGC-AAAGCTCCATAGA-AAGACC---TTTTTTGATCTGGGAGACTGGGATTACCAACAA-T-TCCCGTAACATATAAGTAAGACTAGACATCGGTATGTGTATGATAGTAACGAAATGTAATATTATC--TTGACA--TTTTTTACATAAC-CT-TGAAGATAGTGAAAAAAAACTGAAAAATTCAATATA---------C
>YBR019C_Spar
TATA-GT--TTGTTCTTTCCTTG----ACGATAAA-GTATAGATGTATATTAACTACTTGTTGTTAATATTTCCATGAC----AATAAAACAAAACGAAAT------TAA-AAATCTTGGAGGAGCAAATGTATTTAA-------------AGTG--GCTATGCAACTTTTGCATTTATATGTTTGCTAATAGACCA-AAAACCATTATTTTTTGCGGATGAATTGTCTCAGG---AATAAGG--AAAACCCTAACTGTACTATTACC------ATATGCTGG------CCAATTTGATTCAATAGTTTGAAAGGTATA-TGGGGTCAGTGG----GCTGCTTATTTTGTTCCTGTCATAACCATAAAAAAACAAGTATTTCTTG-TAGAATCCTTATTGTTCGGAGCAGGGCGGCGCGAGGTACATCGATG-TTAAGGGACACAACCGACGAAGACGGCGACGCACGGAGGGGAATATTCCGTCGGAGGGCTGTCGTCCGCTCGGCGGCTTCTGATCCGTACTTCAATATTGC-----AATGAGCAGTTTAGCGTAC-AGTTCAAGTTTCATAGA-AAG-GC---TTTCTTGGGTTAAGATACTGGGGTTATCAACAC-TCT-CA-CAACATATAAGTAAGACTAGATATGGATATGTGTATGATAGGAAAGCGATGTAATATAATA--ATCACAC-TTCTTTGCACCCATCT--TAATATAGTG-AAAAATTCTGAAAAGATCAACATA---------G
monkey2.0/example/YBR020W_upstream.aln.fasta0000644000163100010400000000602410432332204017177 0ustar  alanall>YBR020W_Sbay
ATACGTTTATTTTTTTGATTTCTTTAGTTTTCTTTCTTTAACTTC-AAA-ATTATAAAAGAAAGTGTAGTCACATCATGCTATCTGTCACTATCAC-ATATATAGAT--ATTTCTGATCTTTCTTATATATTATAGAGAGATGCCAATAAACGTGCTACCTCG-AACAAAAGAAGGGGATT-TTCTGTAGGGCTTTCCCTATTTTGTCTTATTGTCCATTACTTCGCAATGTTGAAATACGGATCAGAAGCTGCCGACCGGATGACAGTACTCCGGCGGAAAACTGTCCTCCGTGCGAAGTCGTCTGTGCGGATCACGTCCCTGAT-TACTGAAGCGTCTCGCCCCGCCATACCCCGAACAATGCAAATGCAAG-AACAAATGCCTGTAGT--GGCAGTTATGGTGAACG---TGAAATGACAATTCCTTGCCCCTCCCCAATATACTTT-GTTCCGTGTACAGCACACTGGATAGAACAATGATGGGGTTGCGGTCA-AG-CCTACTCG--TTTTCCGTTTTACTTCTGTAGTGGCTCATG--CAGAAAGTAATGGTTTTCTGTTCCTTTTGCAAACATATAAATATGAAAGTAAGATCGCCTCAATTGTATAGTTTTTCTTTATTCCGTTTGTACTTCTTAGA-TTTGTTATTTCCGGTTTTACTTTGTCTCCAA-------TTATCAAAACATCAATAACAAGTATTCAACATTTG-------TTTATCC--CAAAAAAACAACAACAACATATA
>YBR020W_Scer
---------TTATATTGAATTTT-CAAAAATTCT-TACTTTTTTTTTGGATGGACGCAAAGAAGTTTAATAA--TCATATTACATGGCATTACCACCATATACATATCCATATCTAATCTTACTTATATGTTGT-GGA-AATGTAAAGAGCCCCATTATCTTA-GCCTAAAAAA---ACCTTCTCTTTGGAACTTTCAGTAATACGCTTAACTGCTCATT-----GCTATATTGAAGTACGGATTAGAAGCCGCCGAGCGGGCGACAGCCCTCCGACGGAAGACTCTCCTCCGTGCGTCCTCGTCTTCACCGGTCGCGTTCCTGAAACGCAGATGTGCCTCGCGCCGCACTGCTCCGAACAATAAAGATTCTAC-AA----TACTAGCTTTT--ATGGTTATGAAG-AGGAA--AAATTGGCAGTAACC-----TGGCCCCACAAACCTT-CAAATTAACGAATCAAATTAACAACC------ATAG------GATGATAATGCGATTAGTTTTTT-AGCCTTATTTCTGGGGTAATTAATCAGCGAAGCG--ATGATTTT-TGATCTATTAACAGATATATAAATGGAAAAGCTGCATAA--CCACTT-------------TAA-----CTAATACTTTCAACA-TTTTCAGTTT---GTATTACTTCTTATTCAAA-T---GTCATAAAAGTATCAACAAAAAATTGTTAATATACCTCTATACTTTAACGT-CAAGGAGA--AAAA--AC-TATA
>YBR020W_Smik
G---------TATATTGAATTTTTCAGTTTTTTTTCACTATCTTC-AAG-GTTATGTAAAAAA-TGTCAAGA--TAATATTACATTTCGTTACTATCATACACATACCGATGTCTAGTCTTACTTATATGTTAC-GGGAATTGTTGGTAATCCCAGTCTCCCAGATCAAAAAA---GGTCT-TTCTATGGAGCTTTGC-T-ATATGTTTAGCTGTTCAA--------GATATTGAAATACGGATGAGAAGCCGCCGAACGGACGACAATTCCCCGACGGAACATTCTCCTCCGCGCGGCGTCCTCTACGTTGGTCGCGTCCCTGAA-CATAGGTACGGCTCGCACCACCGTGGTCCGAACTATAATACTGGCATAAAGAGGTACTAATTTCT--ACGGTGATGCCCAACG--CAAAATAAACAGTCCCC-----CGGCCCCACATACCTT-CAAATCGATGCGTAAAACTGGCTAGC------ATAGAATTTTGGTAGCAAAA-TATTAGGTTCTC--ACCTT-TCTCTGTGATAATTCATCACCGAAA--TGATGGTTTA--GGACTATTAGCAAACATATAAATGCAAAAGTCGCAGAG-ATCAATT-------------CAT-TCCATTCGAACCTTTGAGACTAATTA-TAT---TTAGTAC---TAGTTTTCTTTGGAGTTATAGAAATACCAAAAAAAAATAGTCAGTATCTATACATAC---ATCGTTCATCAAGAAAAAAA--AC---TA
>YBR020W_Spar
C---------TATGTTGATCTTTTCAGAATTTTT-CACTATATT--AAGATGGGTGCAAAGAAGTGTGATTA--TTATATTACATCGCTTTCCTATCATACACATATCCATATCTAGTCTTACTTATATGTTGT-GAG-AGTGTTGATAACCCCAGTATCTTA-ACCCAAGAAA---GCCT-TTCTATGAAACTTGAACT-GTACGCTAAACTGCTCATT-----GCAATATTGAAGTACGGATCAGAAGCCGCCGAGCGGACGACAGCCCTCCGACGGAATATTCCCCTCCGTGCGTCGCCGTCTTCGTCGGTTGTGTCCCTTAA-CATCGATGTACCTCGCGCCGCCCTGCTCCGAACAATAAGGATTCTAC-AAGAAATACTTGTTTTTTTATGGTTATG-AC-AGGAACAAAATAAGCAGCCCAC-----TGACCCCATATACCTTTCAAACTATTGAATCAAATTGGCCAGC------ATA------TGGTAATAGTACAGTTAGGGTTTT---CCTTATTCCTGAGACAATTCATCCGCAAAAAATAATGGTTTT-TGGTCTATTAGCAAACATATAAATGCAAAAGTTGCATAG--CCACTT-------------TAAATACATTTGCTCCTCCAAGA-TTTTTAATTT---------CGTTTTGTTTTAT-T---GTCATGGAAATATTAACAACAAGTAGTTAATATACATCTATACTTTATCGT-CAAGGAAAGAACAA--AC-TATA
monkey2.0/example/YML051W_upstream.aln.fasta0000644000163100010400000000240410432332204017206 0ustar  alanall>YML051W_Sbay
-AGCGAGGACTCTTTCAT--TTTGTAATATATATTTATACCCTCCCCATCCCTTCC---ACAATATAATAGTATAATTCTAATATTAATAATA---TCCTATATTTTTTTCACC-TCCGGAGCACTCCCGTCCGAAAAACTGCGAAAATGTTTGCTTGGTATACATAATACTCAAAAAT---TCACTATTTTCT------TGCATGTAGAAACACACTTTACATCATTGTTGGCCACTGCAGACCACTCC--ACCTCCGTATACAAGTCTGACCGCTAGTTTTTCCCGTGCC-TTCC-----CGCC
>YML051W_Scer
ATGGCGCAAGTTTTCCGC--TTTGTAATATATATTTATACCCCTTTCTTCTCTCCCCTG-CAATATAATAGTTTAATTCTAATATTAATAATA---TCCTATATTTTCTTCATTTACCGGCGCACTCTCGCCCGAACGACCTC-AAAATGTCTGCT--ACATTCATAATAACCAAAAGCTCATAACT-TTTTTTTT----TGAACCTGAATATA---TATACATCA--CAT-ATCACTGCTGGTCCTTGCCGACCAGCGTATACAATCTCGATAGTTGGT--TTCCCGTTCT-TTCCACTCCCGTC
>YML051W_Smik
TAGGGACAAAG--CCCGCCTTTTGTAATATATACTTATACCCTCTCCTTCTCTCCCCTG-CAACATAATAGTTAACTCCTAATATTAATAATAATATCCTACAATTTCCTTAGC-ACCGGGGCACTCTCGCCCGAACGACCTC-AAAACGCTTGCT--ACATCCATAATATTCAGAACTACATCACT-TTTTTTTT-----GTACATAAAAATA---TATAC--CA--CAT-GTCACTGCTGATCCTTGCTGACGAGCGTATACAAGCTCGATAGCTGGTCTTTACCGTGCCATTCC-CTGCCGTC
>YML051W_Spar
AGGGGCCAAAGCTCCCGC--TCTGTAAAATATATTTATATCCCTTCCTTCTCTCCCCTG-CAATATAATAGTTTAATTCTAATATTAATAATA---TCCTATATTTTCCTTACC-ACCGGCGCACTCTCGCCCGAACGACCTC-AAAATGCTTGCT--ACATTCATAATAATCAAAAGCTTATAACT-TTTTTTTTCCTTTGTACCTGAATATA---TATACATCT--CAT-GTCACTGCTGGTCCTTGCCGGCCAGCGTATACAACCTCGATAGCTGGT-TTTCCCGTTCT-TCCCACTCCTGTC
monkey2.0/example/YOR120W_upstream.aln.fasta0000644000163100010400000000757010432332204017224 0ustar  alanall>YOR120W_Sbay
-GACGTT--AGAAATATCGAGCGTACAAGTTTATAA-TATGTCTCA-GCTATCTTTCTCC-ATTTTGGACAAG----G-----TACCTCATCTCATCTTAA-TTCTCTATTTTTCACCATAAAATTTT-CAAATTTTTCGTCGCGCGATGATG-TCTCTTAGAAAGTACGTAATGAATATTACGAGTATG----AATATGAAAG----ATACAAAATTGCTATTCT-GAACCTTT--TTCTAGACCGATGAAGAACTTGCATTACAGCAAAAGCGG-ATTGCCAATAGAAAGAGAAGATAGAACGCGTCTACCCCTCTGTTAGTACCCCATCACCGTT----TTTTTTTTACTTTTTTCGCTTACGCCGGTTGAAGCTTCACGG-GTCCTCTTAGAAG--CTGACATAGATG--GCGCACACGAGTGAGCTGTCATATCACCTGTCTCAAGGCTTTTGAAACACT-CACTGGTGTTATTACTCTTCAAA-----AGTGAAAAATTC-CATAGATTTATAGATTTTAACACCGACCGAAAATCGGGGCAAAGTATTCCGGGAAAGGGAAAGGATTGCCAAGGGCGGTCTTTTGT--TATTTCTCAAAGTCTT---GGCAAGATTGGCATGCTTTCACCATAAGGGGGTAA----AACTTGCCCAGCACGTAAACTGAAATGGCGACTTACTAATATGCACACATGGACGTGGCATCGTGTGGGTGGGTGATT-GCAAACGCCCGGTCTTGGTCTCCCTGGTTCTCTCG--AGAGGTAACCGAACCTCTAGGGAAAA-AGTCGTCTGACAATAAAGGAAATTAGATACTGGTATCCCTTTGTTCGTCGTATAAAAGCGGACACATTAAAGCTTGATCATGTCTCTTAGCTTC---TTGCCTGTAGAATACACTGCAACTC----AAGCTAGAAGTTAATATAGC--ACTATCACAAAACT-------TAGAAAA---
>YOR120W_Scer
-GACCTTCTAAAAATAACG-GAGTTTGAATCTCTAAATCTGTTTTAA-CT-TCTTTTTACTA-TTATTTTTAGTC--TT--AGTATCTCATCTCATCTCAA-TT-TCTATATTCCACTATAAAATTTTTCACTCTTTCTG-CGCGCGCCAATG-TCCCC-GCA-ACTACTCAATAGGTAACATGAGAATATTTCAGTTCGTAAGAGAGAAGAGATGAAGTTATT-TGGGCTCTTTG-CTCGAGGTTACAGAAGGGC-CGCATTAGAGTGAATGAGCTG---ATGATATTTCGCCCAGTTCTACATT-TTTTTTTTTTTGGAAGTATGAC--CTCTGTT-AAATTTTTTTTTTTTT---AAATTTCACTTTCTAAAGTC-CCAGAAATCCGCTTGAATG-TCTTACATAT-TGCAATGGATATGCTTGGGTGATCATACTTCCTGGCTTTAGATATTTGAA-A-CTTAACTCTTGTCAACAAACTTCCTATGGA--GTGTATA--------AGAATTGTAAGTTATAACACCGGCGAACAATCGGGGCAGACTATTCCGGGGAAGAACAAGG------AAGGGCGGTCTTTTC-----TCCCTCATTGTCAT---AGCAAGGTCATTTCGCCTTCTC-AGAAAGGGGTAG----AATCAATCTAGCACGCAGATTGCAAACACGGCTTAATAATATGCCTAT-CAGGCATTCACCCGTGTGA-----CGAATCGCACAC------CGCTGCTCTCCTTAATTCCCTAG--AGTAGAAACCGAGCTTTCAGG-AAAAGACT---ACGGCAGTAAA-G-A-------------ATTGCTTTACTGGGCGTATAAAACCGGGAGAATCAAGACATTCTAATGACT-TGATT-CAGGATGA-GAGCTTAATAGG-TGCATCTTAGC-AAGCTAAAA-TTTGGACAGC-TCTCATTACTAAATT---AAGATAGAAAA---
>YOR120W_Smik
-CAAATTAAAAAAATAACGAGCTTTCTAGATTATAATTATCGTCCAAGCCTTCTCTTAT-TAATTCGGATAAGAA---AATTTCACCTCATCTCATCTCTAGTT-T-TAATCTTTTCTATAAAATTTT-CAATCTTTTTG-CGCGCGCTAATG-ACCC--GCTAACTACGTTATAAGTAATA-----ATG-----AAAAGAA-------------GACGCTGTTCT-AACGCTTTGACCCAAGGCAGCGGAATAGC-AGTATTACAGAAAACGAGATGAAAGCTATTT---------GACAAAATTTATATCTCTCTCGACAATACTAATAAAC---------GTTCTTTTCCCT---GAATTTCAATGATGAAAGGC-CCAGAAGTCACCTTTGAAATTCAAACACAT--GAGGTACCTACAATTTAATGATCATATCATCAACTTATTGACTATTGAGAAGC-TCACTGATTTTT-TTGGCTTCCTATAGT--GTGT-------GCA-AGAACTAGATATCATAACACCGACTAGAAATCGGGGCATACTATTCCGGCAAAGAATAAAG------AAGGTCGGTCTTTT---T---CGCTCAATGTCATCATGGCAAGATCGTTTTGATTTCCC-CGAAGGGGGTAG----AACAAAGCCAGCACGTAAAACGCACTCATAACTTGGTAATACTCCGAT-GAGACGTAGACTCGTGTGATCAGATCAGTTACAGAC------TATCGCTCTCCTCGATTCCTTCA--TGTAGTTCCCGAGTTTGAAGA-AAAAAAATATGGAGATA-TAAA-G-A------------TATCGCTTTGCCGGATGAATAAAACTAGAGGAATCAAGGCATTCCAAAGA-TCTTTTTTTCTAGTGAAGTGTTTTATTTA-TGTATCTCAGTGAAGCCAAAA-TTTAGATAGC-TTTTATTGCTAGACT----AAGAAGTAAA---
>YOR120W_Spar
----CTT--CAAAATAAC-AGAGTTCGATTCTATAAATTTGCTTCAAG-TTTCCTTTTACTAATTTTTTCAAGTAATGTATCTTATCTCATCTCATCTCAA-TT-TCAAC-TTTCACTATAAAATTTTTCA-GGTTTCTG-CGCGCGCTAATG-TCCCCCGCATACTGCTTGATAAGTAATATGAGACTATTTCAGTTTGTTAGGGAAAAGAAAAGATGCTGTTCG-GACTCTTTG-CTCAAGATAGCAAAAAGGC-TGTATTAGAAAGAATTAGA--TGAT-GATATCTCGATCAGGTCCAGAT------TTCTTTTGAAAATAT--CAACTCTATTACAATATTTTCTCTTTTCCTAACTTTCACTTACTAAAATA-CTAGAAG-CCCCTTCGATG-TCTAACATGT--GTAATGTATATGATTGAGTGATTATACTGCGAGTCTTAAG-CGCTTGAAGAACTTATCTCTTGTTA----GCTCCCTACGGA--GTGCA-------CA-AGAATTGTAGGTTATAACACCGGCTAAAAATCGGGGCAGACAATTCCGGGGAAGAACAAAG------AAGGGCGGCCTTTTTT---ATCTCTAATTGTCAT---AGCAAGATCGTTTCGCTTTCTCCAGAAGGGGGTAG----AACCAATCTAGCACGCACACTGCAAGCACGGCTTAGTAATATGCCTATTGAGGCATAAACTCGTGTGA-----CGAATTGTAAAC------CGCCGCTCTCCTTAATTCTTTAG--AGTAACAACCGAGCTTTTAGG-AAAAAAATCGT--GACAGTAAG-G-A-------------ATTGCTTTACTGGACGTATAAAACAGGGAGAATCAAGACGTTCCAAAGGCT-TGATCGTCG-TTGA-AGGTTTAGCAGG-TACATTTCAGC-AAACTAAT--TTCAGATAACCTTTCCGTACTAAATC---AAAATAAAAAAAAA
monkey2.0/example/YBR018C_upstream.aln.tree0000644000163100010400000000014510432332204017021 0ustar  alanallt 7
r 6	n 5	l 0
n 5	n 4	l 2
n 4	l 1	l 3
b
4	0.08485	0.17257	
5	0.04130	0.21132	
6	0.182065	0.182065	
monkey2.0/example/YBR019C_upstream.aln.tree0000644000163100010400000000014310432332204017020 0ustar  alanallt 7
r 6	n 5	l 0
n 5	n 4	l 2
n 4	l 1	l 3
b
4	0.08857	0.16250	
5	0.06602	0.19358	
6	0.21246	0.21246	
monkey2.0/example/YBR020W_upstream.aln.tree0000644000163100010400000000014510432332204017036 0ustar  alanallt 7
r 6	n 5	l 0
n 5	n 4	l 2
n 4	l 1	l 3
b
4	0.08296	0.16779	
5	0.07998	0.17453	
6	0.193335	0.193335	
monkey2.0/example/YML051W_upstream.aln.tree0000644000163100010400000000014310432332204017045 0ustar  alanallt 7
r 6	n 5	l 0
n 5	n 4	l 2
n 4	l 1	l 3
b
4	0.03899	0.06740	
5	0.06264	0.07309	
6	0.13631	0.13631	
monkey2.0/example/YOR120W_upstream.aln.tree0000644000163100010400000000014510432332204017054 0ustar  alanallt 7
r 6	n 5	l 0
n 5	n 4	l 2
n 4	l 1	l 3
b
4	0.09787	0.13253	
5	0.12665	0.23934	
6	0.154675	0.154675	
monkey2.0/example/YBR018C_upstream.aln.occ0000644000163100010400000000314510466224604016644 0ustar  alanallNAME	STRAND	POS(0-based)	UG(1-based)	YBR018C_Sbay	p(LR>lr)	POS	YBR018C_Scer	p(LR>lr)	POS	YBR018C_Smik	p(LR>lr)	POS	YBR018C_Spar	p(LR>lr)	p(S^>s)	T	p(T<t|HB)	p(T>t|HKY)	
example/YBR018C_upstream.aln gal4sites.matrix	0.333333	290-311	275	CCGAGGAATTTAACGCT	0.00542566	290-309	CCGCGGAGTTCATTTCG	0.000844666	284-309	CGAGCAACAGATGGTTG	0.00821303	286-303	CGAACAGCGGGGGTCCC	0.00052041	0.00616568	-2.13285	0.00998625	0.249663
example/YBR018C_upstream.aln gal4sites.matrix	0.5	438-458	383	CCCCTGAACCATGTGCG	0.00921064	423-445	GGGGGTAATTTTTCCCC	0.00237353	427-448	CGCATTTTTTCCCCTTG	0.0353334	451-468	AGGTACATACATCCCTC	0.0374886	0.0238562	-3.60199	0.00276398	0.439866
example/YBR018C_upstream.aln gal4sites.matrix	0.333333	504-521	442	CGGAGCACTGTTGGACG	1.37785e-06	504-521	CGGAGCACTGTTGAGCG	5.2603e-07	504-521	CGGACCACTGTTAAGCG	6.92649e-06	504-521	CGGAGCACTGTTGAGCG	5.2603e-07	5.14507e-09	5.75389	0.712073	0.00029392
example/YBR018C_upstream.aln gal4sites.matrix	0.333333	530-548	468	CGGATCATTTTCCTGTT	0.00511583	528-546	ATTAGATATATTTTCTG	0.26698	530-548	CGGATAATGTTCCTGTC	0.0150272	543-560	CTGTCATTTCCCATTTG	0.074529	0.0664249	-7.74165	3.39306e-05	0.912812
example/YBR018C_upstream.aln gal4sites.matrix	0.5	576-593	513	AAGGGGGCCCTTGACCG	0.000638099	572-589	AGGGAAAGGGTCCAAAA	0.057172	572-590	AGGAAAAGGGTCCTAAG	0.00472359	578-595	AGGGTCCAAGAGCCGCT	0.0121489	0.00417077	-4.01732	0.00185818	0.500448
example/YBR018C_upstream.aln gal4sites.matrix	0.666667	595-612	532	CGGACAGCAGTCGACCG	1.15123e-07	595-612	CGGACAACTGTTGACCG	3.29799e-07	595-612	CGGACAACAGTTGACCG	3.5145e-07	595-612	CGGACAACAGTCGACCG	5.23295e-08	1.64032e-10	6.50215	0.816749	0.000112665
monkey2.0/gal4sites.matrix0000644000163100010400000000130210432347252014156 0ustar  alanall0.041666667	0.875	0.041666667	0.041666667
0.041666667	0.041666667	0.875	0.041666667
0.041666667	0.125	0.791666667	0.041666667
0.458333333	0.166666667	0.25	0.125
0.125	0.333333333	0.458333333	0.083333333
0.166666667	0.291666667	0.375	0.166666667
0.625	0.083333333	0.208333333	0.083333333
0.125	0.458333333	0.333333333	0.083333333
0.571428571	0.071428571	0.071428571	0.285714286
0.083333333	0.333333333	0.458333333	0.125
0.083333333	0.208333333	0.083333333	0.625
0.166666667	0.375	0.291666667	0.166666667
0.083333333	0.458333333	0.333333333	0.125
0.125	0.25	0.166666667	0.458333333
0.041666667	0.791666667	0.125	0.041666667
0.041666667	0.875	0.041666667	0.041666667
0.041666667	0.041666667	0.875	0.041666667
monkey2.0/Makefile0000644000163100010400000000034510432335500012470 0ustar  alanallSRC_monkey = src/monkey1.0.cpp
SRC_rmonkey = src/rmonkey1.0.cpp
PROG = monkey rmonkey
CC = g++

all: ${PROG}

monkey : ${SRC_monkey} 
	${CC} -o Monkey ${SRC_monkey}



rmonkey : ${SRC_rmonkey} 
	${CC} -o rMonkey ${SRC_rmonkey}


monkey2.0/scripts/0000755000163100010400000000000010457514562012533 5ustar  alanallmonkey2.0/scripts/clustal2fasta.pl0000755000163100010400000000077010432332204015627 0ustar  alanall#!/usr/bin/perl
die "usage: [.aln file] ...\n" unless (@ARGV);
while (@ARGV) {
my $file=shift (@ARGV);
my %ALN;
open (IN, "<", $file) or die "$file not found";
while (my $line=<IN>) { #print $line;
	next if ($line=~m/CLUSTAL/);
	if ($line =~ m/(\S+)\s+([-A-Z]+)/) {
		#print "$1\t$2\n";
		chomp $line;
		my $sp=$1;
		my $seq=$2;
		$ALN{$sp} .= $seq;
	}
	#elsif ($line) { print $line; }
}
close IN;

open (OUT, ">", $file.".fasta");
for my $sp (keys %ALN) {
	print OUT ">$sp\n$ALN{$sp}\n";
}
close OUT;
}
monkey2.0/scripts/dnd2tree.pl0000755000163100010400000000561010432332204014564 0ustar  alanall#!/usr/bin/perl
die "usage: [.dnd file] ...\n" unless (@ARGV);
while (@ARGV) {
my $file = shift @ARGV;
my $treestring;
open (FILE, "<", $file);
while (my $line = <FILE>) {
	$treestring=$line if ($line =~m/;/);
}
close FILE;

my %tree;
die "couldn't parse $treestring\n" unless (parse_rooted_tree_with_branchlengths ($treestring, \%tree));

print tree_string(\%tree,root(\%tree)), " found\n";

my $treefilename=$file.".tree";
print "wrote $treefilename\n" unless (print_tree_file(\%tree,$treefilename));
}

sub print_tree_file {
	my ($t, $fname) = @_;
	open (OUT, ">", $fname);
	print OUT "t ";
	print OUT $$t{Nnodes};
	print OUT "\n";
	
	for (my $n=$$t{Nnodes}-1;$n>=0;$n--) {
		next if ($$t{$n}{leaf}==1);
		if ($$t{$n}{root}==1) { print OUT "r $n"; }
		else { print OUT "n $n"; }
			
		for my $br ("left","right") { 
			my $chn=$$t{$n}{$br};
			if ($$t{$chn}{leaf}==1) { print OUT "\tl $chn"; }
			else { print OUT "\tn $chn"; }
		}
		print OUT "\n";
	}
	print OUT "b\n";
	for (my $n=0;$n<$$t{Nnodes};$n++) {
		next if ($$t{$n}{leaf}==1);
		print OUT "$n ";			
		for my $br ("left","right") { 
			my $chn=$$t{$n}{$br};
			print OUT "\t".$$t{$chn}{branch};
		}
		print OUT "\n";
	}
	close OUT;		
	return 1;
}

sub parse_rooted_tree_with_branchlengths {

	my ($str, $t) = @_; 
	#print $str, "\n";
	$str =~ s/\s//g;
	my $nn=0;
	while ($str =~ m/(\w+):(\d+\.\d+)/g) {
		#print "leaf $1, branch $2\n";
		$$t{$nn}{"name"}=$1;
		$$t{$nn}{"branch"}=$2;
		$$t{$nn}{"leaf"}=1;		
		$nn++;
	}
	$$t{Nnodes} = 2*$nn - 1;
	#print keys %{$t}, "\n";
	my $i=0;
	#print "$nn leafs found\n";
	while ($nn < $$t{Nnodes}) {
		LOOP: for my $l (0 .. $nn-1) {
			next if ($$t{$l}{anc} ne "");
			for my $r ( 0 .. $nn-1) {
				next if ($r == $l);
				next if ($$t{$r}{anc} ne "");
				my $test = "(".$$t{$l}{"name"}.":".$$t{$l}{"branch"}.",".$$t{$r}{"name"}.":".$$t{$r}{"branch"}.")";
				#print "$test\n";
				
				for my $p (0 .. length($str) - length($test) -1) {
					my $match = substr($str,$p,length($test));
					if ($match eq $test) {
						#print "$test $match\n$str\n";
						$$t{$nn}{"name"} = $test;
						$$t{$nn}{left}=$l; $$t{$nn}{right}=$r;
						$$t{$l}{anc}=$nn; $$t{$r}{anc}=$nn;
						my $brstr = substr($str,$p+length($test),length($str)-length($test)-$p);
						$$t{$nn}{branch}=$1 if ($brstr =~ m/:(\d+\.\d+)/);	
						if ($brstr eq ";") {
							$$t{$nn}{root}=1;	
						}
						#print $$t{$nn}{branch}, "\n";
						$nn++;
						last LOOP;
					}	
				}
			}
		}
		$i++;
		if ($i>100) {
			print "couldn't parse\n"; return 0;
		} 		
	}
	return 1;
}


sub root {
	my ($t) = (@_);
	for my $n ( 0 .. $$t{Nnodes} - 1 ) {
		return $n if ($$t{$n}{root}==1);
	}
}

sub tree_string {
	my ($t, $n) = @_;
	
	return $$t{$n}{"name"} if ($$t{$n}{"leaf"}==1);
	
	return "(".tree_string($t,$$t{$n}{left}).",".tree_string($t,$$t{$n}{right}).")";
}
monkey2.0/scripts/nhx2tree.pl0000755000163100010400000001641010432335654014630 0ustar  alanall#!/usr/bin/perl

#(c) 2006 by Alan Moses 
my $MIN_BRANCH=0.001;

die "usage: [NHX file]\nconverts nhx format to the monkey tree format \n" unless (@ARGV);
TREE: while (@ARGV) {
my $file = shift @ARGV;

my $trstr=""; 
open (FILE,"<", $file) or die "$file not found \n";
while (my $line=<FILE>) { chomp $line; $trstr .= $line; }
close FILE;
$trstr=~s/\s//g;
my %t; 
unless (parse_nhx2($trstr,\%t)) { print STDERR "$file not parsed\n"; next TREE; }
print STDERR "parsed: ".tree_string(\%t,get_root(\%t)), "\n";
my $trif=0; my $branchless=0;
for my $n (keys %{$t{node}}) {
	$trif++ if (scalar(@{$t{child}{$n}})>2);
	next if ($t{root}{$n}==1);
	if (($t{br}{$n} eq "") || ($t{br}{$n}<=0)) { $branchless++; $t{br}{$n}=$MIN_BRANCH; }
}
if ($trif>0) { print STDERR $trif. " multifucations (>2 leafs per node) found ...\n";	next TREE; }
if ($branchless>0) { 
	print STDERR $branchless. " nodes had no branches. (set to $MIN_BRANCH)\n"; 
}

print STDERR $file.".tree written \n" if (print_tree_file(\%t, $file.".tree"));
}

sub print_tree_file {
	my ($t, $fname) = @_;
	open (OUT, ">", $fname);
	print OUT "t ";
	print OUT $$t{Nnodes};
	print OUT "\n";
	
	for (my $n=$$t{Nnodes}-1;$n>=0;$n--) {
		next if ($$t{leaf}{$n}==1);
		if ($$t{root}{$n}==1) { print OUT "r $n"; }
		else { print OUT "n $n"; }
			
		for my $br ("1","0") { 
			my $chn=$$t{child}{$n}[$br];
			if ($$t{leaf}{$chn}==1) { print OUT "\tl $chn"; }
			else { print OUT "\tn $chn"; }
		}
		print OUT "\n";
	}
	print OUT "b\n";
	for (my $n=0;$n<$$t{Nnodes};$n++) {
		next if ($$t{leaf}{$n}==1);
		print OUT "$n ";			
		for my $br ("1","0") { 
			my $chn=$$t{child}{$n}[$br];
			print OUT "\t".$$t{br}{$chn};
		}
		print OUT "\n";
	}
	close OUT;		
	return 1;
}


sub parse_nhx2 {

	my ($str,$t) = @_; my $PMESS=0;
	return 1 unless (($str =~ m/\,/)&&($str =~ m/;/)); #trees must have more than one node.
	warn $str if ($PMESS); 
	$$t{Nnodes}=0;
	#find the leaves - the must have either a '(' or ',' in front and a ')' or ',' behind, and contain no '(', 
							#but essentially everythign else
	my $stuff = '[a-zA-Z0-9\.\=\_\:\s\&\$\[\]\-\{\}\/\%\!\'\+\?]';	
	my %NH;
	#my $stuff2 = '[a-zA-Z0-9\.\=\_\:\s\&\$\[\]\-\{\}\/\,]';					
	while ($str =~m/[\(\,]($stuff+)[\,\)]/g) { #find the leaves		
		my $ln = $$t{Nnodes};
		$$t{leaf}{$ln} = 1;
		create_node($1,$t);
		$NH{$$t{node}{$ln}}=$ln;
		pos($str) = pos($str)-1; # leaf match can overlap by 1
		print "$ln $1: ".$$t{name}{$ln}." ".$$t{br}{$ln}."\n" if ($PMESS);
	}
	#die;
	my $iter=1; my $max_ITER=1000; 
	while ($iter<$max_ITER) {
		
		my $nodesearch=0;
		NODE: for my $n (sort numerically keys %{$$t{node}}) { # find ancestors - nodes that contain this one
			next if (($$t{anc}{$n} ne "")||($$t{root}{$n}==1));
			my $nstr=$$t{node}{$n};
			#print "looking for ancestor for $n: \n".$$t{node}{$n}."\n...\n";
			$nodesearch++;
			$nstr =~ s/\W/\./g; #print "$nstr\n";
			next unless ($str =~ m/[\(\,]$nstr[\,\)]/);
			my $unresolved=1; my @ch = ($n); my $looptime=0;
			while ($unresolved) {
				$looptime++;
				last if ($looptime>10); #maximum number of daughters for each node
				if ($str =~ m/\($nstr\,(.+)/) {
					my $rest=$1;
					print "looking for $n rest=".substr($rest,0,10)."...\n" if ($PMESS);				
					my $nextnode = next_node_right($rest);
					my $nnstr = $nextnode;	$nnstr =~ s/\W/./g;
					print "next node is $nextnode\n" if ($PMESS);
					die "problem with $nextnode\n" unless ($nextnode);
					$nextnode .=$1 if ($str=~m/\($nstr\,$nnstr($stuff*)/);	
					
					$nnstr= $nextnode; $nnstr =~ s/\W/./g;				
					
					my $rn; #check the new sister
					if ($NH{$nextnode}) {
						$rn = $NH{$nextnode};
						print "old node $rn\n" if ($PMESS);
					}	
					else {
						$rn= $$t{Nnodes};
						#$nextnode = $1 if ($str =~m/$nextnode
						create_node($nextnode,$t);
						$NH{$nextnode}=$rn;
						print "new node $rn\n" if ($PMESS);
					}
					push @ch, $rn;
					#check the ancestor
					$nstr = $nstr."\,".$nnstr;
					if ($str =~ m/(\($nstr\)$stuff*)/) {
						$unresolved=0; my $nn; my $newn=$1;
						if ($NH{$newn}) {
							$nn = $NH{$newn};							
						}
						else {	
							$nn= $$t{Nnodes};
							create_node($newn,$t);
							$NH{$newn}=$nn;
						}	
						for my $m (@ch) { $$t{anc}{$m}=$nn; }
						print "$iter : resolved $nn (@ch) $1\n" if ($PMESS);
						if ($$t{node}{$nn}.";" eq $str) {
							print "root found\n" if ($PMESS);
							$$t{root}{$nn}=1;
						}
						
					}
					else { print "$nstr unresolved\n" if ($PMESS); }
				}
				elsif ($str =~ m/\,$nstr[\,\)](.+)/) { $unresolved=0; }
				#elsif ($$t{node}{$n}.";" eq $str) { $$t{root}{$n}=1; $unresolved=0; }
				else { print "node$n $nstr not found\n" if ($PMESS); }
			}			
			
		}
		#die;
		$iter++;
		last if ($nodesearch==0);
	}
	my $rn = get_root($t);
	if ($rn ne "") {
		my $alldef=1;
		for my $n (keys %{$$t{node}}) {
			next if ($n ==$rn);
			if ($$t{anc}{$n} eq "") { 
				$alldef=0;
				print "no anc for $n\n" if ($PMESS);
			}
		}
		if ($alldef) {
			find_children($t);
				#print STDERR "parsed ...\n";
			return 1;
		}		
	}			
	return 0;
}

sub next_node_right {
	my $stuff = '[a-zA-Z0-9\.\=\_\:\s\&\$\[\]\-\{\}\'\/]';
	my ($rstr) = @_;
	my $lb=0; my $rb=0; my $pos;
	while ($rstr =~ m/([\(\)\,])/g) {
		my $brac=$1; #my $ext=$2;
		if ($brac eq '(') { $lb++; }
		elsif ($brac eq ')') {
			$pos=pos($rstr);			
			return substr($rstr,0,$pos-1) if ($lb==0);						 
			$rb++;
			return substr($rstr,0,$pos) if ($lb==$rb);			
		}
		elsif ($brac eq ',') {
			$pos=pos($rstr);
			#warn "middle node in multif ..." if ($lb==0);			
			return substr($rstr,0,$pos-1) if ($lb==0);
		}
	}
	warn "node find failed $pos and $rstr";
	return "";			
}


sub tree_string {
	my ($t,$n)=@_;
	#print $n, "\n";
	if ($$t{leaf}{$n}) {
		#print $$t{name}{$n}."\n";
		return ($$t{br}{$n} ne "") ? $$t{name}{$n}.":".$$t{br}{$n} : $$t{name}{$n};
	}
	else {
		my $str="(";
		for my $ch (@{$$t{child}{$n}}) {
			$str.="," unless (length($str)==1);
			my $temp=tree_string($t,$ch);
			$str.=$temp;
			#print "$n $ch $str ".$$t{br}{$n}."\n";
		}
		return $str.")".$$t{name}{$n}.":".$$t{br}{$n} if ($$t{br}{$n} ne "");
		return $str.")".$$t{name}{$n}.";" if ($$t{root}{$n} ==1 );
		return $str.")".$$t{name}{$n};
	}			
	
}

sub find_children {
	my ($t)=@_;
	#print "findign childern...\n";
	for my $n (keys %{$$t{node}}) { next if ($$t{leaf}{$n}==1);
		for my $m (keys %{$$t{node}}) { 
			push @{$$t{child}{$n}}, $m if ($$t{anc}{$m} ==$n)
		}
		print "$n ".$$t{name}{$n}." ".$$t{anc}{$n}." has no children \n" unless (@{$$t{child}{$n}});
		#for my $c (@{$$t{child}{$n}}) { print $c. " "; }
		#print "\n";
	}
	return 1;		
}

sub get_root {
	my ($t)=@_;
	return "" unless (scalar (keys %{$$t{root}})>0);
	for my $n (keys %{$$t{root}}) {
		return $n if ($$t{root}{$n}==1) ;		
	}
	return "";	
}

sub create_node {
	my ($nstr,$t)=@_;
	my $stuff = '[a-zA-Z0-9\.\=\_\:\s\&\$\[\]\-\{\}\'\/]';
	my $tagstuff='[a-zA-Z0-9\.\=\_\:\s\&\$\-\{\}\/]';
	my $num = $$t{Nnodes};
	$$t{node}{$num} = $nstr;
	$$t{tag}{$num}=$1 if ($nstr=~m/\[($tagstuff+)\]\Z/);	
	$nstr =~ s/\[.+?\]//g;
	#print $nstr."\n";
	$$t{name}{$num} = $1 if ( $nstr =~ m/($stuff+)\:[\d\.]+\Z/);
	$$t{br}{$num} = $1 if ($nstr =~ m/\S+\:([\d\.]+)\Z/);	
	$$t{name}{$num} = $nstr if (($$t{name}{$num} eq "")&&($$t{leaf}{$num}==1));
	#if ($$t{br}{$num} eq "") {
	#	warn "no branch found for $nstr\n" if ($PMESS);
	#}	
	#print "$num ".
	$$t{Nnodes}++;
	return 1;
}

sub numerically { $a <=> $b }
monkey2.0/src/0000755000163100010400000000000010457514562011633 5ustar  alanallmonkey2.0/src/monkey1.0.cpp0000644000163100010400000016010610432332204014044 0ustar  alanallusing namespace std;

#include <iostream.h>
#include <fstream.h>
#include <math.h>
//#include <time.h>
//#include <stdio.h>
#include <stdlib.h>
//#include <map>
#include <string>

const int MAX_ALIGN_SIZE=15; //maximum number of species
const int MAX_ALIGNMENTS = 10000;
//const int MAX_W=20; // biggest motif size
double INITIAL_BRANCH=0.1;
bool PRINT_TO_SCREEN=1;
const bool BOTH_STRANDS=1;

const double MOTIF_PRIOR=0.001;
double thr=0.0;

const double MIN_RATE=0.001;
const int MAX_HITS=10000;
const int EVOL_MODEL=2; // 0= JC bg rate,  1= JC 1/2 bg rate 2=HB selection
const int N_MODELS=5;
double FUDGE=1.0; //correction because the HB model is based on the synonymous rate
double N_BINS=200.0;

struct sequence { // sequence object - 2 copies of the sequence, a string, and as ints
	int length;
	string name;
	string bases;
	int * asints;

	sequence();
	bool integerize();
	void erase();
	bool print(int start, int stop);
	bool is_base(int pos);
	//double Pr(int start, int stop, double * f);
};

struct node {
	bool leaf; bool root;
	struct node * anc;
	struct node * left;
	struct node * right;
	double abr,lbr,rbr;
	//double * seq_post;
	int base; int number;

	node();
	bool print();
};

struct tree { //hold a list of nodes, a pointer to the root and a some substitution matrices
	struct node * nodelist;
	int Nnodes;
	struct node * root;
	double ** psubs[4][4];

	tree();
	bool read_from_file(string treefile);
	bool initialize_pairwise(double dis);
	bool print(struct sequence * seq);
	bool compute_PJC();
	bool compute_PHKY(double * f, double kap);
	string tree_string( struct node * n, struct sequence * seq);
};

struct alignment { // alignment - holds some sequences, positions of the gaps. tree
	int length;
	int Nseqs;
	int gapped;
	int bad;
	int leftend;
	int rightend;
	
	string name;
	struct sequence * seq;
	int * nextgap;
	int * nextbad;
	long int * asints;
	struct tree phyl;
	double treelength;
	
	alignment();
	bool read_bg_freq(string freqfile);
	bool read_fasta(string filename);
	bool print();
	float div();
	
	bool compute_nextgap();
	void integerize();
	
	double * bgfreq;
	double kappa;
	
	bool read_tree(string filename);
	
};


struct motif { // holds a pwm and vector of evol. rates
	int length;
	string name;
	double ** f;
	double pri;
	double * K;
	double br;
	motif * rc;
	double ***** psubs;
	int *** pssm;
	int * minscr;	
	double ** pvalue;
	//double ** score_to_LR;
	//struct tree * phyl;

	motif();
	//bool intialize_from_word(int * word, int start, int stop, double pseudocounts);
	bool initialize_neutral_motif(int w, double pseudocounts);
	bool read_from_file(string mfile);
	bool print();
	double halpern_bruno(int pos);
	struct motif make_rc();
	//bool optimize_by_EM( struct SequenceData * sd,bool update_prior, bool update_matrix);
	double Pr(int s[], int start);
};

double p_subtree2 (struct sequence * seq, struct node * n, struct tree * phy, int base, int pos);
double p_subtree_hb (struct sequence * seq, struct node * n, int base, int sp, int mp, struct motif * m);
double p_subtree_jc (struct sequence * seq, struct node * n,  int base, int pos, double k);

double PJC (int a, int b, double alpha);


int gap_count (int * seq, int start, int stop);
struct sequence remove_gaps (struct sequence * seq, int start, int size);
double pw_pid (struct sequence * seq1,struct sequence * seq2, int start1, int start2, int length);

int main(int argc, char * argv[]) {	
	if (argc<=1) { cout << "usage: [pwm file] [fasta alignment] (options) ... or try -help for more details" <<endl; return 0; }
	int read_alignment_list(struct alignment * ali, string file);
	void help();
	//cout << argc <<endl;
	int ref_num=0,mmod=2,bgmod=0, Nalign;
	
	string freqfile = "?";
	string mfile = "?";
	string treefile = "?";
	string fastafile = "?";
	string outfile = "?";
	
	bool HEURISTIC=1; bool PVALS=1 ; bool LLRS=0; bool LIST=0; bool EVALS=0; 
	double KAPPA = 2.0; bool DAN=0;
	
	struct alignment * ALI = new alignment [MAX_ALIGNMENTS];	
	struct motif * M = new motif [2]; 
	
	int i=1;
	while (i<argc) {
		//cout << argv[i] <<endl;
		if (argv[i][0] == '-') {
			if (argv[i][1] == 'm') {
				if (argv[i+1][0] == 'm') mmod=0;
				else if (argv[i+1][0] == 'J') mmod=1;
				else if (argv[i+1][0] == 'H') mmod=2;
				else if (argv[i+1][0] == 's') mmod=3;
				else { cout << "didn't recognize "<< argv[i+1] <<endl; help(); return 0; }
				i++;
			}	
			else if (argv[i][1] == 'b') {
				if (argv[i+1][0] == 'J') bgmod=0;
				else if (argv[i+1][0] == 'H') bgmod=1;
				else { cout << "didn't recognize "<< argv[i+1]<<endl; help(); return 0; }			
				i++;
			}
			else if (argv[i][1] == 's') { LLRS=1; }
			else if (argv[i][1] == 'c') { thr = atof(argv[i+1]); i++; }
			else if (argv[i][1] == 'k') { KAPPA = atof(argv[i+1]); i++; }
			else if (argv[i][1] == 't') { treefile = argv[i+1];	i++; }
			else if (argv[i][1] == 'f') { freqfile = argv[i+1];	i++; }
			else if (argv[i][1] == 'r') { ref_num = atoi(argv[i+1]); i++; }
			else if (argv[i][1] == 'h') { help(); return 0; }
			else if (argv[i][1] == 'l') { LIST=1; }
			else if (argv[i][1] == 'u') { HEURISTIC=0; }
			else if (argv[i][1] == 'e') { EVALS=1; }
			else if (argv[i][1] == 'o') { outfile=argv[i+1]; i++; }
			else if (argv[i][1] == 'D') { cout << "hey dan"<< endl; DAN=1; i++; }
			else { cout << argv[i] << " is not a recognized option"<<endl; help(); return 0; }
			
		}
		else if (mfile[0] == '?') { 
			mfile = argv[i]; 
			if (!M[1].read_from_file(mfile)) { 
				cerr <<"matrix file not read\n"; 
				return 0; 
			}
		}
		else if (fastafile[0] == '?') { fastafile = argv[i]; }
		else { cout << " too many arguments: "<<argv[i]<<endl; help(); return 0; }
				
		i++;
	}
	
	if (LIST) {
		Nalign = read_alignment_list(ALI, fastafile);
		if (Nalign==0) return 0; 
		else if (PRINT_TO_SCREEN) cerr << Nalign << " alignments read"<<endl;
	}
	else {
		Nalign=1;	
		if (!ALI[0].read_fasta(fastafile)) { cerr << "error reading alignment "<< fastafile<<endl; return 0; }
	}	
		
	int a;		
	//return 0;
	//cout << ALI[0].name <<endl; return 0;
	for (a=0;a<Nalign;a++) {		//cout << ALI[a].name<<endl;
		if (ALI[a].Nseqs>2) {		 			
			string temp;
			if (treefile[0]=='?') temp=ALI[a].name + ".tree";
			else temp = treefile;
			if (!ALI[a].read_tree(temp)) { cerr << "tree file "<<temp<<" not found"<<endl; return 0; }		
		}
		else if (treefile[0]!='?') {
			float k= atof(treefile.c_str());
			//if (PRINT_TO_SCREEN) cout << "using k="<< k<<endl;
			ALI[a].phyl.initialize_pairwise(k);
		}
		if (freqfile[0]!='?') ALI[a].read_bg_freq(freqfile);
		//else ALI[a].print();
	}
	
	//if (ALI[0].Nseqs>2) ALI[0].phyl.print(ALI[0].seq);
	if (ALI[0].Nseqs>2) 
		cerr << "tree interpreted as: " <<ALI[0].phyl.tree_string(ALI[0].phyl.root,ALI[0].seq) << endl;
	M[0].initialize_neutral_motif(M[1].length,0.0); //make a "motif" for the background model

	int p,m,s,b;

	struct motif * RC= new motif [2];
	for (m=0;m<2;m++) {
		RC[m]=M[m].make_rc();
		M[m].rc = &RC[m];
		M[m].pssm = new int ** [N_MODELS];
		M[m].pvalue = new double * [N_MODELS];
		M[m].minscr = new int [N_MODELS];
		//M[m].score_to_LR = new double * [N_MODELS];
		RC[m].pssm = new int ** [N_MODELS];
		RC[m].pvalue = new double * [N_MODELS];
		RC[m].minscr = new int [N_MODELS];
		//RC[m].score_to_LR = new double * [N_MODELS];
	}
	
	for (p=0;p<M[1].length;p++) for (b=0;b<4;b++) M[0].f[p][b]=ALI[0].bgfreq[b];
			
	//double logLR (struct motif * mod, struct alignment * ali, int pos);
	double logLR2 (struct motif * m, struct sequence * seq, struct tree * phy, double * bgf, int pos, int e_mod);
	double skip_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start);
	double keep_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start);
	struct sequence remove_gaps (struct sequence * seq, int start, int size);
	double scr;//double rcscr;
	int ugp=0;
	bool compute_pHB (struct motif * m, struct tree * t, int bg, double * bgf, double kap);
		
	for (a=0;a<Nalign;a++) { 
		if (bgmod==0) ALI[a].phyl.compute_PJC();
		else if (bgmod==1) { 
			ALI[a].kappa=KAPPA;
			ALI[a].phyl.compute_PHKY(ALI[a].bgfreq, ALI[a].kappa);
		}			
	}
	if (mmod==2) {
		if (LIST) {
			if (treefile[0] != '?') {		
				if (!compute_pHB (&M[1],&ALI[0].phyl,bgmod,ALI[0].bgfreq, ALI[0].kappa)) return 0;
				compute_pHB (&RC[1], &ALI[0].phyl,bgmod,ALI[0].bgfreq, ALI[0].kappa);
			}		
			else { cerr << "run seperately"<<endl; return 0; }	
		}
		else { 
			if (!compute_pHB (&M[1],&ALI[0].phyl, bgmod,ALI[0].bgfreq, ALI[0].kappa)) return 0;
			compute_pHB (&RC[1], &ALI[0].phyl, bgmod,ALI[0].bgfreq, ALI[0].kappa);
		}		
	}	
	//return 1;
	void compute_pssm_pvals (struct motif * m, struct motif * bg, struct tree * t, int N, int e_mod, int ref);
	
	if (PVALS) {
		if (PRINT_TO_SCREEN) cerr << "working on model " << mmod << "...";
		compute_pssm_pvals(&M[1], &M[0], &ALI[0].phyl, ALI[0].Nseqs,mmod,ref_num);	
		if (PRINT_TO_SCREEN) cerr<<"done"<<endl;
		if (PRINT_TO_SCREEN) cerr << "working on model " << 4 << "...";
		compute_pssm_pvals(&M[1], &M[0], &ALI[0].phyl, ALI[0].Nseqs,4,ref_num);	
		if (PRINT_TO_SCREEN) cerr<<"done"<<endl;
	}	
	//return 1;	
	if (outfile[0] == '?') {
		if (LIST) outfile=fastafile+".occ";
		else outfile=ALI[0].name+".occ";
	}
	else if (PRINT_TO_SCREEN) cerr << "writing to "<< outfile<<endl;	
	//return 1;
	ofstream out(outfile.c_str()); //open the output file and write the headers
	
	out <<"NAME\tSTRAND\tPOS\tUG\t";
	out << ALI[0].seq[ref_num].name << "\tlog(LR)\t";
	if (EVALS) out << "E(LR>lr)\t";
	if (PVALS) out << "p(LR>lr)\t";
	for (s=0;s<ALI[0].Nseqs;s++) if (s!=ref_num) out << ALI[0].seq[s].name<<"\tlog(LR)\t"; 		
	
	if (EVALS) out << "E(s>S)\t";
	if ((LLRS)&&(PVALS)) out <<"S\tp(s>S)\tpID\tpUG\n";
	else if (LLRS) out <<"S\tpID\tpUG\n";
	else if (PVALS) out <<"p(s>S)\tpID\tpUG\n";
	else out <<"pID\tpUG\n";	
	
	//return 1;
	//int gap_count (int * seq, int start, int stop);
	double pid (struct sequence * seq, int ns, int start, int size);
	double pug (struct sequence * seq, int ns, int start, int size);
	struct sequence optimized_seq(struct sequence * ref, struct sequence * orth, int start, int length, int offset);
	double aligned_motif_pval(struct motif * m, long int * ali, int ns, int pos, int e_mod);
	double gapped_aligned_motif_pval(struct motif * m, struct motif * bg, struct alignment * ali, int ns, int pos, int e_mod, int ref); 
	long int col_to_int (struct sequence * seq, int p, int n);	
	long int gapped_complem_int (long int i, int l);
	void general_variant(long int mn, int ml, int * seq, int N);
	void generate_variant (long int mn, int ml, int * seq);
	
	int total_bp=0;
	
	for (a=0;a<Nalign;a++) {
		ugp=0;
		int * spos = new int [ALI[a].Nseqs];
		for (s=0;s<ALI[a].Nseqs;s++) spos[s]=0;
	for (p=0;p<ALI[a].length-M[1].length;p++) { 
		if (ALI[a].seq[ref_num].is_base(p)) {
			scr=skip_gap_log_LR(&M[1], &M[0], &ALI[a].seq[ref_num], p);
			
			if (scr>thr) {						
				double * sp_scr= new double [ALI[a].Nseqs];				
				
				sp_scr[ref_num]=scr;
				//sp_pos[ref_num]=p;
				int offset=gap_count(ALI[a].seq[ref_num].asints, p, p+M[1].length);
				out << ALI[a].name << " "<<M[1].name << "\t"; //"= "<<scr <<endl;
				
				double ns=1.0; double average=scr;
				//if (PRINT_TO_SCREEN) cout << p<<" "<<ALI[a].nextgap[p]<<endl;
				//if (PRINT_TO_SCREEN) for (s=0;s<ALI[a].Nseqs;s++) ALI[a].seq[s].print(p,p+M[1].length);				
				if (((ALI[a].nextgap[p]<=(p+M[1].length))||(ALI[a].nextbad[p]<=(p+M[1].length)))&&(HEURISTIC==1)) {
					struct sequence * bestalign = new struct sequence [ALI[a].Nseqs];
 					for (s=0;s<ALI[a].Nseqs;s++) {					
 						if (s == ref_num) {							
 							bestalign[s] = remove_gaps(&ALI[a].seq[s],p, M[1].length); 
 							double strand=M[1].Pr(bestalign[s].asints,0)/(M[1].Pr(bestalign[s].asints,0)+M[1].rc->Pr(bestalign[s].asints,0));						
 							out <<strand<<"\t"<< p <<"\t"<<ugp<<"\t"<<bestalign[s].bases<<"\t"<<sp_scr[ref_num]<<"\t"; 						
 						}
 						else {						
 							bestalign[s] = optimized_seq(&ALI[a].seq[ref_num],&ALI[a].seq[s],p, M[1].length, offset);
 							sp_scr[s]=skip_gap_log_LR(&M[1], &M[0], &bestalign[s],0);
 						} 						
 					} 						 					
  					bool usall=1; int pos;
  					for (s=0;s<ALI[a].Nseqs;s++) if (bestalign[s].length!=M[1].length) usall=0;
 					if (usall) { 
						// double aligned_motif_pval(struct motif * m, struct motif * bg, struct alignment * ali, int ns, int pos, int e_mod) 
 						long int * asints = new long int [M[1].length];
  						for (pos=0;pos<M[1].length;pos++) asints[pos] = col_to_int(bestalign, pos, ALI[a].Nseqs);
						if (EVALS) out << ((double)(ALI[a].length - M[1].length))*aligned_motif_pval(&M[1],asints,ALI[a].Nseqs,0,4)<<"\t";
 						if (PVALS) out << aligned_motif_pval(&M[1],asints,ALI[a].Nseqs,0,4)<<"\t";						
 						//cout << bestalign[ref_num].bases<<endl;
						for (s=0;s<ALI[a].Nseqs;s++) if (s!=ref_num) {
						  //cout <<p<<" "<< bestalign[s].bases<<endl;
						  	if (DAN) out << spos[s] <<"\t";
  							out <<bestalign[s].bases<<"\t"<<sp_scr[s]<<"\t";
  							average+=sp_scr[s];
  							ns+=1.0;
  						}
  						int emo; 						
						if (EVALS) out <<((double)(ALI[a].length - M[1].length))*aligned_motif_pval(&M[1],asints,ALI[a].Nseqs,0,mmod)<<"\t";
  						for (emo=0;emo<3;emo++) if (emo==mmod) {							
  							if (LLRS) out << logLR2(&M[1],bestalign,&ALI[a].phyl,M[0].f[0],0,emo) <<"\t";
  							if (PVALS) out <<aligned_motif_pval(&M[1],asints,ALI[a].Nseqs,0,emo)<<"\t";
  						}
						if (mmod==3) {
  							if (LLRS) out <<average/ns<<"\t";
 							if (PVALS) out <<aligned_motif_pval(&M[1],asints,ALI[a].Nseqs,0,3)<<"\t";
						}
  						out <<pid(bestalign,ALI[a].Nseqs,0,M[1].length)<<"\t";
  						out <<pug(ALI[a].seq,ALI[a].Nseqs,p,M[1].length);						
 						out <<endl;
						delete asints;
  					}
  					else {
 						if (EVALS) out<<((double)(ALI[a].length - M[1].length))*gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,4, ref_num)<<"\t";
						if (PVALS) out<<gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,4, ref_num)<<"\t";
						
 						for (s=0;s<ALI[a].Nseqs;s++) {
 							if (s!=ref_num) {
 								if (bestalign[s].length==M[1].length) out <<bestalign[s].bases<<"\t"<<sp_scr[s]<<"\t";
 								else out<< "-\t-\t";
 							}	
 						}	
 						out << "-\t-\t-\t-\t-\t-\t-\t-\t-\t-\t"<<endl;	
						total_bp--;
 					}	
					//for (s=0;s<ALI[a].Nseqs;s++) bestalign[s].erase();
 					//bestalign=0;
					//delete bestalign;				  
				}
				else {
 					int pos;
					sp_scr[ref_num] = keep_gap_log_LR(&M[1],&M[0],&ALI[a].seq[ref_num],p);
 					double strand=M[1].Pr(ALI[a].seq[ref_num].asints,p)/(M[1].Pr(ALI[a].seq[ref_num].asints,p)+M[1].rc->Pr(ALI[a].seq[ref_num].asints,p));
 					out <<strand<<"\t"<< p <<"\t"<<ugp<<"\t"; for (pos=0;pos<M[1].length;pos++) out << ALI[a].seq[ref_num].bases[p+pos];
 					out<<"\t"<<sp_scr[ref_num]<<"\t";
					if (EVALS) out<<((double)(ALI[a].length - M[1].length))*gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,4, ref_num)<<"\t";
					if (PVALS) out<<gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,4, ref_num)<<"\t";
					
 					for (s=0;s<ALI[a].Nseqs;s++) if (s!=ref_num) { 
						//sp_scr[s]=skip_gap_log_LR(&M[1],&M[0],&ALI[a].seq[s],p);
						//sp_scr[s]=log(M[1].Pr(ALI[a].seq[s].asints,p)/M[0].Pr(ALI[a].seq[s].asints,p));
						sp_scr[s]=keep_gap_log_LR(&M[1],&M[0],&ALI[a].seq[s],p);
						if (DAN) out << spos[s] <<"\t";
 						for (pos=0;pos<M[1].length;pos++) out << ALI[a].seq[s].bases[p+pos];
 						out << "\t"<<sp_scr[s] <<"\t";
 						average+=sp_scr[s];
 						ns+=1.0;
 					}	
 					int emo;
					if (EVALS) out <<((double)(ALI[a].length - M[1].length))*gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,mmod, ref_num)<<"\t";
					if (mmod==3) {
						if (LLRS) out <<average/ns<<"\t";
						if (PVALS) out <<gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,3, ref_num)<<"\t";
					}	
 					for (emo=0;emo<3;emo++) if (emo==mmod) {
 						if (LLRS) out << logLR2(&M[1],ALI[a].seq,&ALI[a].phyl,M[0].f[0],p,emo) <<"\t";
 						if (PVALS) out <<gapped_aligned_motif_pval(&M[1],&M[0],&ALI[a],ALI[a].Nseqs,p,emo, ref_num)<<"\t";
 					}
					out <<pid(ALI[a].seq,ALI[a].Nseqs,p,M[1].length)<<"\t";
 					out <<pug(ALI[a].seq,ALI[a].Nseqs,p,M[1].length);					
					out <<endl;						
				}
			}
			ugp++;
			total_bp++;
		}
		for (s=0;s<ALI[a].Nseqs;s++) if (ALI[a].seq[s].is_base(p)) spos[s]++;	
	}
	}	
	//cout <<"log(LR)\n";
	out <<total_bp<<" bp searched"<<endl;
}


void help () {
	//cout << "help contents here"<<endl;
	cout << "you are using monkey v1.0"<<endl<<"usage: [matrix file] [fasta alignment] (options)"<<endl;
	cout << "\toption\texplanation"<<endl;
	cout << "\t-m\tmotif model : mk,JC0.5,HB or simple (default is HB)"<<endl;
	cout << "\t-b\tbackground model : JC,HKY (default is JC)"<<endl;
	cout << "\t-kappa\ttransition transversion rate ratio for HKY model (default is 2.0)"<<endl;
	cout << "\t-tree\ttree file name (default is alignment.tree) or pairwise distance (default is to compute it)"<<endl;
	cout << "\t-list\tsecond argument is a text file that is a list of fasta alignments (instead of a fasta alignment)"<<endl;
	cout << "\t-help\tdisplay this mesage"<<endl;
	cout << "\t-scr\tprint scores (as well as pvalues)"<<endl;
	cout << "\t-freq\tbackground frequency file (default is to compute them from the sequences)"<<endl;
	cout << "\t-ref\treference sequence number (used in heuristics; default is 0, the first one in the file)"<<endl;
	cout << "\t-cut\theuristic score cutoff for reference genome sites to test (default is 0.0)"<<endl;
	cout << "\t-out\t.occ output filename (default is fasta alignment.occ)"<<endl;
	 
	return;
}


int read_alignment_list(struct alignment * ali, string filename) {
	ifstream data(filename.c_str());	
	if (!data) { cout << filename<<" not found" << endl; return 0; }
	int length=-1;
	while (!data.fail()) {
		string line;
		getline(data,line);
		//cout << line << endl;
		length++;
	}
	if (PRINT_TO_SCREEN) cout << length << " alignemnts found" << endl;
	//ali = new alignment [length];
	ifstream data2(filename.c_str());
	string fastafile;
	int a;
	int not_read=0;
	for (a=0;a<length;a++) {
		data2>>fastafile;
		//cout << a << " "<< fastafile<<endl;
		if (!ali[a].read_fasta(fastafile)) not_read++;
	}	
	if ((PRINT_TO_SCREEN)&&(not_read>0)) cout << not_read << " alignemnts not read" << endl;
	return length;
}


 double aligned_motif_pval(struct motif * m, long int * ali, int ns, int pos, int e_mod) {
 	long int complem_int (long int i, int l);
 	int p, scr,rcscr;
	scr=rcscr=0;
	for (p=0;p<m->length;p++) if (ali[p+pos] !=-1) {
		scr+=m->pssm[e_mod][p][ali[p+pos]];
		rcscr+=m->pssm[e_mod][m->length-p-1][complem_int(ali[p+pos],ns)];
	}
	
	return ( (scr>rcscr) ? m->pvalue[e_mod][scr] : m->pvalue[e_mod][rcscr] );	
 }
 
 double gapped_aligned_motif_pval(struct motif * m, struct motif * bg, struct alignment * ali, int ns, int pos, int e_mod, int ref) {
 	long int complem_int (long int i, int l);
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	
	//long int gapped_to_ungapped_int(long int gi, int l);
 	int p, scr,rcscr;
	scr=rcscr=0;
	//int gapscr =0; int rcgapscr=0;
	for (p=0;p<m->length;p++) {
		//long int ugint = gapped_to_ungapped_int(ali->asints[p+pos],ns);
		if (ali->asints[p+pos] !=-1) {
			scr+=m->pssm[e_mod][p][ali->asints[p+pos]];
			rcscr+=m->pssm[e_mod][m->length-p-1][complem_int(ali->asints[p+pos],ns)];
			//scr+=m->pssm[e_mod][p][ali[p+pos]];
			//cout << ali[p+pos]<<" "<<gapped_complem_int(ali[p+pos],ns)<<endl;
			//rcscr+=m->pssm[e_mod][m->length-p-1][complem_int(ali[p+pos],ns)];
		}
		else {
			//gapped=1;
			int b,s; double pm=0.0; double pbg=0.0; double prc=0.0; double llr=0.0; double rcllr=0.0;
			//double pbgrc=0.0;
			for (b=0;b<4;b++) {
				
				//double bcond = p_subtree_jc(ali->seq,ali->phyl.root,b,p+pos,1.0) ;
				double bcond = p_subtree2(ali->seq,ali->phyl.root,&ali->phyl,b,p+pos);
				pbg+=bg->f[p][b]*bcond;
			//	pbgrc+=bg->f[p][C[b]]*bcond;
				if (e_mod==0) { pm += m->f[p][b]*bcond; prc += m->rc->f[p][b]*bcond; }
				else if (e_mod==1) { 
					double mcond = p_subtree_jc(ali->seq,ali->phyl.root,b,p+pos,0.5);
					pm += m->f[p][b]*mcond; prc += m->rc->f[p][b]*mcond;
				}	
				else if (e_mod==2) {
					pm += m->f[p][b]*p_subtree_hb(ali->seq,ali->phyl.root,b,p+pos,p,m);
					prc += m->rc->f[p][b]*p_subtree_hb(ali->seq,ali->phyl.root,b,p+pos,p,m->rc);
				}	
			}
			if (e_mod<3) {
				scr += (int) (N_BINS*log(pm/pbg));
				rcscr += (int) (N_BINS*log(prc/pbg));
			}
			else if (e_mod==3) {				
				for (s=0;s<ns;s++) if (ali->seq[s].is_base(p+pos)) {
					int b = ali->seq[s].asints[p+pos];
					llr+=log(m->f[p][b]/bg->f[p][b]);
					rcllr+=log(m->rc->f[p][b]/bg->f[p][C[b]]);
				}
				scr += (int) (N_BINS*llr/((int) ns));
				rcscr+=(int) (N_BINS*rcllr/((int) ns));
			}
			else if (e_mod==4) {					
				if (ali->seq[ref].is_base(p+pos)) {
					int b = ali->seq[ref].asints[p+pos];
					llr=log(m->f[p][b]/bg->f[p][b]);
					rcllr=log(m->rc->f[p][b]/bg->f[p][C[b]]);
					//rcllr=log(m->f[m->length-p-1][C[b]]/bg->f[p][C[b]]);
				}
				scr += (int) (N_BINS*llr);
				rcscr+=(int) (N_BINS*rcllr);
				//int b;
				//cout << p <<" "<< (b=ali->seq[ref].asints[p+pos])<<" " << (((int) (N_BINS*llr)) - m->minscr[e_mod]) <<" "<<(scr-m->minscr[e_mod])<< endl;
			
				//cout << m->f[p][b] << " "<<m->rc->f[m->length-p-1][C[b]]<<endl;
			}		
			
			scr -= m->minscr[e_mod];			
			rcscr -= m->minscr[e_mod];
		}	
		//cout<< pos <<" " <<e_mod<< " "<<scr<< " "<<m->pvalue[e_mod][scr]<< " "<<rcscr<< " "<<m->pvalue[e_mod][rcscr]<<endl;
	}			
	//cout<< pos <<" " <<e_mod<< " "<<m->minscr[e_mod]<<" "<<scr<< " "<<m->pvalue[e_mod][scr]<< " "<<rcscr<< " "<<m->pvalue[e_mod][rcscr]<<endl;
	delete C;
	return ( (scr>rcscr) ? m->pvalue[e_mod][scr] : m->pvalue[e_mod][rcscr] );	
 }
 
 
 //void compute_pssm_pvals
 
 void compute_pssm_pvals (struct motif * m, struct motif * bg, struct tree * t, int N, int e_mod, int ref) {
 	void generate_variant (long int mn, int ml, int * seq);
	//void general_variant(long int mn, int ml, int * seq, int N);
	long int Nalpha = (long int) pow(4.0,N);
	//long int Nalpha = (long int) pow(5.0,N);
	m->pssm[e_mod] = new int * [m->length];
	
	int p,b,s;
	long int i;
	int min=10000;	
	//int worst=10000;
	double ** bgdist = new double * [m->length];
	//double 
	int best=0;
	for (p=0;p<m->length;p++) {
		//double check=0.0;
		m->pssm[e_mod][p] = new int [Nalpha];
		bgdist[p] = new double [Nalpha];
		for (i=0;i<Nalpha;i++) {
			struct sequence * col = new struct sequence [N];
			int * variant = new int [N];
			generate_variant(i,N,variant);
			//general_variant(i,N,variant,5);
			double llr=0.0;
			//bool gapped=0;
			for (s=0;s<N;s++) {
				if (e_mod==3) llr += log(m->f[p][variant[s]]/bg->f[p][variant[s]]);
				//if (variant[s]==4) gapped=1;
				col[s].asints = new int [1];
				col[s].asints[0]=variant[s];				
				//if (e_mod==4) cout << variant[s];
			}
			if (e_mod==3) m->pssm[e_mod][p][i] = (int) (N_BINS*(llr/((double) N)));
			if (e_mod==4) m->pssm[e_mod][p][i] = (variant[ref] == 5) ? 0 : ((int) (N_BINS*log(m->f[p][variant[ref]]/bg->f[p][variant[ref]])));
			//else {
			//if (e_mod==4) { cout << " |"<<p<<"| "; for (b=0;b<4;b++) cout << m->f[p][b] << " ";}
			double pm=0.0; 
			double pbg=0.0;			
			for (b=0;b<4;b++) {
				double bcond = p_subtree2(col,t->root,t,b,0) ;
				pbg+=bg->f[p][b]*bcond;
				if (e_mod==0) pm += m->f[p][b]*bcond;
				if (e_mod==1) pm += m->f[p][b]*p_subtree_jc(col,t->root,b,0,0.5);
				if (e_mod==2) pm += m->f[p][b]*p_subtree_hb(col,t->root,b,0,p,m);
				//if (e_mod==4) cout <<p_subtree_hb(col,t->root,b,0,p,m) <<" ";				
			}
			bgdist[p][i] = pbg;
			//if (e_mod==4) cout << endl	;
			//check += pbg;			
			if (e_mod<3) m->pssm[e_mod][p][i] = ((int) (N_BINS*log(pm/pbg)));
			//if (e_mod==2) cout << i << " "<<pm<<" "<<pbg<<" "<<m->pssm[e_mod][p][i]<<endl;
			//}	
			if (m->pssm[e_mod][p][i]<min) min=m->pssm[e_mod][p][i];
		}
		//cout << check << "=? 1"<<endl;
	}
	m->minscr[e_mod]=min;
	for (p=0;p<m->length;p++) { 
		int posbest=0;
		for (i=0;i<Nalpha;i++) { 
			m->pssm[e_mod][p][i]-=min;			
			if (m->pssm[e_mod][p][i]>posbest) posbest=m->pssm[e_mod][p][i];
		}
		best+=posbest;
	}	
	//if (e_mod==4) for (p=0;p<m->length;p++) { for (i=0;i<Nalpha;i++) cout << m->pssm[e_mod][p][i]<< " "; cout <<endl; }
	double ** pdf = new double * [m->length];
	for (p=0;p<m->length;p++) pdf[p]= new double [best+1];	
	for (p=0;p<m->length;p++) for (i=0;i<(best+1);i++) pdf[p][i] = -1.0;
	m->pvalue[e_mod] = new double [best+1];
	
	double score_pdf (int scr, int p, long int Na,double ** pbg, int ** pssm, double ** pdf);
	
	m->pvalue[e_mod][best] = score_pdf(best,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf);
	//m->pvalue[e_mod][best] = score_pdf(459,2,Nalpha,bgdist,m->pssm[e_mod],pdf);
	for (s=best-1;s>=0;s--)  {		
		m->pvalue[e_mod][s] = score_pdf(s,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf)
								+ m->pvalue[e_mod][s+1] ;
		//if (e_mod==4) cout << s <<" "<<score_pdf(s,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf)<< " "<<m->pvalue[e_mod][s]<<endl;	
		//cout << s << " "<< pdf[2][s]<<endl;
	}	
	delete pdf;
	delete bgdist;							
 }
 
 double score_pdf (int scr, int p, long int Na, double ** pbg, int ** pssm, double ** pdf) {
 	if (scr<0) return 0.0;
 	if (p==-1) return ((scr==0) ? 1.0 : 0.0);
	//cout <<p<<" "<<scr<<" "<<pdf[p][scr]<<endl;
	if (pdf[p][scr] != -1.0) return pdf[p][scr];
	long int i;
	pdf[p][scr]=0.0;
	for (i=0;i<Na;i++) pdf[p][scr]+=score_pdf(scr-pssm[p][i], p-1, Na, pbg, pssm, pdf) * pbg[p][i];
	//cout <<"p("<< scr<<" after "<<p <<")";
	//cout << "="<< pdf[p][scr]<<endl;
	//delete i;
	return pdf[p][scr];
 }
 
 long int gapped_to_ungapped_int(long int gi, int l) {
 	void general_variant (long int mn, int ml, int * seq, int N);
	long int seq_to_num (int * seq, int start, int ml);
 	int * s = new int [l];
	general_variant(gi,l,s,5);
	int p; bool gapped=0;
	for (p=0;p<l;p++) if (s[p]==4) gapped=1;
	if (gapped) return -1;
	else return seq_to_num(s,0,l);
 }
 
 long int col_to_int (struct sequence * seq, int p, int n) {
 	long int seq_to_num (int * seq, int start, int ml);
 	int s; int * asints = new int [n];	
	for (s=0;s<n;s++) asints[s] = seq[s].asints[p];
	//for (s=0;s<n;s++) cout << asints[s] << " "<<seq[s].bases[p]<< " "<<seq[s].asints[p]<<endl;
	long int i = seq_to_num(asints,0,n);
	delete asints;
	return i;
 }
 
 long int gapped_col_to_int (struct sequence * seq, int p, int n) {
 	long int gapped_seq_to_num(int * seq, int start, int ml);
	int s; int * asints = new int [n];	
	for (s=0;s<n;s++) asints[s] = seq[s].asints[p];
	long int i = gapped_seq_to_num(asints,0,n);
	delete asints;
	return i;
 }
 
 long int complem_int (long int i, int l) {
 	if (i<=-1) return -1;
 	long int seq_to_num (int * seq, int start, int ml);
	void generate_variant (long int mn, int ml, int * seq);
 	int * seq = new int [l];
	generate_variant(i,l,seq);
	//long int * rc = new long int [l];
	int p;
	int * com = new int [4]; com[0]=3; com[1]=2; com[2]=1; com[3]=0;
	//int com[4];  com[0]=3; com[1]=2; com[2]=1; com[3]=0; 
	for (p=0;p<l;p++) seq[p] = com[seq[p]];	
	long int ci = seq_to_num(seq,0,l);
	delete seq;
	delete com;
	return ci; 
 }
 
 long int gapped_complem_int (long int i, int l) {
 	long int gapped_seq_to_num (int * seq, int start, int ml);
	void general_variant (long int mn, int ml, int * seq, int N);
 	int * seq = new int [l];
	general_variant(i,l,seq,5);
	//cout << i << " "<<l<<endl;
	//long int * rc = new long int [l];
	int p;
	int * com = new int [4]; com[0]=3; com[1]=2; com[2]=1; com[3]=0;
	//int com[4];  com[0]=3; com[1]=2; com[2]=1; com[3]=0; 
	///for (p=0;p<l;p++) cout << seq[p] << " "<<endl;
	for (p=0;p<l;p++) seq[p] = ((seq[p]==-1)||(seq[p]==4)) ? 4 : com[seq[p]];	
	//for (p=0;p<l;p++) cout << seq[p] << " "<<endl;
	//long int ci = seq_to_num(seq,0,l);
	long int ci = gapped_seq_to_num(seq,0,l);
	delete seq;
	delete com;
	return ci; 
 }
 
void general_variant (long int mn, int ml, int * seq, int N) {
	int p,b;
	for (p=0;p<ml;p++) {
		b = mn % N;
		seq[p] = b;
		mn /= N;
	}
}
 
long int gapped_seq_to_num (int * seq, int start, int ml) {
	long int mn=0;
	int p;
	for (p=(ml-1);p>=0;p--) { 
		mn *= 5;
		mn += ((seq[start+p]==-1) ? 4 : seq[start+p]);
	}
	return mn;	
}

void generate_variant (long int mn, int ml, int * seq) {
	//int * seq = new int [ml];
	int p,b;
	for (p=0;p<ml;p++) {
		b = mn % 4;
		seq[p] = b;
		mn /= 4;
	}
	//print_seq(seq,ml);
	//return seq;
}


long int seq_to_num (int * seq, int start, int ml) {
	long int mn=0;
	int p;
	for (p=(ml-1);p>=0;p--) {
		mn*=4;
		mn+=seq[start+p];
	}
	return mn;
}


bool compute_pHB (struct motif * m, struct tree * t, int bg, double * bgf, double kap) {
	//cout << m->name<<endl;
	void matrix_exponential (double ** a, double ** b, int n);
	m->psubs = new double **** [m->length];
	int p,a,b,n;
	for (p=0;p<m->length;p++) {
		m->psubs[p] = new double *** [t->Nnodes];
		for (n=0;n<t->Nnodes;n++) m->psubs[p][n] = new double ** [t->Nnodes];
		
		for (n=0;n<t->Nnodes;n++) if (!t->nodelist[n].root) {
			double ** R = new double * [4];	
			m->psubs[p][t->nodelist[n].anc->number][n] = new double * [4];
			for (a=0;a<4;a++) {
				R[a] = new double [4];
				m->psubs[p][t->nodelist[n].anc->number][n][a] = new double [4];
			}
				
			if (bg==0) {
				for (a=0;a<4;a++) {
					R[a][a]=0.0;
					for (b=0;b<4;b++) if (a != b) {					
						R[a][b] = (t->nodelist[n].abr/3.0);
						if (m->f[p][a] != m->f[p][b]) 
							R[a][b] *= (log(m->f[p][b]/m->f[p][a])/(1.0 - (m->f[p][a]/m->f[p][b])));
												
						R[a][a] -= R[a][b];	
					}	
				}
			}	
			else if (bg==1) {
				double ** Q = new double * [4];
				for (a=0;a<4;a++) {
					double tot=0.0;  Q[a] = new double [4];
					for (b=0;b<4;b++) if (a != b) {	
						Q[a][b] = bgf[b];
						if (((a==0)&&(b==2)) || ((a==2)&&(b==0)) || ((a==1)&&(b==3)) || ((a==3)&&(b==1)))
							Q[a][b] *= kap;
						tot+=Q[a][b];	
					}				
					for (b=0;b<4;b++) if (a!=b) Q[a][b] *= (t->nodelist[n].abr/tot);
				}
				for (a=0;a<4;a++) {
					double idrate=0.0;	
					for (b=0;b<4;b++) if (a!=b) {
						double temp = (m->f[p][b]*Q[b][a])/(m->f[p][a]*Q[a][b]);
						//double temp = 0.5;
						if (m->f[p][a] != m->f[p][b]) R[a][b] = Q[a][b] * (log(temp)/(1.0 - (1.0/temp)));
						else R[a][b] = Q[a][b];
						idrate -= R[a][b];
					}
					R[a][a] = idrate;						
				}	
			}						 
			matrix_exponential(R, m->psubs[p][t->nodelist[n].anc->number][n], 4);
//  			void print_matrix (double ** m, int n);
//  			for (a=0;a<4;a++) cout << m->f[p][a] << " ";
//  			cout << " "<<p<<" "<<t->nodelist[n].anc->number <<"->"<<n<<" k="<<t->nodelist[n].abr<<endl;
//  			cout << "rates:"<<endl; print_matrix(R,4);
//  			cout << "probabilities:"<<endl; print_matrix(m->psubs[p][t->nodelist[n].anc->number][n], 4);
// 			
		}		
	}
// 	for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b] = new double * [Nnodes];
// 	for (n=0;n<Nnodes;n++) for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][n] = new double [Nnodes];
// 	for (n=0;n<Nnodes;n++) if (!nodelist[n].root) for (a=0;a<4;a++) for (b=0;b<4;b++) {
// 		if (a != b) psubs[a][b][nodelist[n].anc->number][n] = 0.25 - 0.25*exp(-(4.0/3.0)*nodelist[n].abr);
// 		else psubs[a][b][nodelist[n].anc->number][n] = 0.25 + 0.75*exp(-(4.0/3.0)*nodelist[n].abr);
// 	}
	
	return 1;
}

void print_matrix (double ** m, int n) {
	int a,b;
	for (a=0;a<n;a++) { 
		for (b=0;b<n;b++) cout << m[a][b]<< " ";
		cout << endl;
	}
}		

void matrix_exponential (double ** a, double ** b, int n) { //b= exp(a) a is n x n
	double temp[n][n];
	double temp2[n][n];
	int i,j,k;
	for (i=0;i<n;i++) for (j=0;j<n;j++) {
		b[i][j] = a[i][j]; 
		if (i==j) b[i][j]+=1.0;
		temp[i][j]=0.0;
		for (k=0;k<n;k++) temp[i][j]+=a[i][k]*a[k][j];
	}
	double fact=1.0;	
	double MAX_ITER=200;
	int iter;
	for (iter=2;iter<MAX_ITER;iter++) {
		fact /= ((double) iter);
		for (i=0;i<n;i++) for (j=0;j<n;j++) {
			b[i][j]+= temp[i][j]*fact;
			temp2[i][j]=0.0;
			for (k=0;k<n;k++) temp2[i][j]+=a[i][k]*temp[k][j];
		}
		for (i=0;i<n;i++) for (j=0;j<n;j++) temp[i][j] = temp2[i][j];
	}

}

struct sequence optimized_seq(struct sequence * ref, struct sequence * orth, int start, int length, int offset) {
	int s;
	//int gap_count (int * seq, int start, int stop);
	//struct sequence remove_gaps (struct sequence * seq, int start, int size);
	//double pw_pid (struct sequence * seq1,struct sequence * seq2, int start1, int start2, int length);
	int MAX_GAPS=5;
	struct sequence bestseq;
	int minus,plus,ngaps;	
	ngaps = gap_count(orth->asints,start,start+length);
	
	struct sequence refseq=remove_gaps(ref,start,length);	
	if ((ngaps>MAX_GAPS)||(offset>MAX_GAPS)) { 
	  bestseq.length=0; 
	  return bestseq; 
	}
	
	minus=plus=offset+ngaps;
	int begin=start;
	if (minus==0) bestseq = remove_gaps(orth, start, length); 	
	else {
		//cout << ngaps << " gaps and "<<offset<< " gaps in reference" <<endl;
		int ugbp=0;
		int end = ((start+plus+length)<=orth->length) ? (start+plus) : (orth->length-length) ;	
		while ((begin>0)&&(ugbp<minus)) {
			begin--; 
			if (orth->is_base(begin)) ugbp++;
		}		
			//find the best ungapped local alignment									
		//cout << "looking from " << begin << " to " << end <<" to match "<<refseq.bases<<endl;	
		double bestscr=0.0;
		int bestq,q;
		for (q=begin;q<=end;q++) if (orth->is_base(q)) {
			struct sequence temp = remove_gaps(orth,q,length);
			if (temp.length!=length) return temp;		
			double scr=pw_pid(&refseq,&temp,0,0,length);
			if (scr>bestscr) { bestq=q; bestscr=scr; }
			else if ((scr == bestscr)&& (abs(bestq-start)>abs(q-start))) { bestq=q; bestscr=scr; }
			//cout << q << " " << scr << " " <<temp.bases<<endl;	
		}
			//store the best ungapped wmer
		bestseq = remove_gaps(orth, bestq, length); 
	}
	//delete refseq;
	refseq.erase();
	return bestseq;	
}

double pw_pid (struct sequence * seq1,struct sequence * seq2, int start1, int start2, int length) {
	int p; int id=0;
	for (p=0;p<length;p++) 
		if ((seq1->is_base(p+start1)) && (seq1->asints[p+start1]==seq2->asints[p+start2])) id++;
	//cout << id << " " <<endl;	
	return ((double) id / (double) length);
}

int gap_count (int * seq, int start, int stop) {
	int pos;
	int n=0;
	for (pos=start;pos<stop;pos++) if ((seq[pos]==4)||(seq[pos]==-1)) n++;
	return n;
}

double pid (struct sequence * seq, int ns, int start, int size) {
	int s,p; float id=0;
	for (p=start;p<(size+start);p++) {
		bool same=1;
		for (s=1;s<ns;s++) if (seq[s].asints[p] != seq[0].asints[p]) same=0;
		if (same) id+=1.0;
	}
	return id/((double) size);
}

double pug (struct sequence * seq, int ns, int start, int size) {
	int s,p; float ug=0;
	for (p=start;p<(size+start);p++) {
		bool no_gap=1;
		for (s=0;s<ns;s++) if ((seq[s].asints[p] == 4) || (seq[s].asints[p] == -1)) no_gap=0;
		if (no_gap) ug+=1.0;
	}
	return ug/((double) size);
}


struct sequence remove_gaps (struct sequence * seq, int start, int size) {
	struct sequence ug;
	ug.name = seq->name + "_sub";
	int pos=start;
	while (pos<seq->length) {
		if (seq->is_base(pos)) ug.bases += seq->bases[pos];
		if(ug.bases.size()>=size) pos=seq->length;
		//cout << ug.bases <<endl;
		pos++;
	}
	ug.integerize();
	//ug.print(0,size-1);
	return ug;
}		
	
double skip_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start) {
	int pos=0; int skipped=0; 
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	double llr=0.0; double rcllr=0.0;
	while (pos<M->length) {
		int i=start+pos+skipped;
		if (i >= seq->length) return -1000.0;		
		else if (seq->is_base(i)) {
			llr += log( M->f[pos][seq->asints[i]] / B->f[pos][seq->asints[i]] );
			rcllr += log( M->rc->f[pos][seq->asints[i]] / B->f[pos][C[seq->asints[i]]] );
			pos++;
		}			
		else skipped++;
	}
	delete C;
	if (rcllr>llr) return rcllr;
	else return llr;	
}

double keep_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start) {
	int pos;
	double llr=0.0; double rcllr=0.0;
	for (pos=0;pos<M->length;pos++) {
		int i=start+pos;
		if (i >= seq->length) return -1000.0;		
		else if (seq->is_base(i)) {
			llr += log( M->f[pos][seq->asints[i]] / B->f[pos][seq->asints[i]] );
			rcllr += log( M->rc->f[pos][seq->asints[i]] / B->f[pos][seq->asints[i]] );			
		}			
		//else skipped++;
	}
	
	if (rcllr>llr) return rcllr;
	else return llr;	
}

double div_to_subs_JC (double pid) { // Jukes-cantor parameter as a function of Pid
	if (pid>0.25) {
		double k=-0.75*log((4.0*pid-1.0)/3.0);
		return ((k>MIN_RATE) ? k : MIN_RATE);
	}
	else return -1.0;
}

double p_subtree2 (struct sequence * seq, struct node * n, struct tree * phy, int base, int pos) {
	if (n->leaf) {
		if (seq[n->number].is_base(pos)) return ( (seq[n->number].asints[pos] == base) ? 1.0 : 0.0 );
		else return 1.0;
	}	
	int lb,rb; double Prl=0.0; double Prr=0.0;
	for (lb=0;lb<4;lb++) Prl+= ( p_subtree2(seq,n->left,phy,lb,pos) * phy->psubs[base][lb][n->number][n->left->number] );	
	for (rb=0;rb<4;rb++) Prr+= ( p_subtree2(seq,n->right,phy,rb,pos) * phy->psubs[base][rb][n->number][n->right->number] );
	return Prr*Prl;
}	

double p_subtree_hb (struct sequence * seq, struct node * n, int base, int sp, int mp, struct motif * m) {
	//if (n->leaf) return ( (seq[n->number].asints[sp] == base) ? 1.0 : 0.0 );
	if (n->leaf) { 
		if (seq[n->number].is_base(sp)) return ( (seq[n->number].asints[sp] == base) ? 1.0 : 0.0 );
		else return 1.0;
	}	
	int lb,rb; double Prl=0.0; double Prr=0.0;
	//for (lb=0;lb<4;lb++) Prl+= ( p_subtree_hb(seq,n->left,lb,pos,f,bgf) * PHB(base,lb,n->lbr,f,bgf) );
	//for (rb=0;rb<4;rb++) Prr+= ( p_subtree_hb(seq,n->right,rb,pos,f,bgf) * PHB(base,rb,n->rbr,f,bgf) );
	for (lb=0;lb<4;lb++) Prl+= ( p_subtree_hb(seq,n->left,lb,sp,mp,m) * m->psubs[mp][n->number][n->left->number][base][lb] );
	for (rb=0;rb<4;rb++) Prr+= ( p_subtree_hb(seq,n->right,rb,sp,mp,m) * m->psubs[mp][n->number][n->right->number][base][rb] );
	return Prr*Prl;
}

double p_subtree_jc (struct sequence * seq, struct node * n,  int base, int pos, double k) {
	//if (n->leaf) return ( (seq[n->number].asints[pos] == base) ? 1.0 : 0.0 );
	if (n->leaf) { 
		if (seq[n->number].is_base(pos)) return ( (seq[n->number].asints[pos] == base) ? 1.0 : 0.0 );
		else return 1.0;
	}	
	int lb,rb; double Prl=0.0; double Prr=0.0;
	for (lb=0;lb<4;lb++) Prl+= ( p_subtree_jc(seq,n->left,lb,pos,k) * PJC(base,lb,k*n->lbr) );
	for (rb=0;rb<4;rb++) Prr+= ( p_subtree_jc(seq,n->right,rb,pos,k) * PJC(base,rb,k*n->rbr) );
	return Prr*Prl;
}

double logLR2 (struct motif * m, struct sequence * seq, struct tree * phy, double * bgf, int pos, int e_mod) {
	double pm,prc,pbg,lhbg,lhm,lhrc;
	int p,b;	
	//double bcond[4]; double mcond[4]; double rcond[4];
	double * bcond = new double [4];
	double * mcond = new double [4];
	double * rcond = new double [4];
	
	pm=prc=pbg=1.0;
	for (p=0;p<m->length;p++) {
		for (b=0;b<4;b++) bcond[b]= p_subtree2(seq,phy->root,phy,b,p+pos);
		if (e_mod==0) { for (b=0;b<4;b++) mcond[b]=rcond[b]=bcond[b]; }
		else if (e_mod==1) { for (b=0;b<4;b++) mcond[b]=rcond[b]=p_subtree_jc(seq,phy->root,b,p+pos,0.5); }
		else if (e_mod==2) { for (b=0;b<4;b++) mcond[b]=p_subtree_hb(seq,phy->root,b,p+pos,p,m); }
		lhbg=lhm=lhrc=0.0;
		for (b=0;b<4;b++) {
			lhbg+= bgf[b] * bcond[b];
			lhm+= m->f[p][b] * mcond[b];
			
			//if (e_mod!=0) 
			//cout << m->f[p][b]<<" x "<<  mcond[b]<<"\t";
		}
		pm*=lhm;
		pbg*=lhbg;
		//if (e_mod!=0) { 
		//cout <<"mod"<<e_mod <<" pos"<<p+pos<<" "; for (b=0;b<4;b++) cout <<bcond[b] << " "; cout <<endl;
		//}
		if (BOTH_STRANDS) { 
			if (e_mod==2) for (b=0;b<4;b++) rcond[b]=p_subtree_hb(seq,phy->root,b,p+pos,p,m->rc);
			for (b=0;b<4;b++) lhrc += m->rc->f[p][b] * rcond[b];
			prc *= lhrc;
		}	
	}
	//if (PRINT_TO_SCREEN) {
	//	cout << e_mod<<": pm="<<log(pm) << " (prc="<<log(prc)<<") pbg="<< log(pbg) <<"LR= " ;
	//	if ((BOTH_STRANDS)&&(prc>pm)) cout << log(prc)-log(pbg);
	//	else cout << log(pm)-log(pbg);
	//	cout <<endl;
	//}	
	
	delete rcond; delete mcond; delete bcond;
	if ((BOTH_STRANDS)&&(prc>pm)) return (log(prc)-log(pbg));	
	return (log(pm)-log(pbg));
}

double PJC (int a, int b, double alpha) { //jukes-cantor probability of a to b in alpha//
	//cout << a<<"->"<<b<<endl;
	if (a!=b) return (0.25-0.25*exp((-4.0/3.0)*alpha));
	else return (0.25+0.75*exp((-4.0/3.0)*alpha));

}

sequence::sequence() {
	length=0;
}

//double sequence::Pr(int start, int stop, double * f) {
//	int p; double pr=1.0;
//	for (p=start;p<stop;p++) if ((asints[p]!=4)&&(asints[p]!=-1)) pr *= f[p];
//	return pr;
//}
void sequence::erase() {
  bases.erase();
  delete asints;
}
bool sequence::print(int start, int stop) {
	if (start>=stop) return 0;
	if (stop>length) return 0;
	int p,b; char B[5]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T'; B[4]='-';
	cout << name <<"("<<start<<"-"<<stop-1<<")\t";
	for (p=start;p<stop;p++) {
		if (asints[p]==-1) cout <<"N";
		else cout <<B[asints[p]];
	}
	cout << endl;
	return 1;
}
bool sequence::integerize() {
	length=bases.size();
	int p,b;
	//char B[4];
	asints = new int [length];
	//B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	//for (p=0;p<length;p++) for (b=0;b<4;b++) if (B[b] == bases[p]) asints[p]=b;
	//for (p=0;p<length;p++) cout <<seq[p]<<endl;
	for (p=0;p<length;p++) asints[p] = (
				(bases[p]=='A') ? 0 : (
					(bases[p]=='C') ? 1 : (
						(bases[p]=='G') ? 2 : (
							(bases[p]=='T') ? 3 : (
								(bases[p]=='-') ? 4 : -1
									)
								)
							)
						)
					);

	return 1;
}

bool sequence::is_base(int pos) {
	if ((asints[pos]==4)||(asints[pos]==-1)) return 0;
	else return 1;
}	
	

node::node() {
	leaf=0;
	root=0;
	//seq_post = new double [4];
	//int b;
	//for (b=0;b<4;b++) seq_post[b]=0.0;
	base=4;
	lbr=abr=rbr=INITIAL_BRANCH;

}

bool node::print() {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	int b;
	if (root) cout << "\troot"<<endl;
	else cout <<"\tn"<<anc->number<<" "<<endl;
	cout <<"\t|"<<endl;
	if (leaf) cout << "\tleaf" <<endl;
	else cout <<"n"<< left->number << "\t-\tn" <<right->number<<endl;
	return 1;
}

tree::tree() {
	Nnodes=0;
}


bool tree::print(struct sequence * seq) {
	cout << tree_string(root, seq) << endl;
	//int n=0;
// 	for (n=0;n<Nnodes;n++) {
// 	//while (!done) {
// 		cout << "node " <<n<<endl;
// 		nodelist[n].print();
// 	}
	return 1;
}

string tree::tree_string( struct node * n, struct sequence * seq) {
	if (n->leaf) return seq[n->number].name;
	else return "("+tree_string(n->left,seq) +","+ tree_string(n->right,seq)+")";
}

bool tree::compute_PJC() {
	int a,b,n;
	for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b] = new double * [Nnodes];
	for (n=0;n<Nnodes;n++) for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][n] = new double [Nnodes];
	for (n=0;n<Nnodes;n++) if (!nodelist[n].root) for (a=0;a<4;a++) for (b=0;b<4;b++) {
		if (a != b) psubs[a][b][nodelist[n].anc->number][n] = 0.25 - 0.25*exp(-(4.0/3.0)*nodelist[n].abr);
		else psubs[a][b][nodelist[n].anc->number][n] = 0.25 + 0.75*exp(-(4.0/3.0)*nodelist[n].abr);
	}
	//for (n=0;n<Nnodes;n++) if (!nodelist[n].root) for (a=0;a<4;a++) for (b=0;b<4;b++) {
		//cout << a<<"->"<<b<< " "<<nodelist[n].anc->number<<"<"<<nodelist[n].abr<<">"<< n<<":";
	//	cout<<psubs[a][b][nodelist[n].anc->number][n]<<endl;
	//}	
	return 1;
}
bool tree::compute_PHKY(double * f, double kap) {
	int a,b,n;
	void matrix_exponential (double ** a, double ** b, int n);
	for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b] = new double * [Nnodes];
	for (n=0;n<Nnodes;n++) for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][n] = new double [Nnodes];
	for (n=0;n<Nnodes;n++) if (!nodelist[n].root) {
		double ** R = new double * [4];
		double ** P = new double * [4];		
		for (a=0;a<4;a++) {
			double tot=0.0; double idrate=0.0;
			R[a] = new double [4];
			P[a] = new double [4];
			for (b=0;b<4;b++) if (a != b) {
				R[a][b] = f[b];
				//R[a][b]=1.0;
				if (((a==0)&&(b==2)) || ((a==2)&&(b==0)) || ((a==1)&&(b==3)) || ((a==3)&&(b==1)))
					R[a][b] *= kap;
				tot+=R[a][b];
			}	
			for (b=0;b<4;b++) if (a!=b) {
				R[a][b] *= (nodelist[n].abr/tot);
				idrate -= R[a][b];
			}	
			R[a][a] = idrate;			
		}
		//for (a=0;a<4;a++) { for (b=0;b<4;b++) cout << R[a][b]<< " "; cout <<endl; }
		//cout << endl;
		matrix_exponential(R, P, 4);
		for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][nodelist[n].anc->number][n] = P[a][b];
		//for (a=0;a<4;a++) { for (b=0;b<4;b++) cout << psubs[a][b][nodelist[n].anc->number][n]<< " "; cout <<endl; }
		//cout << endl;
		delete R;
		delete P;
	}	
		
	return 1;	
}

bool tree::initialize_pairwise(double dis) {
	int n;
	Nnodes=3;
	nodelist=new struct node [3];
	nodelist[2].number=2;
	root= &nodelist[2];
	nodelist[2].root=1;
	nodelist[2].left = &nodelist[0];
	nodelist[2].right= &nodelist[1];
	nodelist[2].lbr=nodelist[2].rbr=0.5*dis;
	for (n=0;n<2;n++) {
		nodelist[n].number=n;
		nodelist[n].leaf=1;
		nodelist[n].anc= &nodelist[2];
		nodelist[n].abr=0.5*dis;
	}
	return 1;
}
bool tree::read_from_file(string treefile) {
	ifstream treedata(treefile.c_str());
	if (!treedata) return 0;
	char comm;
	int thisnum,num;
	treedata>>comm;
	if (comm == 't') {
		treedata>>Nnodes;
		//cout << Nnodes <<endl;
		nodelist=new struct node [Nnodes];
	}
	int i=0;
	while (i < Nnodes) {
		treedata>>comm; //print();
		//if (!comm) return 0;
		if ((comm == 'r')||(comm == 'n')) {
			treedata>> thisnum;
			//cout <<thisnum<<endl;
			//cout << i << endl;
			nodelist[thisnum].number=thisnum;
			if (comm=='r') {
				root = &nodelist[thisnum];
				nodelist[thisnum].root=1;
			}
			treedata>>comm;
			if (comm == 'l') {
				treedata>>num;
				nodelist[thisnum].left = &nodelist[num];
				nodelist[num].leaf=1;
				nodelist[num].number=num;
				nodelist[num].anc = &nodelist[thisnum];
				//nodelist[num].lbr=nodelist[num].rbr=-1.0;
				i++;
			}
			else if (comm == 'n') {
				treedata>>num;
				nodelist[thisnum].left = &nodelist[num];
				nodelist[num].anc= &nodelist[thisnum];
			}
			treedata>>comm;
			if (comm == 'l') {
				treedata>>num;
				nodelist[thisnum].right= &nodelist[num];
				nodelist[num].leaf=1;
				nodelist[num].number=num;
				nodelist[num].anc= &nodelist[thisnum];
				i++;
			}
			else if (comm == 'n') {
				treedata>>num;
				nodelist[thisnum].right= &nodelist[num];
				nodelist[num].anc= &nodelist[thisnum];
			}
			i++;
		}
		//else {cout << i ; print();}
	//print();
	}
	treedata>>comm;
	if (comm == 'b') {
		for (i=0;i<Nnodes;i++) if (!nodelist[i].leaf) {
			treedata>>num; if (num!=i) return 0;
			treedata>>nodelist[num].lbr;
			treedata>>nodelist[num].rbr;
			nodelist[num].left->abr=nodelist[num].lbr;
			nodelist[num].right->abr=nodelist[num].rbr;
		}
	}
	//cout << "read "<<endl;
	return 1;

}

alignment::alignment() {
	length=0;
	Nseqs=0;
	gapped=0;
	bad=0;
	seq = new struct sequence [MAX_ALIGN_SIZE];
	bgfreq = new double [4];
}
bool alignment::read_bg_freq(string freqfile) {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	ifstream data(freqfile.c_str());
	if (!data) return 0;
	//else if (PRINT_TO_SCREEN) cerr << "reading "<<freqfile<<" ";
	int b;
	for (b=0;b<4;b++) {
		char base;
		data>>base;
		if (B[b] == base) { 
		  data>>bgfreq[b];  
		  //if (PRINT_TO_SCREEN) cerr<<B[b]<< ": "<<bgfreq[b]<<" "; 
		}
		else { cerr << "coudn't parse " << freqfile <<endl; return 0; }
	}
	// if (PRINT_TO_SCREEN) cerr<<endl;
	char comm;
	data>>comm;
	if (comm == 'f')  data >> FUDGE; 
	else if (comm == 'k') data >> kappa;
	//	if (PRINT_TO_SCREEN) cerr << "ki/ks=" << 1.0/FUDGE<<endl; 
	
	return 1;

}

bool alignment::print() {
	int b; char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	cout << name << " " <<length<<" bp (" <<(length-gapped-bad) << " ungapped)" <<endl;
	for (b=0;b<4;b++) cout << B[b] <<": "<<bgfreq[b]<<" ";
	cout <<endl;
	return 1;
}

void alignment::integerize() {
	int p, b,s;
	asints = new long int [length];
	long int seq_to_num (int * seq, int start, int ml);
	//long int gapped_seq_to_num (int * seq, int start, int ml);
	for (p=0;p<length;p++) {
		bool gapped=0;
		for (s=0;s<Nseqs;s++) if (!seq[s].is_base(p)) gapped=1;
		if (gapped) { asints[p] =-1; }
		else { 
			int * col = new int [Nseqs];
			for (s=0;s<Nseqs;s++) col[s] = seq[s].asints[p];
		//asints[p] = gapped_seq_to_num(col,0,Nseqs);			
			asints[p] = seq_to_num(col,0,Nseqs);
		}
		//else asints[p]=-1;
		//cout << p << " " << asints[p] << " ";
		//for (s=0;s<Nseqs;s++) cout << " "<<seq[s].asints[p];
		//cout << endl;
	}
}

bool alignment::compute_nextgap() {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	int p,b,s;
	int lastgap=length;
	int lastbad=length;
	leftend=0; rightend=-1;
	//double total=0.0;
	nextgap= new int [length];
	nextbad= new int [length];
	double totalbases=0.0;
	gapped=0;
	bad=0;
	for (b=0;b<4;b++) bgfreq[b]=0.0;
	for (p=length-1;p>=0;p--) {
		for (s=0;s<Nseqs;s++) {
			if ((seq[s].bases[p] != 'A') && (seq[s].bases[p] != 'C')
				&& (seq[s].bases[p] != 'G') && (seq[s].bases[p] != 'T')) {
				if (seq[s].bases[p] == '-') lastgap=p; //seq[s].asints[p]=-1;
				else lastbad=p;
			}
			else for (b=0;b<4;b++) if (seq[s].bases[p] == B[b]) {
				bgfreq[b]+=1.0;
				totalbases+=1.0;
			}

		}
		nextgap[p]=lastgap;
		nextbad[p]=lastbad;
		if (p==lastgap) gapped ++;
		else leftend=p;
		if (p==lastbad) bad ++;
		if (rightend==-1) if (gapped != (length-p)) rightend=p;
	}
	for (b=0;b<4;b++) bgfreq[b] /= totalbases;
	//for (p=0;p<length;p++) {
		//for (s=0;s<Nseqs;s++) cout << seq[s].bases[p];
		// cout << p <<" "<< nextgap[p] <<" "<<endl;
	//}
	return 1;
}

float alignment::div() {
	float id=0.0;
	//float ungapped=0.0;
	int p,s1,s2;
	for (p=0;p<length;p++) {
		int thisposdiff=0;
		bool thisgapped=0;
		for (s1=0;s1<Nseqs;s1++) if (seq[s1].asints[p]==4) thisgapped=1;
		if (!thisgapped) {
			//for (s1=0;s1<Nseqs;s1++)
			for (s2=1;s2<Nseqs;s2++)
				if (seq[0].asints[p] != seq[s2].asints[p]) thisposdiff=1;
			//ungapped+=1.0;
			if (!thisposdiff) id+=1.0;
		}
	}
	id /= ((float) (length-gapped-bad));
	return (1.0-id);
}

bool alignment::read_tree(string filename) {
	int p,n,b;
	//if (PRINT_TO_SCREEN) cerr << "reading tree from " <<filename<<"...";
	//phyl= new struct tree [length];
	if (!phyl.read_from_file(filename)) return 0;
	treelength=0.0;
	for (n=0;n<phyl.Nnodes;n++) if (!phyl.nodelist[n].root)
		treelength += phyl.nodelist[n].abr;
	//cout << global_tree.root->lbr<<endl;
	//double eqf[4]; eqf[0]=eqf[1]=eqf[2]=eqf[3]=0.25;
	//double loglh=0.0;
	//for (p=0;p<length;p++) if (nextgap[p]!=p) {
	//for (p=0;p<1;p++) {
		//if (!phyl[p].read_from_file(filename)) return 0;
		//phyl[p].print();
		//ALI[0].phyl[p].print();
		//for (n=0;n<phyl[p].Nnodes;n++) {
		//	if (phyl[p].nodelist[n].leaf) {
		//		phyl[p].nodelist[n].base=seq[n].asints[p];
				//for (b=0;b<4;b++) ALI[0].phyl[p].nodelist[n].seq_post[b]=0.0;
 		//		phyl[p].nodelist[n].seq_post[phyl[p].nodelist[n].base]=1.0;
		//	}
		//	else if ((global_tree.root->lbr==INITIAL_BRANCH)&&(global_tree.root->rbr==INITIAL_BRANCH)) {
		//			for (b=0;b<4;b++) phyl[p].nodelist[n].seq_post[b]=bgf[b];
		//	}
		//	else if (EVOL_MODEL==1) compute_post(phyl[p].nodelist + n, bgf);
		//	else if (EVOL_MODEL==0) compute_post_JC(phyl[p].nodelist + n, bgf);
		//}
	//	loglh+=log(compute_post_JC(phyl[p].root,bgf));
		//compute_post(phyl[p].root->left);
		//compute_post(phyl[p].root->right);
		//loglh+=compute_post(phyl[p].root, bgf);
	//}
	//cout << "initial bg log liklihood " << loglh <<endl;
	//if (EVOL_MODEL==0) {
	//if ((global_tree.root->lbr==INITIAL_BRANCH)&&(global_tree.root->rbr==INITIAL_BRANCH))
	//	optimize_branch_lengths(EPSILON,MAX_ITER);
	if ((phyl.root->lbr!=INITIAL_BRANCH)&&(phyl.root->rbr!=INITIAL_BRANCH))
	  //	if (PRINT_TO_SCREEN) cerr << "branch lengths found" <<endl;
	  //else if (PRINT_TO_SCREEN) cerr << "branch lengths not found" <<endl;
	//}
	//loglh=0.0;
	//for (p=0;p<length;p++) if (nextgap[p]!=p)
	//	loglh+=log(compute_post_JC(phyl[p].root,bgf));
	//cout << "maximum bg liklihood "<<loglh<<endl;
	return 1;
}

bool alignment::read_fasta(string filename) {
	ifstream data(filename.c_str());
	if (!data) return 0;
	//else if (PRINT_TO_SCREEN) cerr << "reading "<<filename<<"...";
	//string seqdata;
	Nseqs=-1;
	bool foundmax=0;
	int p; bool ext=1;
	for (p=filename.size();p>=0;p--) { //get rid of the last extension
		if (ext) { if (filename[p]=='.') ext=0; }
		else name=filename[p] + name;
		//cout << filename[p] <<" "<<ext<<" "<<name<<endl;
	}
	if (ext) name=filename; // if there was no extension
	//name[filename.size()]='\0';
	//name=filename;
	while ((!data.fail())&&(!foundmax)) {
		string line;
		getline(data,line);
		int linelength=line.size();
		//cout << linelength << endl;
		if ( line[0] == '>' ) {
			Nseqs++;
			if (Nseqs<MAX_ALIGN_SIZE) { for (p=1;p<linelength;p++) seq[Nseqs].name += line[p]; }
			else foundmax=1;
			//if (Nseqs>0) if (seq[Nseqs].bases.size()>length)
			//		length=seq[Nseqs].bases.size();
				//seq[Nseqs].name = new char [linelength];
				//for (p=1;p<linelength;p++) {
				//	seq[Nseqs].name+=line[p];
					//cout << p << " " << line[p]<< endl;
				//}
				//seq[Nseqs].name[linelength-1]='\0';
				//seq[Nseqs].name=line.c_str();
				//cout << Nseqs<< " "<<seq[Nseqs].name << endl;
			//}
		}
		else if ((linelength>0)&&(!foundmax)) {
			seq[Nseqs].bases += line;
			//cout << line << endl;
			//for (p=0;p<linelength;p++) if (line[p] == '\n') cout << p <<endl;
		}
	}
	int s;
	if (Nseqs<MAX_ALIGN_SIZE) Nseqs++;
	for (s=0;s<Nseqs;s++) {
		seq[s].integerize();
		if (length < seq[s].length) length=seq[s].length;
	}
	compute_nextgap();
	integerize();
	//P_BG_COMPUTED=0;
	//pbg=new double [length];
	kappa=2.0;
	double div_to_subs_JC (double pid);
	if (Nseqs==2) { // initialise the pairwise tree

		treelength=div_to_subs_JC(1.0-div()); //JC distance
		phyl.initialize_pairwise(treelength);
		//double d=div();
		//double r=(1.0-d)/d;
		//double k=-0.75*log((1.0-(r/3.0))/(1.0-r));
		//phyl.initialize_pairwise(k);
		//int b;
		//phyl= new struct tree [length];
		//for (p=0;p<length;p++) if (nextgap[p]!=p) {
		//	phyl[p].Nnodes=3;
		//	phyl[p].nodelist= new struct node [3];
		//	phyl[p].nodelist[2].number=2;
		//	phyl[p].root= &phyl[p].nodelist[2];
		//	phyl[p].nodelist[2].root=1;
		//	phyl[p].nodelist[2].left = &phyl[p].nodelist[0];
		//	phyl[p].nodelist[2].right= &phyl[p].nodelist[1];
		//	for (b=0;b<4;b++) phyl[p].nodelist[2].seq_post[b]=bgfreq[b];
		//	for (s=0;s<Nseqs;s++) {
		//		phyl[p].nodelist[s].number=s;
		//		phyl[p].nodelist[s].leaf=1;
		//		phyl[p].nodelist[s].anc= &phyl[p].nodelist[2];
		//		phyl[p].nodelist[s].base=seq[s].asints[p];
		//		phyl[p].nodelist[s].seq_post[phyl[p].nodelist[s].base]=1.0;
		//	}
		//}
		//phyl.Nnodes=3;
		//phyl.nodelist=new struct node [3];
		//phyl.nodelist[2].number=2;
		//phyl.root= &phyl.nodelist[2];
		//phyl.nodelist[2].root=1;
		//phyl.nodelist[2].left = &phyl.nodelist[0];
		//phyl.nodelist[2].right= &phyl.nodelist[1];
		//for (s=0;s<Nseqs;s++) {
		//	phyl.nodelist[s].number=s;
		//	phyl.nodelist[s].leaf=1;
		//	phyl.nodelist[s].anc= &phyl.nodelist[2];
		//}
	}
	//if (PRINT_TO_SCREEN) { cerr << "done" <<" %id="<<1.0-div()<<" k="<<treelength<<endl; }
	if (Nseqs>1) return 1;
	else return 0;
}
motif::motif() {
	//name = new char [MAX_NAME_LENGTH];
	//f = new double * [4];
	pri=MOTIF_PRIOR;
}
double motif::Pr(int s[],int start) {
	int pos; double p=1.0;
	for (pos=0;pos<length;pos++) {
		int b=s[start+pos];		
		if ((b!=-1)&&(b!=4)) p *= f[pos][b];
	}	
	return p;
}

struct motif motif::make_rc() {
	struct motif rc;
	int b, pos;
	int revcom[4]; revcom[0]=3; revcom[1]=2; revcom[2]=1; revcom[3]=0;
	rc.length=length;
	rc.pri=pri;
	rc.K = new double [length];
	rc.f = new double * [length];
	for (pos=0;pos<length;pos++) rc.f[pos] = new double [5];
	for (pos=0;pos<length;pos++) {
		for (b=0;b<4;b++) rc.f[pos][b]=f[length-pos-1][revcom[b]];
		rc.K[pos]=K[length-pos-1];
		rc.f[pos][4]=f[length-pos-1][4];
	}
	int a;	
	return rc;
}
bool motif::initialize_neutral_motif(int w, double pseudocounts) {
	int b, pos;
	length=w;
	//cout << "initialising matrix of size "<< length <<endl;
	f = new double * [length];
	K = new double [length];
	for (pos=0;pos<length;pos++) f[pos] = new double [5];
	for (pos=0;pos<length;pos++) {
		for (b=0;b<4;b++) f[pos][b]=0.25*pseudocounts;
		K[pos]=1.0;
		f[pos][4]=1.0;
	}
	return 1;
}


bool motif::print() {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	int b,pos;
	cout << "pos\t";
	for (b=0;b<4;b++) cout << B[b] <<"\t";
	cout << "Kp"<<"\t"<<"Khb"<<endl;
	for (pos=0;pos<length;pos++) {
		cout << pos << "\t";
		for (b=0;b<4;b++) cout << 0.01*((int) (100.0*f[pos][b])) <<"\t";
		cout << K[pos]<< "\t"<<halpern_bruno(pos)<<endl;

	}
	return 1;
}
double motif::halpern_bruno(int pos) {
	double hbrate=0.0;
	int b1,b2;
	for (b1=0;b1<4;b1++) for (b2=0;b2<4;b2++) if (b2!=b1) {
		if (f[pos][b1]==f[pos][b2]) hbrate+=f[pos][b1];
		else hbrate+=(f[pos][b1]*log(f[pos][b2]/f[pos][b1])/(1-(f[pos][b1]/f[pos][b2])));
	}
	return (hbrate/3.0);
}
bool motif::read_from_file(string filename) {
	ifstream data(filename.c_str());
	if (!data) { cout << filename<<" not found" << endl; return 0; }
	length=-1;
	while (!data.fail()) {
		string line;
		getline(data,line);
		//cout << line << endl;
		length++;
	}
	if (PRINT_TO_SCREEN) cout << "motif of length " << length << " found" << endl;
	int b,pos;
	f = new double * [length];
	for (pos=0;pos<length;pos++) f[pos] = new double [5];
	//int pos=0;
	ifstream datatoo(filename.c_str());
	for (pos=0;pos<length;pos++) {		
		for (b=0;b<4;b++) datatoo>>f[pos][b];			
		f[pos][4]=1.0;
	}

	K=new double [length];
	for (pos=0;pos<length;pos++) K[pos]=1.0;
	name=filename;
	return 1;
}
monkey2.0/src/rmonkey1.0.cpp0000644000163100010400000024132710466223644014251 0ustar  alanallusing namespace std;

#include <iostream.h>
#include <fstream.h>
#include <math.h>
//#include <time.h>
//#include <stdio.h>
#include <stdlib.h>
//#include <map>
#include <string>

const int MAX_ALIGN_SIZE=15; //maximum number of species
const int MAX_ALIGNMENTS = 10000;
//const int MAX_W=20; // biggest motif size
double INITIAL_BRANCH=0.1;
bool PRINT_TO_SCREEN=1;
const bool BOTH_STRANDS=1;
int MAX_GAPS=5;
const double MOTIF_PRIOR=0.001;
double thr=0.01;

const double MIN_RATE=0.001;
const int MAX_HITS=10000;
const int EVOL_MODEL=2; // 0= JC bg rate,  1= JC 1/2 bg rate 2=HB selection
const int N_MODELS=100;
double FUDGE=1.0; //correction because the HB model is based on the synonymous rate
double N_BINS=200.0;
double MODEL_ONE_RATE=0.5;
bool HEURISTIC=1;


struct sequence { // sequence object - 2 copies of the sequence, a string, and as ints
	int length;
	string name;
	string bases;
	int * asints;

	sequence();
	//~sequence();
	bool integerize();
	void erase();
	bool print(int start, int stop);
	bool is_base(int pos);
	int next_base(int pos);
	int ungapped_pos(int pos);
	//double Pr(int start, int stop, double * f);
};

struct node {
	bool leaf; bool root;
	struct node * anc;
	struct node * left;
	struct node * right;
	double abr,lbr,rbr;
	//double * seq_post;
	int base; int number; int state;

	node();
	bool print();
	struct node * sister();
};

struct tree { //hold a list of nodes, a pointer to the root and a some substitution matrices
	struct node * nodelist;
	int Nnodes;
	struct node * root;
	double ** psubs[4][4];

	tree();
	bool read_from_file(string treefile);
	bool initialize_pairwise(double dis);
	bool print(struct sequence * seq);
	bool compute_PJC();
	bool compute_PHKY(double * f, double kap);
	string tree_string( struct node * n, struct sequence * seq);
	bool copy_topology(struct tree * t);
};

struct alignment { // alignment - holds some sequences, positions of the gaps. tree
	int length;
	int Nseqs;
	int gapped;
	int bad;
	int leftend;
	int rightend;

	string name;
	struct sequence * seq;
	int * nextgap;
	int * nextbad;
	long int * asints;
	struct tree phyl;
	double treelength;

	alignment();
	bool read_bg_freq(string freqfile);
	bool read_fasta(string filename);
	bool print();
	float div();

	bool compute_nextgap();
	void integerize();

	double * bgfreq;
	double kappa;

	bool read_tree(string filename);

};


struct motif { // holds a pwm and vector of evol. rates
	int length;
	string name;
	double ** f;
	double pri;
	double * K;
	double br;
	motif * rc;
	double ***** psubs;
	int *** pssm;
	int * minscr;
	double ** pvalue;
	//double ** score_to_LR;
	//struct tree * phyl;

	motif();
	//bool intialize_from_word(int * word, int start, int stop, double pseudocounts);
	bool initialize_neutral_motif(int w, double pseudocounts);
	bool read_from_file(string mfile);
	bool print();
	double halpern_bruno(int pos);
	struct motif make_rc();
	//bool optimize_by_EM( struct SequenceData * sd,bool update_prior, bool update_matrix);
	double Pr(int s[], int start);
};


struct hit {
	int * apos;
	int * aend;
	int * ugpos;
	int refs; int leftend; int rightend;
	struct alignment ali;
	double sspval; double evolpval;
	double strand;

	hit();
	bool print();

};

double p_subtree2 (struct sequence * seq, struct node * n, struct tree * phy, int base, int pos);
double p_subtree_hb (struct sequence * seq, struct node * n, int base, int sp, int mp, struct motif * m);
double p_subtree_jc (struct sequence * seq, struct node * n,  int base, int pos, double k);

double PJC (int a, int b, double alpha);


int gap_count (int * seq, int start, int stop);

double pw_pid (struct sequence * seq1,struct sequence * seq2, int start1, int start2, int length);

int main(int argc, char * argv[]) {
	if (argc<=1) { cout << "usage: [pwm file] [fasta alignment] (options) ... or try -help for more details" <<endl; return 0; }
	int read_alignment_list(struct alignment * ali, string file);
	void help();
	//cout << argc <<endl;
	int ref_num=0,mmod=2,bgmod=0, Nalign;

	string freqfile = "?";
	string mfile = "?";
	string treefile = "?";
	string fastafile = "?";
	string outfile = "?";

	bool PVALS=1 ; bool LLRS=0; bool LIST=0; bool EVALS=0;
	bool ANY_SP=0;  int C_MEM=0;
	double KAPPA = 2.0; bool DAN=1; bool BLOCK_MODE=0; int WIDTH = 10;

	struct alignment * ALI = new alignment [MAX_ALIGNMENTS];
	struct motif * M = new motif [2];

	int i=1;
	while (i<argc) {
		//cout << argv[i] <<endl;
		if (argv[i][0] == '-') {
			if (argv[i][1] == 'm') {
				if (argv[i+1][0] == 'm') mmod=0;
				else if (argv[i+1][0] == 'J') mmod=1;
				else if (argv[i+1][0] == 'H') mmod=2;
				else if (argv[i+1][0] == 's') mmod=3;
				else { cout << "didn't recognize "<< argv[i+1] <<endl; help(); return 0; }
				i++;
			}
			else if (argv[i][1] == 'b') {
				if (argv[i+1][0] == 'J') bgmod=0;
				else if (argv[i+1][0] == 'H') bgmod=1;
				else { cout << "didn't recognize "<< argv[i+1]<<endl; help(); return 0; }
				i++;
			}
			else if (argv[i][1] == 's') { LLRS=1; }
			else if (argv[i][1] == 'c') { thr = atof(argv[i+1]); i++; }
			else if (argv[i][1] == 'k') { KAPPA = atof(argv[i+1]); i++; }
			else if (argv[i][1] == 't') { treefile = argv[i+1];	i++; }
			else if (argv[i][1] == 'f') { freqfile = argv[i+1];	i++; }
			else if (argv[i][1] == 'r') { ref_num = atoi(argv[i+1]); i++; }
			else if (argv[i][1] == 'h') { help(); return 0; }
			else if (argv[i][1] == 'l') { LIST=1; }
			else if (argv[i][1] == 'u') { HEURISTIC=0; }
			//else if (argv[i][1] == 'e') { EVALS=1; }
			else if (argv[i][1] == 'q') { PRINT_TO_SCREEN=0; }
			else if (argv[i][1] == 'o') { outfile=argv[i+1]; i++; }
			else if (argv[i][1] == 'D') { cout << "hey dan"<< endl; DAN=1; }
			else if (argv[i][1] == 'R') { MODEL_ONE_RATE = atof(argv[i+1]); i++; }
			//else if (argv[i][1] == 'B') { BLOCK_MODE = 1; }
			else if (argv[i][1] == 'a') { ANY_SP=1; }
			else if (argv[i][1] == 'M') { C_MEM=atoi(argv[i+1]); i++; }
			//else if (argv[i][1] == 'W') { WIDTH=atoi(argv[i+1]); i++; }
			else { cout << argv[i] << " is not a recognized option"<<endl; help(); return 0; }

		}
		else if (mfile[0] == '?') { mfile = argv[i]; }
		else if (fastafile[0] == '?') { fastafile = argv[i]; }
		else { cout << " too many arguments: "<<argv[i]<<endl; help(); return 0; }

		i++;
	}

	if (LIST) {
		Nalign = read_alignment_list(ALI, fastafile);
		if (Nalign==0) return 0;
		else if (PRINT_TO_SCREEN) cerr << Nalign << " alignments read"<<endl;
	}
	else {
		Nalign=1;
		if (!ALI[0].read_fasta(fastafile)) { cerr << "error reading alignment "<< fastafile<<endl; return 0; }
	}

	//if (BLOCK_MODE) {
	//	if (mmod!=1) { cerr << " use -m JC for conserved blocks"<<endl; return 0; }
	//	thr=-1.0;
	//	WIDTH = atoi(mfile.c_str());
	//	M[1].initialize_neutral_motif(WIDTH,0.0);
	//	if (PRINT_TO_SCREEN) cerr << "testing blocks of size "<<WIDTH <<" at "<<MODEL_ONE_RATE<<endl;
	//}
	//else 
	if (!M[1].read_from_file(mfile)) {
		cerr <<"matrix file not read\n";
		return 0;
	}

	int a;
	//return 0;
	//cout << ALI[0].name <<endl; return 0;
	for (a=0;a<Nalign;a++) {		//cout << ALI[a].name<<endl;
		if (ALI[a].Nseqs>2) {
			string temp;
			if (treefile[0]=='?') temp=ALI[a].name + ".tree";
			else temp = treefile;
			if (!ALI[a].read_tree(temp)) { cerr << "tree file "<<temp<<" not found"<<endl; return 0; }
		}
		else if (treefile[0]!='?') {
			float k= atof(treefile.c_str());
			//if (PRINT_TO_SCREEN) cout << "using k="<< k<<endl;
			ALI[a].phyl.initialize_pairwise(k);
		}
		if (freqfile[0]!='?') ALI[a].read_bg_freq(freqfile);
		//else ALI[a].print();
	}

	//if (ALI[0].Nseqs>2) ALI[0].phyl.print(ALI[0].seq);
	if ((ALI[0].Nseqs>2)&&(PRINT_TO_SCREEN))
		cerr << "tree interpreted as: " <<ALI[0].phyl.tree_string(ALI[0].phyl.root,ALI[0].seq) << endl;
	M[0].initialize_neutral_motif(M[1].length,0.0); //make a "motif" for the background model

	int p,m,s,b;

	struct motif * RC= new motif [2];
	for (m=0;m<2;m++) {
		RC[m]=M[m].make_rc();
		M[m].rc = &RC[m];
		M[m].pssm = new int ** [N_MODELS];
		M[m].pvalue = new double * [N_MODELS];
		M[m].minscr = new int [N_MODELS];
		//M[m].score_to_LR = new double * [N_MODELS];
		RC[m].pssm = new int ** [N_MODELS];
		RC[m].pvalue = new double * [N_MODELS];
		RC[m].minscr = new int [N_MODELS];
		//RC[m].score_to_LR = new double * [N_MODELS];
	}
	//if (BLOCK_MODE) for (p=0;p<M[1].length;p++) for (b=0;b<4;b++) M[1].f[p][b]=ALI[0].bgfreq[b];
	for (p=0;p<M[1].length;p++) for (b=0;b<4;b++) M[0].f[p][b]=ALI[0].bgfreq[b];

	//double logLR (struct motif * mod, struct alignment * ali, int pos);
	int remove_gaps (struct sequence * seq, struct sequence * ug, int start, int size);
	double scr;//double rcscr;
	int ugp=0;
	bool compute_pHB (struct motif * m, struct tree * t, int bg, double * bgf, double kap);
	for (a=0;a<Nalign;a++) {
		if (bgmod==0) ALI[a].phyl.compute_PJC();
		else if (bgmod==1) {
			ALI[a].kappa=KAPPA;
			ALI[a].phyl.compute_PHKY(ALI[a].bgfreq, ALI[a].kappa);
		}
	}
	if (mmod==2) {
		if (LIST) {
			if (treefile[0] != '?') {
				if (!compute_pHB (&M[1],&ALI[0].phyl,bgmod,ALI[0].bgfreq, ALI[0].kappa)) return 0;
				compute_pHB (&RC[1], &ALI[0].phyl,bgmod,ALI[0].bgfreq, ALI[0].kappa);
			}
			else { cerr << "run seperately"<<endl; return 0; }
		}
		else {
			if (!compute_pHB (&M[1],&ALI[0].phyl, bgmod,ALI[0].bgfreq, ALI[0].kappa)) return 0;
			compute_pHB (&RC[1], &ALI[0].phyl, bgmod,ALI[0].bgfreq, ALI[0].kappa);
		}
	}
	//return 1;
	void compute_pssm_pvals (struct motif * m, struct motif * bg, struct tree * t, int N, int e_mod, int ref);

	if (PVALS) {
		if (PRINT_TO_SCREEN) cerr << "working on model " << mmod << "...";
		compute_pssm_pvals(&M[1], &M[0], &ALI[0].phyl, ALI[0].Nseqs,mmod,ref_num);
		if (PRINT_TO_SCREEN) cerr<<"done"<<endl;
		if (PRINT_TO_SCREEN) cerr << "working on model " << 4 << "...";
		compute_pssm_pvals(&M[1], &M[0], &ALI[0].phyl, ALI[0].Nseqs,4,ref_num);
		if (PRINT_TO_SCREEN) cerr<<"done"<<endl;
		if (C_MEM<2) {
		if (PRINT_TO_SCREEN) cerr << "working on model " << 5 << "..";
		int r; 
		for (r=0;r<ALI[0].Nseqs;r++) {
			if ((r==ref_num) || (ANY_SP)) {
				compute_pssm_pvals(&M[1], &M[0], &ALI[0].phyl,ALI[0].Nseqs,5+r,r);		
				if (PRINT_TO_SCREEN) cerr <<".";
			}	
		}	
		if (PRINT_TO_SCREEN) cerr<<"done"<<endl;
		}
		if (C_MEM==0) {
			int r; 
		if (PRINT_TO_SCREEN) cerr << "working on model "<<5+ALI[0].Nseqs << "..";
		 
		for (r=0;r<ALI[0].Nseqs;r++) {
			if ((r==ref_num) || (ANY_SP)) {
				compute_pssm_pvals(&M[1], &M[0],&ALI[0].phyl,ALI[0].Nseqs,5+ALI[0].Nseqs+r,r);		
				if (PRINT_TO_SCREEN) cerr <<".";
			}	
		}
		//compute_pssm_pvals(&M[1], &M[0], &ALI[0].phyl,ALI[0].Nseqs,6,ref_num);
		if (PRINT_TO_SCREEN) cerr<<"done"<<endl;
		}


	}
	//return 1;
	if (outfile[0] == '?') {
		if (LIST) outfile=fastafile+".occ";
		else outfile=ALI[0].name+".occ";
	}
	else if (PRINT_TO_SCREEN) cerr << "writing to "<< outfile<<endl;
	//return 1;
	ofstream out(outfile.c_str()); //open the output file and write the headers

	out <<"NAME\tSTRAND\tPOS(0-based)\tUG(1-based)\t";
	out << ALI[0].seq[ref_num].name;
	if (LLRS) out << "\tlog(LR)\t";
	else out << "\tp(LR>lr)\t";
	//if (EVALS) out << "E(LR>lr)\t";
	//if (PVALS) out << "p(LR>lr)\t";
	for (s=0;s<ALI[0].Nseqs;s++) if (s!=ref_num) {
	//	if (DAN) out << "POS\t";
		out <<"POS\t"<< ALI[0].seq[s].name<<"\t";
		if (LLRS) out << "log(LR)\t";
		else out << "p(LR>lr)\t";
	}
	if (LLRS) out << "S^";
	else out << "p(S^>s)";
	out<<"\tT\tp(T<t|HB)\tp(T>t|HKY)\t";//<<endl;
	//return 1;
	//int gap_count (int * seq, int start, int stop);
	double pid (struct sequence * seq, int ns, int start, int size);
	double pug (struct sequence * seq, int ns, int start, int size);
	int optimized_seq(struct sequence * ref, struct sequence * orth, struct sequence * news, int start, int length, int offset);
	double aligned_motif_pval(struct motif * m, long int * ali, int ns, int pos, int e_mod);
	double gapped_aligned_motif_pval(struct motif * m, struct motif * bg, struct alignment * ali, int ns, int pos, int e_mod, int ref);
	double ss_pval(struct motif * m, int * seq, int pos);
	long int col_to_int (struct sequence * seq, int p, int n);
	long int gapped_complem_int (long int i, int l);
	void general_variant(long int mn, int ml, int * seq, int N);
	void generate_variant (long int mn, int ml, int * seq);
	void generate_two_state_tree (int n, struct tree * st);
	
	double logLR2 (struct motif * m, struct sequence * seq, struct tree * phy, double * bgf, int pos, int e_mod);
	double skip_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start);
	double skip_gap_ss_pval (struct motif * M,  struct sequence * seq, int start);
	double keep_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start);
	

	out<<endl;
	//return 1;
	
	int total_bp=0;
	int h;
	for (a=0;a<Nalign;a++) {
		double * ssp = new double [ALI[a].length];
		int * bestsp = new int [ALI[a].length];
		struct hit * hitlist;
		int nabove=0;
		if (ANY_SP) {
			//
			for (p=0;p<ALI[a].length-M[1].length;p++) {
				double bp=1.0; int bs;
				for (s=0;s<ALI[a].Nseqs;s++) {
					double pval = skip_gap_ss_pval(&M[1], &ALI[a].seq[s], p);
					if (pval<bp) { bs=s; bp=pval; }
				}
				ssp[p] = bp; bestsp[p] = bs;
				if (ssp[p]<thr) nabove++;
				//if (ssp[p]<thr) { cout << p << " "<<ssp[p]<<" "<<bestsp[p]<<endl; }
			}
			hitlist = new hit [nabove];
			 for (h=0;h<nabove;h++) {
				hitlist[h].apos = new int [ALI[a].Nseqs];
				//0hitlist[h].refs=-1;
			}
		}
		else {
			for (p=0;p<ALI[a].length-M[1].length;p++) {
				//if (ALI[a].seq[ref_num].is_base(p)) {
					if ((p!=0) && (p % 10000==0)&&(PRINT_TO_SCREEN==1)) cerr << ".";
					ssp[p] = skip_gap_ss_pval(&M[1],  &ALI[a].seq[ref_num], p);
					bestsp[p] = ref_num;
					if (ssp[p]<thr) nabove++;
				//if (ssp[p]<thr)
				 //cout << p <<" "<<ssp[p]<< endl;
			//if (ALI[a].seq[ref_num].is_base(p)) ugp++;
			}
			if (ALI[a].length>10000) cerr << endl;
		//cerr << nabove << " ss hits in " << ALI[a].name <<endl;
			hitlist = new hit [nabove];
			for (h=0;h<nabove;h++) {
				hitlist[h].apos = new int [ALI[a].Nseqs];
				//hitlist[h].refs=ref_num;
			}
		}
		bool find_hits ( struct alignment * a, struct motif * m, struct hit * hlist ,
						double * sspvals, int * sssp, double cut, int * hnum, int l, int r);
		int * HN = new int [1]; HN[0]=0;
		find_hits(&ALI[a],&M[1],hitlist,ssp,bestsp,thr,HN,0,ALI[a].length-M[1].length);


		delete [] bestsp; // = new int [ALI[a].length];


		//bool find_top_hit ( double * ssplist, struct hit * tophit, int left, int right, double cut);

		//find_top_hit(ssp,&hitlist[0],0,ALI[a].length-M[1].length,thr);
		//bool make_hit_alignment ( struct alignment * a, struct motif * m, struct hit * h );
		//make_hit_alignment(&ALI[a],&M[1],&hitlist[0]);
		//cout << " top hit:" <<endl; hitlist[0].print();
		 int ugp=0; int s;
		double ss_pval(struct motif * m, int * seq, int pos);
		 double aligned_motif_pval(struct motif * m, long int * ali, int ns, int pos, int e_mod);
		 double TO_statistic (struct motif * m, struct motif * bg, struct tree * t, struct alignment * ali, int pos, int ref, bool strand);
		double Tstat_pval(struct motif * m, long int * ali, int ns, int pos, bool strand, int e_mod);
		double test_leaf_loss(int s, struct alignment * ali, struct tree * t, int pos, struct motif * m, double * bgf, bool strand);
		double test_model_change_below (struct node * n, struct alignment * ali, struct tree * t, int pos,	struct motif * m, double * bgf, bool strand, int mod_above);
		double score_given_state_tree (struct motif * M, struct motif * B, struct sequence * seq, struct tree * phy, struct tree * st1, struct tree * st2, int pos);
		double state_tree_score (struct tree * st1, struct tree * st2, double br, double dr);
		double TO_statistic_given_state_tree (struct motif * m, struct motif * bg, struct alignment * ali, struct tree * t, struct tree * st1, struct tree * st2, int pos, int ref);
		//return 1;

		for (p=0;p<ALI[a].length-M[1].length;p++) {
			if (ALI[a].seq[ref_num].is_base(p)) ugp++;
			//if (ssp[p]<thr)
			for (h=0;h<HN[0];h++) if (hitlist[h].apos[ref_num] == p) {
				bool good=1;
				out << ALI[a].name<<" "<<M[1].name<<"\t"<<hitlist[h].strand
				<<"\t"<<hitlist[h].apos[ref_num]<<"-"<<hitlist[h].aend[ref_num]
				<<"\t"<<ugp<<"\t"<<hitlist[h].ali.seq[ref_num].bases<<"\t";
				if (hitlist[h].ali.seq[ref_num].length < M[1].length) {
					out << "g\t"; good=0;
				}
				else if (LLRS) { out << skip_gap_log_LR (&M[1],&M[0], &hitlist[h].ali.seq[ref_num],0)<<"\t"; }			
				else { out <<ss_pval(&M[1],hitlist[h].ali.seq[ref_num].asints,0)<<"\t"; }
				for (s=0;s<hitlist[h].ali.Nseqs;s++) if (s!=ref_num) {
					out <<hitlist[h].apos[s]<<"-"<<hitlist[h].aend[s]<<"\t"
					<<hitlist[h].ali.seq[s].bases<<"\t";
					if (hitlist[h].ali.seq[s].length < M[1].length) {
						out << "g\t"; good=0;
					}
					else if (LLRS) { out << skip_gap_log_LR (&M[1],&M[0], &hitlist[h].ali.seq[s],0)<<"\t"; }			
					else { out<<ss_pval(&M[1],hitlist[h].ali.seq[s].asints,0)<<"\t"; }
				}
				//	<<ss_pval(&M[1],hitlist[h].ali.seq[s].asints,0)<<"\t";
				//bool ungapped=1;

//					int pos; for (pos=0;pos<M[1].length;pos++) {
//								cout << hitlist[h].ali.asints[pos] << endl;
//
//							}
//							hitlist[h].print(); return 0;

				if (good) {
					if (LLRS) out << logLR2 (&M[1], hitlist[h].ali.seq, &ALI[0].phyl, ALI[0].bgfreq,0,mmod);
					else out <<aligned_motif_pval(&M[1],hitlist[h].ali.asints,hitlist[h].ali.Nseqs,0,2);
				}	
				else { out << "g"; }
				//cout << hitlist[h].apos[ref_num]<<" "<<TO_statistic(&M[1], &M[0], &ALI[a].phyl, &hitlist[h].ali,0,ref_num);
				//cout << " ";
				//cout <<aligned_motif_pval(&M[1],hitlist[h].ali.asints,hitlist[h].ali.Nseqs,0,2);
				//cout << endl;
				//out << "\t"<< hitlist[h].refs;// " was the reference"<<endl;
				if (good) {

				
				if (hitlist[h].strand<0.5) {
					//if (hitlist[h].refs==ref_num) 
					out << "\t"<<TO_statistic(&M[1], &M[0], &ALI[a].phyl, &hitlist[h].ali,0,hitlist[h].refs,0);

					out <<"\t";
					if (C_MEM<2) {
						out<<Tstat_pval(&M[1],hitlist[h].ali.asints,hitlist[h].ali.Nseqs,0,0,5+hitlist[h].refs);
					}
					else { out <<"nc"; }
					out<< "\t";
					if (C_MEM==0) { 
						out<<Tstat_pval(&M[1],hitlist[h].ali.asints,hitlist[h].ali.Nseqs,0,0,5+ALI[0].Nseqs+hitlist[h].refs);
					}
					else { out << "nc"; }
					//else out <<"\tnr\tnr";
		//			for (s=0;s<hitlist[h].ali.Nseqs;s++) {
		//				out << "\t"<< test_leaf_loss(s,&hitlist[h].ali,&ALI[a].phyl,0,&M[1],M[0].f[p],0);
					//	out << "\t" << test_model_change_below (ALI[a].phyl.nodelist[s].anc, &hitlist[h].ali, &ALI[a].phyl,0,&M[1],M[0].f[0],0,1);
					//}
				}
				else {
					//if (hitlist[h].refs==ref_num) 
					out << "\t"<<TO_statistic(&M[1], &M[0], &ALI[a].phyl, &hitlist[h].ali,0,hitlist[h].refs,1);

					out << "\t";
					if (C_MEM<2) {				
						out<<Tstat_pval(&M[1],hitlist[h].ali.asints,hitlist[h].ali.Nseqs,0,1,5+hitlist[h].refs);
					}
					else { out << "nc"; }
					out<< "\t";
					if	(C_MEM==0) {
						 out << Tstat_pval(&M[1],hitlist[h].ali.asints,hitlist[h].ali.Nseqs,0,1,5+hitlist[h].refs+ALI[0].Nseqs);
					}
					else { 	 out<<"nc"; }
					//else out <<"\tnr\tnr";
		//			for (s=0;s<hitlist[h].ali.Nseqs;s++) {
		//				out << "\t"<< test_leaf_loss(s,&hitlist[h].ali,&ALI[a].phyl,0,&M[1],M[0].f[p],1);
					//	out << "\t" << test_model_change_below (ALI[a].phyl.nodelist[s].anc, &hitlist[h].ali, &ALI[a].phyl,0,&M[1],M[0].f[0],1,1);
					//}
				
				}
				}
				else { out <<"\tg\tg\tg"; }
				out << endl;
			}
		}
		delete [] hitlist; delete [] ssp; // = new double [ALI[a].length];
	}
}

bool find_hits ( struct alignment * a, struct motif * m, struct hit * hlist ,
				double * sspvals, int * sssp, double cut, int * hnum, int l, int r) {

	bool find_top_hit ( double * ssplist, int * sssp, struct hit * tophit, int left, int right, double cut, struct alignment * ali, int w);
	bool make_hit_alignment ( struct alignment * a, struct motif * m, struct hit * h , int rightbound, int leftbound);
	//cout << "searching " << l << " to " << r <<endl;
	if (find_top_hit(sspvals,sssp,&hlist[hnum[0]],l,r,cut,a,m->length)) {
		//cout << "hit "<<hnum[0]<<" "<<l<< "-"<<r<<": "<<hlist[hnum[0]].leftend<< " in "<<hlist[hnum[0]].refs<<endl ;
		make_hit_alignment(a,m,&hlist[hnum[0]],l,r);
		//hlist[hnum[0]].print();
		int hitleft=hlist[hnum[0]].leftend;
		int hitright=hlist[hnum[0]].rightend;
		hnum[0]++;
		
		find_hits(a,m,hlist,sspvals,sssp,cut,hnum,l,hitleft-m->length);
		//find_hits(a,m,hlist,sspvals,sssp,cut,hnum,hitright,r-m->length);
		find_hits(a,m,hlist,sspvals,sssp,cut,hnum,hitright,r);
	}
	//else { cout << "none found"<<endl; }
	return 1;
}

bool make_hit_alignment ( struct alignment * a, struct motif * m, struct hit * h , int leftbound, int rightbound) {
	int s,p;
	int remove_gaps (struct sequence * seq, struct sequence * ug, int start, int size);
	int n_skipped_gaps (struct sequence * seq, int start, int w);
	h->ali.seq = new struct sequence [ a->Nseqs ];
	h->ali.Nseqs=a->Nseqs;	int ngaps_in_ref=remove_gaps (&a->seq[h->refs],&h->ali.seq[h->refs],h->apos[h->refs],m->length);
	h->aend = new int [ a->Nseqs ];
	h->aend[h->refs] = h->apos[h->refs] + m->length + ngaps_in_ref;
	//cout << h->ali.seq[h->refs].bases << endl;
	h->rightend=h->aend[h->refs]; h->leftend=h->apos[h->refs];

	h->strand = m->Pr(h->ali.seq[h->refs].asints,0)/(m->Pr(h->ali.seq[h->refs].asints,0)+m->rc->Pr(h->ali.seq[h->refs].asints,0));
	//cout << " strand = "<< strand << " " <<endl;
	//a->seq[h->refs].print(h->apos[h->refs],h->apos[h->refs]+m->length);
	double skip_gap_ss_pval_given_strand (struct motif * M, struct sequence * seq, int start, bool strand);
	
	for (s=0;s<a->Nseqs;s++) {//if (s!=h->refs) {

		int left = h->apos[h->refs]; 
		int right = (HEURISTIC) ? h->aend[h->refs]-1 : left+1;
		if (left <leftbound) cerr << "warning: "<<left <<" < "<<leftbound<<endl;
		int ugb=0; int bestp=left; double bestpval=1.0;
		//cout << s<<" "<<left <<endl;
		if (HEURISTIC) {
		
		while ((left>leftbound)&&(ugb<(m->length-1)))	{
			if (a->seq[s].is_base(left)) ugb++;	left--;
		}
		}
		//if (s==h->refs) cout << left <<" "<<leftbound ;
		//if (s==h->refs) cout << " "<<h->apos[h->refs] << endl;

		for (p=left;p<right;p++) {
			double ssp = (h->strand >= 0.5) ?
				skip_gap_ss_pval_given_strand (m, &a->seq[s], p, 1)
								:
				skip_gap_ss_pval_given_strand (m, &a->seq[s], p, 0);
			//if (s==h->refs) cout << p << " " << ssp <<endl;
			if (ssp<=bestpval) { //use the rightmost 
				int n = n_skipped_gaps(&a->seq[s],p,m->length);
				if ((p + n)<=rightbound) { bestp=p; bestpval=ssp; }
			}
		}
		h->apos[s] = h->aend[s] = bestp;
		h->aend[s] += m->length;
		int nsk = n_skipped_gaps(&a->seq[s],bestp,m->length);
		if ((nsk<50)&&(HEURISTIC)) h->aend[s] += remove_gaps (&a->seq[s],&h->ali.seq[s],bestp,m->length);
		else remove_gaps (&a->seq[s],&h->ali.seq[s],bestp,m->length);
		
		if (h->aend[s] > h->rightend) h->rightend=h->aend[s];
		if (h->apos[s] < h->leftend) h->leftend=h->apos[s];
		//if (h->refs==s) cout << "* ";
		//cout << h->apos[s]<<"-"<<h->aend[s]<< " "<< h->ali.seq[s].bases << " " << h->ali.seq[s].name << " "<<leftbound<< " "<<rightbound<<endl;

	}
	//cout <<h->sspval<<endl;
	h->ali.length=m->length;
	h->ali.integerize();

	//h->ali.print();
	return 1;
}

bool find_top_hit ( double * ssplist, int * sssp, struct hit * tophit, int left, int right, double cut, struct alignment * ali, int w) {
	int n_skipped_gaps (struct sequence * seq, int start, int w);
	int p; double bestpval=1.0; int bestp=-1;
//cout <<	"searching " << left << " to "<<right<<endl;
	for (p=left;p<right;p++) if ((ssplist[p]<bestpval) &&
	(ali->seq[sssp[p]].is_base(p))) { //don't start with a gap
		if ((n_skipped_gaps(&ali->seq[sssp[p]],p,w)+p)<right) {
			bestp=p; bestpval=ssplist[p];
			//cout << p << " : "<< ssplist[p] << " " <<sssp[p]<<endl;
		}
		
	}
	if (bestpval>=cut) return 0;
	//if (bestp==-1) return 0;
	tophit->refs = sssp[bestp];

	tophit->apos[tophit->refs]=bestp; tophit->sspval=bestpval;
	//cout << left<<"-"<<right<<" "<<tophit->apos[tophit->refs]<<" "<<bestpval<<" "<<cut<<endl;
	tophit->leftend=bestp; tophit->rightend=bestp;
	return 1;

}

void help () {
	//cout << "help contents here"<<endl;
	cout << "you are using Rmonkey v1.0"<<endl<<"usage: [matrix file] [fasta alignment] (options)"<<endl;
	cout << "\toption\texplanation"<<endl;
	//cout << "\t-m\tmotif model : mk,JC,HB or simple (default is HB)"<<endl;
	//cout << "\t-R\trate for JC motif model (default is 0.5)"<<endl;
	cout << "\t-b\tbackground model : JC,HKY (default is JC)"<<endl;
	cout << "\t-kappa\ttransition transversion rate ratio for HKY model (default is 2.0)"<<endl;
	cout << "\t-tree\ttree file name (default is alignment.tree) or pairwise distance (default is to compute it)"<<endl;
	cout << "\t-list\tsecond argument is a text file that is a list of fasta alignments (instead of a fasta alignment)"<<endl;
	cout << "\t-help\tdisplay this mesage"<<endl;
	cout << "\t-scr\tprint log LR scores (instead of pvalues)"<<endl;
	cout << "\t-freq\tbackground frequency file (default is to compute them from the sequences)"<<endl;
	cout << "\t-ref\treference sequence number (used in heuristics; default is 0, the first one in the file)"<<endl;
	cout << "\t-any\tallow match in any sequence to be the reference"<<endl;	
	cout << "\t-cut\tpvalue cutoff (default is 0.01)"<<endl;
	cout << "\t-out\t.occ output filename (default is fasta alignment.occ)"<<endl;
	cout << "\t-M\t 0,1,2 (save memory by not calculating some p-values 0=default, calculate all)"<<endl;
	cout << "\t-quiet\tprint less messages"<<endl;
	//cout << "\t-n\tmin number of sequences that have a match

	return;
}


int read_alignment_list(struct alignment * ali, string filename) {
	ifstream data(filename.c_str());
	if (!data) { cout << filename<<" not found" << endl; return 0; }
	int length=-1;
	while (!data.fail()) {
		string line;
		getline(data,line);
		//cout << line << endl;
		length++;
	}
	if (PRINT_TO_SCREEN) cout << length << " alignemnts found" << endl;
	//ali = new alignment [length];
	ifstream data2(filename.c_str());
	string fastafile;
	int a;
	int not_read=0;
	for (a=0;a<length;a++) {
		data2>>fastafile;
		//cout << a << " "<< fastafile<<endl;
		if (!ali[a].read_fasta(fastafile)) not_read++;
	}
	if ((PRINT_TO_SCREEN)&&(not_read>0)) cout << not_read << " alignemnts not read" << endl;
	return length;
}


double aligned_motif_pval(struct motif * m, long int * ali, int ns, int pos, int e_mod) {
 	long int complem_int (long int i, int l);
 	int p, scr,rcscr;
	scr=rcscr=0;
	for (p=0;p<m->length;p++) if (ali[p+pos] !=-1) {
		scr+=m->pssm[e_mod][p][ali[p+pos]];
		rcscr+=m->pssm[e_mod][m->length-p-1][complem_int(ali[p+pos],ns)];
		//cout << p << " "<< ali[p+pos]<<" "<<e_mod << " "<<m->pssm[e_mod][p][ali[p+pos]]<< " "<<m->pssm[e_mod][m->length-p-1][complem_int(ali[p+pos],ns)]<<endl;
	}

	return ( (scr>rcscr) ? m->pvalue[e_mod][scr] : m->pvalue[e_mod][rcscr] );
 }

double Tstat_pval(struct motif * m, long int * ali, int ns, int pos, bool strand,int e_mod) {
 	long int complem_int (long int i, int l);
 	int p, scr,rcscr;
	scr=rcscr=0;
	for (p=0;p<m->length;p++) if (ali[p+pos] !=-1) {
		if (strand==1) scr+=m->pssm[e_mod][p][ali[p+pos]];
		else scr+=m->pssm[e_mod][m->length-p-1][complem_int(ali[p+pos],ns)];
		//cout << p << " "<< ali[p+pos]<<" "<<m->pssm[5][p][ali[p+pos]]<< " "<<m->pssm[5][m->length-p-1][complem_int(ali[p+pos],ns)]<<endl;
	}
	return m->pvalue[e_mod][scr];
	//cout << scr << " " <<m->pvalue[5][scr]<< " "<<rcscr<<" "<<m->pvalue[5][rcscr]<<endl;
	//return ( (scr>rcscr) ? m->pvalue[5][scr] : m->pvalue[5][rcscr] );
 }

// double ss_pval_given_strand(struct motif * m, int * seq, int pos, bool strand) {
 //	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
//	int p, scr,rcscr;
//	scr=rcscr=0;
//	for (p=0;p<m->length;p++) if (seq[p+pos] !=-1) {
//		scr += m->pssm[4][p][seq[p+pos]];
//		rcscr+=m->pssm[4][m->length-p-1][C[seq[p+pos]]];
//	}
//	delete [] C;
//	//cout << scr <<" "<< rcscr<< " "<<m->pvalue[4][scr]<<" "<<m->pvalue[4][rcscr]<<endl;
//	return ( (scr>rcscr) ? m->pvalue[4][scr] : m->pvalue[4][rcscr] );
// }

bool strand_w_min_pval (struct motif * m, int * seq, int pos) {
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	int p, scr,rcscr;
	scr=rcscr=0;
	for (p=0;p<m->length;p++) if (seq[p+pos] !=-1) {
		scr += m->pssm[4][p][seq[p+pos]];
		rcscr+=m->pssm[4][m->length-p-1][C[seq[p+pos]]];
	}
	delete [] C;
	//cout << scr <<" "<< rcscr<< " "<<m->pvalue[4][scr]<<" "<<m->pvalue[4][rcscr]<<endl;
	return ( (scr>rcscr) ? 1 : 0 );
}

 double ss_pval(struct motif * m, int * seq, int pos) {
 	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	int p, scr,rcscr;
	scr=rcscr=0;
	for (p=0;p<m->length;p++) if (seq[p+pos] !=-1) {
		scr += m->pssm[4][p][seq[p+pos]];
		rcscr+=m->pssm[4][m->length-p-1][C[seq[p+pos]]];
	}
	delete [] C;
	//cout << scr <<" "<< rcscr<< " "<<m->pvalue[4][scr]<<" "<<m->pvalue[4][rcscr]<<endl;
	return ( (scr>rcscr) ? m->pvalue[4][scr] : m->pvalue[4][rcscr] );
 }

 double gapped_aligned_motif_pval(struct motif * m, struct motif * bg, struct alignment * ali, int ns, int pos, int e_mod, int ref) {
 	long int complem_int (long int i, int l);
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;

	//long int gapped_to_ungapped_int(long int gi, int l);
 	int p, scr,rcscr;
	scr=rcscr=0;
	//int gapscr =0; int rcgapscr=0;
	for (p=0;p<m->length;p++) {
		//long int ugint = gapped_to_ungapped_int(ali->asints[p+pos],ns);
		 if (ali->asints[p+pos] !=-1) {
		 	if (e_mod!=4) {
				scr+=m->pssm[e_mod][p][ali->asints[p+pos]];
				rcscr+=m->pssm[e_mod][m->length-p-1][complem_int(ali->asints[p+pos],ns)];
			}
			else {
				scr+=m->pssm[e_mod][p][ali->seq[ref].asints[p+pos]];
				rcscr+=m->pssm[e_mod][m->length-p-1][C[ali->seq[ref].asints[p+pos]]];
			}
			//scr+=m->pssm[e_mod][p][ali[p+pos]];
			//cout << ali[p+pos]<<" "<<gapped_complem_int(ali[p+pos],ns)<<endl;
			//rcscr+=m->pssm[e_mod][m->length-p-1][complem_int(ali[p+pos],ns)];
			//if (e_mod==4) cout << p << " "<< e_mod << " "<<m->pssm[e_mod][p][ali->asints[p+pos]]<< " "<<m->pssm[e_mod][m->length-p-1][complem_int(ali->asints[p+pos],ns)]<<endl;
		}
		else {
			//gapped=1;
			int b,s; double pm=0.0; double pbg=0.0; double prc=0.0; double llr=0.0; double rcllr=0.0;
			//double pbgrc=0.0;
			for (b=0;b<4;b++) {
				//double bcond = p_subtree_jc(ali->seq,ali->phyl.root,b,p+pos,1.0) ;
				double bcond = p_subtree2(ali->seq,ali->phyl.root,&ali->phyl,b,p+pos);
				pbg+=bg->f[p][b]*bcond;
			//	pbgrc+=bg->f[p][C[b]]*bcond;
				if (e_mod==0) { pm += m->f[p][b]*bcond; prc += m->rc->f[p][b]*bcond; }
				else if (e_mod==1) {
					double mcond = p_subtree_jc(ali->seq,ali->phyl.root,b,p+pos,MODEL_ONE_RATE);
					pm += m->f[p][b]*mcond; prc += m->rc->f[p][b]*mcond;
				}
				else if (e_mod==2) {
					pm += m->f[p][b]*p_subtree_hb(ali->seq,ali->phyl.root,b,p+pos,p,m);
					prc += m->rc->f[p][b]*p_subtree_hb(ali->seq,ali->phyl.root,b,p+pos,p,m->rc);
				}
			}
			if (e_mod<3) {
				scr += (int) (N_BINS*log(pm/pbg));
				rcscr += (int) (N_BINS*log(prc/pbg));
			}
			else if (e_mod==3) {
				for (s=0;s<ns;s++) if (ali->seq[s].is_base(p+pos)) {
					int b = ali->seq[s].asints[p+pos];
					llr+=log(m->f[p][b]/bg->f[p][b]);
					rcllr+=log(m->rc->f[p][b]/bg->f[p][C[b]]);
				}
				scr += (int) (N_BINS*llr/((int) ns));
				rcscr+=(int) (N_BINS*rcllr/((int) ns));
			}
			else if (e_mod==4) {
				if (ali->seq[ref].is_base(p+pos)) {
					int b = ali->seq[ref].asints[p+pos];
					llr=log(m->f[p][b]/bg->f[p][b]);
					rcllr=log(m->rc->f[p][b]/bg->f[p][C[b]]);
					//rcllr=log(m->f[m->length-p-1][C[b]]/bg->f[p][C[b]]);
				}
				//else cout << "this position if gapped"<<endl;
				scr += (int) (N_BINS*llr);
				rcscr+=(int) (N_BINS*rcllr);
				//int b;
				//cout << p <<" "<< ali->seq[ref].asints[p+pos]<<" " <<llr<<" "<< (((int) (N_BINS*llr)) - m->minscr[e_mod]) <<" "<<scr<< endl;

				//cout << m->f[p][b] << " "<<m->rc->f[m->length-p-1][C[b]]<<endl;
			}

			scr -= m->minscr[e_mod];
			rcscr -= m->minscr[e_mod];
		}
		//cout<< pos <<" " <<e_mod<< " "<<scr<< " "<<m->pvalue[e_mod][scr]<< " "<<rcscr<< " "<<m->pvalue[e_mod][rcscr]<<endl;
	}
	//cout<< pos <<" " <<e_mod<< " "<<m->minscr[e_mod]<<" "<<scr<< " "<<m->pvalue[e_mod][scr]<< " "<<rcscr<< " "<<m->pvalue[e_mod][rcscr]<<endl;
	delete [] C;
	return ( (scr>rcscr) ? m->pvalue[e_mod][scr] : m->pvalue[e_mod][rcscr] );
 }



double TO_statistic (struct motif * m, struct motif * bg, struct tree * t, struct alignment * ali,
					int pos, int ref, bool strand) {
	//int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	double llh=0.0; double llhrc=0.0; int p; int b;
	for (p=0;p<m->length;p++) {
		 double pm=0.0; double pbg=0.0; double prc=0.0;
		for (b=0;b<4;b++) {
			if (strand) pm += m->f[p][b]*p_subtree_hb(ali->seq,t->root,b,pos+p,p,m);
			else pm += m->rc->f[p][b]*p_subtree_hb(ali->seq,t->root,b,pos+p,p,m->rc);
			pbg += bg->f[p][b]*p_subtree2(ali->seq,t->root,t,b,pos+p);
		}
		llh+=log(pm/pbg); //llhrc += log(prc/pbg);
		//cout << p<<" "<<pm << " "<<prc<<" "<<pbg<<endl;
	}
	//cout << llh<<" " <<llhrc<<endl;
	//return 1.0;

	double sscorr=0.0;
	double marginalized_px_hb (struct node * n, int p, int base, struct motif * m, struct node * x, int * seq);
	double marginalized_px_bg (struct node * n, int base, struct tree * t, struct node * x, int * seq);
	double p_above (int base, struct node * n, struct alignment * ali, struct tree * t, int pos, double * bgf);
	double  p_above_hb (struct sequence * seq, struct node * n, int base, int sp, int mp, struct motif * m);
	//struct node * sis = t->nodelist[ref].sister(); int c;
	for (p=0;p<m->length;p++) {
		int * col = new int [ali->Nseqs]; int s;
		for (s=0;s<ali->Nseqs;s++) col[s] = ali->seq[s].asints[p+pos];
		double prefgbg = 0.0;
		//int refb = ali->seq[ref].asints[p+pos];
		//cout <<p<<" " << refb << " " <<endl;

		for (b=0;b<4;b++) //{
			prefgbg += bg->f[p][b] * marginalized_px_bg (t->root, b,t, &t->nodelist[ref], col);
			//double sispr=0.0;
			//= p_subtree2(ali->seq, sis, t,
		//prefgbg += p_above(refb,&t->nodelist[ref],ali,t,pos+p,bg->f[p]);
		//}
		//cout << p << " "<<prefgbg<<endl;
		double prefgm=0.0;
		if (strand) 	{
			for (b=0;b<4;b++)
				prefgm += m->f[p][b] * marginalized_px_hb (t->root, p,b,m, &t->nodelist[ref], col);
				//prefgm  += p_above_hb (ali->seq, &t->nodelist[ref],refb, pos+p, p, m) ;
			//}
		}
		else {
			for (b=0;b<4;b++)
				prefgm += m->rc->f[p][b] * marginalized_px_hb (t->root, p,b,m->rc, &t->nodelist[ref], col);
				//prefgm += p_above_hb (ali->seq, &t->nodelist[ref],refb,pos+p,p,m->rc);
		}
		//cout << prefgm << " "<< prefgbg<<" "<<p<<endl;
		sscorr += log (prefgm/prefgbg);
	}
	//cout << llh <<" "<<llhrc <<" "<<sscorr<<endl;
	//if (llh>llhrc) return (llh - sscorr);
	//else 
	return (llh - sscorr);
	//double ssllr= log(m->Pr(ali->seq[ref].asints,0)/bg->Pr(ali->seq[ref].asints,0));
	//cout << " "<<sscorr <<" "<<ssllr<<endl;
	//return (llh - sscorr);
//
////
}


double p_above (int base, struct node * n, struct alignment * ali, struct tree * t,
					int pos, double * bgf) {
	if (n->root==1) return bgf[base];
	int b,c;
	double pr=0.0;
	struct node * sis = n->sister();
	for (b=0;b<4;b++)  {
		double ps=0.0;
		for (c=0;c<4;c++) ps+= (p_subtree2(ali->seq,sis,t,c,pos)
								* t->psubs[b][c][n->anc->number][sis->number]);
		pr += (ps * p_above(b,n->anc,ali,t,pos,bgf)
				* t->psubs[b][base][n->anc->number][n->number]);
	}
	//cout << base<<" "<<pr<<" "<<n->number<<" "<<sis->number<<endl;
	return pr;
}



//void internal_node_posteriors (struct node * n, double * post, struct alignment * ali, int pos,double * bgf) {
//	int b; double tot=0.0; double lh=0.0;
//	for (b=0;b<4;b++) post[b] = p_above(b,n,ali,pos,bgf)*p_subtree2(ali->seq, n, &ali->phyl, b, pos);
//	for (b=0;b<4;b++) { tot+=post[b]; lh += post[b]; }
//	for (b=0;b<4;b++) post[b] /= tot;
//	//cout <<lh;
//}


double  p_above_hb (struct sequence * seq, struct node * n, int base, int sp, int mp, struct motif * m) {
	if (n->root==1) return m->f[mp][base];
	int b,c; double pr=0.0;
	struct node * sis = n->sister();
	for (b=0;b<4;b++) {
		double ps = 0.0;
		for (c=0;c<4;c++) ps+= (p_subtree_hb(seq,sis,c,sp,mp,m) *
								m->psubs[mp][n->anc->number][sis->number][b][c] );
		pr+= (ps * p_above_hb(seq,n->anc,b,sp,mp,m)
				* m->psubs[mp][n->anc->number][n->number][b][base]);
	}
	return pr;
}

void internal_node_post_HB (struct node * n, double * post, struct alignment * ali, int pos, struct motif * m, int mp) {
	int b; double tot=0.0; double lh=0.0;
	for (b=0;b<4;b++) post[b] = p_above_hb(ali->seq,n,b,pos,mp,m) * p_subtree_hb(ali->seq,n,b,pos,mp,m);
	for (b=0;b<4;b++) { tot+=post[b]; lh += post[b]; }
	for (b=0;b<4;b++) post[b] /= tot;
	//cout <<lh;
}

double marginalized_px_hb (struct node * n, int p, int base, struct motif * m, struct node * x, int * seq) {
	bool is_parent (struct node * n, struct node * m);
	if (x->anc->number == n->number)
		return m->psubs[p][n->number][x->number][base][seq[x->number]];
	int b; double pr=0.0;
	if (is_parent(x,n->left)) for (b=0;b<4;b++)
		pr+= m->psubs[p][n->number][n->left->number][base][b] * marginalized_px_hb(n->left,p,b,m,x,seq);
	else if (is_parent(x,n->right)) for (b=0;b<4;b++)
		pr+= m->psubs[p][n->number][n->right->number][base][b]* marginalized_px_hb(n->right,p,b,m,x,seq);
	return pr;
}


double marginalized_px_bg (struct node * n, int base, struct tree * t, struct node * x, int * seq) {
	bool is_parent (struct node * n, struct node * m);
	if (x->anc->number == n->number)
		return t->psubs[base][seq[x->number]][n->number][x->number];
	int b; double pr=0.0;
	if (is_parent(x,n->right)) for (b=0;b<4;b++)
		pr+= t->psubs[base][b][n->number][n->right->number]*marginalized_px_bg(n->right,b,t,x,seq);
	else if (is_parent(x,n->left)) for (b=0;b<4;b++)
		pr+= t->psubs[base][b][n->number][n->left->number]*marginalized_px_bg(n->left,b,t,x,seq);
	return pr;
}

bool is_parent (struct node * n, struct node * m) { //is m a parent of n
	if (m->leaf) return 0;
	if (m->root) return 1;
	if (n->root) return 0;
	if (n->anc->number == m->number) return 1;
	if (is_parent(n,m->left)) return 1;
	if (is_parent(n,m->right)) return 1;
	return 0;
}


 void compute_pssm_pvals (struct motif * m, struct motif * bg, struct tree * t, int N, int e_mod, int ref) {
 	void generate_variant (long int mn, int ml, int * seq);
	//void general_variant(long int mn, int ml, int * seq, int N);
	double marginalized_px_hb (struct node * n, int p, int base, struct motif * m, struct node * x, int * seq);
	double marginalized_px_bg (struct node * n, int base, struct tree * t, struct node * x, int * seq);

	long int Nalpha = (long int) pow(4.0,N);
	if (e_mod==4) Nalpha=4;
	//long int Nalpha = (long int) pow(5.0,N);
	m->pssm[e_mod] = new int * [m->length];

	int p,b,s;
	long int i;
	int min=10000;
	//int worst=10000;
	double ** bgdist = new double * [m->length];
	double * bgpr = new double [Nalpha];
	//double
	int best=0;
	//cout <<endl;
	if (e_mod!=4) {
	for (p=0;p<m->length;p++) {
		//double check=0.0;
		m->pssm[e_mod][p] = new int [Nalpha];
		bgdist[p] = new double [Nalpha];
		for (i=0;i<Nalpha;i++) {
			struct sequence * col = new struct sequence [N];
			int * variant = new int [N];
			generate_variant(i,N,variant);
			//general_variant(i,N,variant,5);
			double llr=0.0;
			//bool gapped=0;

			for (s=0;s<N;s++) {
				if (e_mod==3) llr += log(m->f[p][variant[s]]/bg->f[p][variant[s]]);
				//if (variant[s]==4) gapped=1;
				col[s].asints = new int [1];
				col[s].length=1;
				col[s].asints[0]=variant[s];
				//if ((p==0)&&(e_mod==4)) cout << variant[s];
			}
			if (e_mod==3) m->pssm[e_mod][p][i] = (int) (N_BINS*(llr/((double) N)));
			//if (e_mod==4) m->pssm[e_mod][p][i] = (variant[ref] == 5) ? 0 : ((int) (N_BINS*log(m->f[p][variant[ref]]/bg->f[p][variant[ref]])));
			//else {
			//if (e_mod==2) { cout << " |"<<p<<"| "; for (b=0;b<4;b++) cout << m->f[p][b] << " ";}
			double pm=0.0;
			double pbg=0.0;
			double sscorr=0.0; // used only in e_mod==5
			//for (b=0;b<4;b++) {

				//calculate bg and bgdist
			if (p==0) {
				for (b=0;b<4;b++) pbg += bg->f[p][b]*p_subtree2(col,t->root,t,b,0);
				//if (e_mod==4) cout << "b"<<pbg<<" p"<<pow(bg->f[p][b], N)<<endl;
			}

			if (p==0) bgdist[p][i] = bgpr[i] = pbg;
			else bgdist[p][i] = pbg = bgpr[i];

				//calculate the motif probabilities
			if (e_mod==0) for (b=0;b<4;b++) pm += m->f[p][b]*p_subtree2(col,t->root,t,b,0);
			else if (e_mod==1) for (b=0;b<4;b++) pm += m->f[p][b]*p_subtree_jc(col,t->root,b,0,MODEL_ONE_RATE);
			else if (e_mod==2) for (b=0;b<4;b++) pm += m->f[p][b]*p_subtree_hb(col,t->root,b,0,p,m);
			else if (e_mod>=5) { //conditioned on reference models
				double prefgbg=0.0; double prefgm=0.0;
				for (b=0;b<4;b++) {
					pm += m->f[p][b]*p_subtree_hb(col,t->root,b,0,p,m);
					prefgbg += bg->f[p][b] * marginalized_px_bg (t->root, b,t, &t->nodelist[ref], variant);
					prefgm += m->f[p][b] * marginalized_px_hb (t->root, p,b,m, &t->nodelist[ref], variant);
				}
				sscorr+=log(prefgm/prefgbg);
				if ((e_mod>=5)&&(e_mod<(5+N))) bgdist[p][i]=pm; //use HB as bgdist
			}
			//for (b=0;b<4;b++) if (e_mod==2) cout <<t->nodelist[b].abr<< " "<<p_subtree_hb(col,t->root,b,0,p,m) <<" ";
			//cout <<endl;
			//}
			//if (e_mod==2) cout << p << " "<<pm <<" "<< pbg << endl	;
			//check += pbg;
			if (e_mod<3) m->pssm[e_mod][p][i] = ((int) (N_BINS*log(pm/pbg)));
			if ((e_mod>=5)&&(e_mod<(5+N))) m->pssm[e_mod][p][i] = ((int) (-1.0*N_BINS*(log(pm/pbg) - sscorr)));
			//reverse the sign for <=T
			if (e_mod>=(5+N)) m->pssm[e_mod][p][i] = ((int) (N_BINS*(log(pm/pbg) - sscorr)));
				
			//if ((p==5)&&(e_mod==2))
			//	cout <<" "<< i <<" "<< p<<" "<<variant[ref]<< " "<<pm<<" "<<bgdist[p][i]<<" "<<m->pssm[e_mod][p][i]<<endl;

			if (m->pssm[e_mod][p][i]<min) min=m->pssm[e_mod][p][i];
		}
		//cout << check << "=? 1"<<endl;
	}
	}
	else {
		for (p=0;p<m->length;p++) {
			m->pssm[e_mod][p] = new int [Nalpha];
			bgdist[p] = new double [Nalpha];
			for (i=0;i<Nalpha;i++) {
				bgdist[p][i] = bg->f[p][i];
				m->pssm[e_mod][p][i] = ((int) (N_BINS*log(m->f[p][i]/bg->f[p][i])));
				//cout << m->pssm[e_mod][p][i] << " ";
				if (m->pssm[e_mod][p][i]<min) min=m->pssm[e_mod][p][i];
			}
			//cout << endl;
		}

	}

	m->minscr[e_mod]=min;
	for (p=0;p<m->length;p++) {
		int posbest=0;
		for (i=0;i<Nalpha;i++) {
			m->pssm[e_mod][p][i]-=min;
			if (m->pssm[e_mod][p][i]>posbest) posbest=m->pssm[e_mod][p][i];
		}
		best+=posbest;
	}
	//if (e_mod==2) for (p=0;p<m->length;p++) { for (i=0;i<Nalpha;i++) cout << m->pssm[e_mod][p][i]<< " "; cout <<endl; }
	double ** pdf = new double * [m->length];
	for (p=0;p<m->length;p++) pdf[p]= new double [best+1];
	for (p=0;p<m->length;p++) for (i=0;i<(best+1);i++) pdf[p][i] = -1.0;
	m->pvalue[e_mod] = new double [best+1];

	double score_pdf (int scr, int p, long int Na,double ** pbg, int ** pssm, double ** pdf);

	m->pvalue[e_mod][best] = score_pdf(best,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf);
	//m->pvalue[e_mod][best] = score_pdf(459,2,Nalpha,bgdist,m->pssm[e_mod],pdf);
	for (s=best-1;s>=0;s--)  {
		m->pvalue[e_mod][s] = score_pdf(s,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf)
								+ m->pvalue[e_mod][s+1] ;
		//if ((e_mod==9)&&(score_pdf(s,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf)!=0.0)) cout << s <<" "<<score_pdf(s,m->length-1,Nalpha,bgdist,m->pssm[e_mod],pdf)<< " "<<m->pvalue[e_mod][s]<<endl;
		//cout << s << " "<< pdf[2][s]<<endl;
	}
	delete [] pdf;
	delete [] bgpr;
	delete [] bgdist;
 }

 double score_pdf (int scr, int p, long int Na, double ** pbg, int ** pssm, double ** pdf) {
 	if (scr<0) return 0.0;
 	if (p==-1) return ((scr==0) ? 1.0 : 0.0);
	//cout <<p<<" "<<scr<<" "<<pdf[p][scr]<<endl;
	if (pdf[p][scr] != -1.0) return pdf[p][scr];
	long int i;
	pdf[p][scr]=0.0;
	for (i=0;i<Na;i++) pdf[p][scr]+=score_pdf(scr-pssm[p][i], p-1, Na, pbg, pssm, pdf) * pbg[p][i];
	//cout <<"p("<< scr<<" after "<<p <<")";
	//cout << "="<< pdf[p][scr]<<endl;
	//delete i;
	return pdf[p][scr];
 }

 long int gapped_to_ungapped_int(long int gi, int l) {
 	void general_variant (long int mn, int ml, int * seq, int N);
	long int seq_to_num (int * seq, int start, int ml);
 	int * s = new int [l];
	general_variant(gi,l,s,5);
	int p; bool gapped=0;
	for (p=0;p<l;p++) if (s[p]==4) gapped=1;
	if (gapped) return -1;
	else return seq_to_num(s,0,l);
 }

 long int col_to_int (struct sequence * seq, int p, int n) {
 	long int seq_to_num (int * seq, int start, int ml);
 	int s; int * asints = new int [n];
	for (s=0;s<n;s++) asints[s] = seq[s].asints[p];
	//for (s=0;s<n;s++) cout << asints[s] << " "<<seq[s].bases[p]<< " "<<seq[s].asints[p]<<endl;
	long int i = seq_to_num(asints,0,n);
	delete [] asints;
	return i;
 }

 long int gapped_col_to_int (struct sequence * seq, int p, int n) {
 	long int gapped_seq_to_num(int * seq, int start, int ml);
	int s; int * asints = new int [n];
	for (s=0;s<n;s++) asints[s] = seq[s].asints[p];
	long int i = gapped_seq_to_num(asints,0,n);
	delete [] asints;
	return i;
 }

 long int complem_int (long int i, int l) {
 	if (i<=-1) return -1;
 	long int seq_to_num (int * seq, int start, int ml);
	void generate_variant (long int mn, int ml, int * seq);
 	int * seq = new int [l];
	generate_variant(i,l,seq);
	//long int * rc = new long int [l];
	int p;
	int * com = new int [4]; com[0]=3; com[1]=2; com[2]=1; com[3]=0;
	//int com[4];  com[0]=3; com[1]=2; com[2]=1; com[3]=0;
	for (p=0;p<l;p++) seq[p] = com[seq[p]];
	long int ci = seq_to_num(seq,0,l);
	delete [] seq;
	delete [] com;
	return ci;
 }

 long int gapped_complem_int (long int i, int l) {
 	long int gapped_seq_to_num (int * seq, int start, int ml);
	void general_variant (long int mn, int ml, int * seq, int N);
 	int * seq = new int [l];
	general_variant(i,l,seq,5);
	//cout << i << " "<<l<<endl;
	//long int * rc = new long int [l];
	int p;
	int * com = new int [4]; com[0]=3; com[1]=2; com[2]=1; com[3]=0;
	//int com[4];  com[0]=3; com[1]=2; com[2]=1; com[3]=0;
	///for (p=0;p<l;p++) cout << seq[p] << " "<<endl;
	for (p=0;p<l;p++) seq[p] = ((seq[p]==-1)||(seq[p]==4)) ? 4 : com[seq[p]];
	//for (p=0;p<l;p++) cout << seq[p] << " "<<endl;
	//long int ci = seq_to_num(seq,0,l);
	long int ci = gapped_seq_to_num(seq,0,l);
	delete [] seq;
	delete [] com;
	return ci;
 }

void general_variant (long int mn, int ml, int * seq, int N) {
	int p,b;
	for (p=0;p<ml;p++) {
		b = mn % N;
		seq[p] = b;
		mn /= N;
	}
}

long int gapped_seq_to_num (int * seq, int start, int ml) {
	long int mn=0;
	int p;
	for (p=(ml-1);p>=0;p--) {
		mn *= 5;
		mn += ((seq[start+p]==-1) ? 4 : seq[start+p]);
	}
	return mn;
}

void generate_variant (long int mn, int ml, int * seq) {
	//int * seq = new int [ml];
	int p,b;
	for (p=0;p<ml;p++) {
		b = mn % 4;
		seq[p] = b;
		mn /= 4;
	}
	//print_seq(seq,ml);
	//return seq;
}


long int seq_to_num (int * seq, int start, int ml) {
	long int mn=0;
	int p;
	for (p=(ml-1);p>=0;p--) {
		mn*=4;
		mn+=seq[start+p];
	}
	return mn;
}


bool compute_pHB (struct motif * m, struct tree * t, int bg, double * bgf, double kap) {
	//cout << m->name<<endl;
	void matrix_exponential (double ** a, double ** b, int n);
	m->psubs = new double **** [m->length];
	int p,a,b,n;
	for (p=0;p<m->length;p++) {
		m->psubs[p] = new double *** [t->Nnodes];
		for (n=0;n<t->Nnodes;n++) m->psubs[p][n] = new double ** [t->Nnodes];

		for (n=0;n<t->Nnodes;n++) if (!t->nodelist[n].root) {
			double ** R = new double * [4];
			m->psubs[p][t->nodelist[n].anc->number][n] = new double * [4];
			for (a=0;a<4;a++) {
				R[a] = new double [4];
				m->psubs[p][t->nodelist[n].anc->number][n][a] = new double [4];
			}

			if (bg==0) {
				for (a=0;a<4;a++) {
					R[a][a]=0.0;
					for (b=0;b<4;b++) if (a != b) {
						R[a][b] = (t->nodelist[n].abr/3.0);
						if (m->f[p][a] != m->f[p][b])
							R[a][b] *= (log(m->f[p][b]/m->f[p][a])/(1.0 - (m->f[p][a]/m->f[p][b])));

						R[a][a] -= R[a][b];
					}
				}
			}
			else if (bg==1) {
				double ** Q = new double * [4];
				double tot=0.0;
				for (a=0;a<4;a++) {
					  Q[a] = new double [4];
					for (b=0;b<4;b++) if (a != b) {
						Q[a][b] = bgf[b];
						if (((a==0)&&(b==2)) || ((a==2)&&(b==0)) || ((a==1)&&(b==3)) || ((a==3)&&(b==1)))
							Q[a][b] *= kap;
						tot+= (bgf[a]*Q[a][b]);
					}
				}
				for (a=0;a<4;a++) for (b=0;b<4;b++)
					if (a!=b) Q[a][b] *= (t->nodelist[n].abr/tot);
				for (a=0;a<4;a++) {
					double idrate=0.0;
					for (b=0;b<4;b++) if (a!=b) {
						double temp = (m->f[p][b]*Q[b][a])/(m->f[p][a]*Q[a][b]);
						//double temp = 0.5;
						if (fabs(1.0-(1.0/temp))>1e-100) R[a][b] = Q[a][b] * (log(temp)/(1.0 - (1.0/temp)));
						else R[a][b] = Q[a][b];
						idrate -= R[a][b];
					}
					R[a][a] = idrate;
				}
			}
			matrix_exponential(R, m->psubs[p][t->nodelist[n].anc->number][n], 4);
//  			void print_matrix (double ** m, int n);
//  			for (a=0;a<4;a++) cout << m->f[p][a] << " ";
//  			cout << " "<<p<<" "<<t->nodelist[n].anc->number <<"->"<<n<<" k="<<t->nodelist[n].abr<<endl;
//  			cout << "rates:"<<endl; print_matrix(R,4);
//  			cout << "probabilities:"<<endl; print_matrix(m->psubs[p][t->nodelist[n].anc->number][n], 4);
//
		}
	}
// 	for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b] = new double * [Nnodes];
// 	for (n=0;n<Nnodes;n++) for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][n] = new double [Nnodes];
// 	for (n=0;n<Nnodes;n++) if (!nodelist[n].root) for (a=0;a<4;a++) for (b=0;b<4;b++) {
// 		if (a != b) psubs[a][b][nodelist[n].anc->number][n] = 0.25 - 0.25*exp(-(4.0/3.0)*nodelist[n].abr);
// 		else psubs[a][b][nodelist[n].anc->number][n] = 0.25 + 0.75*exp(-(4.0/3.0)*nodelist[n].abr);
// 	}

	return 1;
}

void print_matrix (double ** m, int n) {
	int a,b;
	for (a=0;a<n;a++) {
		for (b=0;b<n;b++) cout << m[a][b]<< " ";
		cout << endl;
	}
}

void matrix_exponential (double ** a, double ** b, int n) { //b= exp(a) a is n x n
	double temp[n][n];
	double temp2[n][n];
	int i,j,k;
	for (i=0;i<n;i++) for (j=0;j<n;j++) {
		b[i][j] = a[i][j];
		if (i==j) b[i][j]+=1.0;
		temp[i][j]=0.0;
		for (k=0;k<n;k++) temp[i][j]+=a[i][k]*a[k][j];
	}
	double fact=1.0;
	double MAX_ITER=200;
	int iter;
	for (iter=2;iter<MAX_ITER;iter++) {
		fact /= ((double) iter);
		for (i=0;i<n;i++) for (j=0;j<n;j++) {
			b[i][j]+= temp[i][j]*fact;
			temp2[i][j]=0.0;
			for (k=0;k<n;k++) temp2[i][j]+=a[i][k]*temp[k][j];
		}
		for (i=0;i<n;i++) for (j=0;j<n;j++) temp[i][j] = temp2[i][j];
	}

}

int optimized_seq(struct sequence * ref, struct sequence * orth, struct sequence * news,int start, int length, int offset) {

	//					****						//
	//this is substantially different from monkey1.0//
	//mix and match with caution					//
	//					****						//

	int s;
	//int  (int * seq, int start, int stop);
	int remove_gaps (struct sequence * seq, struct sequence * ug, int start, int size);
	//struct sequence remove_gaps (struct sequence * seq, int start, int size);
	//double pw_pid (struct sequence * seq1,struct sequence * seq2, int start1, int start2, int length);

	//struct sequence bestseq;
	int minus,plus,ngaps;
	ngaps = gap_count(orth->asints,start,start+length);

	if ((ngaps>MAX_GAPS)||(offset>MAX_GAPS)) {
	  news->length=0;
	  return -1;
	}

	struct sequence refseq;
	remove_gaps(ref,&refseq,0,length);
	//cout <<start<< " "<< ref->bases << endl;
	minus=plus=offset+ngaps;
	int begin=start;
	if (minus==0) { remove_gaps(orth, news, start, length); refseq.erase(); return start; }
	else {
		//cout << ngaps << " gaps and "<<offset<< " gaps in reference" <<endl;
		int ugbp=0;
		int end = ((start+plus+length)<=orth->length) ? (start+plus) : (orth->length-length) ;
		while ((begin>0)&&(ugbp<minus)) {
			begin--;
			if (orth->is_base(begin)) ugbp++;
		}
			//find the best ungapped local alignment
		//cout << ref->name << endl;
		//cout << "looking from " << begin << " to " << end <<" to match "<<refseq.bases<< " in "<<orth->print(begin,end+length);
		double bestscr=0.0;
		int bestq=-1; int q;
		for (q=begin;q<=end;q++) if (orth->is_base(q)) {
			struct sequence temp;
			remove_gaps(orth,&temp,q,length);
			if (temp.length==length) {
				//remove_gaps(orth,news,q,length);
				//return;
				double scr=pw_pid(&refseq,&temp,0,0,length);
				if (scr>bestscr) { bestq=q; bestscr=scr; }
				else if ((scr == bestscr)&&(abs(bestq-start)>abs(q-start))) { bestq=q; bestscr=scr; }
			}
			//cout << q << " " << scr << " " <<temp.bases<<endl;
		}
		refseq.erase();
		if (bestq==-1) return -1;
		else {	//store the best ungapped wmer
			remove_gaps(orth, news,bestq, length);
			//cout <<news->bases<< " chosen"<<endl;
			return bestq;
		}

		//cout << news->bases << " is the sequence"<<endl;
	}
	//delete refseq;
	//cout << "erasing ";

	//cout << " done"<<endl;


}

double pw_pid (struct sequence * seq1,struct sequence * seq2, int start1, int start2, int length) {
	int p; int id=0;
	for (p=0;p<length;p++)
		if ((seq1->is_base(p+start1)) && (seq1->asints[p+start1]==seq2->asints[p+start2])) id++;
	//cout << id << " " <<endl;
	return ((double) id / (double) length);
}

int gap_count (int * seq, int start, int stop) {
	int pos;
	int n=0;
	for (pos=start;pos<stop;pos++) if ((seq[pos]==4)||(seq[pos]==-1)) n++;
	return n;
}

double pid (struct sequence * seq, int ns, int start, int size) {
	int s,p; float id=0;
	for (p=start;p<(size+start);p++) {
		bool same=1;
		for (s=1;s<ns;s++) if (seq[s].asints[p] != seq[0].asints[p]) same=0;
		if (same) id+=1.0;
	}
	return id/((double) size);
}

double pug (struct sequence * seq, int ns, int start, int size) {
	int s,p; float ug=0;
	for (p=start;p<(size+start);p++) {
		bool no_gap=1;
		for (s=0;s<ns;s++) if ((seq[s].asints[p] == 4) || (seq[s].asints[p] == -1)) no_gap=0;
		if (no_gap) ug+=1.0;
	}
	return ug/((double) size);
}


int remove_gaps (struct sequence * seq, struct sequence * ug, int start, int size) {
	//struct sequence ug;
	ug->name = seq->name + "_sub"; ug->bases.erase();
	int pos=start; int skipped=0;
	//cout << pos << " "<<ug->name<<endl;
	while (pos<seq->length) {
		if (seq->is_base(pos)) ug->bases += seq->bases[pos];
		else skipped++;
		if ((skipped > 50)||(ug->bases.size()>=size)) pos=seq->length;
		//cout << pos <<" "<<ug->bases <<endl;
		pos++;
	}
	//cout << " out of the loop"<<endl;

	ug->integerize();
	//ug->print(0,ug->length);
	return skipped;
}

double skip_gap_ss_pval_given_strand (struct motif * M, struct sequence * seq, int start, bool strand) {
	int pos=0; int skipped=0;
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	int scr=0; int rcscr=0;
	while (pos<M->length) {
		int i=start+pos+skipped;
		if ((i >= seq->length)||(skipped>50)) return 100.0;
		else if (seq->is_base(i)) {
			if (strand) scr += M->pssm[4][pos][seq->asints[i]];
			else scr+=M->pssm[4][M->length-pos-1][C[seq->asints[i]]];
			pos++;
		}
		else skipped++;
	}
	delete [] C;
	//cout << scr << " "<<rcscr<< " "<<M->pvalue[4][scr]<<" "<<M->pvalue[4][rcscr]<<endl;
	//if (rcscr>scr) return M->pvalue[4][rcscr];
	//else
	return M->pvalue[4][scr];
}

int n_skipped_gaps (struct sequence * seq, int start, int w) {
	int pos=0; int skipped=0;
	while (pos<w) {
		int i=start+pos+skipped;
		if ((i >= seq->length)||(skipped>50)) return skipped;
		if (seq->is_base(i)) pos++;
		else skipped++;
	}
	return skipped;
}

double skip_gap_ss_pval (struct motif * M,  struct sequence * seq, int start) {
	int pos=0; int skipped=0;
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	int scr=0; int rcscr=0;
	while (pos<M->length) {
		int i=start+pos+skipped;
		if ((i >= seq->length)||(skipped>50)) return 100.0;
		else if (seq->is_base(i)) {
			scr += M->pssm[4][pos][seq->asints[i]];
			rcscr+=M->pssm[4][M->length-pos-1][C[seq->asints[i]]];
			pos++;
		}
		else skipped++;
	}
	delete [] C;
	//cout << scr << " "<<rcscr<< " "<<M->pvalue[4][scr]<<" "<<M->pvalue[4][rcscr]<<endl;
	if (rcscr>scr) return M->pvalue[4][rcscr];
	else return M->pvalue[4][scr];
}

double skip_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start) {
	int pos=0; int skipped=0;
	int * C= new int [4]; C[0]=3; C[1]=2; C[2]=1; C[3]=0;
	double llr=0.0; double rcllr=0.0;
	while (pos<M->length) {
		int i=start+pos+skipped;
		if ((i >= seq->length)||(skipped>50)) return -1000.0;
		else if (seq->is_base(i)) {
			llr += log( M->f[pos][seq->asints[i]] / B->f[pos][seq->asints[i]] );
			rcllr += log( M->rc->f[pos][seq->asints[i]] / B->f[pos][C[seq->asints[i]]] );
			pos++;
		}
		else skipped++;
	}
	delete [] C;
	if (rcllr>llr) return rcllr;
	else return llr;
}

double keep_gap_log_LR (struct motif * M, struct motif * B, struct sequence * seq, int start) {
	int pos;
	double llr=0.0; double rcllr=0.0;
	for (pos=0;pos<M->length;pos++) {
		int i=start+pos;
		if (i >= seq->length) return -1000.0;
		else if (seq->is_base(i)) {
			llr += log( M->f[pos][seq->asints[i]] / B->f[pos][seq->asints[i]] );
			rcllr += log( M->rc->f[pos][seq->asints[i]] / B->f[pos][seq->asints[i]] );
		}
		//else skipped++;
	}

	if (rcllr>llr) return rcllr;
	else return llr;
}

double div_to_subs_JC (double pid) { // Jukes-cantor parameter as a function of Pid
	if (pid>0.25) {
		double k=-0.75*log((4.0*pid-1.0)/3.0);
		return ((k>MIN_RATE) ? k : MIN_RATE);
	}
	else return -1.0;
}

double p_subtree2 (struct sequence * seq, struct node * n, struct tree * phy, int base, int pos) {
	if (n->leaf) {
		if (seq[n->number].is_base(pos)) return ( (seq[n->number].asints[pos] == base) ? 1.0 : 0.0 );
		else return 1.0;
	}
	int lb,rb; double Prl=0.0; double Prr=0.0;
	for (lb=0;lb<4;lb++) Prl+= ( p_subtree2(seq,n->left,phy,lb,pos) * phy->psubs[base][lb][n->number][n->left->number] );
	for (rb=0;rb<4;rb++) Prr+= ( p_subtree2(seq,n->right,phy,rb,pos) * phy->psubs[base][rb][n->number][n->right->number] );
	return Prr*Prl;
}

double p_subtree_hb (struct sequence * seq, struct node * n, int base, int sp, int mp, struct motif * m) {
	//if (n->leaf) return ( (seq[n->number].asints[sp] == base) ? 1.0 : 0.0 );
	if (n->leaf) {
		if (seq[n->number].is_base(sp)) return ( (seq[n->number].asints[sp] == base) ? 1.0 : 0.0 );
		else return 1.0;
	}
	int lb,rb; double Prl=0.0; double Prr=0.0;
	//for (lb=0;lb<4;lb++) Prl+= ( p_subtree_hb(seq,n->left,lb,pos,f,bgf) * PHB(base,lb,n->lbr,f,bgf) );
	//for (rb=0;rb<4;rb++) Prr+= ( p_subtree_hb(seq,n->right,rb,pos,f,bgf) * PHB(base,rb,n->rbr,f,bgf) );
	for (lb=0;lb<4;lb++) Prl+= ( p_subtree_hb(seq,n->left,lb,sp,mp,m) * m->psubs[mp][n->number][n->left->number][base][lb] );
	for (rb=0;rb<4;rb++) Prr+= ( p_subtree_hb(seq,n->right,rb,sp,mp,m) * m->psubs[mp][n->number][n->right->number][base][rb] );
	return Prr*Prl;
}

double p_subtree_jc (struct sequence * seq, struct node * n,  int base, int pos, double k) {
	//if (n->leaf) return ( (seq[n->number].asints[pos] == base) ? 1.0 : 0.0 );
	if (n->leaf) {
		if (seq[n->number].is_base(pos)) return ( (seq[n->number].asints[pos] == base) ? 1.0 : 0.0 );
		else return 1.0;
	}
	int lb,rb; double Prl=0.0; double Prr=0.0;
	for (lb=0;lb<4;lb++) Prl+= ( p_subtree_jc(seq,n->left,lb,pos,k) * PJC(base,lb,k*n->lbr) );
	for (rb=0;rb<4;rb++) Prr+= ( p_subtree_jc(seq,n->right,rb,pos,k) * PJC(base,rb,k*n->rbr) );
	return Prr*Prl;
}

double logLR2 (struct motif * m, struct sequence * seq, struct tree * phy, double * bgf, int pos, int e_mod) {
	double pm,prc,pbg,lhbg,lhm,lhrc;
	int p,b;
	//double bcond[4]; double mcond[4]; double rcond[4];

	int n;
	pm=prc=pbg=1.0;

	if (e_mod==3) {
		int s,ns=0; for (n=0;n<phy->Nnodes;n++) if (phy->nodelist[n].leaf) ns++;
		lhbg=lhm=lhrc=1.0;
		for (p=0;p<m->length;p++) for (s=0;s<ns;s++) if (seq[s].is_base(p+pos)) {
			pm *= m->f[p][seq[s].asints[p+pos]];
			pbg *= bgf[seq[s].asints[p+pos]];
			if (BOTH_STRANDS) prc *=  m->rc->f[p][seq[s].asints[p+pos]];
		}
		if ((BOTH_STRANDS)&&(prc>pm)) return ((log(prc/pbg))/((double) ns));
		return ((log(pm/pbg))/((double) ns));
	}
	else {
	double * bcond = new double [4];
	double * mcond = new double [4];
	double * rcond = new double [4];
	for (p=0;p<m->length;p++) {
		for (b=0;b<4;b++) bcond[b]= p_subtree2(seq,phy->root,phy,b,p+pos);
		if (e_mod==0) { for (b=0;b<4;b++) mcond[b]=rcond[b]=bcond[b]; }
		else if (e_mod==1) { for (b=0;b<4;b++) mcond[b]=rcond[b]=p_subtree_jc(seq,phy->root,b,p+pos,MODEL_ONE_RATE); }
		else if (e_mod==2) { for (b=0;b<4;b++) mcond[b]=p_subtree_hb(seq,phy->root,b,p+pos,p,m); }
		lhbg=lhm=lhrc=0.0;
		for (b=0;b<4;b++) {
			lhbg+= bgf[b] * bcond[b];
			lhm+= m->f[p][b] * mcond[b];
		}

		pm*=lhm;
		pbg*=lhbg;
		//if (e_mod!=0) {
		//cout <<"mod"<<e_mod <<" pos"<<p+pos<<" "; for (b=0;b<4;b++) cout <<bcond[b] << " "; cout <<endl;
		//}
		if (BOTH_STRANDS) {
			if (e_mod==2) for (b=0;b<4;b++) rcond[b]=p_subtree_hb(seq,phy->root,b,p+pos,p,m->rc);
			for (b=0;b<4;b++) lhrc += m->rc->f[p][b] * rcond[b];
			prc *= lhrc;
		}
	}
	delete [] rcond; delete [] mcond; delete [] bcond;
	if ((BOTH_STRANDS)&&(prc>pm)) return (log(prc)-log(pbg));
	return (log(pm)-log(pbg));
	}
	//if (PRINT_TO_SCREEN) {
	//	cout << e_mod<<": pm="<<log(pm) << " (prc="<<log(prc)<<") pbg="<< log(pbg) <<"LR= " ;
	//	if ((BOTH_STRANDS)&&(prc>pm)) cout << log(prc)-log(pbg);
	//	else cout << log(pm)-log(pbg);
	//	cout <<endl;
	//}


}

double PJC (int a, int b, double alpha) { //jukes-cantor probability of a to b in alpha//
	//cout << a<<"->"<<b<<endl;
	if (a!=b) return (0.25-0.25*exp((-4.0/3.0)*alpha));
	else return (0.25+0.75*exp((-4.0/3.0)*alpha));

}

hit::hit() {
	sspval=evolpval=1.0;
	ali.Nseqs=0;
}

bool hit::print() {
	cout << apos[refs] <<"("<<leftend<<"-"<<rightend << ") "<<sspval<<endl;
	if (ali.Nseqs>0) {
		int s;
		for (s=0;s<ali.Nseqs;s++)
			cout << apos[s] << " "<<ali.seq[s].bases<< " "<<ali.seq[s].name<<endl;
	}
	return 1;
}

sequence::sequence() {
	length=0;
}

//sequence::~sequence() {
	//cout << "deleting "<<name<<" "<<length<<"bp ...";
	//name.erase();
  	//bases.erase();

  	//delete [] asints;
  	//cout<<"done"<<endl;
	//erase();
	//if (length>0) delete [] asints;
//}
//double sequence::Pr(int start, int stop, double * f) {
//	int p; double pr=1.0;
//	for (p=start;p<stop;p++) if ((asints[p]!=4)&&(asints[p]!=-1)) pr *= f[p];
//	return pr;
//}
void sequence::erase() {
	name.erase();
  bases.erase();
 // cerr << " about to erase...";
  delete [] asints;
 // cerr << " done"<<endl;
  return;
}
bool sequence::print(int start, int stop) {
	if (start>=stop) return 0;
	if (stop>length) return 0;
	int p,b; char B[5]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T'; B[4]='-';
	cout << name <<"("<<start<<"-"<<stop-1<<")\t";
	for (p=start;p<stop;p++) {
		if (asints[p]==-1) cout <<"N";
		else cout <<B[asints[p]];
	}
	cout << endl;
	return 1;
}
bool sequence::integerize() {
	length=bases.size();
	int p,b;
	//char B[4];
	//cout << length << endl;
	asints = new int [length];
	//cout << "declared " <<endl;
	//B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	//for (p=0;p<length;p++) for (b=0;b<4;b++) if (B[b] == bases[p]) asints[p]=b;
	//for (p=0;p<length;p++) cout <<bases[p]<<endl;
	//cout << bases <<endl;
	for (p=0;p<length;p++) asints[p] = (
				(bases[p]=='A') ? 0 : (
					(bases[p]=='C') ? 1 : (
						(bases[p]=='G') ? 2 : (
							(bases[p]=='T') ? 3 : (
								(bases[p]=='-') ? 4 : -1
									)
								)
							)
						)
					);

	return 1;
}

bool sequence::is_base(int pos) {
	if ((pos<0)||(pos>=length)) return 0;
	if ((asints[pos]==4)||(asints[pos]==-1)) return 0;
	else return 1;
}
int sequence::next_base(int pos) {
	while (pos<length) {
		if (is_base(pos)) return pos;
		pos++;
	}
	return -1;
}
int sequence::ungapped_pos(int pos) {
	int p; int ugp=pos;
	for (p=0;p<pos;p++) if (asints[p]==4) ugp--;
	//for (p=0;p<pos;p++) cout << asints[p] <<" ";
	return ugp;
}

node::node() {
	leaf=0;
	root=0;
	//seq_post = new double [4];
	//int b;
	//for (b=0;b<4;b++) seq_post[b]=0.0;
	base=4;
	lbr=abr=rbr=INITIAL_BRANCH;

}

struct node * node::sister() {
	return (anc->left->number == number) ? anc->right : anc->left;
}

bool node::print() {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	int b;
	if (root) cout << "\troot"<<endl;
	else cout <<"\tn"<<anc->number<<" "<<endl;
	cout <<"\t|"<<endl;
	if (leaf) cout << "\tleaf" <<endl;
	else cout <<"n"<< left->number << "\t-\tn" <<right->number<<endl;
	return 1;
}

tree::tree() {
	Nnodes=0;
}


bool tree::print(struct sequence * seq) {
	cout << tree_string(root, seq) << endl;
	//int n=0;
// 	for (n=0;n<Nnodes;n++) {
// 	//while (!done) {
// 		cout << "node " <<n<<endl;
// 		nodelist[n].print();
// 	}
	return 1;
}


bool tree::copy_topology(struct tree * t) {
	int i;
	Nnodes=t->Nnodes;
	nodelist = new struct node [Nnodes];
	for (i=0;i<Nnodes;i++) {
		nodelist[i].number = t->nodelist[i].number;
		if (t->nodelist[i].leaf) { nodelist[i].leaf=1; nodelist[i].abr = t->nodelist[i].abr; }
		else {
			nodelist[i].left = &nodelist[t->nodelist[i].left->number];
			nodelist[i].lbr = t->nodelist[i].lbr;
			nodelist[i].right = &nodelist[t->nodelist[i].right->number];
			nodelist[i].rbr = t->nodelist[i].rbr;
			nodelist[i].leaf=0;
			nodelist[i].abr = t->nodelist[i].abr;
		}
		if (t->nodelist[i].root) { root = &nodelist[i]; nodelist[i].root=1; }
		else { nodelist[i].anc = &nodelist[t->nodelist[i].anc->number]; nodelist[i].root=0;}
	}
	return 1;
}


string tree::tree_string( struct node * n, struct sequence * seq) {
	if (n->leaf) return seq[n->number].name;
	else return "("+tree_string(n->left,seq) +","+ tree_string(n->right,seq)+")";
}

bool tree::compute_PJC() {
	int a,b,n;
	for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b] = new double * [Nnodes];
	for (n=0;n<Nnodes;n++) for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][n] = new double [Nnodes];
	for (n=0;n<Nnodes;n++) if (!nodelist[n].root) for (a=0;a<4;a++) for (b=0;b<4;b++) {
		if (a != b) psubs[a][b][nodelist[n].anc->number][n] = 0.25 - 0.25*exp(-(4.0/3.0)*nodelist[n].abr);
		else psubs[a][b][nodelist[n].anc->number][n] = 0.25 + 0.75*exp(-(4.0/3.0)*nodelist[n].abr);
	}
	//for (n=0;n<Nnodes;n++) if (!nodelist[n].root) for (a=0;a<4;a++) for (b=0;b<4;b++) {
		//cout << a<<"->"<<b<< " "<<nodelist[n].anc->number<<"<"<<nodelist[n].abr<<">"<< n<<":";
	//	cout<<psubs[a][b][nodelist[n].anc->number][n]<<endl;
	//}
	return 1;
}
bool tree::compute_PHKY(double * f, double kap) {
	int a,b,n;
	void matrix_exponential (double ** a, double ** b, int n);
	for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b] = new double * [Nnodes];
	for (n=0;n<Nnodes;n++) for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][n] = new double [Nnodes];
	for (n=0;n<Nnodes;n++) if (!nodelist[n].root) {
		double ** R = new double * [4];
		double ** P = new double * [4];
		double tot=0.0;
		for (a=0;a<4;a++) {

			R[a] = new double [4];
			P[a] = new double [4];
			for (b=0;b<4;b++) if (a != b) {
				R[a][b] = f[b];
				//R[a][b]=1.0;
				if (((a==0)&&(b==2)) || ((a==2)&&(b==0)) || ((a==1)&&(b==3)) || ((a==3)&&(b==1)))
					R[a][b] *= kap;
				tot+=(f[a]*R[a][b]);
			}
		}
		for (a=0;a<4;a++) {
			double idrate=0.0;
			for (b=0;b<4;b++) if (a!=b) {
				R[a][b] *= (nodelist[n].abr/tot);
				idrate -= R[a][b];
			}
			R[a][a] = idrate;
		}
		//for (a=0;a<4;a++) { for (b=0;b<4;b++) cout << R[a][b]<< " "; cout <<endl; }
		//cout << endl;
		matrix_exponential(R, P, 4);
		for (a=0;a<4;a++) for (b=0;b<4;b++) psubs[a][b][nodelist[n].anc->number][n] = P[a][b];
		//for (a=0;a<4;a++) { for (b=0;b<4;b++) cout << psubs[a][b][nodelist[n].anc->number][n]<< " "; cout <<endl; }
		//cout << endl;
		delete [] R;
		delete [] P;
	}

	return 1;
}

bool tree::initialize_pairwise(double dis) {
	int n;
	Nnodes=3;
	nodelist=new struct node [3];
	nodelist[2].number=2;
	root= &nodelist[2];
	nodelist[2].root=1;
	nodelist[2].left = &nodelist[0];
	nodelist[2].right= &nodelist[1];
	nodelist[2].lbr=nodelist[2].rbr=0.5*dis;
	for (n=0;n<2;n++) {
		nodelist[n].number=n;
		nodelist[n].leaf=1;
		nodelist[n].anc= &nodelist[2];
		nodelist[n].abr=0.5*dis;
	}
	return 1;
}
bool tree::read_from_file(string treefile) {
	//cout << " reading tree file"<<endl;
	ifstream treedata(treefile.c_str());
	if (!treedata) return 0;
	char comm;
	int thisnum,num;
	treedata>>comm;
	if (comm == 't') {
		treedata>>Nnodes;
		//cout << Nnodes <<endl;
		nodelist=new struct node [Nnodes];
	}
	int i=0;
	while (i < Nnodes) {
		treedata>>comm; //print();
		//if (!comm) return 0;
		if ((comm == 'r')||(comm == 'n')) {
			treedata>> thisnum;
			//cout <<thisnum<<endl;
			//cout << i << endl;
			nodelist[thisnum].number=thisnum;
			if (comm=='r') {
				root = &nodelist[thisnum];
				nodelist[thisnum].root=1;
			}
			treedata>>comm;
			if (comm == 'l') {
				treedata>>num;
				nodelist[thisnum].left = &nodelist[num];
				nodelist[num].leaf=1;
				nodelist[num].number=num;
				nodelist[num].anc = &nodelist[thisnum];
				//nodelist[num].lbr=nodelist[num].rbr=-1.0;
				i++;
			}
			else if (comm == 'n') {
				treedata>>num;
				nodelist[thisnum].left = &nodelist[num];
				nodelist[num].anc= &nodelist[thisnum];
			}
			treedata>>comm;
			if (comm == 'l') {
				treedata>>num;
				nodelist[thisnum].right= &nodelist[num];
				nodelist[num].leaf=1;
				nodelist[num].number=num;
				nodelist[num].anc= &nodelist[thisnum];
				i++;
			}
			else if (comm == 'n') {
				treedata>>num;
				nodelist[thisnum].right= &nodelist[num];
				nodelist[num].anc= &nodelist[thisnum];
			}
			i++;
		}
		//else {cout << i ; print();}
	//print();
	}
	treedata>>comm;
	if (comm == 'b') {
		for (i=0;i<Nnodes;i++) if (!nodelist[i].leaf) {
			treedata>>num; if (num!=i) return 0;
			treedata>>nodelist[num].lbr;
			treedata>>nodelist[num].rbr;
			nodelist[num].left->abr=nodelist[num].lbr;
			nodelist[num].right->abr=nodelist[num].rbr;
		}
	}
	//cout << "read "<<endl;
	return 1;

}

alignment::alignment() {
	length=0;
	Nseqs=0;
	gapped=0;
	bad=0;
	seq = new struct sequence [MAX_ALIGN_SIZE];
	bgfreq = new double [4];
}
bool alignment::read_bg_freq(string freqfile) {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	ifstream data(freqfile.c_str());
	if (!data) return 0;
	//else if (PRINT_TO_SCREEN) cerr << "reading "<<freqfile<<" ";
	int b;
	for (b=0;b<4;b++) {
		char base;
		data>>base;
		if (B[b] == base) {
		  data>>bgfreq[b];
		  //if (PRINT_TO_SCREEN) cerr<<B[b]<< ": "<<bgfreq[b]<<" ";
		}
		else { cerr << "coudn't parse " << freqfile <<endl; return 0; }
	}
	// if (PRINT_TO_SCREEN) cerr<<endl;
	char comm;
	data>>comm;
	if (comm == 'f')  data >> FUDGE;
	else if (comm == 'k') data >> kappa;
	//	if (PRINT_TO_SCREEN) cerr << "ki/ks=" << 1.0/FUDGE<<endl;

	return 1;

}

bool alignment::print() {
	int b; char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	cout << name << " " <<length<<" bp (" <<(length-gapped-bad) << " ungapped)" <<endl;
	for (b=0;b<4;b++) cout << B[b] <<": "<<bgfreq[b]<<" ";
	cout <<endl;
	return 1;
}

void alignment::integerize() {
	int p, b,s;
	asints = new long int [length];
	long int seq_to_num (int * seq, int start, int ml);
	//long int gapped_seq_to_num (int * seq, int start, int ml);
	for (p=0;p<length;p++) {
		bool gapped=0;
		for (s=0;s<Nseqs;s++) if (!seq[s].is_base(p)) gapped=1;
		if (gapped) { asints[p] =-1; }
		else {
			int * col = new int [Nseqs];
			for (s=0;s<Nseqs;s++) col[s] = seq[s].asints[p];
		//asints[p] = gapped_seq_to_num(col,0,Nseqs);
			asints[p] = seq_to_num(col,0,Nseqs);
		}
		//else asints[p]=-1;
		//cout << p << " " << asints[p] << " ";
		//for (s=0;s<Nseqs;s++) cout << " "<<seq[s].asints[p];
		//cout << endl;
	}
}

bool alignment::compute_nextgap() {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	int p,b,s;
	int lastgap=length;
	int lastbad=length;
	leftend=0; rightend=-1;
	//double total=0.0;
	nextgap= new int [length];
	nextbad= new int [length];
	double totalbases=0.0;
	gapped=0;
	bad=0;
	for (b=0;b<4;b++) bgfreq[b]=0.0;
	for (p=length-1;p>=0;p--) {
		for (s=0;s<Nseqs;s++) {
			if ((seq[s].bases[p] != 'A') && (seq[s].bases[p] != 'C')
				&& (seq[s].bases[p] != 'G') && (seq[s].bases[p] != 'T')) {
				if (seq[s].bases[p] == '-') lastgap=p; //seq[s].asints[p]=-1;
				else lastbad=p;
			}
			else for (b=0;b<4;b++) if (seq[s].bases[p] == B[b]) {
				bgfreq[b]+=1.0;
				totalbases+=1.0;
			}

		}
		nextgap[p]=lastgap;
		nextbad[p]=lastbad;
		if (p==lastgap) gapped ++;
		else leftend=p;
		if (p==lastbad) bad ++;
		if (rightend==-1) if (gapped != (length-p)) rightend=p;
	}
	for (b=0;b<4;b++) bgfreq[b] /= totalbases;
	//for (p=0;p<length;p++) {
		//for (s=0;s<Nseqs;s++) cout << seq[s].bases[p];
		// cout << p <<" "<< nextgap[p] <<" "<<endl;
	//}
	return 1;
}

float alignment::div() {
	float id=0.0;
	//float ungapped=0.0;
	int p,s1,s2;
	for (p=0;p<length;p++) {
		int thisposdiff=0;
		bool thisgapped=0;
		for (s1=0;s1<Nseqs;s1++) if (seq[s1].asints[p]==4) thisgapped=1;
		if (!thisgapped) {
			//for (s1=0;s1<Nseqs;s1++)
			for (s2=1;s2<Nseqs;s2++)
				if (seq[0].asints[p] != seq[s2].asints[p]) thisposdiff=1;
			//ungapped+=1.0;
			if (!thisposdiff) id+=1.0;
		}
	}
	id /= ((float) (length-gapped-bad));
	return (1.0-id);
}

bool alignment::read_tree(string filename) {
	int p,n,b;
	//if (PRINT_TO_SCREEN) cerr << "reading tree from " <<filename<<"...";
	//phyl= new struct tree [length];
	if (!phyl.read_from_file(filename)) return 0;
	treelength=0.0;
	for (n=0;n<phyl.Nnodes;n++) if (!phyl.nodelist[n].root)
		treelength += phyl.nodelist[n].abr;
	//cout << global_tree.root->lbr<<endl;
	//double eqf[4]; eqf[0]=eqf[1]=eqf[2]=eqf[3]=0.25;
	//double loglh=0.0;
	//for (p=0;p<length;p++) if (nextgap[p]!=p) {
	//for (p=0;p<1;p++) {
		//if (!phyl[p].read_from_file(filename)) return 0;
		//phyl[p].print();
		//ALI[0].phyl[p].print();
		//for (n=0;n<phyl[p].Nnodes;n++) {
		//	if (phyl[p].nodelist[n].leaf) {
		//		phyl[p].nodelist[n].base=seq[n].asints[p];
				//for (b=0;b<4;b++) ALI[0].phyl[p].nodelist[n].seq_post[b]=0.0;
 		//		phyl[p].nodelist[n].seq_post[phyl[p].nodelist[n].base]=1.0;
		//	}
		//	else if ((global_tree.root->lbr==INITIAL_BRANCH)&&(global_tree.root->rbr==INITIAL_BRANCH)) {
		//			for (b=0;b<4;b++) phyl[p].nodelist[n].seq_post[b]=bgf[b];
		//	}
		//	else if (EVOL_MODEL==1) compute_post(phyl[p].nodelist + n, bgf);
		//	else if (EVOL_MODEL==0) compute_post_JC(phyl[p].nodelist + n, bgf);
		//}
	//	loglh+=log(compute_post_JC(phyl[p].root,bgf));
		//compute_post(phyl[p].root->left);
		//compute_post(phyl[p].root->right);
		//loglh+=compute_post(phyl[p].root, bgf);
	//}
	//cout << "initial bg log liklihood " << loglh <<endl;
	//if (EVOL_MODEL==0) {
	//if ((global_tree.root->lbr==INITIAL_BRANCH)&&(global_tree.root->rbr==INITIAL_BRANCH))
	//	optimize_branch_lengths(EPSILON,MAX_ITER);
	if ((phyl.root->lbr!=INITIAL_BRANCH)&&(phyl.root->rbr!=INITIAL_BRANCH))
	  //	if (PRINT_TO_SCREEN) cerr << "branch lengths found" <<endl;
	  //else if (PRINT_TO_SCREEN) cerr << "branch lengths not found" <<endl;
	//}
	//loglh=0.0;
	//for (p=0;p<length;p++) if (nextgap[p]!=p)
	//	loglh+=log(compute_post_JC(phyl[p].root,bgf));
	//cout << "maximum bg liklihood "<<loglh<<endl;
	return 1;
}

bool alignment::read_fasta(string filename) {
	ifstream data(filename.c_str());
	if (!data) return 0;
	//else if (PRINT_TO_SCREEN) cerr << "reading "<<filename<<"...";
	//string seqdata;
	Nseqs=-1;
	bool foundmax=0;
	int p; bool ext=1;
	for (p=filename.size();p>=0;p--) { //get rid of the last extension
		if (ext) { if (filename[p]=='.') ext=0; }
		else name=filename[p] + name;
		//cout << filename[p] <<" "<<ext<<" "<<name<<endl;
	}
	if (ext) name=filename; // if there was no extension
	//name[filename.size()]='\0';
	//name=filename;
	while ((!data.fail())&&(!foundmax)) {
		string line;
		getline(data,line);
		int linelength=line.size();
		//cout << linelength << endl;
		if ( line[0] == '>' ) {
			Nseqs++;
			if (Nseqs<MAX_ALIGN_SIZE) { for (p=1;p<linelength;p++) seq[Nseqs].name += line[p]; }
			else foundmax=1;
			//if (Nseqs>0) if (seq[Nseqs].bases.size()>length)
			//		length=seq[Nseqs].bases.size();
				//seq[Nseqs].name = new char [linelength];
				//for (p=1;p<linelength;p++) {
				//	seq[Nseqs].name+=line[p];
					//cout << p << " " << line[p]<< endl;
				//}
				//seq[Nseqs].name[linelength-1]='\0';
				//seq[Nseqs].name=line.c_str();
				//cout << Nseqs<< " "<<seq[Nseqs].name << endl;
			//}
		}
		else if ((linelength>0)&&(!foundmax)) {
			seq[Nseqs].bases += line;
			//cout << line << endl;
			//for (p=0;p<linelength;p++) if (line[p] == '\n') cout << p <<endl;
		}
	}
	int s;
	if (Nseqs<MAX_ALIGN_SIZE) Nseqs++;
	for (s=0;s<Nseqs;s++) {
		seq[s].integerize();
		if (length < seq[s].length) length=seq[s].length;
	}
	compute_nextgap();
	integerize();
	//P_BG_COMPUTED=0;
	//pbg=new double [length];
	kappa=2.0;
	double div_to_subs_JC (double pid);
	if (Nseqs==2) { // initialise the pairwise tree

		treelength=div_to_subs_JC(1.0-div()); //JC distance
		phyl.initialize_pairwise(treelength);
		//double d=div();
		//double r=(1.0-d)/d;
		//double k=-0.75*log((1.0-(r/3.0))/(1.0-r));
		//phyl.initialize_pairwise(k);
		//int b;
		//phyl= new struct tree [length];
		//for (p=0;p<length;p++) if (nextgap[p]!=p) {
		//	phyl[p].Nnodes=3;
		//	phyl[p].nodelist= new struct node [3];
		//	phyl[p].nodelist[2].number=2;
		//	phyl[p].root= &phyl[p].nodelist[2];
		//	phyl[p].nodelist[2].root=1;
		//	phyl[p].nodelist[2].left = &phyl[p].nodelist[0];
		//	phyl[p].nodelist[2].right= &phyl[p].nodelist[1];
		//	for (b=0;b<4;b++) phyl[p].nodelist[2].seq_post[b]=bgfreq[b];
		//	for (s=0;s<Nseqs;s++) {
		//		phyl[p].nodelist[s].number=s;
		//		phyl[p].nodelist[s].leaf=1;
		//		phyl[p].nodelist[s].anc= &phyl[p].nodelist[2];
		//		phyl[p].nodelist[s].base=seq[s].asints[p];
		//		phyl[p].nodelist[s].seq_post[phyl[p].nodelist[s].base]=1.0;
		//	}
		//}
		//phyl.Nnodes=3;
		//phyl.nodelist=new struct node [3];
		//phyl.nodelist[2].number=2;
		//phyl.root= &phyl.nodelist[2];
		//phyl.nodelist[2].root=1;
		//phyl.nodelist[2].left = &phyl.nodelist[0];
		//phyl.nodelist[2].right= &phyl.nodelist[1];
		//for (s=0;s<Nseqs;s++) {
		//	phyl.nodelist[s].number=s;
		//	phyl.nodelist[s].leaf=1;
		//	phyl.nodelist[s].anc= &phyl.nodelist[2];
		//}
	}
	//if (PRINT_TO_SCREEN) { cerr << "done" <<" %id="<<1.0-div()<<" k="<<treelength<<endl; }
	if (Nseqs>1) return 1;
	else return 0;
}
motif::motif() {
	//name = new char [MAX_NAME_LENGTH];
	//f = new double * [4];
	pri=MOTIF_PRIOR;
}
double motif::Pr(int s[],int start) {
	int pos; double p=1.0;
	for (pos=0;pos<length;pos++) {
		int b=s[start+pos];
		if ((b!=-1)&&(b!=4)) p *= f[pos][b];
	}
	return p;
}

struct motif motif::make_rc() {
	struct motif rc;
	int b, pos;
	int revcom[4]; revcom[0]=3; revcom[1]=2; revcom[2]=1; revcom[3]=0;
	rc.length=length;
	rc.pri=pri;
	rc.K = new double [length];
	rc.f = new double * [length];
	for (pos=0;pos<length;pos++) rc.f[pos] = new double [5];
	for (pos=0;pos<length;pos++) {
		for (b=0;b<4;b++) rc.f[pos][b]=f[length-pos-1][revcom[b]];
		rc.K[pos]=K[length-pos-1];
		rc.f[pos][4]=f[length-pos-1][4];
	}
	int a;
	return rc;
}
bool motif::initialize_neutral_motif(int w, double pseudocounts) {
	int b, pos;
	length=w;
	//cout << "initialising matrix of size "<< length <<endl;
	f = new double * [length];
	K = new double [length];
	for (pos=0;pos<length;pos++) f[pos] = new double [5];
	for (pos=0;pos<length;pos++) {
		for (b=0;b<4;b++) f[pos][b]=0.25*pseudocounts;
		K[pos]=1.0;
		f[pos][4]=1.0;
	}
	return 1;
}


bool motif::print() {
	char B[4]; B[0]='A'; B[1]='C'; B[2]='G'; B[3]='T';
	int b,pos;
	cout << "pos\t";
	for (b=0;b<4;b++) cout << B[b] <<"\t";
	cout << "Kp"<<"\t"<<"Khb"<<endl;
	for (pos=0;pos<length;pos++) {
		cout << pos << "\t";
		for (b=0;b<4;b++) cout << 0.01*((int) (100.0*f[pos][b])) <<"\t";
		cout << K[pos]<< "\t"<<halpern_bruno(pos)<<endl;

	}
	return 1;
}
double motif::halpern_bruno(int pos) {
	double hbrate=0.0;
	int b1,b2;
	for (b1=0;b1<4;b1++) for (b2=0;b2<4;b2++) if (b2!=b1) {
		if (f[pos][b1]==f[pos][b2]) hbrate+=f[pos][b1];
		else hbrate+=(f[pos][b1]*log(f[pos][b2]/f[pos][b1])/(1-(f[pos][b1]/f[pos][b2])));
	}
	return (hbrate/3.0);
}
bool motif::read_from_file(string filename) {
	ifstream data(filename.c_str());
	if (!data) { cout << filename<<" not found" << endl; return 0; }
	length=-1;
	while (!data.fail()) {
		string line;
		getline(data,line);
		//cout << line << endl;
		length++;
	}
	if (PRINT_TO_SCREEN) cout << "motif of length " << length << " found" << endl;
	int b,pos;
	f = new double * [length];
	for (pos=0;pos<length;pos++) f[pos] = new double [5];
	//int pos=0;
	ifstream datatoo(filename.c_str());
	for (pos=0;pos<length;pos++) {
		for (b=0;b<4;b++) datatoo>>f[pos][b];
		f[pos][4]=1.0;
	}

	K=new double [length];
	for (pos=0;pos<length;pos++) K[pos]=1.0;
	name=filename;
	return 1;
}
monkey2.0/monkey-doc.txt0000644000163100010400000003246610466225440013657 0ustar  alanallMonkey (doc file updated May 2006)

version 2.0 
copyright 2004-2006 by Alan Moses. This software is provided "as is" without
warranty of any kind.  The author assumes no responsibility for the results it
produces or conclusions based thereupon.  It is distributed free of charge for
academic use only.  Permission to copy and use it is granted free of charge
provided that no fee is charged and this copyright notice is not removed.



Contents

I	Introduction
II	Getting started with monkey
III	The matrix file
IV	The alignment file
V	The tree file
VI	The .occ output file
VII 	More options
VIII	rMonkey


I introduction

This program is designed to search aligned non-coding DNA for matches to a known 
motif matrix.   It runs on multiple alignments (2 or more sequences) of non-coding DNA
and provides a variety of information of about the putative motif instances that
it encounters.  

The basic operation of the program is to scan a single 'reference' genome (which can be 
specified by the user) ignoring gaps.  For all regions of the 'reference' genome that may
be instances of the motif, the corresponding region of the alignment is examined using
a model of motif evolution.  The sequences from each genome corresponding to this
region are reported, along with the various scores and p-values that have been computed.  
Monkey uses a statistical model to report a pvalue that is an estimate of the probability
of the observed match in random genomic alignments.


II Getting started with monkey

To install monkey on linux 

type 

tar -xvf monkey2.0.tar

this will create a directory called monkey2.0 .  Change to that directory, and
compile the program.

$ cd monkey2.0
$ make all

Monkey can now be run by typing ./Monkey

Monkey requires two command line arguments: the names of the matrix file and the fasta 
alignment file. The formats for these will be discussed below. The simplest
run of monkey might look like  

$ ./Monkey gal4sites.matrix example/YBR018C_upstream.aln.fasta 

which scans the alignment example/YBR018C_upstream.aln.fasta for sequences
matching the matrix in gal4sites.matrix.  The results should appear in a new
file called example/YBR018C_upstream.aln.occ.

Monkey can be run on windows under cygwin (http://www.cygwin.com/). Once the
make and g++ compiler have been installed, just follow the instructions for
linux.


III The matrix file

Monkey must read in the frequency matrix representing the sequence specificty of the motif.
The format for this 'matrix' file is simply four columns of numbers, with each row
representing the frequencies of A, C, G and T for one position in the motif.
Frequencies should be greater than 0 and sum to one. An example matrix file is 
provided with the distribution.


IV The alignment file

Monkey reads in a 'multiple fasta alignment.'  The format for file is as a fasta file, with 
gapped positions represented by '-' and nucleotides represented as A, C, G and
T. Other characters, such as 'a','c','g','t' or 'N' may be included, but
columns containing these will be ignored. Examples of these files are included with the 
distribution and the Monkey distribution includes a perl script (clustal2fasta.pl) to convert 
clustalw format files to multiple fasta. Typing 

$ perl scripts/clustal2fasta.pl filename.aln

will create a multiple fasta file called filename.aln.fasta.


V The tree file
 
If the fasta alignment file contains more than 2 sequences a tree file is required to 
specify the tree topology and branchlengths. This can be estimated using one of many
availble phylogentic analysis packages (such as paml).  The tree file is structured as 
follows

t 7
r 6     n 5     l 0
n 5     n 4     l 2
n 4     l 1     l 3
b
4       0.09787 0.13253
5       0.12665 0.23934
6       0.154675        0.154675


The first section of the file, beginning with t, specifies the tree topology. t is followed by
the total number of nodes in the tree. r stands for the root, n for interior nodes and l for 
leaves. The line 

r 6     n 5     l 0

means that the root is node 6 and it joins the internal node 5 with the leaf 0.  The numbers of
the leaves MUST correspond to the order of the sequences in the corresponding multi-fasta 
alignment file, 0 being the first sequence, 1 the next and so on.  The numbers of the other nodes
must be self-consistent, but are arbitrary.

The second second section of the tree file starts with b and lists the branch lengths for the 
interior nodes.  The line 

4       0.09787 0.13253

means that the left branch of node 4 (to leaf 1) is of length 0.09787 substitutions per site, and
the right branch (to leaf 3) is of length 0.13253.  Phylogenetic trees can be computed by many 
freely available programs (such as the PAML package.)


The distribution contains a perl script (dnd2tree.pl) to convert 'Newick' format files to the tree 
format used for monkey.  This script requires that the species occur in left
to right order corresponding to the top to bottom order of the fasta file.  In
addition, branchlengths must be specified in 'X.X' format, (such that 0 must be
written as 0.0 or 0.000), the tree be strictly bifurcating (including the
root) and does not support internal node names.

Monkey2.0 includes a parsing script 'nhx2tree.pl' to convert newick files to 
the format used for monkey.  This script should handle the full specification of 
the extended newick format (http://www.genetics.wustl.edu/eddy/forester/NHX.html) 
and will produce warnings if the tree is not strictly bifurcating, or does not 
contain branchlengths (requirements for monkey and Rmonkey).  As with dnd2tree.pl,
the left to right order of species names in the tree should agree with the top
to bottom order in the multi-fasta file.  

 
If the alignment contains more than 2 sequences,  Monkey will by default 
look for a file with the same name as the multiple fasta file, except with the extension '.tree'. 
It is possible to specify a tree file using the -tree option. For example, typing
 
$ ./Monkey gal4sites.matrix example/YBR018C_upstream.aln.fasta -tree all.tree 

uses all.tree instead of example/YBR018C_upstream.aln.tree

For pairwise alignments, no tree is required and Monkey will by default calculate the distance between 
the two sequences.  If -tree is used it should be followed by a number - monkey will
take this to be the distance between the sequences (in substituions per site.)



VI The .occ output file

Monkey will produce a file with the same name as the multiple fasta alignment file that 
was given as input, except the extension will now be .occ instead of .fasta.  This 
file is in tab delimited format.  Each row in this file represents a putative instance of
the motif. It is possible to specify a different filename using the -out option.  For example, 

$ ./Monkey gal4sites.matrix example/YBR018C_upstream.aln.fasta -out file.occ

will write the output to file.occ.

The first 7 columns in the .occ file contain 1) the name of the fasta file and
matrix file, 2) a probabilistic estimate of which strand the binding site is on,
where 1 is the forward strand, 0 is the reverse, 3) the position in the
alignment, 4) the position in the ungapped version of the reference sequence,
5) the segment of reference sequence, 6) the single sequence liklihood ratio
score and 7) the p-value associated with that score.  

The next group of columns contain the sequences and single sequence liklihood
ratio scores for all of the other sequences in the alignment.  

Finally, the last three columns contain the p-value associated with the score
under the evolutionary model, the % identitiy of the subsequences in the
alignment, and the % ungapped.  If the percent ungapped is high, monkey will
heuristically modify the alignment and calculate a p-value.  If the percent
gaps is large, monkey will print out '-'s instead of the sequences and
p-values.


VII More options

-help
To see all the options currently available type

$ ./Monkey -help

-m 
Monkey allows the use of several scores based on different assumptions about
the evolution of the motif. The default is to use the Halpern-Bruno model.

-b
This specifies the background substitution model.  The Default is to use the
Jukes-Cantor (JC) model, but the HKY model is also available.  If the HKY
model is used (by adding -b HKY to the command line) the transition
transversion rate ratio, or kappa, should also be specified using the -kappa
option.  The default is kappa=2.0 . 

-list
If the same tree is to be used for searches of many alignment files, there is
a considerable speed advantage to scanning them all at once.  Monkey allows
the user to provide a file containing a list of alignments. 

$ ./Monkey gal4sites.matrix eg.list -tree all.tree -list 

will scan all of the alignments listed in eg.list using the tree found in
all.tree, and write the output to a single file called eg.list.occ. 

-scr
Prints out the evolutionary score as well as the associated p-value.

-ref
specifies which sequence in the alignment file is to be treated as the
reference.  This is done by adding -ref n, where n is the position in the file
of the desired reference sequence, starting with 0.

-cut
This option can be used to reduce the total amount of hits produced in a scan
of many alignments.  This is the minimum score that a binding sites must have
in the reference genome before it is considered for evolutionary analysis.
The default is to use a score of zero.

-freq 
To specify a background frequency file use -freq followed by the filename
This file has the following form.

A 0.3 
C 0.2
G 0.2
T 0.3

This means that average frequencies in the intergenic region for A, C, G and T are 30%, 20%
20% and 30% respectively.  If a file is not specified monkey will calculate
these frequencies from the sequences given.


VIII rMonkey


Rmonkey is another program to detect matches to specificity matrices is multiple 
sequence alignments.  


Differences between rMonkey and the original monkey.  

rMonkey uses a different heuristic to identify binding sites in alignments (see 
Moses et al 2006).  Breifly, rather than adjusting the alignment based on the 
number of gaps in the cloumns that correspond to the binding site, rmonkey 
selects an reference match, and then choses the best overlapping matches to 
the matrix in each other sequence.  This makes rmonkey more 'agressive' at 
assigning orthology to binding sites, and it can often consider non-orthologous 
sequences to be aligned binding sites,  Unlike the original monkey, in which the 
'reference' species was held constant for a given run, rmonkey includes an option 
(specified by adding -any to the command line) which will allow the 'reference' 
sequence for matches to be in any species.

In addition to the p-value assusiated with the evolutionary generalization of the 
likelihood ratio score (as in monkey), rmonkey computes p-values associated with 
a conditional likelihood ratio, which assumes a match in one 'reference' species 
has already been observed (see Moses et al, 2006).   These p-values can be take 
a long time and use a lot of memory to calculate; they can be disabled using the 
-M option.

A final difference between rmonkey and monkey is that cutoffs in rmonkey are in 
p-values rather than likelihood scores, and p-values are returned for each sequence 
by default.  Single sequence p-values are computed by a different method, and will
differ slightly (but be more accurate) than those produced by the original monkey.  


Inputs

The input files to rmonkey are the same as for the original monkey, but some of 
the options have slightly different behaviours in rmonkey.

-scr 
produces scores for each sequence and the evolutionary log likelihood ratio instead
of p-values

-cut
specifies a single sequence p-value cutoff for a match to be included in 
evolutionary analysis


Outputs

The rmonkey output file is very similar to that of the original monkey.  rMonkey 
produces a tab delimted file with one row corresponding to each match to the matrix.  
The first 5 columns contain the names, strand,  aligned and ungapped positions in 
the reference sequence, and the sequence for the match in the reference.  Note that 
positions in the aligned sequences are now given as 'start'-'stop'.  The next 
column gives the single species p-value for the match in the reference. The next 
sets of columns give positions, sequence and single species p-values for each of 
the orthologous matches.  The next column gives the p-value associated with the 
likelihood ratio score for a conserved match.  This is the p-value computed in 
the original monkey.  If a sequence contains too many gaps, rmonkey reports a 'g'
instead of the '-' character used in the original monkey.

Finally, the last three columns give the conditional likelihood ratio (T statistic,
see Moses et al. 2006), and associated p-values under either the HB null hypothesis
(a test for lack of conservation) or the HKY null hypothesis (a more conservative 
test for conservation).


When should you use rMonkey?

-If you are trying to find all the binding sites that might be conserved (and 
are willing to include some that are actually not conserved)

-If you are trying to be conservative in identifying non-conserved bindng sites.

When should you not use rMonkey?

-If you are doing an evolutionary analysis of the bases in the binding sites - 
in rmonkey the alignment is modified in a way that depends on the matrix.

-If you are interested in a region that contains many partially overlapping 
binding sites, and want to preserve the orthology relationships between them.