blast_ui.api.inc 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868
  1. <?php
  2. /**
  3. * @file
  4. * Contains more generally applicable functions as well as some meant to help developers
  5. * Plug-in to the BLAST UI functionality
  6. */
  7. /**
  8. * Get a specific BlastDB.
  9. *
  10. * @param $identifiers
  11. * An array of identifiers used to determine which BLAST DB to retrieve.
  12. *
  13. * @return
  14. * A fully-loaded BLAST DB Node
  15. */
  16. function get_blast_database($identifiers) {
  17. $node = FALSE;
  18. if (isset($identifiers['nid'])) {
  19. $node = node_load($identifiers['nid']);
  20. }
  21. elseif (isset($identifiers['name'])) {
  22. $nid = db_query('SELECT nid FROM {blastdb} WHERE name=:name', array(':name' => $identifiers['name']))->fetchField();
  23. $node = node_load($nid);
  24. } elseif (isset($identifiers['path'])) {
  25. $nid = db_query('SELECT nid FROM {blastdb} WHERE path LIKE :path', array(':path' => db_like($identifiers['path']) . '%'))->fetchField();
  26. $node = node_load($nid);
  27. }
  28. return $node;
  29. }
  30. /**
  31. * Returns a list BLAST DATABASE options
  32. *
  33. * @param $type
  34. * The type of BLAST dabases to restrict the list to (ie: n: nucleotide or p: protein)
  35. *
  36. * @return
  37. * An array where the nid is the key and the value is the human-readable name of the option
  38. */
  39. function get_blast_database_options($type) {
  40. global $user;
  41. // Get all BlastDB nodes
  42. $nodes = get_blast_database_nodes();
  43. // Support obsolete database type n/p
  44. $obs_type = '';
  45. if ($type == 'protein') {
  46. $obs_type = 'p';
  47. }
  48. else {
  49. $obs_type = 'n';
  50. }
  51. $options = array();
  52. foreach ($nodes as $node) {
  53. if ( isset($node) && isset($node->db_dbtype) ) {
  54. if ( ($node->db_dbtype == $type) OR ($node->db_dbtype == $obs_type) ) {
  55. $options[$node->nid] = $node->db_name;
  56. }
  57. }
  58. }
  59. // Sort alphabetically
  60. asort($options);
  61. return $options;
  62. }
  63. function get_blast_database_nodes() {
  64. // Use the Entity API to get a list of BLAST Nodes to load
  65. // We use this function in order respect node access control so that
  66. // administrators can use this module in combination with a node access module
  67. // of their choice to limit access to specific BLAST databases.
  68. $query = new EntityFieldQuery();
  69. $query->entityCondition('entity_type', 'node')
  70. // Restrict to BLASTDB nodes.
  71. ->entityCondition('bundle', 'blastdb')
  72. // Restrict to Published nodes.
  73. ->propertyCondition('status', 1)
  74. // Restrict to nodes the current user has permission to view.
  75. ->addTag('node_access');
  76. $entities = $query->execute();
  77. // Get all BlastDB nodes
  78. return node_load_multiple(array_keys($entities['node']));
  79. }
  80. /**
  81. * Retrieve all the information for a blast job in a standardized node-like format.
  82. *
  83. * @param $job_id
  84. * The non-encoded tripal job_id.
  85. * @retun
  86. * An object describing the blast job.
  87. */
  88. function get_BLAST_job($job_id) {
  89. $blastjob = db_query('SELECT * FROM blastjob WHERE job_id=:id', array(':id' => $job_id))->fetchObject();
  90. if (!$blastjob) {
  91. return false;
  92. }
  93. $tripal_job = tripal_get_job($job_id);
  94. $job = new stdClass();
  95. $job->job_id = $job_id;
  96. $job->program = $blastjob->blast_program;
  97. $job->options = unserialize($blastjob->options);
  98. $job->date_submitted = $tripal_job->submit_date;
  99. $job->date_started = $tripal_job->start_time;
  100. $job->date_completed = $tripal_job->end_time;
  101. // TARGET BLAST DATABASE.
  102. // If a provided blast database was used then load details.
  103. if ($blastjob->target_blastdb ) {
  104. $job->blastdb = get_blast_database(array('nid' => $blastjob->target_blastdb));
  105. }
  106. // Otherwise the user uploaded their own database so provide what information we can.
  107. else {
  108. $job->blastdb = new stdClass();
  109. $job->blastdb->db_name = 'User Uploaded';
  110. $job->blastdb->db_path = $blastjob->target_file;
  111. $job->blastdb->linkout = new stdClass();
  112. $job->blastdb->linkout->none = TRUE;
  113. if ($job->program == 'blastp' OR $job->program == 'tblastn') {
  114. $job->blastdb->db_dbtype = 'protein';
  115. }
  116. else {
  117. $job->blastdb->db_dbtype = 'nucleotide';
  118. }
  119. }
  120. // FILES.
  121. $job->files = new stdClass();
  122. $job->files->query = $blastjob->query_file;
  123. $job->files->target = $blastjob->target_file;
  124. $job->files->result = new stdClass();
  125. $job->files->result->archive = $blastjob->result_filestub . '.asn';
  126. $job->files->result->xml = $blastjob->result_filestub . '.xml';
  127. $job->files->result->tsv = $blastjob->result_filestub . '.tsv';
  128. $job->files->result->html = $blastjob->result_filestub . '.html';
  129. $job->files->result->gff = $blastjob->result_filestub . '.gff';
  130. return $job;
  131. }
  132. /**
  133. * Run BLAST (should be called from the command-line)
  134. *
  135. * @param $program
  136. * Which BLAST program to run (ie: 'blastn', 'tblastn', tblastx', 'blastp','blastx')
  137. * @param $query
  138. * The full path and filename of the query FASTA file
  139. * @param $database
  140. * The full path and filename prefix (excluding .nhr, .nin, .nsq, etc.)
  141. * @param $output_filestub
  142. * The filename (not including path) to give the results. Should not include file type suffix
  143. * @param $options
  144. * An array of additional option where the key is the name of the option used by
  145. * BLAST (ie: 'num_alignments') and the value is relates to this particular
  146. * BLAST job (ie: 250)
  147. */
  148. function run_BLAST_tripal_job($program, $query, $database, $output_filestub, $options, $job_id = NULL) {
  149. $output_file = $output_filestub . '.asn';
  150. $output_file_xml = $output_filestub . '.xml';
  151. $output_file_tsv = $output_filestub . '.tsv';
  152. $output_file_html = $output_filestub . '.html';
  153. $output_file_gff = $output_filestub . '.gff';
  154. print "\nExecuting $program\n\n";
  155. print "Query: $query\n";
  156. print "Database: $database\n";
  157. print "Results File: $output_file\n";
  158. print "Options:\n";
  159. // Allow administrators to use an absolute path for these commands.
  160. // Defaults to using $PATH.
  161. $blast_path = variable_get('blast_path', '');
  162. $blast_threads = variable_get('blast_threads', 1);
  163. // Strip the extension off the BLAST target
  164. $database = preg_replace("/(.*)\.[pn]\w\w$/", '$1', $database);
  165. // The BLAST executeable.
  166. $program = $blast_path . $program;
  167. if (!file_exists($program)) {
  168. tripal_report_error(
  169. 'blast_ui',
  170. TRIPAL_ERROR,
  171. "Unable to find the BLAST executable (ie: /urs/bin/blastn). This can be changed in the admin settings; you supplied: @command",
  172. array('@command' => $program),
  173. array('print' => TRUE)
  174. );
  175. return FALSE;
  176. }
  177. // The blast db formatter executable.
  178. $blast_formatter_command = $blast_path . 'blast_formatter';
  179. if (!file_exists($blast_formatter_command)) {
  180. tripal_report_error(
  181. 'blast_ui',
  182. TRIPAL_ERROR,
  183. "Unable to find the BLAST Formatter executable (ie: /urs/bin/blast_formatter). This can be changed in the admin settings; you supplied: @command",
  184. array('@command' => $blast_formatter_command),
  185. array('print' => TRUE)
  186. );
  187. return FALSE;
  188. }
  189. // Note: all variables are escaped (adds single quotes around their values) for security reasons.
  190. $blast_cmd = escapeshellarg($program) . ' -query ' . escapeshellarg($query) . ' -db ' . escapeshellarg($database) . ' -out ' . escapeshellarg($output_file) . ' -outfmt=11';
  191. if (!empty($options)) {
  192. foreach ($options as $opt => $val) {
  193. $val = trim($val);
  194. if (!empty($val)) {
  195. print "\t$opt: $val\n";
  196. // We want to escape all the option values since they were supplied via
  197. // user input. These values should also have been checked in the
  198. // advanced form _validate functions but this adds an extra layer of
  199. // protection.
  200. $blast_cmd .= ' -' . escapeshellarg($opt) . ' ' . escapeshellarg($val);
  201. }
  202. }
  203. }
  204. // Setting the value of threads by admin page
  205. $blast_cmd .= ' -num_threads ' . escapeshellarg($blast_threads);
  206. print "\nExecuting the following BLAST command:\n" . $blast_cmd . "\n";
  207. system($blast_cmd);
  208. if (!file_exists($output_file)) {
  209. tripal_report_error(
  210. 'blast_ui',
  211. TRIPAL_ERROR,
  212. "BLAST did not complete successfully as is implied by the lack of output file (%file). The command run was @command",
  213. array('%file' => $output_file, '@command' => $blast_cmd),
  214. array('print' => TRUE)
  215. );
  216. return FALSE;
  217. }
  218. print "\nGenerating additional download formats...\n";
  219. print "\tXML\n";
  220. $format_cmd = escapeshellarg($blast_formatter_command) . ' -archive ' . escapeshellarg($output_file) . ' -outfmt 5 -out ' . escapeshellarg($output_file_xml);
  221. print "\t\tExecuting $format_cmd\n\n";
  222. system($format_cmd);
  223. if (!file_exists($output_file_xml)) {
  224. tripal_report_error(
  225. 'blast_ui',
  226. TRIPAL_ERROR,
  227. "Unable to convert BLAST ASN.1 archive to XML (%archive => %file).",
  228. array('%archive' => $output_file, '%file' => $output_file_xml),
  229. array('print' => TRUE)
  230. );
  231. }
  232. print "\tTab-delimited\n";
  233. $format_cmd = escapeshellarg($blast_formatter_command) . ' -archive ' . escapeshellarg($output_file) . ' -outfmt 7 -out ' . escapeshellarg($output_file_tsv);
  234. print "\t\tExecuting $format_cmd\n\n";
  235. system($format_cmd);
  236. if (!file_exists($output_file_tsv)) {
  237. tripal_report_error(
  238. 'blast_ui',
  239. TRIPAL_WARNING,
  240. "Unable to convert BLAST ASN.1 archive to Tabular Output (%archive => %file).",
  241. array('%archive' => $output_file, '%file' => $output_file_tsv),
  242. array('print' => TRUE)
  243. );
  244. }
  245. print "\tGFF\n";
  246. convert_tsv2gff3($output_file_tsv,$output_file_gff);
  247. if (!file_exists($output_file_gff)) {
  248. tripal_report_error(
  249. 'blast_ui',
  250. TRIPAL_WARNING,
  251. "Unable to convert BLAST Tabular Output to GFF Output (%archive => %file).",
  252. array('%archive' => $output_file, '%file' => $output_file_gff),
  253. array('print' => TRUE)
  254. );
  255. }
  256. print "\tHTML (includes alignments)\n";
  257. $format_cmd = escapeshellarg($blast_formatter_command) . ' -archive ' . escapeshellarg($output_file) . ' -outfmt 0 -out ' . escapeshellarg($output_file_html) . ' -html';
  258. print "\t\tExecuting $format_cmd\n\n";
  259. system($format_cmd);
  260. if (!file_exists($output_file_tsv)) {
  261. tripal_report_error(
  262. 'blast_ui',
  263. TRIPAL_WARNING,
  264. "Unable to convert BLAST ASN.1 archive to HTML Output (%archive => %file).",
  265. array('%archive' => $output_file, '%file' => $output_file_html),
  266. array('print' => TRUE)
  267. );
  268. }
  269. print "\nDone!\n";
  270. }
  271. /**
  272. * FASTA validating parser
  273. *
  274. * A sequence in FASTA format begins with a single-line description, followed
  275. * by lines of sequence data.The description line is distinguished from the
  276. * sequence data by a greater-than (">") symbol in the first column. The word
  277. * following the ">" symbol is the identifier of the sequence, and the rest of
  278. * the line is the description (both are optional). There should be no space
  279. * between the ">" and the first letter of the identifier. The sequence ends
  280. * if another line starting with a ">" appears which indicates the start of
  281. * another sequence.
  282. *
  283. * @param $type
  284. * The type of sequence to be validated (ie: either nucleotide or protein).
  285. * @param $sequence
  286. * A string of characters to be validated.
  287. *
  288. * @return
  289. * Return a boolean. 1 if the sequence does not pass the format valifation stage and 0 otherwise.
  290. *
  291. */
  292. function validate_fasta_sequence($type, $sequence) {
  293. //Includes IUPAC codes.
  294. $fastaSeqRegEx = ($type == 'nucleotide')
  295. ? '/^[ATCGNUKMBVSWDYRHatcgnukmbvswdyrh\[\/\]\s\n\r]*$/'
  296. : '/^[ABCDEFGHIKLMNPQRSTUVWYZXabcdefghiklmnpqrstuvwyzx\*\-\s\n\r]*$/';
  297. $defRegEx = '/^>.*(\\n|\\r)(.*)$/sm';
  298. if (preg_match($defRegEx, $sequence, $matches)) {
  299. if (isset($matches[2]) && $matches[2] != '' && preg_match($fastaSeqRegEx, $matches[2])) {
  300. return true;
  301. }
  302. }
  303. else if ($sequence != '' && preg_match($defRegEx, $sequence)) {
  304. return true;
  305. }
  306. return false;
  307. }
  308. /**
  309. * Retrieve the regex to capture the Link-out Accession from the Hit Def.
  310. *
  311. * @param $nid
  312. * The node ID of the BLAST database the hit is from.
  313. * @param $options
  314. * An array of options that can be passed to this function. Supported
  315. * options include:
  316. * -
  317. *
  318. * @return
  319. * A PHP regex for use with preg_match to cature the Link-out Accession.
  320. */
  321. function get_blastdb_linkout_regex($node, $options = array()) {
  322. if (empty($node->linkout->regex)) {
  323. switch ($node->linkout->regex_type) {
  324. case 'default':
  325. $regex = '/^(\S+).*/';
  326. break;
  327. case 'genbank':
  328. $regex = '/^gb\|([^\|])*\|.*/';
  329. break;
  330. case 'embl':
  331. $regex = '/^embl\|([^\|])*\|.*/';
  332. break;
  333. case 'swissprot':
  334. $regex = '/^sp\|([^\|])*\|.*/';
  335. break;
  336. }
  337. }
  338. else {
  339. $regex = $node->linkout->regex;
  340. }
  341. return $regex;
  342. }
  343. /**
  344. * Return a list of recent blast jobs to be displayed to the user.
  345. *
  346. * @param $programs
  347. * An array of blast programs you want jobs to be displayed for (ie: blastn, blastx, tblastn, blastp)
  348. *
  349. * @return
  350. * An array of recent jobs.
  351. */
  352. function get_recent_blast_jobs($programs = array()) {
  353. $filter_jobs = !empty($programs);
  354. // Retrieve any recent jobs from the session variable.
  355. if (isset($_SESSION['blast_jobs'])) {
  356. $jobs = array();
  357. foreach ($_SESSION['blast_jobs'] as $job_secret) {
  358. $add = TRUE;
  359. $job_id = blast_ui_reveal_secret($job_secret);
  360. if ($job = get_BLAST_job($job_id)) {
  361. // @TODO: Check that the results are still available.
  362. // This is meant to replace the arbitrary only show jobs executed less than 48 hrs ago.
  363. // Remove jobs from the list that are not of the correct program.
  364. if ($filter_jobs AND !in_array($job->program, $programs)) {
  365. $add = FALSE;
  366. }
  367. if ($add) {
  368. $job->query_summary = format_query_headers($job->files->query);
  369. $jobs[] = $job;
  370. }
  371. }
  372. }
  373. return $jobs;
  374. }
  375. else {
  376. return array();
  377. }
  378. }
  379. /**
  380. * Retrieve the number of recent jobs.
  381. */
  382. function get_number_of_recent_jobs() {
  383. if (isset($_SESSION['blast_jobs'])) {
  384. return sizeof($_SESSION['blast_jobs']);
  385. }
  386. return 0;
  387. }
  388. /**
  389. * Summarize a fasta file based on it's headers.
  390. *
  391. * @param $file
  392. * The full path to the FASTA file.
  393. *
  394. * @return
  395. * A string describing the number of sequences and often including the first query header.
  396. */
  397. function format_query_headers($file) {
  398. $headers = array();
  399. exec('grep ">" ' . escapeshellarg($file), $headers);
  400. // Easiest case: if there is only one query header then show it.
  401. if (sizeof($headers) == 1 AND isset($headers[0])) {
  402. return ltrim($headers[0], '>');
  403. }
  404. // If we have at least one header then show that along with the count of queries.
  405. elseif (isset($headers[0])) {
  406. return sizeof($headers) . ' queries including "' . ltrim($headers[0], '>') . '"';
  407. }
  408. // If they provided a sequence with no header.
  409. elseif (empty($headers)) {
  410. return 'Unnamed Query';
  411. }
  412. // At the very least show the count of queries.
  413. else {
  414. return sizeof($headers) . ' queries';
  415. }
  416. }
  417. /**
  418. * Sort recent blast jobs by the date they were submitted.
  419. * Ascending is in order of submission.
  420. *
  421. * THIS FUNCTION SHOULD BY USED BY USORT().
  422. */
  423. function sort_blast_jobs_by_date_submitted_asc($a, $b) {
  424. return ($a->date_submitted - $b->date_submitted);
  425. }
  426. /**
  427. * Sort recent blast jobs by the date they were submitted.
  428. * Descending is most recent first.
  429. *
  430. * THIS FUNCTION SHOULD BY USED BY USORT().
  431. */
  432. function sort_blast_jobs_by_date_submitted_desc($a, $b) {
  433. return ($b->date_submitted - $a->date_submitted);
  434. }
  435. /**
  436. * Generate an image of HSPs for a given hit.
  437. *
  438. * history:
  439. * 09/23/10 Carson created
  440. * 04/16/12 eksc adapted into POPcorn code
  441. * 03/12/15 deepak Adapted code into Tripal BLAST
  442. * 10/23/15 lacey Fixed deepak's code to be suitable for Drupal.
  443. *
  444. * @param $acc
  445. * target name
  446. * @param $name
  447. * query name, false if none
  448. * @param $tsize
  449. * target size
  450. * @param $qsize
  451. * query size
  452. * @param $hits
  453. * each hit represented in URL as: targetstart_targetend_hspstart_hspend;
  454. * @param $score
  455. * score for each hit
  456. *
  457. * @returm
  458. * A base64 encoded image representing the hit information.
  459. */
  460. function generate_blast_hit_image($acc = '', $scores, $hits, $tsize, $qsize, $name, $hit_name) {
  461. $tok = strtok($hits, ";");
  462. $b_hits = Array();
  463. while ($tok !== false) {
  464. $b_hits[] = $tok;
  465. $tok = strtok(";");
  466. }
  467. // extract score information from score param
  468. $tokscr = strtok($scores, ";");
  469. $b_scores = Array();
  470. while ($tokscr !== false) {
  471. $b_scores[] = $tokscr;
  472. $tokscr = strtok(";");
  473. }
  474. // image measurements
  475. $height = 200 + (count($b_hits) * 16);
  476. $width = 520;
  477. $img = imagecreatetruecolor($width, $height);
  478. $white = imagecolorallocate($img, 255, 255, 255);
  479. $black = imagecolorallocate($img, 0, 0, 0);
  480. $darkgray = imagecolorallocate($img, 100, 100, 100);
  481. $strong = imagecolorallocatealpha($img, 202, 0, 0, 15);
  482. $moderate = imagecolorallocatealpha($img, 204, 102, 0, 20);
  483. $present = imagecolorallocatealpha($img, 204, 204, 0, 35);
  484. $weak = imagecolorallocatealpha($img, 102, 204, 0, 50);
  485. $gray = imagecolorallocate($img, 190, 190, 190);
  486. $lightgray = $white; //imagecolorallocate($img, 230, 230, 230);
  487. imagefill($img, 0, 0, $lightgray);
  488. // Target coordinates
  489. $maxlength = 300;
  490. $t_length = ($tsize > $qsize)
  491. ? $maxlength : $maxlength - 50;
  492. $q_length = ($qsize > $tsize)
  493. ? $maxlength : $maxlength - 50;
  494. $tnormal = $t_length / $tsize;
  495. $qnormal = $q_length / $qsize;
  496. $t_ystart = 30;
  497. $t_yend = $t_ystart + 20;
  498. $t_xstart = 50;
  499. $t_xend = $t_xstart + $t_length;
  500. $t_center = $t_xstart + ($t_length / 2);
  501. // Target labels
  502. $warn = '"'. $hit_name . '"';
  503. imagestring($img, 5, $t_xstart, $t_ystart-20, $acc.$warn, $black);
  504. imagestring($img, 3, 5, $t_ystart+2, "Target", $black);
  505. // Draw bar representing target
  506. imagefilledrectangle($img, $t_xstart, $t_ystart, $t_xend, $t_yend, $gray);
  507. imagerectangle($img, $t_xstart, $t_ystart, $t_xend, $t_yend, $darkgray);
  508. // query coordinates
  509. $q_maxheight = 250;
  510. $q_ystart = $t_yend + 100;
  511. $q_yend = $q_ystart + 20;
  512. $q_xstart = $t_center - $q_length / 2;
  513. $q_xend = $q_xstart + $q_length;
  514. $q_center = ($q_xend + $q_xstart) / 2;
  515. $q_xwidth = $q_xend - $q_xstart;
  516. // Query labels
  517. imagestring($img, 5, $q_xstart, $q_yend+2, $name, $black);
  518. imagestring($img, 3, $q_xstart, $q_ystart+2, 'Query', $black);
  519. // Draw bar representing query
  520. imagefilledrectangle($img, $q_xstart, $q_ystart, $q_xend, $q_yend, $gray);
  521. imagerectangle($img ,$q_xstart, $q_ystart, $q_xend, $q_yend, $darkgray);
  522. // HSP bars will start here
  523. $hsp_bary = $q_yend + 20;
  524. // Draw solids for HSP alignments
  525. for ($ii=count($b_hits)-1; $ii>=0; $ii--) {
  526. // alignment
  527. $cur_hit = $b_hits[$ii];
  528. $cur_score = intval($b_scores[$ii]);
  529. // set color according to score
  530. $cur_color = $darkgray;
  531. if ($cur_score > 200) {
  532. $cur_color = $strong;
  533. }
  534. else if ($cur_score > 80 && $cur_score <= 200) {
  535. $cur_color = $moderate;
  536. }
  537. else if ($cur_score > 50 && $cur_score <= 80) {
  538. $cur_color = $present;
  539. }
  540. else if ($cur_score > 40 && $cur_score <= 50) {
  541. $cur_color = $weak;
  542. }
  543. $t_start = $tnormal * intval(strtok($cur_hit, "_")) + $t_xstart;
  544. $t_end = $tnormal * intval(strtok("_")) + $t_xstart;
  545. $q_start = $qnormal * intval(strtok("_")) + $q_xstart;
  546. $q_end = $qnormal * intval(strtok("_")) + $q_xstart;
  547. $hit1_array = array($t_start, $t_yend, $t_end, $t_yend, $q_end,
  548. $q_ystart, $q_start, $q_ystart);
  549. // HSP coords
  550. imagefilledpolygon($img, $hit1_array, 4, $cur_color);
  551. }//each hit
  552. // Draw lines over fills for HSP alignments
  553. for ($ii=0; $ii<count($b_hits); $ii++) {
  554. // alignment
  555. $cur_hit = $b_hits[$ii];
  556. $t_start = $tnormal * intval(strtok($cur_hit, "_")) + $t_xstart;
  557. $t_end = $tnormal * intval(strtok("_")) + $t_xstart;
  558. $q_start = $qnormal * intval(strtok("_")) + $q_xstart;
  559. $q_end = $qnormal * intval(strtok("_")) + $q_xstart;
  560. $hit1_array = array($t_start, $t_yend, $t_end, $t_yend, $q_end, $q_ystart,
  561. $q_start, $q_ystart,);
  562. imagerectangle($img, $t_start, $t_ystart, $t_end, $t_yend, $black);
  563. imagerectangle($img, $q_start, $q_ystart, $q_end, $q_yend, $black);
  564. imagepolygon ($img, $hit1_array, 4, $black);
  565. // show HSP
  566. imagestring($img, 3, 2, $hsp_bary, ($acc ."HSP" . ($ii + 1)), $black);
  567. $cur_score = intval($b_scores[$ii]);
  568. // set color according to score
  569. $cur_color = $darkgray;
  570. if ($cur_score > 200) {
  571. $cur_color = $strong;
  572. }
  573. else if ($cur_score > 80 && $cur_score <= 200) {
  574. $cur_color = $moderate;
  575. }
  576. else if ($cur_score > 50 && $cur_score <= 80) {
  577. $cur_color = $present;
  578. }
  579. else if ($cur_score > 40 && $cur_score <= 50) {
  580. $cur_color = $weak;
  581. }
  582. imagefilledrectangle($img, $q_start, $hsp_bary, $q_end, $hsp_bary+10, $cur_color);
  583. $hsp_bary += 15;
  584. }//each hit
  585. // Draw the key
  586. $xchart = 390;
  587. $ychart = 10;
  588. $fontsize = 4;
  589. $yinc = 20;
  590. $ywidth = 7;
  591. $xinc = 10;
  592. imagestring($img, 5, $xchart, $ychart - 5, "Bit Scores", $black);
  593. imagestring($img, $fontsize, $xchart + $yinc + $xinc,$ychart + ($yinc * 1) + $ywidth, ">= 200" , $black);
  594. imagestring($img, $fontsize, $xchart + $yinc + $xinc,$ychart + ($yinc * 2) + $ywidth, "80 - 200" , $black);
  595. imagestring($img, $fontsize, $xchart + $yinc + $xinc,$ychart + ($yinc * 3) + $ywidth, "50 - 80" , $black);
  596. imagestring($img, $fontsize, $xchart + $yinc + $xinc,$ychart + ($yinc * 4) + $ywidth, "40 - 50" , $black);
  597. imagestring($img, $fontsize, $xchart + $yinc + $xinc,$ychart + ($yinc * 5) + $ywidth, "< 40" , $black);
  598. imagefilledRectangle($img, $xchart, $ychart + ($yinc * 1) + $xinc, $xchart + $yinc, $ychart + ($yinc * 2), $strong);
  599. imagefilledRectangle($img, $xchart, $ychart + ($yinc * 2) + $xinc, $xchart + $yinc, $ychart + ($yinc * 3), $moderate);
  600. imagefilledRectangle($img, $xchart, $ychart + ($yinc * 3) + $xinc, $xchart + $yinc, $ychart + ($yinc * 4), $present);
  601. imagefilledRectangle($img, $xchart, $ychart + ($yinc * 4) + $xinc, $xchart + $yinc, $ychart + ($yinc * 5), $weak);
  602. imagefilledRectangle($img, $xchart, $ychart + ($yinc * 5) + $xinc, $xchart + $yinc, $ychart + ($yinc * 6), $darkgray);
  603. // Now, we have a completed image resource and need to change it to an actual image
  604. // that can be displayed. This is done using imagepng() but unfortuatly that function
  605. // either saves the image to a file or outputs it directly to the screen. Thus, we use
  606. // the following code to capture it and base64 encode it.
  607. ob_start(); // Start buffering the output
  608. imagepng($img, null, 0, PNG_NO_FILTER);
  609. $b64_img = base64_encode(ob_get_contents()); // Get what we've just outputted and base64 it
  610. imagedestroy($img);
  611. ob_end_clean();
  612. return $b64_img;
  613. }
  614. /**
  615. * Convert tsv blast output to gff output file.
  616. *
  617. * Created by Sofia Robb
  618. * 09/15/2016
  619. * counter bugfix 10/27/2016
  620. *
  621. * The subject (hit) will be the source feature.
  622. * The query will be the target.
  623. *
  624. * @todo: find a more efficient way since currently the first loop stores all the blast
  625. * results into an array and then the second loop prints them.
  626. *
  627. * @param $blast_tsv
  628. * The name of the blast tsv output file.
  629. * @param $blast_gff
  630. * The name of the blast gff output file.
  631. */
  632. function convert_tsv2gff3($blast_tsv,$blast_gff){
  633. // Open a new file for writting the gff.
  634. $gff = fopen($blast_gff,"w");
  635. fwrite($gff,"##gff-version 3\n");
  636. // Open the TSV file to read from.
  637. $tsv = fopen($blast_tsv, "r") or die("Unable to open tsv file!");
  638. // For each line in the TSV file...
  639. // Need to go thru each line of tsv to find the first and last hsp of a hit.
  640. $last_s = NULL;
  641. $hsp = NULL;
  642. $HitResult=array();
  643. while(!feof($tsv)) {
  644. $line = fgets($tsv);
  645. $line = rtrim($line);
  646. // Skip the line if it's empty.
  647. if (preg_match('/^#/',$line) or preg_match('/^\s*$/',$line)){
  648. continue;
  649. }
  650. ## for keeping track of new queries and hits
  651. // Each line has the following parts:
  652. // 0: query id,
  653. // 1: subject id,
  654. // 2: % identity,
  655. // 3: alignment length,
  656. // 4: mismatches,
  657. // 5: gap opens,
  658. // 6: q. start,
  659. // 7: q. end,
  660. // 8: s. start,
  661. // 9: s. end,
  662. // 10: evalue,
  663. // 11: bit score
  664. $parts = preg_split('/\t/', $line);
  665. // Assign the important parts of the line to readable variables.
  666. $s = $parts[1];
  667. $q = $parts[0];
  668. $ss = $parts[8];
  669. $se = $parts[9];
  670. $qs = $parts[6];
  671. $qe = $parts[7];
  672. $e = $parts[10];
  673. // if this is a new hit print the last and
  674. // empty the $HitResult array and
  675. // reset hsp counter
  676. if ($last_s != NULL and $s != $last_s ) {
  677. printGFF_parent_children($gff,$HitResult);
  678. $HitResult = array();
  679. $hsp=0;
  680. }
  681. // every line is a new hsp
  682. $hsp++;
  683. // determine query strand to use in match_part line, no need to store, just print
  684. $q_strand = '+';
  685. if ($qs > $qe) {
  686. list($qs,$qe) = array($qe,$qs);
  687. $q_strand = '-';
  688. }
  689. // determine subject (hit) strand to use in match line, needs to be stored
  690. $HitResult["$s,$q"]['strand']='+';
  691. list($start,$end) = array($ss,$se);
  692. if($ss > $se) {
  693. list($start,$end) = array($se,$ss);
  694. $HitResult["$s,$q"]['strand']='-';
  695. }
  696. // store smallest start
  697. if (!array_key_exists('SS',$HitResult["$s,$q"]) or $ss < $HitResult["$s,$q"]['SS']) {
  698. $HitResult["$s,$q"]['SS'] = $ss;
  699. }
  700. // store largest end
  701. if (!array_key_exists('SE',$HitResult["$s,$q"]) or $se > $HitResult["$s,$q"]['SE']) {
  702. $HitResult["$s,$q"]['SE'] = $se;
  703. }
  704. // store best evalue
  705. if (!array_key_exists('E',$HitResult["$s,$q"]) or $e < $HitResult["$s,$q"]['E']) {
  706. $HitResult["$s,$q"]['E'] = $e;
  707. }
  708. // generate the match_part line for each hsp
  709. $HitResult["$s,$q"]['HSPs'][] = join("\t", array($s, "BLASTRESULT" , "match_part" , $start , $end , $e , $HitResult["$s,$q"]['strand'] , '.' , "ID=$s.$q.$hsp;Parent=$s.$q;Target=$q $qs $qe $q_strand"));
  710. $last_s = $s;
  711. } // end tsv file while
  712. // print hit and hsp for the last hit
  713. printGFF_parent_children($gff,$HitResult);
  714. // Close the files.
  715. fclose($tsv);
  716. fclose($gff);
  717. }
  718. /**
  719. * printGFF_parent_children
  720. * prints the GFF parent feature and all of its children features
  721. *
  722. *
  723. * @param $blast_feature_array
  724. * an array of the all the child features which is used to generate the smallest and largest coordinates for the parent
  725. *
  726. *
  727. */
  728. function printGFF_parent_children ($gff,$blast_feature_array){
  729. foreach ($blast_feature_array as $sq => $value ) {
  730. list ($s,$q) = preg_split('/,/' , $sq);
  731. $evalue = $blast_feature_array["$s,$q"]['E'];
  732. $parent = join ("\t", array($s, "BLASTRESULT" , "match" , $blast_feature_array["$s,$q"]['SS'] , $blast_feature_array["$s,$q"]['SE'] , $blast_feature_array["$s,$q"]['E'] , $blast_feature_array["$s,$q"]['strand'] , '.' , "ID=$s.$q;Name=$q($evalue)")) . "\n";
  733. $child = join ("\n",$blast_feature_array["$s,$q"]['HSPs']) . "\n";
  734. fwrite($gff,$parent);
  735. fwrite($gff,$child);
  736. }
  737. }