Fast kmer counting table algorithm using perfect hash function: C++ pseudo-code integration into R using Rcpp API

Abstract

Counting kmers (substrings of length k in DNA sequence data) is an essential component of many methods in bioinformatics, including data preprocessing for de novo assembly, repeat detection, and sequencing coverage estimation. We proposed a simple algorithm to calculate the kmer count using perfect hash table implemented in C++ and using Rcpp API to be able exported into R. The pdf version is available at: Fast kmer counting table algorithm using perfect hash function: C++ pseudo-code integration into R using Rcpp API

References

1.            Deorowicz, S., A. Debudaj-Grabysz, and S. Grabowski, Disk-based k-mer counting on a PC. BMC bioinformatics, 2013. 14(1): p. 1.

2.            Melsted, P. and J.K. Pritchard, Efficient counting of k-mers in DNA sequences using a bloom filter. BMC bioinformatics, 2011. 12(1): p. 1.

3.            Zhang, Q., et al., These are not the k-mers you are looking for: efficient online k-mer counting using a probabilistic data structure. PloS one, 2014. 9(7): p. e101271.

4.            Pages, H., et al., String objects representing biological sequences, and matching algorithms. R package version, 2009. 2(2).

Rcpp: the Shell sort code with gaps function exemples

/*
 * The method starts by sorting pairs of elements far apart from each other, then progressively reducing the gap between elements to be compared.
 * The running time of Shellsort is heavily dependent on the gap sequence it uses.
 * Many gap function are evaluable see Link.
 * source file:Rcpp: the Shell sort code with gaps function exemples

 * 
 * The best profile was effectuated when we use sort_shell_sed function.
 * 
 *                                 MACHERKI M E.
 *                                 17/02/2016
 * 
 * 
 */

#include <Rcpp.h>
using namespace Rcpp;



//The following code is the typical body of the shell sort algorithm:
inline void shell ( NumericVector unsorted, int size, NumericVector Gaps)
{
  int gn=Gaps.size();
  int j,gap ;
  
  for (int t = gn-1; t > -1; t--)
  {
    gap=Gaps[t];
    for (int i = gap; i < size; ++i)
    {
      double temp = unsorted[i];
      for (j = i; j >= gap && temp < unsorted[j - gap]; j -= gap)
      {
        unsorted[j] = unsorted[j - gap];
      }
      unsorted[j] = temp;
    } 
  }
}
//******************gaps functions*************************

//Many gaps function can be used 
/*Shell sort function using slicing by two gaps (Shell, 1959)
 Gaps= N/2^k={N/2,N/4,N/8...1}
*/
inline NumericVector gapshell( int size ) {
  NumericVector result;
  for(int i=size/2;i>0;i/=2)
    result.push_back(i);
  return rev(result);}
/*Method of Frank & Lazarus, 1960
 Gaps=2*(N/2^k+1)+1
*/
inline NumericVector gapFrank( int N ) {
  int tmp=10;
  NumericVector result;
  for(int k=2;tmp>1;k++){
    tmp=2*(N/pow(2,k))+1;
    result.push_back(tmp);}
  return rev(result);}

/*Method of Hibbard, 1963
 Gaps=2^k-1
 */
inline NumericVector gapHibbard( int N ) {
  NumericVector result;
  int tmp=0;
  for(int k=1;tmp<N;k++){
    tmp=pow(2,k)-1;
    result.push_back(tmp);}
  return result;}

/*Method of Pratt, 1971
 Gaps=(3^k-1)/2
 */
inline NumericVector gappratt( int N ) {
  int n=N/3;
  NumericVector result;
  int tmp=0;
  for(int k=1;tmp<n;k++){
    tmp=(pow(3,k)-1)/2;
    result.push_back(tmp);}
  return result;}
/*Method of Sedgewick, 1986
 Gaps=4^k+3*2^(k-1)+1
 */
inline NumericVector gapsed( int N ) {
  NumericVector result;
  result.push_back(1);
    int tmp=0;
  for(int k=1;;k++){
    tmp=pow(4,k)+(3*pow(2,k-1))+1;
    if(tmp>N)break;
    result.push_back(tmp);}
  return result;}
/*Method of Tokuda, 1992
 Gaps=(9^k-4^k)/(5*4^(k-1))
*/

inline NumericVector gaptocuda( int N ) {
  NumericVector result;
  int tmp=0;
  for(int k=1;;k++){
    tmp=(pow(9,k)-pow(4,k))/(5*pow(4,k-1));
    if(tmp>N)break;
    result.push_back(tmp);}//we can use tmp+1
  result[0]=1;
  return result;}
/*Method of Tokuda, 1992:generic
 Gaps=h[i]=(h[i-1])*2.25+1; h[1]=1,1 base indexing
*/

inline NumericVector gaptocuda_emp( int N ) {
  NumericVector result;
  result.push_back(1);
  int tmp=1;
  for(int k=1;;k++){
    tmp=2.25*tmp+1;
    if(tmp>N)break;
    result.push_back(tmp);}//we can use tmp+1
  return result;}





//using empirically derived gaps (Ciura, 2001) N<2000

NumericVector gapCiura( int sz ) { 
  NumericVector result;
  result.push_back(1);
  result.push_back(4);
  result.push_back(10);
  result.push_back(23);
  result.push_back(57);
  result.push_back(132);
  result.push_back(301);
  result.push_back(701);
  result.push_back(1705);
  return result;}

/******************sorting functions*************************
************************************************************/







// [[Rcpp::export]]
NumericVector sort_shell_gap(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gapshell(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_frunk(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gapFrank(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_hibbard(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gapHibbard(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_patt(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gappratt(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_sed(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gapsed(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_tacuda(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gaptocuda(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_tacudaemp(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector gaps=gaptocuda_emp(N);
  NumericVector result=clone(unsorted);
  shell(result,N,gaps);
  return result;}
// [[Rcpp::export]]
NumericVector sort_shell_Ciura(NumericVector unsorted) { //run
  int N =unsorted.size();
  NumericVector Ciura=gapCiura(0);
  NumericVector result=clone(unsorted);
  shell(result,N,Ciura);
  return result;}
/***R
##trial mode
x<-runif(1000)
base<-sort(x)
a<-sort_shell_gap(x)
b<-sort_shell_frunk(x)
c<-sort_shell_hibbard(x)
d<-sort_shell_patt(x)
e<-sort_shell_sed(x)
f<-sort_shell_tacuda(x)
g<-sort_shell_tacudaemp(x)
h<-sort_shell_Ciura(x)
all.equal.list(base,a,b,c,e,f,g,h)
## timing :
x<-runif(1000000)
system.time(base<-sort(x))
system.time(a<-sort_shell_gap(x))
system.time(b<-sort_shell_frunk(x))
system.time(c<-sort_shell_hibbard(x))
system.time(d<-sort_shell_patt(x))
system.time(e<-sort_shell_sed(x))
system.time(f<-sort_shell_tacuda(x))
system.time(g<-sort_shell_tacudaemp(x))
system.time(h<-sort_shell_Ciura(x))

*/
  

Rcpp: reading file functions codes

/*
All these function are useful in Reading file.read_file_delime function read the file from start to end.'Tol' is an option to make string into lower case.
                              MACHERKI M E.
                              16/02/2016


Download source file:link

*/


#include <fstream>
#include <sstream>
#include <string>
#include <algorithm>
#include<cctype>
#include<Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
CharacterVector read_file(std::string path){
std::ifstream t(path.c_str());//connect with file 
std::stringstream ss;
ss<<t.rdbuf();                // scan file
return ss.str();
}
//[[Rcpp::export]]
CharacterVector read_file1(std::string path){
std::ifstream in(path.c_str());
std::string contents;
in.seekg(0,std::ios::end);
contents.resize(in.tellg());
in.seekg(0,std::ios::beg);
in.read(&contents[0],contents.size());
in.close();
return contents;
}
// function to read file from start to end pointer
//[[Rcpp::export]]
CharacterVector read_file_delime(std::string path,int start,int end,int TOL){
std::ifstream in(path.c_str());
std::string contents;
in.seekg(0,std::ios::end);
contents.resize(in.tellg());
in.seekg(start,std::ios::beg);
in.read(&contents[0],end);
in.close();
contents.erase(std::remove(contents.begin(),contents.end(),'\n'),contents.end());
//remove lines delimiter
if(TOL>0){//force to lower case transformation
std::transform(contents.begin(),contents.end(),contents.begin(),tolower);
}
return contents.substr(start-1,end-start);
}
/* read a file from start  to end*/
/* set the file path before using */

/*** R
###*******************choose file to read************************####
con<-utils::choose.files()
system.time(FR<-readLines(con))     #vector as the split(x,'\n') 
system.time(Fcpp1<-read_file(con))  #contents Line delimiter
system.time(Fcpp2<-read_file1(con)) #contents Line delimiter
#one string without \n and in lower case
system.time(Fcpp3<-read_file_delime(con,20,100,1)) 

*/ 

Merge sort: Rcpp integration code

/*sorting function using the merge sort algorithm
*/

#include <Rcpp.h>
using namespace Rcpp;

void intercal(int p, int q, int r, NumericVector v,  NumericVector w)
{
  int i, j, k;
  i = p;
  j = q;
  k = 0;
  while (i < q && j < r) {
    if (v[i] < v[j]) {
      w[k] = v[i];
      i++;
    }
    else {
      w[k] = v[j];
      j++;
    }
    k++;
  }
  while (i < q) {
    w[k] = v[i];
    i++;
    k++;
  }
  while (j < r) {
    w[k] = v[j];
    j++;
    k++;
  }
  for (i = p; i < r; i++)
    v[i] = w[i-p];
}

void mergesort(int p, int r, NumericVector v, NumericVector aux)
{
  int q;
  if (p < r - 1) {
    q = (p + r) / 2;
    mergesort(p, q, v,aux);
    mergesort(q, r, v,aux);
    intercal(p, q, r, v,aux);
  }
}

// [[Rcpp::export]]
NumericVector sort_merge(NumericVector vetor) {
  Rcpp::NumericVector res = Rcpp::clone(vetor);
  Rcpp::NumericVector aux = Rcpp::clone(vetor);
  int n = res.size();
  mergesort(0,n,res,aux);
  return res;}

Analyse d'une séquence de protéine //Macherki M E

La question qui m’intéresse et intéresse beaucoup des personnes concernant une séquence codante pour une protéine : quels sont les manipulations à faire à cette séquence ?

La première étape consiste à un alignement dans la base de donné (BLAST, protSWISS…) pour savoir quel sont les organismes qui l’utilisent. Une information générale est élaborée à partir des publications sur la structure, le rôle et les propriétés physico-chimiques (PUBMED).

Supposant qu’on va essayer avec notre ordinateur pour déterminer les caractéristiques d’une séquence aléatoire. Pour les propriétés physico chimique, je propose le package seqinr pour le logiciel R.

Par exemple pour la séquence de référence pour le protéine 16S R d'Aggregatibacter sp. Clone_MB3_C38 :

Le fichier FASTA est définit par :

attr(,"name")

[1] "513_3635"

attr(,"Annot")

[1] ">513_3635 | Aggregatibacter sp. |  HOT_513 | Clone_MB3_C38 | DQ003635 | Phylotype"

attr(,"class")

[1] "SeqFastadna"

>   AAstat(translate(a))###C'est une fonction qui donne des statistiques utiles pour notre protéine

$Compo

 *  A  C  D  E  F  G  H  I  K  L  M  N  P  Q  R  S  T  V  W  Y

31 32  8 18 22  6 48 12 14 11 41 10 20 27 19 34 36 21 37  9  8

$Prop

$Prop$Tiny

[1] 0.21261

$Prop$Small

[1] 0.3621701

$Prop$Aliphatic

[1] 0.1348974

$Prop$Aromatic

[1] 0.05131965

$Prop$Non.polar

[1] 0.3519062

$Prop$Polar

[1] 0.2829912

$Prop$Charged

[1] 0.1422287

$Prop$Basic

[1] 0.08357771

$Prop$Acidic

[1] 0.05865103

$Pi

[1] 8.774443

Pour la visualisation 3D, je propose le logiciel SPDBV_4.10_PC. Il est très simple à manipuler avec des séquence d’ADN.

SWISSMODEL-- séquence de nucléotide---Tools-----> surface

Il existe plusieurs autres manipulations à réaliser avec ce deux logiciels, il faut juste les essayer.

Application de filtre linaire dans l’analyse de génome//Macherki M E

L’application de filtre linaire pour de série temporelle est très applicable. Dans le cas des acides nucléiques, il est très simple d’appliqué un tel genre d’analyse en utilisant un logicielle simple d’analyse tel que R. Le vecteur contenant les positions d’un oligo dans une séquence est la base de l’analyse en appliquant la fonction diff qui permettre de déterminer alors la distance entre deux oligo consécutif soit disant un dérivative de premier dégrée. En utilisant la fonction lm entre la somme cumulatif (position) de dérivatif et le vecteur normal, nous pouvons exprimer de  façon globale la fréquence de l’oligo étudié. Le graphique de résiduelle permet de schématiser la structure de chromosome avec détaille (origine et  terminus pour une bactérie par exemple E.coli K12 en utilisant l'oligo ‘GGG’  fig).

Il est simple de déterminer la valeur originale sur la séquence par juste une étape de retour en arrière. L’intérêt de filtre linéaire est plus clair si en va effectuer un alignement des séquences. Notant que la fréquence d’un oligo de la séquence à alignée est bien connu, en appliquant un retour en arrière    selon le coefficient de corrélation de séquence cible, nous pouvons  déterminer un intervalle ou la séquence sera exister. En utilisant plusieurs oligos (exemple pour n=5), nous somme capable de  prédire l’emplacement avec une exactitude parfaite.
Macherki M E:R course and applications in genomic 2014

La descrimination de la longeur de CDS //Macherki M E

L’analyse de la phylogénie nous permettre de déterminer les relations entre les groupes d’individus qui sont discriminatoire et  différents. Dans ce contexte, il est commun d’utiliser l’ADN dans la comparaison étant qu’il représente le support de l’information génétique. Avec tout les critiques qu’il entoure, la taille de génome était toujours une grandeur de base dans cette analyse. Le pourcentage en GC reste le paramètre dépendant puis qu’il est stable au nivaux de génome. Autrement, nous avons essayé d’utiliser l’ADNc dans la comparaison et non l’ADN génomique pour bien cibler la partie transcrite et non hérédité. Nous avons utilisé une méthode de permutation de Monte Carlo (n=10000) pour déterminer le moyenne  de Log de longueur exprimé en paire de base (table I)

Macherki M E:R course and applications in genomic 2014

Tentative en biologie computationnelle//Macherki M E

L’étude  de l’ADN permet de révéler plusieurs informations permettant l’identification des différents compartiments dans le  génome et entre des génomes distincts. L’information la plus commun est la fréquence d’un oligo au sein d’une séquence. Par contre cet information est non explicatif .C’est juste un indice calculé qui vari énormément au sein de la séquence. En utilisant la méthode de "DNA walk", cette fréquence suit par conséquence une loi normale puisque cette fréquence change d’un emplacement à un autre dans le génome. Nous étudions la distance entre les oligo présent dans la séquence. Nous considérons alors que les bases symbolisent un enchaînement de pièces de dominos mis l’une après l’autre. Nous avons déterminé la loi de probabilité pour ce variable discret notamment géométrique  et une corrélation entre la fréquence et la distance. Pour des buts comparatifs, nous avons essayé de stabiliser nos statistique en utilisant la méthode de karlin S. et al.et comparer les séquences des génomique après une normalisation à des prédit issu de l’approximation géométrique.

En fin, nous avons essayé de comparer  la structure eucaryote (Chr X de l’homme) et celle de procaryote (E.coli k12)en utilisant notre approximation d’indépendance telle que l’indice de rho avec deux oligos à la fois . Nous avons signalé un effet d'usage des  codons  chez E.coli (aire segmentée dans la FIG) avec des zone sur exprimées et sous exprimées consécutives tel qu’il signaler ultérieurement avec karlin et al.

Dans le chromosome X, l’effet d'usage des  codons est omis et nous reportons que les effets d’abondance( aire bleu dans la FIG indique la répression GC). Cette remarque permet de nous conclure que notre tentative peut permettre de distinguer entre une zone transcrit et non transcrit se qui permet de distinguer l’ADN eucaryote de celle procaryote facilement.

Macherki M E:R course and applications in genomic 2014