Consider the possibility of strand=".". Resolves #2158

Author: pd3

csq.c

@@ -965,7 +965,7 @@ int shifted_del_synonymous(args_t *args, splice_t *splice, uint32_t ex_beg, uint
     if ( tr->strand==STRAND_FWD && splice->vcf.pos >= ex_beg + 3 ) return 0;
 
 #if XDBG
-    fprintf(stderr,"shifted_del_synonymous: %d-%d  %s\n",ex_beg,ex_end, tr->strand==STRAND_FWD?"fwd":"rev");
+    fprintf(stderr,"shifted_del_synonymous: %d-%d  %s\n",ex_beg,ex_end, tr->strand==STRAND_FWD?"fwd":(tr->strand==STRAND_REV?"rev":"unk"));
     fprintf(stderr,"   %d  ..  %s > %s\n",splice->vcf.pos+1,splice->vcf.ref,splice->vcf.alt);
 #endif
 
@@ -998,7 +998,7 @@ int shifted_del_synonymous(args_t *args, splice_t *splice, uint32_t ex_beg, uint
         while ( ptr_vcf[i] && ptr_vcf[i]==ptr_ref[i] ) i++;
         if ( ptr_vcf[i] ) return 0;       // the deleted sequence cannot be replaced
     }
-    else
+    else if ( tr->strand==STRAND_FWD )
     {
         // STRAND_FWD
         int32_t vcf_block_beg = splice->vcf.pos + ref_len - 2*ndel;        // the position of the first base of the ref block that could potentially replace the deletion
@@ -1271,13 +1271,13 @@ fprintf(stderr,"mnp: %s>%s .. ex=%d,%d  beg,end=%d,%d  tbeg,tend=%d,%d  check_ut
     {
         if ( splice->check_region_beg ) splice->csq |= CSQ_SPLICE_REGION;
         if ( splice->tr->strand==STRAND_FWD ) { if ( splice->check_start ) splice->csq |= CSQ_START_LOST; }
-        else { if ( splice->check_stop ) splice->csq |= CSQ_STOP_LOST; }
+        else if ( splice->tr->strand==STRAND_REV ) { if ( splice->check_stop ) splice->csq |= CSQ_STOP_LOST; }
     }
     if ( splice->ref_end > ex_end - 3 )
     {
         if ( splice->check_region_end ) splice->csq |= CSQ_SPLICE_REGION;
         if ( splice->tr->strand==STRAND_REV ) { if ( splice->check_start ) splice->csq |= CSQ_START_LOST; }
-        else { if ( splice->check_stop ) splice->csq |= CSQ_STOP_LOST; }
+        else if ( splice->tr->strand==STRAND_FWD ) { if ( splice->check_stop ) splice->csq |= CSQ_STOP_LOST; }
     }
     if ( splice->set_refalt )
     {
@@ -1338,17 +1338,17 @@ int hap_init(args_t *args, hap_node_t *parent, hap_node_t *child, gf_cds_t *cds,
     if ( !(tr->trim & TRIM_5PRIME) )
     {
         if ( tr->strand==STRAND_FWD ) { if ( child->icds==0 ) splice.check_start = 1; }
-        else { if ( child->icds==tr->ncds-1 ) splice.check_start = 1; }
+        else if ( tr->strand==STRAND_REV ) { if ( child->icds==tr->ncds-1 ) splice.check_start = 1; }
     }
     if ( !(tr->trim & TRIM_3PRIME) )
     {
         if ( tr->strand==STRAND_FWD ) { if ( child->icds==tr->ncds-1 ) splice.check_stop = 1; }
-        else { if ( child->icds==0 ) splice.check_stop = 1; }
+        else if ( tr->strand==STRAND_REV ) { if ( child->icds==0 ) splice.check_stop = 1; }
     }
     if ( splice.check_start )   // do not check starts in incomplete CDS, defined as not starting with M
     {
         if ( tr->strand==STRAND_FWD ) { if ( dna2aa(TSCRIPT_AUX(tr)->ref+N_REF_PAD+cds->beg-tr->beg) != 'M' ) splice.check_start = 0; }
-        else { if ( cdna2aa(TSCRIPT_AUX(tr)->ref+N_REF_PAD+cds->beg-tr->beg+cds->len-3) != 'M' ) splice.check_start = 0; }
+        else if ( tr->strand==STRAND_REV ) { if ( cdna2aa(TSCRIPT_AUX(tr)->ref+N_REF_PAD+cds->beg-tr->beg+cds->len-3) != 'M' ) splice.check_start = 0; }
     }
     if ( child->icds!=0 ) splice.check_region_beg = 1;
     if ( child->icds!=tr->ncds-1 ) splice.check_region_end = 1;
@@ -1586,7 +1586,7 @@ fprintf(stderr,"\ntranslate: %d %d %d  fill=%d  seq.l=%d\n",sbeg,rbeg,rend,fill,
             }
         }
     }
-    else    // STRAND_REV
+    else if ( strand==STRAND_REV )
     {
         // right padding - number of bases to take from ref
         npad = (seq.m - (sbeg + seq.l)) % 3;
@@ -1673,6 +1673,7 @@ fprintf(stderr,"\ntranslate: %d %d %d  fill=%d  seq.l=%d\n",sbeg,rbeg,rend,fill,
             }
         }
     }
+    else error("Should not happen: %d\n", strand);
     kputc_(0,tseq); tseq->l--;
 #if DBG
  fprintf(stderr,"    tseq: %s\n", tseq->s);
@@ -1858,7 +1859,7 @@ void kput_vcsq(args_t *args, vcsq_t *csq, kstring_t *str)
     kputs(gf_type2gff_string(csq->biotype), str);
 
     if ( CSQ_PRN_STRAND(csq->type) || csq->vstr.l )
-        kputs(csq->strand==STRAND_FWD ? "|+" : "|-", str);
+        kputs(csq->strand==STRAND_FWD ? "|+" : (csq->strand==STRAND_REV ? "|-" : "|."), str);
 
     if ( csq->vstr.l )
         kputs(csq->vstr.s, str);
@@ -1882,6 +1883,7 @@ void hap_add_csq(args_t *args, hap_t *hap, hap_node_t *node, int tlen, int ibeg,
 {
     int i;
     gf_tscript_t *tr = hap->tr;
+    assert( tr->strand==STRAND_FWD || tr->strand==STRAND_REV );
     int ref_node = tr->strand==STRAND_FWD ? ibeg : iend;
     int icsq = node->ncsq_list++;
     hts_expand0(csq_t,node->ncsq_list,node->mcsq_list,node->csq_list);
@@ -2177,7 +2179,7 @@ void hap_finalize(args_t *args, hap_t *hap)
                 indel = 0;
             }
         }
-        else
+        else if ( tr->strand==STRAND_REV )
         {
             i = istack + 1, ibeg = -1;
             while ( --i > 0 )

gff.c

@@ -53,9 +53,9 @@ typedef struct
     uint32_t beg;
     uint32_t end;
     uint32_t trid;
-    uint32_t strand:1;  // STRAND_REV,STRAND_FWD
+    uint32_t strand:2;  // STRAND_{REV,FWD,UNK}
     uint32_t phase:2;   // 0, 1, 2, or 3 for unknown
-    uint32_t iseq:29;
+    uint32_t iseq:28;
 }
 ftr_t;
 
@@ -474,7 +474,7 @@ static void gff_parse_exon(gff_t *gff, const char *line, ftr_t *ftr)
     // associate with transcript id
     gff_id_register(&gff->tscript_ids, aux->parent, aux->parent_end, &ftr->trid);
 
-    if ( ftr->strand==-1 && gff->verbosity > 0 )
+    if ( ftr->strand==STRAND_UNK && gff->verbosity > 0 )
     {
         if ( !gff->warned.unknown_strand || gff->verbosity > 1 )
             fprintf(stderr,"Warning: Ignoring GFF feature with unknown strand .. %s\n",line);
@@ -568,6 +568,7 @@ static int gff_parse_line(gff_t *gff, char *line, ftr_t *ftr)
     ftr->strand = -1;
     if ( *ss == '+' ) ftr->strand = STRAND_FWD;
     else if ( *ss == '-' ) ftr->strand = STRAND_REV;
+    else ftr->strand = STRAND_UNK;
     ss += 2;
 
     // 8th column: phase (codon offset)
@@ -757,7 +758,7 @@ static void tscript_init_cds(gff_t *gff)
             }
             if ( !tscript_ok ) continue;    // skip this transcript
         }
-        else
+        else if ( tr->strand==STRAND_REV )
         {
             if ( tr->cds[tr->ncds-1]->phase != CDS_PHASE_UNKN )
             {
@@ -820,6 +821,8 @@ static void tscript_init_cds(gff_t *gff)
             }
             if ( !tscript_ok ) continue;    // skip this transcript
         }
+        else
+            continue;   // unknown strand
 
         // set len. At the same check that CDS within a transcript do not overlap
         len = 0;
@@ -868,7 +871,7 @@ static void tscript_init_cds(gff_t *gff)
                     i--;
                 }
             }
-            else
+            else if ( tr->strand==STRAND_REV )
             {
                 i = 0;
                 while ( i<tr->ncds && len%3 )
@@ -910,7 +913,7 @@ static int gff_dump(gff_t *gff, const char *fname)
         gf_gene_t *gene = (gf_gene_t*) kh_val(gff->init.gid2gene, k);
         char *gene_id = gff->init.gene_ids.str[gene->id];
         str.l = 0;
-        ksprintf(&str,"%s\t.\tgene\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tID=%s;Name=%s;used=%d\n",gff->init.seq[gene->iseq],gene->beg+1,gene->end+1,gene->strand==STRAND_FWD?'+':'-',gene_id,gene->name,gene->used);
+        ksprintf(&str,"%s\t.\tgene\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tID=%s;Name=%s;used=%d\n",gff->init.seq[gene->iseq],gene->beg+1,gene->end+1,gene->strand==STRAND_FWD?'+':(gene->strand==STRAND_REV?'-':'.'),gene_id,gene->name,gene->used);
         if ( bgzf_write(out, str.s, str.l) != str.l ) error("Error writing %s: %s\n", fname, strerror(errno));
     }
 
@@ -921,7 +924,7 @@ static int gff_dump(gff_t *gff, const char *fname)
         char *gene_id =  gff->init.gene_ids.str[tr->gene->id];
         const char *type = tr->type==GF_PROTEIN_CODING ? "mRNA" : gf_type2gff_string(tr->type);
         str.l = 0;
-        ksprintf(&str,"%s\t.\t%s\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tID=%s;Parent=%s;biotype=%s;used=%d\n",itr->seq,type,itr->beg+1,itr->end+1,tr->strand==STRAND_FWD?'+':'-',gff->tscript_ids.str[tr->id],gene_id,gf_type2gff_string(tr->type),tr->used);
+        ksprintf(&str,"%s\t.\t%s\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tID=%s;Parent=%s;biotype=%s;used=%d\n",itr->seq,type,itr->beg+1,itr->end+1,tr->strand==STRAND_FWD?'+':(tr->strand==STRAND_REV?'-':'.'),gff->tscript_ids.str[tr->id],gene_id,gf_type2gff_string(tr->type),tr->used);
         if ( bgzf_write(out, str.s, str.l) != str.l ) error("Error writing %s: %s\n", fname, strerror(errno));
     }
     regitr_destroy(itr);
@@ -932,7 +935,7 @@ static int gff_dump(gff_t *gff, const char *fname)
         gf_cds_t *cds = regitr_payload(itr,gf_cds_t*);
         gf_tscript_t *tr = cds->tr;
         str.l = 0;
-        ksprintf(&str,"%s\t.\tCDS\t%"PRIu32"\t%"PRIu32"\t.\t%c\t%c\tParent=%s\n",itr->seq,cds->beg+1,cds->beg+cds->len,tr->strand==STRAND_FWD?'+':'-',cds->phase==3?'.':cds->phase+(int)'0',gff->tscript_ids.str[tr->id]);
+        ksprintf(&str,"%s\t.\tCDS\t%"PRIu32"\t%"PRIu32"\t.\t%c\t%c\tParent=%s\n",itr->seq,cds->beg+1,cds->beg+cds->len,tr->strand==STRAND_FWD?'+':(tr->strand==STRAND_REV?'-':'.'),cds->phase==3?'.':cds->phase+(int)'0',gff->tscript_ids.str[tr->id]);
         if ( bgzf_write(out, str.s, str.l) != str.l ) error("Error writing %s: %s\n", fname, strerror(errno));
     }
     regitr_destroy(itr);
@@ -943,7 +946,7 @@ static int gff_dump(gff_t *gff, const char *fname)
         gf_utr_t *utr = regitr_payload(itr,gf_utr_t*);
         gf_tscript_t *tr = utr->tr;
         str.l = 0;
-        ksprintf(&str,"%s\t.\t%s_prime_UTR\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tParent=%s\n",itr->seq,utr->which==prime3?"three":"five",utr->beg+1,utr->end+1,tr->strand==STRAND_FWD?'+':'-',gff->tscript_ids.str[tr->id]);
+        ksprintf(&str,"%s\t.\t%s_prime_UTR\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tParent=%s\n",itr->seq,utr->which==prime3?"three":"five",utr->beg+1,utr->end+1,tr->strand==STRAND_FWD?'+':(tr->strand==STRAND_REV?'-':'.'),gff->tscript_ids.str[tr->id]);
         if ( bgzf_write(out, str.s, str.l) != str.l ) error("Error writing %s: %s\n", fname, strerror(errno));
     }
     regitr_destroy(itr);
@@ -954,7 +957,7 @@ static int gff_dump(gff_t *gff, const char *fname)
         gf_exon_t *exon = regitr_payload(itr,gf_exon_t*);
         gf_tscript_t *tr = exon->tr;
         str.l = 0;
-        ksprintf(&str,"%s\t.\texon\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tParent=%s\n",itr->seq,exon->beg+1,exon->end+1,tr->strand==STRAND_FWD?'+':'-',gff->tscript_ids.str[tr->id]);
+        ksprintf(&str,"%s\t.\texon\t%"PRIu32"\t%"PRIu32"\t.\t%c\t.\tParent=%s\n",itr->seq,exon->beg+1,exon->end+1,tr->strand==STRAND_FWD?'+':(tr->strand==STRAND_REV?'-':'.'),gff->tscript_ids.str[tr->id]);
         if ( bgzf_write(out, str.s, str.l) != str.l ) error("Error writing %s: %s\n", fname, strerror(errno));
     }
     regitr_destroy(itr);

gff.h

@@ -150,6 +150,7 @@
 
 #define STRAND_REV 0
 #define STRAND_FWD 1
+#define STRAND_UNK 2
 
 #define TRIM_NONE   0
 #define TRIM_5PRIME 1
@@ -273,9 +274,9 @@ struct gf_tscript_t_
 {
     uint32_t id;        // transcript id
     uint32_t beg,end;   // transcript's beg and end coordinate (ref strand, 0-based, inclusive)
-    uint32_t strand:1,  // STRAND_REV or STRAND_FWD
+    uint32_t strand:2,  // STRAND_REV,FWD,UNK
              used:1,    // does it have any exons, UTRs, CDS?
-             ncds:30,   // number of exons
+             ncds:29,   // number of exons
              mcds;
     gf_cds_t **cds;     // ordered list of exons
     uint32_t trim:2,    // complete, 5' or 3' trimmed, see TRIM_* types