1 /*
2 
3   VSEARCH: a versatile open source tool for metagenomics
4 
5   Copyright (C) 2014-2021, Torbjorn Rognes, Frederic Mahe and Tomas Flouri
6   All rights reserved.
7 
8   Contact: Torbjorn Rognes <torognes@ifi.uio.no>,
9   Department of Informatics, University of Oslo,
10   PO Box 1080 Blindern, NO-0316 Oslo, Norway
11 
12   This software is dual-licensed and available under a choice
13   of one of two licenses, either under the terms of the GNU
14   General Public License version 3 or the BSD 2-Clause License.
15 
16 
17   GNU General Public License version 3
18 
19   This program is free software: you can redistribute it and/or modify
20   it under the terms of the GNU General Public License as published by
21   the Free Software Foundation, either version 3 of the License, or
22   (at your option) any later version.
23 
24   This program is distributed in the hope that it will be useful,
25   but WITHOUT ANY WARRANTY; without even the implied warranty of
26   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27   GNU General Public License for more details.
28 
29   You should have received a copy of the GNU General Public License
30   along with this program.  If not, see <http://www.gnu.org/licenses/>.
31 
32 
33   The BSD 2-Clause License
34 
35   Redistribution and use in source and binary forms, with or without
36   modification, are permitted provided that the following conditions
37   are met:
38 
39   1. Redistributions of source code must retain the above copyright
40   notice, this list of conditions and the following disclaimer.
41 
42   2. Redistributions in binary form must reproduce the above copyright
43   notice, this list of conditions and the following disclaimer in the
44   documentation and/or other materials provided with the distribution.
45 
46   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
47   "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
48   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
49   FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
50   COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
51   INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
52   BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
53   LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
54   CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
55   LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
56   ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
57   POSSIBILITY OF SUCH DAMAGE.
58 
59 */
60 
61 #include "vsearch.h"
62 
results_show_fastapairs_one(FILE * fp,struct hit * hp,char * query_head,char * qsequence,int64_t qseqlen,char * rc)63 void results_show_fastapairs_one(FILE * fp,
64                                  struct hit * hp,
65                                  char * query_head,
66                                  char * qsequence,
67                                  int64_t qseqlen,
68                                  char * rc)
69 {
70   /* http://www.drive5.com/usearch/manual/fastapairs.html */
71 
72   if (hp)
73     {
74       char * qrow = align_getrow(hp->strand ? rc : qsequence,
75                                  hp->nwalignment,
76                                  hp->nwalignmentlength,
77                                  0);
78       fasta_print_general(fp,
79                           nullptr,
80                           qrow + hp->trim_q_left + hp->trim_t_left,
81                           hp->internal_alignmentlength,
82                           query_head,
83                           strlen(query_head),
84                           0,
85                           0,
86                           -1.0,
87                           -1,
88                           -1,
89                           nullptr,
90                           0.0);
91       xfree(qrow);
92 
93       char * trow = align_getrow(db_getsequence(hp->target),
94                                  hp->nwalignment,
95                                  hp->nwalignmentlength,
96                                  1);
97       fasta_print_general(fp,
98                           nullptr,
99                           trow + hp->trim_q_left + hp->trim_t_left,
100                           hp->internal_alignmentlength,
101                           db_getheader(hp->target),
102                           db_getheaderlen(hp->target),
103                           0,
104                           0,
105                           -1.0,
106                           -1,
107                           -1,
108                           nullptr,
109                           0.0);
110       xfree(trow);
111 
112       fprintf(fp, "\n");
113     }
114 }
115 
116 
results_show_blast6out_one(FILE * fp,struct hit * hp,char * query_head,char * qsequence,int64_t qseqlen,char * rc)117 void results_show_blast6out_one(FILE * fp,
118                                 struct hit * hp,
119                                 char * query_head,
120                                 char * qsequence,
121                                 int64_t qseqlen,
122                                 char * rc)
123 {
124 
125   /*
126      http://www.drive5.com/usearch/manual/blast6out.html
127 
128      query label
129      target label
130      percent identity
131      alignment length
132      number of mismatches
133      number of gap opens
134      1-based position of start in query
135      1-based position of end in query
136      1-based position of start in target
137      1-based position of end in target
138      E-value
139      bit score
140 
141      Note that USEARCH shows 13 fields when there is no hit,
142      but only 12 when there is a hit. Fixed in VSEARCH.
143   */
144 
145   if (hp)
146     {
147       int qstart, qend;
148 
149       if (hp->strand)
150         {
151           /* minus strand */
152           qstart = qseqlen;
153           qend = 1;
154         }
155       else
156         {
157           /* plus strand */
158           qstart = 1;
159           qend = qseqlen;
160         }
161 
162       fprintf(fp,
163               "%s\t%s\t%.1f\t%d\t%d\t%d\t%d\t%d\t%d\t%" PRIu64 "\t%d\t%d\n",
164               query_head,
165               db_getheader(hp->target),
166               hp->id,
167               hp->internal_alignmentlength,
168               hp->mismatches,
169               hp->internal_gaps,
170               qstart,
171               qend,
172               1,
173               db_getsequencelen(hp->target),
174               -1,
175               0);
176     }
177   else
178     {
179         fprintf(fp, "%s\t*\t0.0\t0\t0\t0\t0\t0\t0\t0\t-1\t0\n", query_head);
180     }
181 }
182 
results_show_uc_one(FILE * fp,struct hit * hp,char * query_head,char * qsequence,int64_t qseqlen,char * rc,int clusterno)183 void results_show_uc_one(FILE * fp,
184                          struct hit * hp,
185                          char * query_head,
186                          char * qsequence,
187                          int64_t qseqlen,
188                          char * rc,
189                          int clusterno)
190 {
191   /*
192     http://www.drive5.com/usearch/manual/ucout.html
193 
194     Columns:
195     H/N
196     cluster no (0-based) (target sequence no)
197     sequence length (query)
198     percent identity
199     strand: + or -
200     0
201     0
202     compressed alignment, e.g. 9I92M14D, or "=" if perfect alignment
203     query label
204     target label
205   */
206 
207   if (hp)
208     {
209       bool perfect;
210 
211       if (opt_cluster_fast)
212         {
213           /* cluster_fast */
214           /* use = for identical sequences ignoring terminal gaps */
215           perfect = (hp->matches == hp->internal_alignmentlength);
216         }
217       else
218         {
219           /* cluster_size, cluster_smallmem, cluster_unoise */
220           /* usearch_global, search_exact, allpairs_global */
221           /* use = for strictly identical sequences */
222           perfect = (hp->matches == hp->nwalignmentlength);
223         }
224 
225       fprintf(fp,
226               "H\t%d\t%" PRId64 "\t%.1f\t%c\t0\t0\t%s\t%s\t%s\n",
227               clusterno,
228               qseqlen,
229               hp->id,
230               hp->strand ? '-' : '+',
231               perfect ? "=" : hp->nwalignment,
232               query_head,
233               db_getheader(hp->target));
234     }
235   else {
236     fprintf(fp, "N\t*\t*\t*\t.\t*\t*\t*\t%s\t*\n", query_head);
237 
238         }
239 }
240 
results_show_userout_one(FILE * fp,struct hit * hp,char * query_head,char * qsequence,int64_t qseqlen,char * rc)241 void results_show_userout_one(FILE * fp, struct hit * hp,
242                               char * query_head,
243                               char * qsequence, int64_t qseqlen,
244                               char * rc)
245 {
246 
247   /*
248     http://drive5.com/usearch/manual/userout.html
249     qlo, qhi, tlo, thi and raw are given more meaningful values here
250   */
251 
252   for (int c = 0; c < userfields_requested_count; c++)
253     {
254       if (c) {
255         fprintf(fp, "\t");
256 
257         }
258 
259       int field = userfields_requested[c];
260 
261       char * tsequence = nullptr;
262       int64_t tseqlen = 0;
263       char * t_head = nullptr;
264 
265       if (hp)
266         {
267           tsequence = db_getsequence(hp->target);
268           tseqlen = db_getsequencelen(hp->target);
269           t_head = db_getheader(hp->target);
270         }
271 
272       char * qrow;
273       char * trow;
274 
275       switch (field)
276         {
277         case 0: /* query */
278           fprintf(fp, "%s", query_head);
279           break;
280         case 1: /* target */
281           fprintf(fp, "%s", hp ? t_head : "*");
282           break;
283         case 2: /* evalue */
284           fprintf(fp, "-1");
285           break;
286         case 3: /* id */
287           fprintf(fp, "%.1f", hp ? hp->id : 0.0);
288           break;
289         case 4: /* pctpv */
290           fprintf(fp, "%.1f", (hp && (hp->internal_alignmentlength > 0)) ? 100.0 * hp->matches / hp->internal_alignmentlength : 0.0);
291           break;
292         case 5: /* pctgaps */
293           fprintf(fp, "%.1f", (hp && (hp->internal_alignmentlength > 0)) ? 100.0 * hp->internal_indels / hp->internal_alignmentlength : 0.0);
294           break;
295         case 6: /* pairs */
296           fprintf(fp, "%d", hp ? hp->matches + hp->mismatches : 0);
297           break;
298         case 7: /* gaps */
299           fprintf(fp, "%d", hp ? hp->internal_indels : 0);
300           break;
301         case 8: /* qlo */
302           fprintf(fp, "%" PRId64, hp ? (hp->strand ? qseqlen : 1) : 0);
303           break;
304         case 9: /* qhi */
305           fprintf(fp, "%" PRId64, hp ? (hp->strand ? 1 : qseqlen) : 0);
306           break;
307         case 10: /* tlo */
308           fprintf(fp, "%d", hp ? 1 : 0);
309           break;
310         case 11: /* thi */
311           fprintf(fp, "%" PRId64, tseqlen);
312           break;
313         case 12: /* pv */
314           fprintf(fp, "%d", hp ? hp->matches : 0);
315           break;
316         case 13: /* ql */
317           fprintf(fp, "%" PRId64, qseqlen);
318           break;
319         case 14: /* tl */
320           fprintf(fp, "%" PRId64, hp ? tseqlen : 0);
321           break;
322         case 15: /* qs */
323           fprintf(fp, "%" PRId64, qseqlen);
324           break;
325         case 16: /* ts */
326           fprintf(fp, "%" PRId64, hp ? tseqlen : 0);
327           break;
328         case 17: /* alnlen */
329           fprintf(fp, "%d", hp ? hp->internal_alignmentlength : 0);
330           break;
331         case 18: /* opens */
332           fprintf(fp, "%d", hp ? hp->internal_gaps : 0);
333           break;
334         case 19: /* exts */
335           fprintf(fp, "%d", hp ? hp->internal_indels - hp->internal_gaps : 0);
336           break;
337         case 20: /* raw */
338           fprintf(fp, "%d", hp ? hp->nwscore : 0);
339           break;
340         case 21: /* bits */
341           fprintf(fp, "%d", 0);
342           break;
343         case 22: /* aln */
344           if (hp) {
345             align_fprint_uncompressed_alignment(fp, hp->nwalignment);
346 
347         }
348           break;
349         case 23: /* caln */
350           if (hp) {
351             fprintf(fp, "%s", hp->nwalignment);
352 
353         }
354           break;
355         case 24: /* qstrand */
356           if (hp) {
357             fprintf(fp, "%c", hp->strand ? '-' : '+');
358 
359         }
360           break;
361         case 25: /* tstrand */
362           if (hp) {
363             fprintf(fp, "%c", '+');
364 
365         }
366           break;
367         case 26: /* qrow */
368           if (hp)
369             {
370               qrow = align_getrow(hp->strand ? rc : qsequence,
371                                   hp->nwalignment,
372                                   hp->nwalignmentlength,
373                                   0);
374               fprintf(fp, "%.*s",
375                       (int)(hp->internal_alignmentlength),
376                       qrow + hp->trim_q_left + hp->trim_t_left);
377               xfree(qrow);
378             }
379           break;
380         case 27: /* trow */
381           if (hp)
382             {
383               trow = align_getrow(tsequence,
384                                   hp->nwalignment,
385                                   hp->nwalignmentlength,
386                                   1);
387               fprintf(fp, "%.*s",
388                       (int)(hp->internal_alignmentlength),
389                       trow + hp->trim_q_left + hp->trim_t_left);
390               xfree(trow);
391             }
392           break;
393         case 28: /* qframe */
394           fprintf(fp, "+0");
395           break;
396         case 29: /* tframe */
397           fprintf(fp, "+0");
398           break;
399         case 30: /* mism */
400           fprintf(fp, "%d", hp ? hp->mismatches : 0);
401           break;
402         case 31: /* ids */
403           fprintf(fp, "%d", hp ? hp->matches : 0);
404           break;
405         case 32: /* qcov */
406           fprintf(fp, "%.1f",
407                   hp ? 100.0 * (hp->matches + hp->mismatches) / qseqlen : 0.0);
408           break;
409         case 33: /* tcov */
410           fprintf(fp, "%.1f",
411                   hp ? 100.0 * (hp->matches + hp->mismatches) / tseqlen : 0.0);
412           break;
413         case 34: /* id0 */
414           fprintf(fp, "%.1f", hp ? hp->id0 : 0.0);
415           break;
416         case 35: /* id1 */
417           fprintf(fp, "%.1f", hp ? hp->id1 : 0.0);
418           break;
419         case 36: /* id2 */
420           fprintf(fp, "%.1f", hp ? hp->id2 : 0.0);
421           break;
422         case 37: /* id3 */
423           fprintf(fp, "%.1f", hp ? hp->id3 : 0.0);
424           break;
425         case 38: /* id4 */
426           fprintf(fp, "%.1f", hp ? hp->id4 : 0.0);
427           break;
428 
429           /* new internal alignment coordinates */
430 
431         case 39: /* qilo */
432           fprintf(fp, "%d", hp ? hp->trim_q_left + 1 : 0);
433           break;
434         case 40: /* qihi */
435           fprintf(fp, "%" PRId64, hp ? qseqlen - hp->trim_q_right : 0);
436           break;
437         case 41: /* tilo */
438           fprintf(fp, "%d", hp ? hp->trim_t_left + 1 : 0);
439           break;
440         case 42: /* tihi */
441           fprintf(fp, "%" PRId64, hp ? tseqlen - hp->trim_t_right : 0);
442           break;
443         }
444     }
445   fprintf(fp, "\n");
446 }
447 
results_show_alnout(FILE * fp,struct hit * hits,int hitcount,char * query_head,char * qsequence,int64_t qseqlen,char * rc)448 void results_show_alnout(FILE * fp,
449                          struct hit * hits,
450                          int hitcount,
451                          char * query_head,
452                          char * qsequence,
453                          int64_t qseqlen,
454                          char * rc)
455 {
456   /* http://drive5.com/usearch/manual/alnout.html */
457 
458   if (hitcount)
459     {
460       fprintf(fp, "\n");
461 
462       fprintf(fp,"Query >%s\n", query_head);
463       fprintf(fp," %%Id   TLen  Target\n");
464 
465       double top_hit_id = hits[0].id;
466 
467       for(int t = 0; t < hitcount; t++)
468         {
469           struct hit * hp = hits + t;
470 
471           if (opt_top_hits_only && (hp->id < top_hit_id)) {
472             break;
473 
474         }
475 
476           fprintf(fp,"%3.0f%% %6" PRIu64 "  %s\n",
477                   hp->id,
478                   db_getsequencelen(hp->target),
479                   db_getheader(hp->target));
480         }
481 
482       for(int t = 0; t < hitcount; t++)
483         {
484           struct hit * hp = hits + t;
485 
486           if (opt_top_hits_only && (hp->id < top_hit_id)) {
487             break;
488 
489         }
490 
491           fprintf(fp,"\n");
492 
493 
494           char * dseq = db_getsequence(hp->target);
495           int64_t dseqlen = db_getsequencelen(hp->target);
496 
497           int qlenlen = snprintf(nullptr, 0, "%" PRId64, qseqlen);
498           int tlenlen = snprintf(nullptr, 0, "%" PRId64, dseqlen);
499           int numwidth = MAX(qlenlen, tlenlen);
500 
501           fprintf(fp," Query %*" PRId64 "nt >%s\n", numwidth,
502                   qseqlen, query_head);
503           fprintf(fp,"Target %*" PRId64 "nt >%s\n", numwidth,
504                   dseqlen, db_getheader(hp->target));
505 
506           int rowlen = opt_rowlen == 0 ? qseqlen+dseqlen : opt_rowlen;
507 
508           align_show(fp,
509                      qsequence,
510                      qseqlen,
511                      hp->trim_q_left,
512                      "Qry",
513                      dseq,
514                      dseqlen,
515                      hp->trim_t_left,
516                      "Tgt",
517                      hp->nwalignment + hp->trim_aln_left,
518                      strlen(hp->nwalignment)
519                      - hp->trim_aln_left - hp->trim_aln_right,
520                      numwidth,
521                      3,
522                      rowlen,
523                      hp->strand);
524 
525           fprintf(fp, "\n%d cols, %d ids (%3.1f%%), %d gaps (%3.1f%%)\n",
526                   hp->internal_alignmentlength,
527                   hp->matches,
528                   hp->id,
529                   hp->internal_indels,
530                   hp->internal_alignmentlength > 0 ?
531                   100.0 * hp->internal_indels / hp->internal_alignmentlength :
532                   0.0);
533 
534 #if 0
535           fprintf(fp, "%d kmers, %d score, %d gap opens. %s %s %d %d %d %d %d\n",
536                   hp->count, hp->nwscore, hp->nwgaps,
537                   hp->accepted ? "accepted" : "not accepted",
538                   hp->nwalignment, hp->nwalignmentlength,
539                   hp->trim_q_left, hp->trim_q_right,
540                   hp->trim_t_left, hp->trim_t_right
541                   );
542 #endif
543         }
544     }
545   else if (opt_output_no_hits)
546     {
547       fprintf(fp, "\n");
548       fprintf(fp,"Query >%s\n", query_head);
549       fprintf(fp,"No hits\n");
550     }
551 }
552 
nucleotide_equal(char a,char b)553 bool inline nucleotide_equal(char a, char b)
554 {
555   return chrmap_4bit[(int)a] == chrmap_4bit[(int)b];
556 }
557 
build_sam_strings(char * alignment,char * queryseq,char * targetseq,xstring * cigar,xstring * md)558 void build_sam_strings(char * alignment,
559                        char * queryseq,
560                        char * targetseq,
561                        xstring * cigar,
562                        xstring * md)
563 {
564   /*
565     convert cigar to sam format:
566     add "1" to operations without run length
567     flip direction of indels in cigar string
568 
569     build MD-string with substitutions
570   */
571 
572   cigar->empty();
573   md->empty();
574 
575   char * p = alignment;
576   char * e = p + strlen(p);
577 
578   int qpos = 0;
579   int tpos = 0;
580 
581   int matched = 0;
582   bool flag = false; /* 1: MD string ends with a number */
583 
584   while(p < e)
585     {
586       int run = 1;
587       int scanned = 0;
588       sscanf(p, "%d%n", & run, & scanned);
589       p += scanned;
590       char op = *p++;
591 
592       switch (op)
593         {
594         case 'M':
595           cigar->add_d(run);
596           cigar->add_c('M');
597 
598           for(int i=0; i<run; i++)
599             {
600               if (nucleotide_equal(queryseq[qpos], targetseq[tpos])) {
601                 matched++;
602               } else
603                 {
604                   if (!flag)
605                     {
606                       md->add_d(matched);
607                       matched = 0;
608                       flag = true;
609                     }
610 
611                   md->add_c(targetseq[tpos]);
612                   flag = false;
613                 }
614               qpos++;
615               tpos++;
616             }
617 
618           break;
619 
620         case 'D':
621           cigar->add_d(run);
622           cigar->add_c('I');
623           qpos += run;
624           break;
625 
626         case 'I':
627           cigar->add_d(run);
628           cigar->add_c('D');
629 
630           if (!flag)
631             {
632               md->add_d(matched);
633               matched = 0;
634               flag = true;
635             }
636 
637           md->add_c('^');
638           for(int i=0; i<run; i++) {
639             md->add_c(targetseq[tpos++]);
640 
641         }
642           flag = false;
643           break;
644         }
645     }
646 
647   if (!flag)
648     {
649       md->add_d(matched);
650       matched = 0;
651       flag = true;
652     }
653 }
654 
results_show_samheader(FILE * fp,char * cmdline,char * dbname)655 void results_show_samheader(FILE * fp,
656                             char * cmdline,
657                             char * dbname)
658 {
659   if (opt_samout && opt_samheader)
660     {
661       fprintf(fp, "@HD\tVN:1.0\tSO:unsorted\tGO:query\n");
662 
663       for(uint64_t i=0; i<db_getsequencecount(); i++)
664         {
665           char md5hex[LEN_HEX_DIG_MD5];
666           get_hex_seq_digest_md5(md5hex,
667                                  db_getsequence(i),
668                                  db_getsequencelen(i));
669           fprintf(fp,
670                   "@SQ\tSN:%s\tLN:%" PRIu64 "\tM5:%s\tUR:file:%s\n",
671                   db_getheader(i),
672                   db_getsequencelen(i),
673                   md5hex,
674                   dbname);
675         }
676 
677       fprintf(fp,
678               "@PG\tID:%s\tVN:%s\tCL:%s\n",
679               PROG_NAME,
680               PROG_VERSION,
681               cmdline);
682     }
683 }
684 
results_show_samout(FILE * fp,struct hit * hits,int hitcount,char * query_head,char * qsequence,int64_t qseqlen,char * rc)685 void results_show_samout(FILE * fp,
686                          struct hit * hits,
687                          int hitcount,
688                          char * query_head,
689                          char * qsequence,
690                          int64_t qseqlen,
691                          char * rc)
692 {
693   /*
694      SAM format output
695 
696      http://samtools.github.io/hts-specs/SAMv1.pdf
697      http://www.drive5.com/usearch/manual/sam_files.html
698      http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml#sam-output
699      http://davetang.org/muse/2011/01/28/perl-and-sam/
700 
701       1: qname, query template name
702       2: flag, bitwise flag (12 bits)
703          (0x004=unmapped, 0x010=rev strand, 0x100 sec. alignment)
704       3: rname, reference sequence name
705       4: pos, 1-based leftmost mapping position (1)
706       5: mapq, mapping quality (255)
707       6: cigar, cigar string (MID)
708       7: rnext, ref name of next/paired read (*)
709       8: pnest, position of next/paired read (0)
710       9: tlen, obs template length (target length)
711      10: seq, segment of sequence
712      11: qual, ascii of phred based quality+33 (*)
713      12: optional tags (tag:type:value)
714 
715      Optional tags AS, XN, XM, XO, XG, NM, MD and YT used in usearch8.
716 
717      Usearch8:
718 
719      AS:i:? alignment score (i.e percent identity)
720      XN:i:? next best alignment score (always 0?)
721      XM:i:? number of mismatches
722      XO:i:? number of gap opens (excluding terminal gaps)
723      XG:i:? number of gap extensions (excluding terminal gaps)
724      NM:i:? edit distance (sum of XM and XG)
725      MD:Z:? variant string
726      YT:Z:UU string representing alignment type
727 
728   */
729 
730   if (hitcount > 0)
731     {
732       double top_hit_id = hits[0].id;
733 
734       for(int t = 0; t < hitcount; t++)
735         {
736           struct hit * hp = hits + t;
737 
738           if (opt_top_hits_only && (hp->id < top_hit_id)) {
739             break;
740 
741         }
742 
743           /*
744 
745           */
746 
747           xstring cigar;
748           xstring md;
749 
750           build_sam_strings(hp->nwalignment,
751                             hp->strand ? rc : qsequence,
752                             db_getsequence(hp->target),
753                             & cigar,
754                             & md);
755 
756           fprintf(fp,
757                   "%s\t%u\t%s\t%" PRIu64
758                   "\t%u\t%s\t%s\t%" PRIu64
759                   "\t%" PRIu64
760                   "\t%s\t%s\t"
761                   "AS:i:%.0f\tXN:i:%d\tXM:i:%d\tXO:i:%d\t"
762                   "XG:i:%d\tNM:i:%d\tMD:Z:%s\tYT:Z:%s\n",
763                   query_head,
764                   0x10 * hp->strand | (t>0 ? 0x100 : 0),
765                   db_getheader(hp->target),
766                   (uint64_t) 1,
767                   255,
768                   cigar.get_string(),
769                   "*",
770                   (uint64_t) 0,
771                   (uint64_t) 0,
772                   hp->strand ? rc : qsequence,
773                   "*",
774                   hp->id,
775                   0,
776                   hp->mismatches,
777                   hp->internal_gaps,
778                   hp->internal_indels,
779                   hp->mismatches + hp->internal_indels,
780                   md.get_string(),
781                   "UU");
782         }
783     }
784   else if (opt_output_no_hits)
785     {
786       fprintf(fp,
787               "%s\t%u\t%s\t%" PRIu64 "\t%u\t%s\t%s\t%" PRIu64 "\t%" PRIu64 "\t%s\t%s\n",
788               query_head,
789               0x04,
790               "*",
791               (uint64_t) 0,
792               255,
793               "*",
794               "*",
795               (uint64_t) 0,
796               (uint64_t) 0,
797               qsequence,
798               "*");
799     }
800 }
801