Matías Neiff's Blog

Con estos pequeños archivos se puede hacer un analisis de kmer rápido y simple en linux. Para compilar solamente ejecutar:
gcc -o kmer kmer.c

Luego para analizar la ocurrencia de 41mer ejecutar:
cat fastQfile.fq | parser | ./kmer 41 > salida.txt

La salida será algo así:

AAATAATGTACGGGGGAGATGCATGAATTTTTCGGATAAAA 190
AATAATGTACGGGGGAGATGCATGAATTTTTCGGATAAAAT 189
ATAATGTACGGGGGAGATGCATGAATTTTTCGGATAAAATC 190
TAATGTACGGGGGAGATGCATGAATTTTTCGGATAAAATCT 193
AATGTACGGGGGAGATGCATGAATTTTTCGGATAAAATCTG 193
ATGTACGGGGGAGATGCATGAATTTTTCGGATAAAATCTGA 193

En la primer columna está el mer y en la segunda la frecuencia. Los archivos son:

–parser–

#!/bin/bash
##parser from fastQ to kmer
if [ $# -eq 0 ]; then
    echo "usage: parser file.fq > file.txt";
    exit 1;
fi
#total of lines
lines=$(cat $1 | wc -l);
for i in $(seq 2 4 $lines); do
    sed $i'q;d' $1
done
exit 0;

–kmer.c–

/**
 * Author:    Matias Neiff and Tatiana Ponce
 * Created:   21.08.2020
 *
 * Software to evaluate the occurrence of kmers in a dna sequence. Use it good, use it for good.
 * otput: kmer [tab] frecuency
 * 
 * This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
 * This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
 * You should have received a copy of the GNU General Public License along with this program. If not, see <https://www.gnu.org/licenses/>.
 **/

#include <stdio.h>
#include <stdlib.h>
#include <errno.h>
#include <string.h>
#include <unistd.h>
// this value has to be sufficient big to host the entire sequence
#define MAX_SEQ_LEN 24845291000
// this is the read buffer
#define MAX_READ_LEN 100 * 1024
# 24845291

//compare to array string and return true if are identical
int compareArrays(char a[], char b[], int n)
{
  for (int i=0; i<n; ++i) if (a[i] != b[i]) return 0;
  return 1;
}

int usage() {
        printf("Usage: cat file | kmer k\n");
	exit(1);
}

int main(int argc, char **argv) {
    int k;
    char* p;
    errno = 0;
    long conv;
    // at least need two arguments
    if (argc == 2) {
	conv = strtol(argv[1], &p, 10);
    } else {
	usage();
    }

    //check if K is integer
    if (errno != 0 || *p != '\0' || conv > 9999 || conv < 1) {
	usage();
    } else {
	k = conv;
    }
    // allocate memory for the three big strings
    char* seq;
    char* buff;
    char* prevkmer;
    seq = (char*) malloc(MAX_SEQ_LEN*sizeof(char));
    buff = (char*) malloc(MAX_READ_LEN*sizeof(char));
    prevkmer = (char*) malloc(MAX_SEQ_LEN*sizeof(char));
    if (seq == NULL) {
        printf("Memory not allocated.\n");
        exit(0);
    }
    // we need to store the bytes readed trouhgt the pipe, and the total of bytes readed
    ssize_t bytesRead=0;
    ssize_t totalRead=0;
    bytesRead = read( STDIN_FILENO, buff, MAX_READ_LEN );
    while (bytesRead != 0) {
	buff[bytesRead]='\0';
	strcat(seq,buff);
	totalRead = totalRead + bytesRead;
	bytesRead = read( STDIN_FILENO, buff, MAX_READ_LEN );
    }
    // add the null to the end to complete the string format of C
    seq[totalRead]='\0';
    char kmer[k];
    char window[k];
    int seqlen = totalRead;
    int freq;
    prevkmer[0]='\0';
    //start the analysis
    for(int index=0;index<seqlen-k;index++) {
	// select kmer
	for (int i=0;i<k;i++) {
	    kmer[i]=seq[i+index];
	}
       // prevkmer store already proccessed kmer to avoid procces twice
	kmer[k]='\0';
	// do not procces if was already proccess with this kmer or if have an return car (incomplete mer)
	if (!strstr(prevkmer,kmer) && !strstr(kmer,"\n")) {
	// store as previous kmer to not proccess again
    	strcat(prevkmer,kmer);
	freq=0;
	// search frecuency of this kmer
	for (int x=0;x<seqlen-k;x++) {
	    //fill the window
	    for (int w=0;w<k;w++) {
		window[w]=seq[x+w];
	    }
	    // if kmer equal to window
	    if(compareArrays(kmer,window,k)) { freq++;x=x+k;}	    
	}
	if (freq>1) printf("%s\t%d\n",kmer,freq);
	}
    } //end for
    free(seq);
    free(buff);
  return 0;
} // end function: main