1 2%{ 3#include "dnaalign.h" 4#include "hitlist.h" 5 6%} 7 8 9struct DnaMatchPara 10DPRunImpl * dpri 11DnaMatrix * mat 12DnaStartEnd * dse 13int gap 14int ext 15 16%{ 17#include "dnamatcher.h" 18 19 20void show_help_DnaMatchPara(FILE * ofp) 21{ 22 fprintf(ofp,"Dna Matching Parameters\n"); 23 fprintf(ofp," -dm_match [5] match score\n"); 24 fprintf(ofp," -dm_mismatch [-4] mismatch score\n"); 25 fprintf(ofp," -dm_gappen [5] gap open penalty\n"); 26 fprintf(ofp," -dm_extpen [1] gap extension penalty\n"); 27 28 show_help_DPRunImpl(ofp); 29} 30 31DnaMatchPara * new_DnaMatchPara_from_argv(int * argc,char ** argv) 32{ 33 DnaMatchPara * out; 34 int match = 5; 35 int mismatch = -10; 36 37 out = DnaMatchPara_alloc(); 38 39 strip_out_integer_argument(argc,argv,"dm_match",&match); 40 strip_out_integer_argument(argc,argv,"dm_mismatch",&mismatch); 41 42 assert(mismatch < 0 ); 43 44 out->dpri = new_DPRunImpl_from_argv(argc,argv); 45 out->mat = identity_DnaMatrix(match,mismatch); 46 out->dse = DnaStartEnd_from_policy("local"); 47 48 out->gap = 30; 49 out->ext = 20; 50 51 strip_out_integer_argument(argc,argv,"dm_gappen",&out->gap); 52 strip_out_integer_argument(argc,argv,"dm_extpen",&out->ext); 53 54 55 return out; 56} 57 58 59HitList * HitList_from_Sequence_SequenceSet_DNA(Sequence * query,SequenceSet * set,DnaMatchPara * p) 60{ 61 int i; 62 HitList * out; 63 HitPair * pair; 64 HitAln * aln; 65 66 AlnBlock * forward; 67 AlnBlock * reverse; 68 69 Sequence * rev; 70 71 char buffer[512]; 72 73 out = HitList_alloc_std(); 74 75 for(i=0;i<set->len;i++) { 76 77 rev = reverse_complement_Sequence(set->set[i]); 78 79 ckfree(rev->name); 80 sprintf(buffer,"%s.reverse",set->set[i]->name); 81 rev->name = stringalloc(buffer); 82 83 pair = HitPair_alloc_std(); 84 85 aln = HitAln_alloc(); 86 87 forward = make_align_dnaalign(query,set->set[i],p->mat,p->dse,-p->gap,-p->ext,-p->gap,-p->ext,p->dpri); 88 reverse = make_align_dnaalign(query,rev,p->mat,p->dse,-p->gap,-p->ext,-p->gap,-p->ext,p->dpri); 89 90 if( forward->score > reverse->score ) { 91 pair->query = hard_link_Sequence(query); 92 pair->target = hard_link_Sequence(set->set[i]); 93 aln->alb = hard_link_AlnBlock(forward); 94 } else { 95 pair->query = hard_link_Sequence(query); 96 pair->target = hard_link_Sequence(rev); 97 aln->alb = hard_link_AlnBlock(reverse); 98 } 99 100 101 102 add_HitPair(pair,aln); 103 104 aln->raw_score = pair->raw_score = aln->alb->score; 105 add_HitList(out,pair); 106 107 108 free_AlnBlock(forward); 109 free_AlnBlock(reverse); 110 111 free_Sequence(rev); 112 113 } 114 115 return out; 116} 117