How to convert Color Space to Fasta (biological sequence)

Quem já se deparou com esse problema deve saber a dificuldade em encontrar um programa que faça essa conversão. Bom, depois de testar um monte de coisas, andei estudando os fontes do "Color Space Mapping tool" da AB (http://www.appliedbiosystems.com/absite/us/en/home/support/software-community/free-ab-software.html) e descobri uma regra de conversão simples:
A0=>'A'
A1=>'C'
A2=>'G'
A3=>'T'
C0=>'C'
C1=>'A'
C2=>'T'
C3=>'G'
G0=>'G'
G1=>'T'
G2=>'A'
G3=>'C'
T0=>'T'
T1=>'G'
T2=>'C'
T3=>'A'

Com essa informação, fiz o programa abaixo (GNU C):

#include 
#include 
#include 

#define MAX_LINE 500

int main(int argc, char *argv[]) {
FILE *fin,*fout;
char buffer[MAX_LINE+1],seq[MAX_LINE];
char p0,p1,base,valid;
int i;

if( argc < 3 ) {
      printf("use: cs2fas  \n");
exit(0);
}

if( (fin = fopen(argv[1],"r")) == NULL ) {
printf("%s not found.",argv[1]);
exit(0);
}

if( (fout = fopen(argv[2], "w")) == NULL ) {
printf("%s can't be created.",argv[2]);
fclose(fin);
exit(0);
}

while( fgets(buffer, MAX_LINE, fin) != NULL ) {
if( buffer[0] == '#' )
continue;
if( buffer[0] == '>' ) {
         for(i=1 ; buffer[i] != '\0' ; i++) {
             if( buffer[i] == '\n' || buffer[i] == '\r')
                 buffer[i] = '\0';
         }
fprintf(fout,"%s\n",buffer);
} else {
// try convert cs 2 fas
memset(seq,'\0',MAX_LINE);
valid = 1;
if( buffer[0] >= 'a' && buffer[0] <= 'z' ) {
            buffer[0] = buffer[0] - 32;
         }
         for(i=1 ; buffer[i] != '\0' && valid ; i++) {
            p0 = buffer[i-1];
            p1 = buffer[i];
            base = 'X';
            if( p0 == 'A' ) {
               switch( p1 ) {
                   case '0': base = 'A'; break;
                   case '1': base = 'C'; break;
                   case '2': base = 'G'; break;
                   case '3': base = 'T'; break;
                   default:
                     valid = 0;
               }

            } else if( p0 == 'C' ) {
               switch( p1 ) {
                   case '0': base = 'C'; break;
                   case '1': base = 'A'; break;
                   case '2': base = 'T'; break;
                   case '3': base = 'G'; break;
                   default:
                     valid = 0;
               }

            } else if( p0 == 'G' ) {
               switch( p1 ) {
                   case '0': base = 'G'; break;
                   case '1': base = 'T'; break;
                   case '2': base = 'A'; break;
                   case '3': base = 'C'; break;
                   default:
                     valid = 0;
               }

            } else if( p0 == 'T' ) {
               switch( p1 ) {
                   case '0': base = 'T'; break;
                   case '1': base = 'G'; break;
                   case '2': base = 'C'; break;
                   case '3': base = 'A'; break;
                   default:
                     valid = 0;
               }

            } else {
              // invalid char
              valid = 0;
            }
            if( valid ) {
               buffer[i] = base;
               seq[i-1] = base;
            }
         }
         fprintf(fout,"%s\n",seq);
      }
   }
   fclose(fin);
   fclose(fout);
   return 0;
}


O programa ainda está em teste, assim que fizer mais analises eu coloco aqui.

Comentários

Postagens mais visitadas deste blog

Jellyfish script

Conversão do encode do MariaDB para atender o moodle 3.8

O GBParsy é uma biblioteca para realizar o parser de arquivos GenBank para o Python