adding matlab code

1 files changed, 52 insertions, 0 deletions

diff --git a/src/matlab/quikr.m b/src/matlab/quikr.m
new file mode 100644
index 0000000..d24596e
--- /dev/null
+++ b/src/matlab/quikr.m
@@ -0,0 +1,52 @@
+function xstar = quikr(inputfasta)

+%xstar=quikr('path/to/input/fasta/file') returns the estimated frequencies

+%of bacteria present when given an input FASTA file of amplicon (454)

+%reads. A k-mer based, L1-regularized, sparsity promoting algorithm is

+%utilized. This is the default implementation of Quikr, here lambda=10,000, k-mer

+%size is 6, and we use the trainset7_112011 database from RDP to estimate

+%the present species. In practice, reconstruction is accurate only down to

+%the genus level (not species or strain).

+%The program ./count-linux or ./count-osx needs to be placed in the same 

+%directory as this script. In Linux, you might need to: chmod a+rx count.

+%Also the inputfasta needs to be a standard FASTA file. The input FASTA file must have reads all in the same

+%orientation, so as a preprocessing step, make sure everything is in the

+%forward (+) orientation.

+if nargin>1

+    error('too many input arguments');

+end

+if (isunix && not(ismac))   %If linux, use ./count-linux

+    [status, counts]=unix(['./count-linux -r 6 -1 -u ' inputfasta]); %count the 6-mers in the fasta file, in the forward direction, return the counts without labels

+    if status ~= 0

+        error('./count-linux failed: ensure count-linux is an executable. Try chmod a+rx count-linux. Be sure matlab/octave is in the same directory as count-linux');

+    end

+elseif ismac %If mac, then use ./count-osx

+    [status, counts]=unix(['./count-osx -r 6 -1 -u ' inputfasta]); %count the 6-mers in the fasta file, in the forward direction, return the counts without labels

+    if status ~= 0

+        error('./count-osx failed: ensure count-linux is an executable. Try chmod a+rx count-osx. Be sure matlab/octave is in the same directory as count-osx');

+    end

+elseif ispc

+    error('Windows is not yet supported');

+end    

+counts=textscan(counts,'%f'); %convert into floats

+counts=counts{:}; %make into a vector

+counts=counts/sum(counts); %form a probability vector from the counts

+lambda=10000; %this is the default lambda value (see the paper), and is used in the formation of trainset7_112011N6Aaux.mat

+yaux=[0;lambda*counts]; % form the sample vector

+load('trainset7_112011N6Aaux.mat') %load the 6-mer sensing matrix...I need to rename this so the variable is just Aaux.

+if not(exist('Aaux2')) %make sure Aaux2 loaded

+    error('The default training database trainset7_112011N6Aaux.mat failed to load. Be sure it exists in the Quikr folder or current directory');

+end

+%based on RDP's trainset 7 from 11/2011

+isOctave = exist('OCTAVE_VERSION') ~= 0; %check to see if running Octave or Matlab

+if isOctave %if octave, need to manually read in sparse-matrix format

+    n=length(Aaux2.jc)-1;

+    ii=repelems(1:n,[1:n; diff(Aaux2.jc)]);

+    Aaux=sparse(Aaux2.ir+1,ii,Aaux2.data);

+else

+    Aaux=Aaux2; %fix the names

+    clear Aaux2;

+end %Matlab automatically reads it in correctly, just need to rename it

+warning off

+xstar=lsqnonneg(Aaux,yaux); %perform the nonnegative least squares

+warning on

+xstar=xstar/sum(xstar); %transform the output into a probability vector. Note this vector is on the same basis as Aaux, so the entries correspond to sequences in trainset7_112011.fa
\ No newline at end of file