Skip to content

Commit

Permalink
fix: getBlast handles spaces in path
Browse files Browse the repository at this point in the history 
also more informative error message
  • Loading branch information
@edkerk
edkerk committed Jun 26, 2024
1 parent 06ecc56 commit 4549f29
Show file tree
Hide file tree
Showing 4 changed files with 147 additions and 163 deletions.
142 changes: 69 additions & 73 deletions doc/external/getBlast.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -177,86 +177,82 @@ <h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" sr
0099
0100 <span class="comment">%Create a database for the new organism and blast each of the refFastaFiles</span>
0101 <span class="comment">%against it</span>
0102 [status, ~]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'makeblastdb'</span> binEnd]) <span class="string">'&quot; -in '</span> fastaFile{1} <span class="string">' -out &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -dbtype prot'</span>]);
0102 [status, message]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'makeblastdb'</span> binEnd]) <span class="string">'&quot; -in &quot;'</span> fastaFile{1} <span class="string">'&quot; -out &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -dbtype prot'</span>]);
0103 <span class="keyword">if</span> developMode
0104 blastReport.dbHashes.phr{numel(blastReport.dbHashes.phr)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.phr'</span>));
0105 blastReport.dbHashes.pot{numel(blastReport.dbHashes.pot)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.pot'</span>));
0106 blastReport.dbHashes.psq{numel(blastReport.dbHashes.psq)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.psq'</span>));
0107 blastReport.dbHashes.pto{numel(blastReport.dbHashes.pto)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.pto'</span>));
0108 <span class="keyword">end</span>
0109 <span class="keyword">if</span> status~=0
0110 EM=[<span class="string">'makeblastdb did not run successfully, error: '</span>, num2str(status)];
0111 dispEM(EM,true);
0112 <span class="keyword">end</span>
0113
0114 <span class="keyword">for</span> i=1:numel(refFastaFiles)
0115 <span class="keyword">if</span> ~hideVerbose
0116 fprintf([<span class="string">'BLASTing &quot;'</span> modelIDs{i} <span class="string">'&quot; against &quot;'</span> organismID{1} <span class="string">'&quot;..\n'</span>]);
0117 <span class="keyword">end</span>
0118 [status, ~]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'blastp'</span> binEnd]) <span class="string">'&quot; -query '</span> refFastaFiles{i} <span class="string">' -out &quot;'</span> outFile <span class="string">'_'</span> num2str(i) <span class="string">'&quot; -db &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -evalue 10e-5 -outfmt &quot;10 qseqid sseqid evalue pident length bitscore ppos&quot; -num_threads &quot;'</span> cores <span class="string">'&quot;'</span>]);
0119 <span class="keyword">if</span> developMode
0120 blastReport.blastTxtOutput{numel(blastReport.blastTxtOutput)+1}=importdata([outFile <span class="string">'_'</span> num2str(i)]);
0121 <span class="keyword">end</span>
0122 <span class="keyword">if</span> status~=0
0123 EM=[<span class="string">'blastp did not run successfully, error: '</span>, num2str(status)];
0124 dispEM(EM,true);
0125 <span class="keyword">end</span>
0126 <span class="keyword">end</span>
0127 delete([tmpDB filesep <span class="string">'tmpDB*'</span>]);
0128
0129 <span class="comment">%Then create a database for each of the reference organisms and blast the</span>
0130 <span class="comment">%new organism against them</span>
0131 <span class="keyword">for</span> i=1:numel(refFastaFiles)
0132 <span class="keyword">if</span> ~hideVerbose
0133 fprintf([<span class="string">'BLASTing &quot;'</span> organismID{1} <span class="string">'&quot; against &quot;'</span> modelIDs{i} <span class="string">'&quot;..\n'</span>]);
0134 <span class="keyword">end</span>
0135 [status, ~]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'makeblastdb'</span> binEnd]) <span class="string">'&quot; -in '</span> refFastaFiles{i} <span class="string">' -out &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -dbtype prot'</span>]);
0136 <span class="keyword">if</span> status~=0
0137 EM=[<span class="string">'makeblastdb did not run successfully, error: '</span>, num2str(status)];
0138 dispEM(EM,true);
0139 <span class="keyword">end</span>
0140 [status, ~]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'blastp'</span> binEnd]) <span class="string">'&quot; -query '</span> fastaFile{1} <span class="string">' -out &quot;'</span> outFile <span class="string">'_r'</span> num2str(i) <span class="string">'&quot; -db &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -evalue 10e-5 -outfmt &quot;10 qseqid sseqid evalue pident length bitscore ppos&quot; -num_threads &quot;'</span> cores <span class="string">'&quot;'</span>]);
0141 <span class="keyword">if</span> developMode
0142 blastReport.dbHashes.phr{numel(blastReport.dbHashes.phr)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.phr'</span>));
0143 blastReport.dbHashes.pot{numel(blastReport.dbHashes.pot)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.pot'</span>));
0144 blastReport.dbHashes.psq{numel(blastReport.dbHashes.psq)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.psq'</span>));
0145 blastReport.dbHashes.pto{numel(blastReport.dbHashes.pto)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.pto'</span>));
0146 blastReport.blastTxtOutput{numel(blastReport.blastTxtOutput)+1}=importdata([outFile <span class="string">'_r'</span> num2str(i)]);
0147 <span class="keyword">end</span>
0148 <span class="keyword">if</span> status~=0
0149 EM=[<span class="string">'blastp did not run successfully, error: '</span>, num2str(status)];
0150 dispEM(EM,true);
0151 <span class="keyword">end</span>
0152 delete([tmpDB filesep <span class="string">'tmpDB*'</span>]);
0153 <span class="keyword">end</span>
0154
0155 <span class="comment">%Done with the BLAST, do the parsing of the text files</span>
0156 <span class="keyword">for</span> i=1:numel(refFastaFiles)*2
0157 tempStruct=[];
0158 <span class="keyword">if</span> i&lt;=numel(refFastaFiles)
0159 tempStruct.fromId=modelIDs{i};
0160 tempStruct.toId=organismID{1};
0161 A=readtable([outFile <span class="string">'_'</span> num2str(i)],<span class="string">'Delimiter'</span>,<span class="string">','</span>,<span class="string">'Format'</span>,<span class="string">'%s%s%f%f%f%f%f'</span>);
0162 <span class="keyword">else</span>
0163 tempStruct.fromId=organismID{1};
0164 tempStruct.toId=modelIDs{i-numel(refFastaFiles)};
0165 A=readtable([outFile <span class="string">'_r'</span> num2str(i-numel(refFastaFiles))],<span class="string">'Delimiter'</span>,<span class="string">','</span>,<span class="string">'Format'</span>,<span class="string">'%s%s%f%f%f%f%f'</span>);
0166 <span class="keyword">end</span>
0167 tempStruct.fromGenes=A{:,1};
0168 tempStruct.toGenes=A{:,2};
0169 tempStruct.evalue=table2array(A(:,3));
0170 tempStruct.identity=table2array(A(:,4));
0171 tempStruct.aligLen=table2array(A(:,5));
0172 tempStruct.bitscore=table2array(A(:,6));
0173 tempStruct.ppos=table2array(A(:,7));
0174 blastStructure=[blastStructure tempStruct];
0175 <span class="keyword">end</span>
0176
0177 <span class="comment">%Remove the old tempfiles</span>
0178 delete([outFile <span class="string">'*'</span>]);
0179 <span class="comment">%Remove the temp fasta files</span>
0180 delete(files{:})
0181 <span class="keyword">end</span></pre></div>
0110 error(<span class="string">'makeblastdb did not run successfully, error:\n%s'</span>,strip(message))
0111 <span class="keyword">end</span>
0112
0113 <span class="keyword">for</span> i=1:numel(refFastaFiles)
0114 <span class="keyword">if</span> ~hideVerbose
0115 fprintf([<span class="string">'BLASTing &quot;'</span> modelIDs{i} <span class="string">'&quot; against &quot;'</span> organismID{1} <span class="string">'&quot;..\n'</span>]);
0116 <span class="keyword">end</span>
0117 [status, message]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'blastp'</span> binEnd]) <span class="string">'&quot; -query &quot;'</span> refFastaFiles{i} <span class="string">'&quot; -out &quot;'</span> outFile <span class="string">'_'</span> num2str(i) <span class="string">'&quot; -db &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -evalue 10e-5 -outfmt &quot;10 qseqid sseqid evalue pident length bitscore ppos&quot; -num_threads &quot;'</span> cores <span class="string">'&quot;'</span>]);
0118 <span class="keyword">if</span> developMode
0119 blastReport.blastTxtOutput{numel(blastReport.blastTxtOutput)+1}=importdata([outFile <span class="string">'_'</span> num2str(i)]);
0120 <span class="keyword">end</span>
0121 <span class="keyword">if</span> status~=0
0122 error(<span class="string">'blastp did not run successfully, error:\n%s'</span>,strip(message))
0123 <span class="keyword">end</span>
0124 <span class="keyword">end</span>
0125 delete([tmpDB filesep <span class="string">'tmpDB*'</span>]);
0126
0127 <span class="comment">%Then create a database for each of the reference organisms and blast the</span>
0128 <span class="comment">%new organism against them</span>
0129 <span class="keyword">for</span> i=1:numel(refFastaFiles)
0130 <span class="keyword">if</span> ~hideVerbose
0131 fprintf([<span class="string">'BLASTing &quot;'</span> organismID{1} <span class="string">'&quot; against &quot;'</span> modelIDs{i} <span class="string">'&quot;..\n'</span>]);
0132 <span class="keyword">end</span>
0133 [status, message]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'makeblastdb'</span> binEnd]) <span class="string">'&quot; -in &quot;'</span> refFastaFiles{i} <span class="string">'&quot; -out &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -dbtype prot'</span>]);
0134 <span class="keyword">if</span> status~=0
0135 error(<span class="string">'makeblastdb did not run successfully, error:\n%s'</span>,strip(message))
0136 <span class="keyword">end</span>
0137 [status, message]=system([<span class="string">'&quot;'</span> fullfile(ravenPath,<span class="string">'software'</span>,<span class="string">'blast+'</span>,[<span class="string">'blastp'</span> binEnd]) <span class="string">'&quot; -query &quot;'</span> fastaFile{1} <span class="string">'&quot; -out &quot;'</span> outFile <span class="string">'_r'</span> num2str(i) <span class="string">'&quot; -db &quot;'</span> fullfile(tmpDB, <span class="string">'tmpDB'</span>) <span class="string">'&quot; -evalue 10e-5 -outfmt &quot;10 qseqid sseqid evalue pident length bitscore ppos&quot; -num_threads &quot;'</span> cores <span class="string">'&quot;'</span>]);
0138 <span class="keyword">if</span> developMode
0139 blastReport.dbHashes.phr{numel(blastReport.dbHashes.phr)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.phr'</span>));
0140 blastReport.dbHashes.pot{numel(blastReport.dbHashes.pot)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.pot'</span>));
0141 blastReport.dbHashes.psq{numel(blastReport.dbHashes.psq)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.psq'</span>));
0142 blastReport.dbHashes.pto{numel(blastReport.dbHashes.pto)+1}=getMD5Hash(fullfile(tmpDB, <span class="string">'tmpDB.pto'</span>));
0143 blastReport.blastTxtOutput{numel(blastReport.blastTxtOutput)+1}=importdata([outFile <span class="string">'_r'</span> num2str(i)]);
0144 <span class="keyword">end</span>
0145 <span class="keyword">if</span> status~=0
0146 error(<span class="string">'blastp did not run successfully, error:\n%s'</span>,strip(message))
0147 <span class="keyword">end</span>
0148 delete([tmpDB filesep <span class="string">'tmpDB*'</span>]);
0149 <span class="keyword">end</span>
0150
0151 <span class="comment">%Done with the BLAST, do the parsing of the text files</span>
0152 <span class="keyword">for</span> i=1:numel(refFastaFiles)*2
0153 tempStruct=[];
0154 <span class="keyword">if</span> i&lt;=numel(refFastaFiles)
0155 tempStruct.fromId=modelIDs{i};
0156 tempStruct.toId=organismID{1};
0157 A=readtable([outFile <span class="string">'_'</span> num2str(i)],<span class="string">'Delimiter'</span>,<span class="string">','</span>,<span class="string">'Format'</span>,<span class="string">'%s%s%f%f%f%f%f'</span>);
0158 <span class="keyword">else</span>
0159 tempStruct.fromId=organismID{1};
0160 tempStruct.toId=modelIDs{i-numel(refFastaFiles)};
0161 A=readtable([outFile <span class="string">'_r'</span> num2str(i-numel(refFastaFiles))],<span class="string">'Delimiter'</span>,<span class="string">','</span>,<span class="string">'Format'</span>,<span class="string">'%s%s%f%f%f%f%f'</span>);
0162 <span class="keyword">end</span>
0163 tempStruct.fromGenes=A{:,1};
0164 tempStruct.toGenes=A{:,2};
0165 tempStruct.evalue=table2array(A(:,3));
0166 tempStruct.identity=table2array(A(:,4));
0167 tempStruct.aligLen=table2array(A(:,5));
0168 tempStruct.bitscore=table2array(A(:,6));
0169 tempStruct.ppos=table2array(A(:,7));
0170 blastStructure=[blastStructure tempStruct];
0171 <span class="keyword">end</span>
0172
0173 <span class="comment">%Remove the old tempfiles</span>
0174 delete([outFile <span class="string">'*'</span>]);
0175 <span class="comment">%Remove the temp fasta files</span>
0176 delete(files{:})
0177 <span class="keyword">end</span></pre></div>
<hr><address>Generated by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" title="Matlab Documentation in HTML">m2html</a></strong> &copy; 2005</address>
</body>
</html>
Loading

0 comments on commit 4549f29

Please sign in to comment.