620a9846b9df834954237a587f149d26c3b6dff5 braney Fri Nov 6 08:49:58 2015 -0800 sometimes Complete Genomics bams have records with no alignment in them. Ignore them #16339 diff --git src/hg/hgTracks/bamTrack.c src/hg/hgTracks/bamTrack.c index a00e929..b00d2d7 100644 --- src/hg/hgTracks/bamTrack.c +++ src/hg/hgTracks/bamTrack.c @@ -26,31 +26,31 @@ struct bamTrackData { struct track *tg; struct hash *pairHash; int minAliQual; char *colorMode; char *grayMode; char *userTag; int aliQualShadeMin; int aliQualShadeMax; int baseQualShadeMin; int baseQualShadeMax; }; -struct psl *pslFromBam(const bam1_t *bam) +static struct psl *pslFromBam(const bam1_t *bam) /* Translate BAM's numeric CIGAR encoding into PSL sufficient for cds.c (just coords, * no scoring info) */ { const bam1_core_t *core = &bam->core; struct psl *psl; AllocVar(psl); boolean isRc = (core->flag & BAM_FREVERSE); psl->strand[0] = isRc ? '-' : '+'; psl->qName = cloneString(bam1_qname(bam)); psl->tName = cloneString(chromName); unsigned blockCount = 0; unsigned *blockSizes, *qStarts, *tStarts; AllocArray(blockSizes, core->n_cigar); AllocArray(qStarts, core->n_cigar); AllocArray(tStarts, core->n_cigar); @@ -81,30 +81,37 @@ case 'D': // deleted from query case 'N': // long deletion from query (intron as opposed to small del) tPos += n; break; case 'S': // skipped query bases at beginning or end ("soft clipping") qPos += n; qLength += n; break; case 'H': // skipped query bases not stored in record's query sequence ("hard clipping") case 'P': // P="silent deletion from padded reference sequence" -- ignore these. break; default: errAbort("pslFromBam: unrecognized CIGAR op %c -- update me", op); } } + +if (blockCount == 0) + { + // sometimes BAM's have alignments with no alignment + return NULL; // leaks allocated PSL. + } + psl->tSize = hChromSize(database, chromName); psl->tStart = tStarts[0]; psl->tEnd = tStarts[blockCount-1] + blockSizes[blockCount-1]; psl->qSize = qLength; psl->qStart = qStarts[0]; psl->qEnd = qStarts[blockCount-1] + blockSizes[blockCount-1]; if (isRc) reverseIntRange(&psl->qStart, &psl->qEnd, psl->qSize); psl->blockCount = blockCount; psl->blockSizes = blockSizes; psl->qStarts = qStarts; psl->tStarts = tStarts; return psl; } @@ -207,40 +214,44 @@ int wordCount = chopString(s, ",", row, ArraySize(row)); if ((wordCount != 3) || (!isdigit(row[0][0]) || !isdigit(row[1][0]) || !isdigit(row[2][0]))) return FALSE; *retR = atoi(row[0]); *retG = atoi(row[1]); *retB = atoi(row[2]); return TRUE; } struct linkedFeatures *bamToLf(const bam1_t *bam, void *data) /* Translate a BAM record into a linkedFeatures item. */ { struct bamTrackData *btd = (struct bamTrackData *)data; const bam1_core_t *core = &bam->core; struct linkedFeatures *lf; +struct psl *original = pslFromBam(bam); +if (original == NULL) + return NULL; + AllocVar(lf); lf->score = core->qual; lf->name = cloneString(bam1_qname(bam)); lf->orientation = (core->flag & BAM_FREVERSE) ? -1 : 1; int length; lf->components = sfFromNumericCigar(bam, &length); lf->start = lf->tallStart = core->pos; lf->end = lf->tallEnd = core->pos + length; lf->extra = bamGetQuerySequence(bam, FALSE); // cds.c reverses if psl != NULL -lf->original = pslFromBam(bam); +lf->original = original; int clippedQLen; bamGetSoftClipping(bam, NULL, NULL, &clippedQLen); if (sameString(btd->colorMode, BAM_COLOR_MODE_GRAY) && sameString(btd->grayMode, BAM_GRAY_MODE_ALI_QUAL)) { lf->grayIx = shadeTransform(btd->aliQualShadeMin, btd->aliQualShadeMax, core->qual); } else if (sameString(btd->colorMode, BAM_COLOR_MODE_GRAY) && sameString(btd->grayMode, BAM_GRAY_MODE_BASE_QUAL)) { UBYTE *quals = bamGetQueryQuals(bam, TRUE); lf->components = expandSfQuals(lf->components, quals, lf->orientation, clippedQLen, btd->baseQualShadeMin, btd->baseQualShadeMax); lf->grayIx = maxShade - 3; } @@ -286,32 +297,35 @@ if (core->flag & BAM_FUNMAP) return FALSE; if (core->qual < btd->minAliQual) return FALSE; return TRUE; } int addBam(const bam1_t *bam, void *data) /* bam_fetch() calls this on each bam alignment retrieved. Translate each bam * into a linkedFeatures item, and add it to tg->items. */ { struct bamTrackData *btd = (struct bamTrackData *)data; if (!passesFilters(bam, btd)) return 0; struct linkedFeatures *lf = bamToLf(bam, data); +if (lf) + { struct track *tg = btd->tg; slAddHead(&(tg->items), lf); + } return 0; } static struct linkedFeatures *lfStub(int startEnd, int orientation) /* Make a linkedFeatures for a zero-length item (so that an arrow will be drawn * toward the given coord. */ { struct linkedFeatures *lf; AllocVar(lf); lf->name = cloneString("stub"); lf->orientation = orientation; struct simpleFeature *sf; AllocVar(sf); sf->start = sf->end = lf->start = lf->end = lf->tallStart = lf->tallEnd = startEnd; lf->components = sf;