1 #include "dyna.h"
2 #include "genestats.h"
3 #include "geneparser4.h"
4 #include "cdnawise10.h"
5 #include "version.h"
6 #include "geneutil.h"
7 #include "geneoutput.h"
8
9 char * program_name = "cdnawise";
10
11 char * codon_file = "codon.table";
12
show_version(FILE * ofp)13 void show_version(FILE * ofp)
14 {
15 fprintf(ofp,"%s\nVersion: %s\nReleased: %s\nCompiled: %s\n",program_name,VERSION_NUMBER,RELEASE_DAY,COMPILE_DATE);
16
17
18 fprintf(ofp,"\nThis program is freely distributed under a Gnu Public License\n");
19 fprintf(ofp,"The source code is copyright (c) EMBL (2001) and others\n");
20 fprintf(ofp,"There is no warranty, implied or otherwise on the performance of this program\n");
21 fprintf(ofp,"For more information read the GNULICENSE file in the distribution\n\n");
22 fprintf(ofp,"Credits: Ewan Birney <birney@sanger.ac.uk> wrote the core code.\n");
23 exit(63);
24 }
25
26
show_help(FILE * ofp)27 void show_help(FILE * ofp)
28 {
29
30 fprintf(ofp,"%s (%s)\n",program_name,VERSION_NUMBER);
31 fprintf(ofp,"cdnawise <cdna-file> <genomic-dna-file> in fasta format\n");
32 fprintf(ofp,"cDna sequence options\n");
33 fprintf(ofp," -s start position in dna\n");
34 fprintf(ofp," -t end position in dna\n");
35 fprintf(ofp,"Genomic Dna sequence options\n");
36 fprintf(ofp," -u start position in dna\n");
37 fprintf(ofp," -v end position in dna\n");
38 fprintf(ofp,"Match algorithm options\n");
39 fprintf(ofp," -g gap open penalty\n");
40 fprintf(ofp," -e gap ext penalty\n");
41
42 show_help_GeneModelParam(ofp);
43 show_help_GeneOutputPara(ofp);
44
45 show_help_DPRunImpl(ofp);
46 show_standard_options(ofp);
47 }
48
49
50
main(int argc,char ** argv)51 int main(int argc,char ** argv)
52 {
53 Sequence * cdna;
54 Sequence * gen;
55 Sequence * active_gen;
56 Sequence * active_cdna;
57
58 int i;
59 int dstart = -1;
60 int dend = -1;
61
62 int cstart = -1;
63 int cend = -1;
64
65 int gapopen = 12;
66 int gapext = 2;
67
68 CodonTable * ct = NULL;
69 CodonMatrixScore * cm = NULL;
70 RandomCodon * rndcodon = NULL;
71 RandomCodonScore * rndcodonscore = NULL;
72 DnaMatrix * dm = NULL;
73
74 DPRunImpl * dpri = NULL;
75
76 GeneModel * gm;
77 GeneModelParam * gmp;
78 GeneStats * gs;
79 GeneParser21 * gp21;
80 GeneParser21Score * gp21s;
81 GeneParser4Score * gp;
82
83 GeneOutputData data;
84 GeneOutputPara * geneout;
85
86
87 ComplexSequenceEvalSet * cdna_cses;
88 ComplexSequenceEvalSet * gen_cses;
89
90 ComplexSequence * cs_cdna;
91 ComplexSequence * cs_gen;
92
93 Genomic * gent;
94 GenomicRegion * gr;
95
96 CompMat * cmat;
97 CompProb * cprob;
98 char * matfile = "blosum62.bla";
99 Protein * trans;
100
101 PackAln * pal;
102 AlnBlock * alb;
103
104 FILE * ofp = stdout;
105
106 dpri = new_DPRunImpl_from_argv(&argc,argv);
107 gmp = new_GeneModelParam_from_argv(&argc,argv);
108
109 geneout = new_GeneOutputPara_from_argv(&argc,argv);
110
111 strip_out_integer_argument(&argc,argv,"u",&dstart);
112 strip_out_integer_argument(&argc,argv,"v",&dend);
113
114 strip_out_integer_argument(&argc,argv,"s",&cstart);
115 strip_out_integer_argument(&argc,argv,"t",&cend);
116
117 strip_out_integer_argument(&argc,argv,"g",&gapopen);
118 strip_out_integer_argument(&argc,argv,"e",&gapext);
119
120
121 strip_out_standard_options(&argc,argv,show_help,show_version);
122
123
124 ct = read_CodonTable_file(codon_file);
125
126 cmat = read_Blast_file_CompMat(matfile);
127 cprob = CompProb_from_halfbit(cmat);
128 cm = naive_CodonMatrixScore_from_prob(ct,cprob);
129
130 gm = GeneModel_from_GeneModelParam(gmp);
131
132 cdna = read_fasta_file_Sequence(argv[1]);
133 gen = read_fasta_file_Sequence(argv[2]);
134
135 if( dstart != -1 || dend != -1 ) {
136 if( dstart == -1 ) {
137 dstart = 1;
138 }
139 if( dend == -1 ) {
140 dend = gen->len;
141 }
142 active_gen = magic_trunc_Sequence(gen,dstart,dend);
143 } else {
144 active_gen = hard_link_Sequence(gen);
145 }
146
147 if( cstart != -1 || cend != -1 ) {
148 if( cstart == -1 ) {
149 cstart = 1;
150 }
151 if( cend == -1 ) {
152 cend = gen->len;
153 }
154 active_cdna = magic_trunc_Sequence(gen,cstart,cend);
155 } else {
156 active_cdna = hard_link_Sequence(gen);
157 }
158
159
160
161 rndcodon = RandomCodon_from_raw_CodonFrequency(gm->codon,ct);
162 fold_in_RandomModelDNA_into_RandomCodon(rndcodon,gm->rnd);
163
164 rndcodonscore = RandomCodonScore_from_RandomCodon(rndcodon);
165
166 assert(active_cdna);
167 assert(active_gen);
168
169 cdna_cses = default_cDNA_ComplexSequenceEvalSet();
170 gen_cses = new_ComplexSequenceEvalSet_from_GeneModel(gm);
171
172 cs_cdna = new_ComplexSequence(active_cdna,cdna_cses);
173 cs_gen = new_ComplexSequence(active_gen,gen_cses);
174
175 gp21 = std_GeneParser21();
176 GeneParser21_fold_in_RandomModelDNA(gp21,gm->rnd);
177 gp21s = GeneParser21Score_from_GeneParser21(gp21);
178 gp = GeneParser4Score_from_GeneParser21Score(gp21s);
179
180 dm = identity_DnaMatrix(Probability2Score(halfbit2Probability(1)),Probability2Score(halfbit2Probability(-1)));
181
182 assert(cs_cdna);
183 assert(cs_gen);
184 assert(gp);
185 assert(rndcodonscore);
186 assert(dm);
187 assert(dpri);
188
189 /* show_CodonMatrixScore(cm,ct,ofp);*/
190 /*
191 pal = PackAln_bestmemory_CdnaWise10(cs_cdna,cs_gen,gp,cm,rndcodonscore,dm,
192 Probability2Score(halfbit2Probability(-gapopen)),
193 Probability2Score(halfbit2Probability(-gapext)),
194 Probability2Score(halfbit2Probability(-5)),
195 +50,
196 NULL,
197 dpri);
198 */
199 pal = PackAln_bestmemory_CdnaWise10(cs_cdna,cs_gen,gp,cm,rndcodonscore,dm,
200 0,
201 0,
202 Probability2Score(halfbit2Probability(-5)),
203 +50,
204 NULL,
205 dpri);
206
207
208 alb = convert_PackAln_to_AlnBlock_CdnaWise10(pal);
209
210 gent = Genomic_from_Sequence(gen);
211 assert(gent);
212
213 gr = new_GenomicRegion(gent);
214 assert(gr);
215
216
217 add_Genes_to_GenomicRegion_GeneWise(gr,active_gen->offset,active_gen->end,alb,cdna->name,0,NULL);
218
219
220 data.pal = pal;
221 data.alb = alb;
222 data.gr = gr;
223 data.gen = gent;
224 data.ct = ct;
225
226 show_GeneOutput(&data,geneout,stdout);
227
228 return 0;
229 }
230
231