adding matlab code

1 files changed, 52 insertions, 0 deletions

diff --git a/src/matlab/Example.m b/src/matlab/Example.m
new file mode 100644
index 0000000..d2a3192
--- /dev/null
+++ b/src/matlab/Example.m
@@ -0,0 +1,52 @@
+% Here is an example of re-training the quikr method, performing the reconstruction, and outputing the results to a CSV file appropriate for Excel.

+mat=quikrTrain('TaxCollectortrainset7_112011.fa',6); %retrain using the TaxCollector version of trainset7_112011.fa

+xstar=quikrCustomTrained(mat,'testfastafile',6,10000); %perform the reconstruction

+quikrWriteOutput('/path/to/output/director','filenameforoutput',full(xstar),'/path/to/TaxCollectortrainset7_112011.fa',2) %write the output in CSV format, the last argument (2) corresponds to the phylum level. See the file quikrWriteOutput.m for more details.

+

+

+

+

+%This is an example of how to run Quikr

+

+%If you want to use the default training-set (that is, trainset7_112011.fa

+%from RPD), then execute the following

+

+%make sure Matlab/Octave is in the Quikr director

+cd /path/to/Quikr

+

+fastafilename='/path/to/your/fastafile.fasta'; %full path name to your data file

+xstar=quikr(fastafilename); %this will give the predicted reconstruction frequencies

+

+

+%xstar will be on the same basis as trainset7_112011.fa, so to get the

+%sequences that are predicted to be present in your sample:

+[headers,~]=fastaread('trainset7_112011.fa'); %read in the training database

+%note fastaread is not by default included in Octave, the fastaread.m file

+%is included in the Quikr download directory however, and is directly

+%compatible with Matlab.

+nonzeroentries=find(xstar); %get the indicies of the sequences quikr predicts are in your sample

+proportionscell=num2cell(xstar(nonzeroentries)); %convert the concentrations into a cell array

+namescell=headers(nonzeroentries); %Get the names of the sequences

+namesandproportions={namescell{:}; proportionscell{:}}; %This cell array contains the (unsorted) names of the reconstructed sequences and their concentrations (in the first and second columns respectively)

+%so to find which sequence is the most abundant in your mixture:

+[val,ind]=max(xstar(nonzeroentries)); %get the maximum value and it's position

+namesandproportions{1:2,ind} %note that this does not imply this specific strain or species is in your sample, just that phylum/class/order/family/genus this species belongs to is in your sample.

+

+

+% If you would like to use a custom training database, follow the following

+% steps

+fastafilename='/path/to/your/fastafile.fasta'; %full path name to your data file

+trainingdatabasefilename='/path/to/your/trainingdatabase.fasta'; %full path to the FASTA file you wish to use as a training database

+k=6; %pick a k-mer size

+trainingmatrix=quikrTrain(trainingdatabasefilename,k); %this will return the training database

+%then to do the reconstuction

+lambda=10000; %pick a lambda (larger lambda -> theoretically predicted concentrations are closer to actual concentrations), this depends on k-mer size picked, also size and condition of the TrainingMatrix

+xstar=quikrCustomTrained(trainingmatrix,fastafilename,k,lambda); %get the predicted reconstruction frequencies

+

+%again Xstar is on the same basis as the TrainingMatrix, so to get the

+%sequences that are predicted to be present in your sample:

+[headers,~]=fastaread(trainingdatabasefilename); %read in the training database

+nonzeroentries=find(xstar); %get the indicies of the sequences quikr predicts are in your sample

+proportionscell=num2cell(xstar(nonzeroentries)); %convert the concentrations into a cell array

+namescell=headers(nonzeroentries); %Get the names of the sequences

+namesandproportions={namescell{:}; proportionscell{:}}; %This cell array contains the (unsorted) names of the reconstructed sequences and their concentrations (in the first and second columns respectively)