The first code:Rcpp

sourceCpp is the important function in Rcpp. The Rcpp code is written into file with'.cpp' extension. The cppFunction can be useful but we can not  easily combine many functions in the same time. This is an example to use the two possibility:

• ccpFunction

R code:

Rcpp::cppFunction(
'int fibonacci(const int x) {
if (x == 0) return(0);
if (x == 1) return(1);
return (fibonacci(x - 1)) + fibonacci(x - 2);
}')

fibonacci(20)

• sourceCpp

I used the example used in the Rstudio C++ source file.

Get the source file: firstrun.cpp and save it in the same R directory.

R code:

Rcpp::sourceCpp('firstrun.cpp')

Rcpp:the first step

Rcpp uses the C++ compiler to produce a sharing library connectable with R environment. The fist step is to install  the Rtools program provided in CRAN site. The second step is to install the Rcpp package.It is important to use a good  editor such Rstudio which integrates the C++ code directly.

Why I use Rcpp ?

Rcpp is an R add-on package which facilitates extending R with C++ functions. It is being used for anything from small and quickly constructed add-on functions written either to fluidly experiment with something new or to accelerate computing by replacing an R function with its C++ equivalent to large-scale bindings for existing libraries, or as a building block in entirely new research computing environments. While still relatively new as a project, Rcpp has already become widely deployed among users and developers in the R community. Rcpp is now the most popular language extension for the R system and used by over 500 CRAN packages as well as ten BioConductor packages. This site aims to provide a solid introduction to Rcpp and provides many codes of the most common algorithms.

Numerical coding of DNA: Rcpp code

/*
In general, DNA sequences are  stored as character data into different formats
of files like FASTA,FNA... ALso numeric coding of DNA in base 4 exists  such 
'A'=0,'G'=1,'C'=2 and 'T'=3. But why I use this kind of coding? 
The first step of sequence alignment algorithms such BLAST or FASTA consists on 
stepwise splitting of the DNA sequence into oligos of 11 bases. This operation 
influences  the execution  running time of the program.
This is an example of codes  ;the first use string manipulation and the second 
use numeric coding in base 5(to overcame  N in the sequence).
PDF file link:source file


                     MACHERKI M.E
                     14/02/2016
 */
 
 /*
 This source file  uses C++ code.You have to install Rcpp package to run it.
 A FASTA file  path is required to test this code in R.Please check the bottom 
 of the page.
 */ 
 
#include <fstream>
#include<Rcpp.h>
#include<string.h>
#include <sstream>
#include<cctype>
#include <algorithm>
using namespace Rcpp;

/*
First of all, this is a function to read FASTA file into R
-----!!!!!!!!!!!!! one sequence only!!!!!!!!!!!!!!-----------
*/
//[[Rcpp::export]]
inline CharacterVector read_fasta(std::string file ){
CharacterVector records;
std::ifstream in(file.c_str());
in.get();// remove first >
std::string rec ;
while(getline(in,rec,'>')){
int newLineLoc=rec.find('\n');
std::string header =rec.substr(0,newLineLoc);
std::string sequence= rec.substr(newLineLoc+1,rec.length()-newLineLoc-2);
sequence.erase(std::remove(sequence.begin(),sequence.end(),'\n'),sequence.end());
std::transform(sequence.begin(),sequence.end(),sequence.begin(),tolower);
records[header]=sequence;
}
return records;
}
/*
String splitter function into N size oligos (Only for lower case string).

*/
// [[Rcpp::export]]
 inline CharacterVector SPLIT(const std::string  strings,const int n){
    int lim=strings.length()-n+1 ;//size of sub string vector
    CharacterVector resultat(lim) ;
          for ( int j=0; j<lim ;j++) resultat[j] = strings.substr(j,n) ;
return resultat ;
} 
/*
String splitter function into N size oligos using numeric coding.
I used a simple hashing function which produce an unique numeric 
identifier for each oligo.It is possible to use integer but I use
double to accept a large N.
(See R help:strtoi or just type ?strtoi in console)

*/

inline double  natoi( const std::string A,const int n ){
//Used as input of numeric_split(see next function)
  int np=n-1;
  IntegerVector temp(n,4);
  char toc;
  for( int i=0;i<n;i++){
    toc=A[i];
if(toc=='a') temp[i]= 0;
if(toc=='t') temp[i]= 3;
if(toc=='c') temp[i]= 1;
if(toc=='g') temp[i]= 2;
}
double  result=0;
 for( int i=0;i<n;i++) result+= pow(5,(np-i))*temp[i];
 return result;
}

// [[Rcpp::export]]
inline NumericVector  numeric_split(const std::string A,const int n)
    {
  const int np=n-1;
  const int Lim=A.length()-np;//-- number of subsequence 
  NumericVector res(Lim);     //-- result vector
  double p=0;                 //-- temp result accumulator
  res[0]=p=natoi(A,n);        //-- first result 
  const double  N=pow(5,np);  //-- cont used in recursion formula
  unsigned short int x,y;     //-- temp to allocate forward and backward values
  char fo,ba;                 //--forward and backward chars
  for( int i=1;i<Lim;i++)
     {
       fo=A[i-1];ba=A[i+np];
       x=y=4;                 //-- empty case:n,p,m... 
       if(fo=='a')x=0;        //-- char detector 
       if(fo=='t')x=3;        // 
       if(fo=='g')x=2;        // 
       if(fo=='c')x=1;
       if(ba=='a')y=0;
       if(ba=='t')y=3;
       if(ba=='g')y=2;
       if(ba=='c')y=1;
     res[i]=p= 5*p-5*N*x+y;   //-- scalar formula 
     }
  return res;
  }
/*
R verification :
I use a file named k12.fna in my computer (full E.coli k12 genome ). 
You have just to change k12.fna into your file name.
*/


/***R
x<-read_fasta("k12.fna")  ##change me to your file name 
system.time(charsplit<-SPLIT(x,11))
#user  system elapsed 
# 5.38    0.02    5.39 
system.time(numsplit<-numeric_split(x,11))
#user  system elapsed 
#0.08    0.01    0.09    ## more then 50 time faster
  
  */

Source file example for some sorting algorithms written in C++ via RCPP to R implimentation

#include <Rcpp.h>
# include<algorithm>
using namespace Rcpp;


/*
 * This file contains some algorithms of sorting written in C++ & implimented in RCPP .
 * All these algorithms are not performed as well as sort inner algorithm 
 * used in R except of the sort_stl function .
 *   

 * get pdf: source      
 *                 MACHERKI M.E 
 *                  12/02/2016 
 */


// [[Rcpp::export]]
NumericVector sort_insert(NumericVector A) {
  int i,j;
  double x;
  for(i=1; i<A.size();i++){
    x=A[i];j=i-1;
    while((j>=0) & (A[j]>x)){ //without using swapping
      A[j+1]=A[j];
      j=j-1;
    }
  A[j+1]=x;
  
  }
  return A;
}
 
 // [[Rcpp::export]]
 NumericVector sort_bubble(NumericVector A) { 
   int n=A.size();
   int swap=1;
   double tmp;
   for(;;){
     if(swap==0)break;
     swap=0;
     for(int i=1;i<n;i++){
       if(A[i-1]>A[i]){
         tmp=A[i];A[i]=A[i-1];A[i-1]=tmp;
         swap=1;
       }
     }
   }
   return A;}

// [[Rcpp::export]]
NumericVector sort_selection(NumericVector A) { 
  int i,j;
  double tmp;
  for(j=0;j<A.size();j++){
    int imin=j;
    for(i=j+1;i<A.size();i++){
      if(A[i]<A[imin])imin=i;
    }
  if(imin!=j){
   tmp=A[j];A[j]= A[imin];A[imin]=tmp;
  }
  
  }
  return A;}
  
 

   // [[Rcpp::export]]
 NumericVector gaps( int sz ) { 
   NumericVector result;
   result.push_back(701);
   result.push_back(301);
   result.push_back(132);
   result.push_back(57);
   result.push_back(23);
   result.push_back(10);
   result.push_back(4);
   result.push_back(1);
   return result;}

// [[Rcpp::export]]
NumericVector sort_shell(NumericVector A) { 
  double temp=0;
  int N =A.size();
  int i,j,gap;
  NumericVector tacuda=gaps(0);
  int n=tacuda.size();
  for(int k=n-1;k>0;k--){
    gap=tacuda[k];
  for( i=gap;i<N;i++){
    temp=A[i];
    for( j=i;(j>=gap)&(A[j-gap]>temp);j-=gap){
      A[j]=A[j-gap];
    }
    A[j]=temp;
    
    }
    
  }
  return A;}


  
  
// This function is faster then R sort function !
// [[Rcpp::export]]
NumericVector sort_stl(NumericVector A) { 
  
  NumericVector tmp=clone(A);
  std::sort(tmp.begin(),tmp.end());
  return tmp;
}
  
  
  // [[Rcpp::export]]
 void Quick(NumericVector arr,int left,int right) { 
  int i=left;int j=right;
  double tmp;
  double pivot=arr[(left+right)/2];
  /*partition*/
  while(i<=j){
    while(arr[i]<pivot)
      i++;
    while(arr[j]>pivot)
      j--;
      if(i<=j){
        tmp=arr[i];
        arr[i]=arr[j];
        arr[j]=tmp;
        i++;
        j--;
        
      }
    }
  /*recursion*/
  if(left<j)
    Quick(arr,left,j);
  if(i<right)
    Quick(arr,i,right);
}
    
    
    
    /*
     * This function is really comparable with the inner sort used in R
     * even the code is sample , the algorithm is good.
     */
    
    
    // [[Rcpp::export]]
    NumericVector sort_quick(NumericVector A) {
      NumericVector tmp=clone(A);
    Quick(tmp,0,tmp.size()-1); //just starting from the middle of the vector
      return tmp;
    }
      
  
  
  
  


// R call

/*** R

###not rapid methods, adapted to small sample size
x<-rnorm(100)
sort_insert(x)
sort (x)
sort_selection(x)
sort_bubble(x)
identical(sort_stl(x),sort (x))## the purpose of C++ integration  
identical(sort_quick(x),sort (x))
identical(sort_insert(x),sort (x))
identical(sort_bubble(x),sort (x))
identical(sort_selection(x),sort (x))
identical(sort_shell(x),sort (x))
x<-rnorm(1000000)    ##large numeric vector
system.time(sort(x))
system.time(sort_quick(x)) ## as R sort function
system.time(sort_stl(x)) ## two time faster the R sort!!!    


  */


   

Question

see link

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.