1 #include "statwise10.h"
2 #include "version.h"
3 #include "genestats.h"
4 #include "geneutil.h"
5
6
7 char * program_name = "statwise";
8
9 char * codon_table = "codon.table";
10
show_version(FILE * ofp)11 void show_version(FILE * ofp)
12 {
13 fprintf(ofp,"%s\nVersion: %s\nReleased: %s\nCompiled: %s\n",program_name,VERSION_NUMBER,RELEASE_DAY,COMPILE_DATE);
14 fprintf(ofp,"\nThis program is freely distributed under a Gnu Public License\n");
15 fprintf(ofp,"The source code is copyright (c) EMBL and others\n");
16 fprintf(ofp,"There is no warranty, implied or otherwise on the performance of this program\n");
17 fprintf(ofp,"For more information read the GNULICENSE file in the distribution\n\n");
18 fprintf(ofp,"Credits: Ewan Birney <birney@ebi.ac.uk>\n");
19 exit(63);
20 }
21
show_help(FILE * ofp)22 void show_help(FILE * ofp)
23 {
24 fprintf(ofp,"%s sequence-as-fasta\n",program_name);
25
26 show_help_GeneModelParam(ofp);
27
28 show_help_DPRunImpl(ofp);
29
30 show_standard_options(ofp);
31 }
32
33 CodonTable * ct;
34
id(int i)35 double id(int i)
36 {
37 return (double)i;
38 }
39
main(int argc,char ** argv)40 int main(int argc,char ** argv)
41 {
42 int i;
43
44 DPRunImpl * dpri = NULL;
45 GeneModelParam * gmp = NULL;
46 GeneModel * gm = NULL;
47
48 Sequence * seq;
49
50 RandomCodon * rc;
51 RandomModelDNA * rmd;
52 RandomCodonScore * rcs;
53
54
55 ComplexSequenceEval * splice5;
56 ComplexSequenceEval * splice3;
57 ComplexSequenceEvalSet * cses;
58 ComplexSequence * cseq;
59
60
61 SyExonScore * exonscore;
62
63 PackAln * pal;
64 AlnBlock * alb;
65
66 Genomic * genomic;
67 GenomicRegion * gr;
68 Protein * trans;
69
70 dpri = new_DPRunImpl_from_argv(&argc,argv);
71 if( dpri == NULL ) {
72 fatal("Unable to build DPRun implementation. Bad arguments");
73 }
74
75 gmp = new_GeneModelParam_from_argv(&argc,argv);
76
77 ct= read_CodonTable_file("codon.table");
78
79 strip_out_standard_options(&argc,argv,show_help,show_version);
80 if( argc != 2 ) {
81 show_help(stdout);
82 exit(12);
83 }
84
85
86 if((gm=GeneModel_from_GeneModelParam(gmp)) == NULL ) {
87 fatal("Could not build gene model");
88 }
89
90
91 seq = read_fasta_file_Sequence(argv[1]);
92
93 assert(seq);
94
95 cses = new_ComplexSequenceEvalSet_from_GeneModel(gm);
96
97 cseq = new_ComplexSequence(seq,cses);
98
99 rc = flat_RandomCodon(ct);
100 rmd = RandomModelDNA_std();
101
102 fold_in_RandomModelDNA_into_RandomCodon(rc,rmd);
103 rcs = RandomCodonScore_from_RandomCodon(rc);
104
105 exonscore = SyExonScore_flat_model(200,250,0.1,0.1);
106 /*
107 for(i=0;i<cseq->length;i++) {
108 fprintf(stdout,"%d PairSeq is %d score %d\n",i,CSEQ_PAIR_PAIRBASE(cseq,i),nonc_score->base[CSEQ_PAIR_PAIRBASE(cseq,i)]);
109 }
110 exit(0);
111 */
112 /*
113 show_RandomCodonScore(rcs,stdout);
114
115
116 for(i=3;i<seq->len;i++) {
117 fprintf(stdout,"seq %d is %c with score %d\n",i,aminoacid_from_seq(ct,seq->seq+i-2),rcs->codon[CSEQ_GENOMIC_CODON(cseq,i)]);
118 }
119
120 exit(0);
121 */
122
123 pal = PackAln_bestmemory_StatWise10(exonscore,cseq,rcs,Probability2Score(1.0/10.0),Probability2Score(1.0/10.0),NULL,dpri);
124 alb = convert_PackAln_to_AlnBlock_StatWise10(pal);
125
126 mapped_ascii_AlnBlock(alb,id,1,stdout);
127
128 genomic = Genomic_from_Sequence(seq);
129 gr = new_GenomicRegion(genomic);
130
131 add_Genes_to_GenomicRegion_GeneWise(gr,1,seq->len,alb,"bollocks",0,NULL);
132
133
134 for(i=0;i<gr->len;i++) {
135 if( gr->gene[i]->ispseudo == TRUE ) {
136 fprintf(stdout,"#Gene %d is a pseudo gene - no translation possible\n",i);
137 } else {
138 trans = get_Protein_from_Translation(gr->gene[i]->transcript[0]->translation[0],ct);
139 write_fasta_Sequence(trans->baseseq,stdout);
140 }
141 }
142
143
144
145 return 0;
146 }
147
148
149
150
151