c864e74a423acbe711e6be084573aa11b687c1ce kent Tue Aug 6 17:47:29 2013 -0700 Making it work ok with /dev/null as fastq output. Fixing an integer overflow bug. diff --git src/utils/fastqStatsAndSubsample/fastqStatsAndSubsample.c src/utils/fastqStatsAndSubsample/fastqStatsAndSubsample.c index 4b23fdd..049919e 100644 --- src/utils/fastqStatsAndSubsample/fastqStatsAndSubsample.c +++ src/utils/fastqStatsAndSubsample/fastqStatsAndSubsample.c @@ -131,59 +131,63 @@ for (i=0; i<arraySize; ++i) total += array[i]; return total; } long long sumIntArray(int *array, int arraySize) /* Return sum of all items in array */ { long long total = 0; int i; for (i=0; i<arraySize; ++i) total += array[i]; return total; } -void printAveDoubleArray(FILE *f, char *label, double *a, int *totalAtPos, int aSize) +void printAveDoubleArray(FILE *f, char *label, double *a, long long *totalAtPos, int aSize) /* Print a[i]/counts[i] for all elements in array */ { fprintf(f, "%s ", label); int i; for (i=0; i<aSize; ++i) fprintf(f, "%g,", a[i]/totalAtPos[i]); fprintf(f, "\n"); } -void printAveIntArray(FILE *f, char *label, int *a, int *totalAtPos, int aSize) +void printAveIntArray(FILE *f, char *label, int *a, long long *totalAtPos, int aSize) /* Print a[i]/totalAtPos[i] for all elements in array */ { fprintf(f, "%s ", label); int i; for (i=0; i<aSize; ++i) fprintf(f, "%g,", ((double)a[i])/totalAtPos[i]); fprintf(f, "\n"); } boolean oneFastqRecord(struct lineFile *lf, FILE *f, boolean copy, boolean firstTime) /* Read next fastq record from LF, and optionally copy it to f. Return FALSE at end of file * Do a _little_ error checking on record while we're at it. The format has already been * validated on the client side fairly thoroughly. */ { char *line; int lineSize; +/* Treat NULL file same as non-copy, so only have one condition to check on . */ +if (f == NULL) + copy = FALSE; + /* Deal with initial line starting with '@' */ if (!lineFileNext(lf, &line, &lineSize)) return FALSE; if (line[0] != '@') errAbort("Expecting line starting with '@' got %s line %d of %s", line, lf->lineIx, lf->fileName); if (copy) mustWrite(f, line, lineSize); /* Deal with line containing sequence. */ if (!lineFileNext(lf, &line, &lineSize)) errAbort("%s truncated in middle of record", lf->fileName); /* Get size and add it to stats */ int seqSize = lineSize-1; @@ -344,71 +348,83 @@ boolean doIt = randomizer[i]; if (!maybeCopyFastqRecord(lf, f, doIt, &seqSize)) internalErr(); if (doIt) basesInSample += seqSize; } freez(&randomizer); lineFileClose(&lf); return basesInSample; } void fastqStatsAndSubsample(char *inFastq, char *outStats, char *outFastq) /* fastqStatsAndSubsample - Go through a fastq file doing sanity checks and collecting * statistics, and also producing a smaller fastq out of a sample of the data. */ { +/* Temporary file if any */ +FILE *smallF = NULL; + +/* Make this work without making input. */ +if (sameString("/dev/null", outFastq)) + outFastq = NULL; + + +/* Open up temp output file. This one will be for the initial scaling. We'll do + * a second round of scaling as well. */ +char smallFastqName[PATH_LEN] = ""; +if (outFastq != NULL) + { /* Split up outFastq path, so we can make a temp file in the same dir. */ char outDir[PATH_LEN]; + if (outFastq) splitPath(outFastq, outDir, NULL, NULL); -/* Open up temp output file. This one will be for the initial scaling. We'll do - * a second round of scaling as well. */ -char smallFastqName[PATH_LEN]; safef(smallFastqName, PATH_LEN, "%sfastqSubsampleXXXXXX", outDir); int smallFd = mkstemp(smallFastqName); -FILE *smallF = fdopen(smallFd, "w"); + smallF = fdopen(smallFd, "w"); + } /* Scan through input, collecting stats, validating, and creating a subset file. */ int downStep = calcInitialReduction(inFastq, sampleSize); struct lineFile *lf = lineFileOpen(inFastq, FALSE); boolean done = FALSE; int readsCopied = 0, totalReads = 0; long long basesInSample = 0; boolean firstTime = TRUE; while (!done) { int hotPosInCycle = rand()%downStep; int cycle; for (cycle=0; cycle<downStep; ++cycle) { boolean hotPos = (cycle == hotPosInCycle); if (!oneFastqRecord(lf, smallF, hotPos, firstTime)) { done = TRUE; break; } if (hotPos) ++readsCopied; firstTime = FALSE; ++totalReads; } } lineFileClose(&lf); carefulClose(&smallF); boolean justUseSmall = FALSE; -if (readsCopied < sampleSize) +if (outFastq != NULL && readsCopied < sampleSize) { /* Our sample isn't big enough. This could have two causes - a bug in * our estimation, maybe caused by read sizes growing a bunch in the time * since this code was written - or a file that is actually smaller smaller * than sampleSize. */ if (sampleSize <= totalReads) { remove(smallFastqName); errAbort("Internal error: bad estimate %d for downStep on %s", downStep, inFastq); } else { if (smallOk) { justUseSmall = TRUE; @@ -420,43 +436,46 @@ remove(smallFastqName); errAbort("SampleSize is set to %d reads, but there are only %d reads in %s", sampleSize, totalReads, inFastq); } } } char *qualType = "solexa"; int qualZero = 64; if (minQual <= 58) { qualType = "sanger"; qualZero = 33; } +if (outFastq != NULL) + { if (justUseSmall) { mustRename(smallFastqName, outFastq); sampleSize = readsCopied; basesInSample= sumReadBases; } else { FILE *f = mustOpen(outFastq, "w"); basesInSample = reduceFastqSample(smallFastqName, f, readsCopied, sampleSize); carefulClose(&f); remove(smallFastqName); } + } FILE *f = mustOpen(outStats, "w"); int posCount = maxReadBases; fprintf(f, "readCount %d\n", totalReads); fprintf(f, "baseCount %lld\n", sumReadBases); fprintf(f, "sampleCount %d\n", sampleSize); fprintf(f, "basesInSample %lld\n", basesInSample); fprintf(f, "readSizeMean %g\n", (double)sumReadBases/totalReads); if (minReadBases != maxReadBases) fprintf(f, "readSizeStd %g\n", calcStdFromSums(sumReadBases, sumSquaredReadBases, totalReads)); else fprintf(f, "readSizeStd 0\n"); fprintf(f, "readSizeMin %d\n", minReadBases); fprintf(f, "readSizeMax %d\n", maxReadBases); double qSum = sumDoubleArray(sumQuals, maxReadBases); @@ -474,38 +493,40 @@ /* Compute overall total nucleotide stats from count arrays. */ long long aSum = sumIntArray(aCount, maxReadBases); long long cSum = sumIntArray(cCount, maxReadBases); long long gSum = sumIntArray(gCount, maxReadBases); long long tSum = sumIntArray(tCount, maxReadBases); long long nSum = sumIntArray(nCount, maxReadBases); fprintf(f, "atRatio %g\n", (double)(aSum + tSum)/(aSum + cSum + gSum + tSum)); fprintf(f, "aRatio %g\n", (double)aSum/sumReadBases); fprintf(f, "cRatio %g\n", (double)cSum/sumReadBases); fprintf(f, "gRatio %g\n", (double)gSum/sumReadBases); fprintf(f, "tRatio %g\n", (double)tSum/sumReadBases); fprintf(f, "nRatio %g\n", (double)nSum/sumReadBases); /* Now deal with array outputs. First make up count of all bases we've seen. */ fprintf(f, "posCount %d\n", posCount); -int totalAtPos[posCount]; +long long totalAtPos[posCount]; int pos; for (pos=0; pos<posCount; ++pos) totalAtPos[pos] = aCount[pos] + cCount[pos] + gCount[pos] + tCount[pos] + nCount[pos]; /* Offset quality by scale */ for (pos=0; pos<posCount; ++pos) sumQuals[pos] -= totalAtPos[pos] * qualZero; +#ifdef SOON +#endif /* SOON */ printAveDoubleArray(f, "qualPos", sumQuals, totalAtPos, posCount); printAveIntArray(f, "aAtPos", aCount, totalAtPos, posCount); printAveIntArray(f, "cAtPos", cCount, totalAtPos, posCount); printAveIntArray(f, "gAtPos", gCount, totalAtPos, posCount); printAveIntArray(f, "tAtPos", tCount, totalAtPos, posCount); printAveIntArray(f, "nAtPos", nCount, totalAtPos, posCount); } int main(int argc, char *argv[]) /* Process command line. */ { optionInit(&argc, argv, options); if (argc != 4)