Script To Calculate Word Frequency For Many Sequences

[thom_otis]

Hello!

Who can help me with a writing a script?

1. For each file in the set (archive, probably):

Open input nucleotide sequence (.embl format is desirable, but fasta format is possible);

Calculate the observable frequency for all possible words based on input word size (e.g. number of all possible words for word size = 3 (triplet) is 64). Obs. Frequency = Obs Count/Total count [Сompseq is work in a similar way, but not support a batch processing].

Shift the reading frame by 1 nucleotide and repeat the previous step (number of shifts = word size - 1). The file with frequencies for each file in the set save as name_of_the_file(id of a contig).dic(+word size). It looks sth like this (word size = 2):



2. Make a summary file, something like this (example for word size = 2)

[Richard Finney]

example usage: cat hg20/chr1.fa | ./kmer3

-bash-4.1$ cat kmer3.c

/*
vi kmer3.c ; gcc -Wall -O3  kmer3.c -o kmer3  -lm 
*/

#include <stdio.h>
#include <string.h>
#include <stdio.h>
#include <ctype.h>
#include <errno.h>
#include <assert.h>

long int x3[4][4][4];

int acgt_index(char c)
{
    if (c=='A') return 0;
    else if (c=='C') return 1;
    else if (c=='G') return 2;
    else if (c=='T') return 3;
    return -1;
}


void wrapup3(void)
{
   int i,j,k; 

   for (i=0;i<4;i++)
   for (j=0;j<4;j++)
   for (k=0;k<4;k++)
   {
       if (i==0) {printf("A");} else if (i==1) {printf("C");} else if (i==2) {printf("G");}  else if (i==3) {printf("T"); }
       if (j==0) {printf("A");} else if (j==1) {printf("C");} else if (j==2) {printf("G");} else if (j==3) {printf("T"); }
       if (k==0) {printf("A");} else if (k==1) {printf("C");} else if (k==2) {printf("G");} else if (k==3) {printf("T"); }
       printf(" %ld\n",x3[i][j][k]);
   }
}

int is[4];

static int kmer(void)
{
   int on_head = 0;
   int c;
   long int cnt = (long int)0;

    while (1)
    {
         c = getchar();
         cnt++;
         if (c==-1) break;
         if (c=='>') {on_head=1; continue; }
         if ((c=='\r')|| (c=='\n')) {on_head=0; continue; }
         if (on_head == 1) continue;
         is[0] = is[1];
         is[1] = is[2];
         is[2] = acgt_index(toupper(c));
         if (is[0] == -1) continue;
         if (is[1] == -1) continue;
         if (is[2] == -1) continue;
         x3[is[0]][is[1]][is[2]]++;
   }
   return 0;
}

int main(int argc, char *argv[])
{
    kmer();
    wrapup3();
    return 0;
}

[thom_otis]

How to run this script?

example usage: cat hg20/chr1.fa | ./kmer3

-bash-4.1$ cat kmer3.c

/*
vi kmer3.c ; gcc -Wall -O3  kmer3.c -o kmer3  -lm 
*/

#include <stdio.h>
#include <string.h>
#include <stdio.h>
#include <ctype.h>
#include <errno.h>
#include <assert.h>

long int x3[4][4][4];

int acgt_index(char c)
{
    if (c=='A') return 0;
    else if (c=='C') return 1;
    else if (c=='G') return 2;
    else if (c=='T') return 3;
    return -1;
}


void wrapup3(void)
{
   int i,j,k; 

   for (i=0;i<4;i++)
   for (j=0;j<4;j++)
   for (k=0;k<4;k++)
   {
       if (i==0) {printf("A");} else if (i==1) {printf("C");} else if (i==2) {printf("G");}  else if (i==3) {printf("T"); }
       if (j==0) {printf("A");} else if (j==1) {printf("C");} else if (j==2) {printf("G");} else if (j==3) {printf("T"); }
       if (k==0) {printf("A");} else if (k==1) {printf("C");} else if (k==2) {printf("G");} else if (k==3) {printf("T"); }
       printf(" %ld\n",x3[i][j][k]);
   }
}

int is[4];

static int kmer(void)
{
   int on_head = 0;
   int c;
   long int cnt = (long int)0;

    while (1)
    {
         c = getchar();
         cnt++;
         if (c==-1) break;
         if (c=='>') {on_head=1; continue; }
         if ((c=='\r')|| (c=='\n')) {on_head=0; continue; }
         if (on_head == 1) continue;
         is[0] = is[1];
         is[1] = is[2];
         is[2] = acgt_index(toupper(c));
         if (is[0] == -1) continue;
         if (is[1] == -1) continue;
         if (is[2] == -1) continue;
         x3[is[0]][is[1]][is[2]]++;
   }
   return 0;
}

int main(int argc, char *argv[])
{
    kmer();
    wrapup3();
    return 0;
}

[Richard Finney]

Create source code using text editor, past in code.
Name the file kmer3.c

Compile with
gcc -Wall -O3 kmer3.c -o kmer3 -lm

Program takes standard input.
Example ...
cat hg20/chr1.fa | ./kmer3

fastas only unless you hack it.