4898794edd81be5285ea6e544acbedeaeb31bf78 max Tue Nov 23 08:10:57 2021 -0800 Fixing pointers to README file for license in all source code files. refs #27614 diff --git src/hg/snp/snpValid/snpValid.c src/hg/snp/snpValid/snpValid.c index be4af3d..fb586db 100644 --- src/hg/snp/snpValid/snpValid.c +++ src/hg/snp/snpValid/snpValid.c @@ -1,920 +1,920 @@ /* snpValid - Test validity of snpMap, dbSnpRsXx, and affy120KDetails, affy10KDetails. */ /* Copyright (C) 2013 The Regents of the University of California - * See README in this or parent directory for licensing information. */ + * See kent/LICENSE or http://genome.ucsc.edu/license/ for licensing information. */ #include "common.h" #include "memalloc.h" #include "linefile.h" #include "hash.h" #include "portable.h" #include "options.h" #include "errCatch.h" #include "ra.h" #include "fa.h" #include "nib.h" #include "hdb.h" /* we have defined it our self, use only needed fields, so we don't need this. #include "snp.h" */ #include "axt.h" #include <time.h> /* Command line variables. */ char *db = NULL; /* DB from command line */ char *flankPath = NULL; /* path to flank files *.seq.gz */ char *chr = NULL; /* chr name from command line e.g "chr1" */ int threshold = 70; /* minimum score for match */ int Verbose = 0; /* set to 1 to see more detailed output (used upper case because var name collision in namespace) */ int maxFlank = 25; /* maximum length of each flanking sequence, left and right. */ void usage() /* Explain usage and exit. */ { errAbort( "snpValid - Validate snp alignments\n" "usage:\n" " snpValid db flankPath\n" " db is the database, e.g. hg16.\n" " flankPath is the path to the flank files *.seq.gz.\n" "options:\n" " -chr=name - Use chrom name (e.g chr1) to limit validation, default all chromosomes.\n" " -threshold=number - Use threshold as minimum score 0-99, default %u.\n" " -verbose=number - Use verbose=1 to see all output mismatches, default %u.\n" " -maxFlank=number - Use maxFlank as maximum length of each flank, default %u.\n" , threshold, Verbose, maxFlank); } static struct optionSpec options[] = { { "chr" , OPTION_STRING }, { "threshold", OPTION_INT }, { "verbose" , OPTION_INT }, { "maxFlank" , OPTION_INT }, { NULL , 0 }, }; /* ========================================================================== */ #define BASE_SNP_EX_NUM 21 /* next avail number we can use */ enum snpExceptionType { snpExNotExact, snpExMismatch, snpExNoFlanks, snpExWrongStrand, snpExCount }; static char *snpExDesc[] = { "Not an exact match for locType=exact and size=1", "Mismatch below threshold.", "No flank data.", "Wrong Strand reported." }; FILE *exf[snpExCount]; void openExOuts(char *db) { int i; char fname[256]; for (i=0;i<snpExCount;i++) { safef(fname,sizeof(fname),"%ssnpException.%d.bed",db,i+BASE_SNP_EX_NUM); exf[i]=mustOpen(fname,"w"); fprintf(exf[i],"# exceptionId: %d\n",i+BASE_SNP_EX_NUM); fprintf(exf[i],"# query: %s\n",snpExDesc[i]); } } void closeExOuts() { int i; for (i=0;i<snpExCount;i++) { carefulClose(&exf[i]); } } /* --- save memory by defining just the fields needed from flank file ---- */ struct flank /* Information from snp */ { struct flank *next; /* Next in singly linked list. */ char *rsId; /* reference snp (rs) identifier */ char *leftFlank; /* leftFlank */ char *rightFlank; /* rightFlank */ char *observed; /* observed */ }; int slFlankCmp(const void *va, const void *vb) /* Compare two slNames. */ { const struct flank *a = *((struct flank **)va); const struct flank *b = *((struct flank **)vb); return strcmp(a->rsId, b->rsId); } struct flank *flankLoad(char **row) /* Load a flank from row fetched from linefile * Dispose of this with flankFree(). */ { struct flank *ret; AllocVar(ret); ret->rsId = cloneString(row[0]); ret->leftFlank = cloneString(row[1]); ret->rightFlank = cloneString(row[2]); ret->observed = cloneString(row[3]); return ret; } void flankFree(struct flank **pEl) /* Free a single dynamically allocated flank such as created * with flankLoad(). */ { struct flank *el; if ((el = *pEl) == NULL) return; freeMem(el->rsId); freeMem(el->leftFlank); freeMem(el->rightFlank); freeMem(el->observed); freez(pEl); } void flankFreeList(struct flank **pList) /* Free a list of dynamically allocated snp's */ { struct flank *el, *next; for (el = *pList; el != NULL; el = next) { next = el->next; flankFree(&el); } *pList = NULL; } struct flank *readFlank(char *chrom) /* Slurp in the flank rows for one chrom */ { char fileName[256]=""; struct lineFile *lf=NULL; struct flank *list=NULL, *el; char *row[4]; int rowCount=4; char *ch = &chrom[3]; /* example flankPath: /cluster/bluearc/snp/hg16/build122/seq/ds_ch%s.xml.contig.seq.gz */ safef(fileName,sizeof(fileName),"%s/ds_ch%s.xml.contig.seq.gz",flankPath,ch); lf=lineFileMayOpen(fileName, TRUE); if (lf == NULL) { safef(fileName,sizeof(fileName),"%s/ds_ch%s.xml.contig.seq",flankPath,ch); lf=lineFileMayOpen(fileName, TRUE); if (lf == NULL) { return NULL; } } while (lineFileNextRowTab(lf, row, rowCount)) { el = flankLoad(row); slAddHead(&list,el); } slReverse(&list); /* could possibly skip if it made much difference in speed. */ lineFileClose(&lf); return list; } /* --- save memory by defining just the fields needed from snp ---- */ struct snp /* Information from snp */ { struct snp *next; /* Next in singly linked list. */ char *chrom; /* chrom */ int chromStart; /* start */ int chromEnd; /* end */ char *name; /* name */ char *strand; /* strand */ char *observed; /* observed variants (usually slash-separated list) */ char *locType; /* location Type */ }; struct snp *snpLoad(char **row) /* Load a snp from row fetched with select * from snp * from database. Dispose of this with snpFree(). */ { struct snp *ret; AllocVar(ret); ret->chrom = cloneString(row[0]); ret->chromStart = atoi(row[1]); ret->chromEnd = atoi(row[2]); ret->name = cloneString(row[3]); ret->strand = cloneString(row[4]); ret->observed = cloneString(row[5]); ret->locType = cloneString(row[6]); return ret; } void snpFree(struct snp **pEl) /* Free a single dynamically allocated snp such as created * with snpLoad(). */ { struct snp *el; if ((el = *pEl) == NULL) return; freeMem(el->chrom); freeMem(el->name); freeMem(el->strand); freeMem(el->observed); freeMem(el->locType); freez(pEl); } void snpFreeList(struct snp **pList) /* Free a list of dynamically allocated snp's */ { struct snp *el, *next; for (el = *pList; el != NULL; el = next) { next = el->next; snpFree(&el); } *pList = NULL; } struct snp *readSnp(char *chrom) /* Slurp in the snp rows for one chrom */ { struct snp *list=NULL, *el; char query[512]; struct sqlConnection *conn = hAllocConn(db); struct sqlResult *sr; char **row; sqlSafef(query, sizeof(query), "select chrom, chromStart, chromEnd, name, strand, observed, locType " "from snp where chrom='%s' order by name", chrom); sr = sqlGetResult(conn, query); while ((row = sqlNextRow(sr)) != NULL) { el = snpLoad(row); slAddHead(&list,el); } sqlFreeResult(&sr); hFreeConn(&conn); slReverse(&list); /* could possibly skip if it made much difference in speed. */ return list; } /* --- ----------------------------------------------------------- ---- */ int mystrcmp(char *s1, char *s2) /* case insensitive compare */ { if ((*s1!=0) && (toupper(*s1)==toupper(*s2))) while ((*s1!=0) && (toupper(*++s1)==toupper(*++s2))) ; return (toupper(*s1) - toupper(*s2)); } boolean allNs(char *s1) /* check for all Ns */ { if ((*s1!=0) && (toupper(*s1)=='N')) while ((*s1!=0) && (toupper(*++s1)=='N')) ; return (*s1 == 0); } void stripDashes(char *s1) /* remove all Dashes from string inplace */ { char *s2 = s1; while (*s1!=0) { if (*s1 != '-') { *s2++ = *s1; } ++s1; } *s2 = 0; } int misses(char *dna1, char *dna2) /* compare two short dna strands of equal length, return number of mismatches */ { int len1 = strlen(dna1); int len2 = strlen(dna2); int misscnt = 0; int i = 0; int len = len1; if (len2 < len) len = len2; for(i=0;i<len;++i) { if (toupper(dna1[i]) != toupper(dna2[i])) { ++misscnt; } } misscnt += abs(len2-len1); return misscnt; } char *checkAndFetchNib(struct dnaSeq *chromSeq, struct snp *snp, int lf, int ls) /* fetch nib return string to be freed later. Reports if segment exceends ends of chromosome. */ { if ((snp->chromStart - lf) < 0) { printf("snp error (chromStart - offset) < 0 : %s %s %u %u %d \n", snp->name, snp->chrom, snp->chromStart, snp->chromEnd, lf ); return NULL; } if ((snp->chromStart - lf + ls) > chromSeq->size) { printf("bed error chromStart - leftFlank + seqlen > chromSeq->size : %s %s %u %u chr-size: %u \n", snp->name, snp->chrom, snp->chromStart, snp->chromEnd, chromSeq->size ); return NULL; } return cloneStringZ(chromSeq->dna + snp->chromStart - lf, ls); } void writeExOut(enum snpExceptionType i, struct snp *snp) { fprintf(exf[i],"%s\t%d\t%d\t%s\t%d\t%s\t%s\n", snp->chrom, snp->chromStart, snp->chromEnd, snp->name, i+BASE_SNP_EX_NUM, snp->strand, snp->locType ); } /* ---------------------------------------------------------------- */ void snpValid() /* Test snp for one assembly. */ { struct axtScoreScheme *simpleDnaScheme = NULL; int match = 0; /* good match of minimal acceptable quality */ int mismatch = 0; /* unacceptable match quality */ int missing = 0; /* unable to find rsId in dbSnpRs/affy */ int goodrc = 0; /* matches after reverse-complement */ int assemblyDash = 0; /* assembly context is just a single dash - (complex cases) */ int gapNib = 0; /* nib returns n's, we are in the gap */ int strandMismatch = 0; /* reported strand differs from what we found */ int totalMatch = 0; int totalMismatch = 0; int totalMissing = 0; int totalGoodrc = 0; int totalAssemblyDash = 0; int totalGapNib = 0; int totalStrandMismatch = 0; void *next = NULL; char *id = NULL; char *seq = NULL; int matchScore = 100; int misMatchScore = 100; int gapOpenPenalty = 400; int gapExtendPenalty = 5; /* was 50, reducing for new snp version. */ int noDna = 0; int snpRows = 0; int goodExact = 0; int badExact = 0; struct slName *cns=NULL,*cn=NULL; if (!hDbIsActive(db)) { printf("Currently no support for db %s\n", db); return; } (void) hOrganism(db); // ignore returned organism name simpleDnaScheme = axtScoreSchemeSimpleDna(matchScore, misMatchScore, gapOpenPenalty, gapExtendPenalty); if (!cns) { printf("testDb: hAllChromNames returned empty list \n"); return; } openExOuts(db); /* open separate files for each snp class of Exception of interest */ printf("maxFlank = %d \n",maxFlank); printf("threshold = %d \n",threshold); cns = hAllChromNames(db); for (cn = cns; cn != NULL; cn = cn->next) { struct dnaSeq *chromSeq = NULL; struct snp *snps = NULL; struct snp *snp = NULL; struct flank *flanks = NULL; struct flank *flank = NULL; if (chr != NULL) if (!sameWord(chr,cn->name)) continue; uglyf("reading flanks %s \n",cn->name); flanks = readFlank(cn->name); if (flanks == NULL) { printf("readFlank returned NULL for chrom %s.\n",cn->name); continue; } printf("readFlank done for chrom %s \n",cn->name); printf("slCount(flanks)=%d for chrom %s \n",slCount(flanks),cn->name); slSort(&flanks, slFlankCmp); printf("slSort done for chrom %s \n",cn->name); uglyf("beginning chrom %s \n",cn->name); chromSeq = hLoadChrom(db, cn->name); printf("chrom %s : size (%u) \n",cn->name,chromSeq->size); snps = readSnp(cn->name); printf("read %s.snp done for chrom=%s \n",db,cn->name); flank = flanks; printf("=========================================================\n"); for (snp = snps; snp != NULL; snp = snp->next) { int cmp = -1; char *nibDna=NULL; ++snpRows; if (Verbose) { printf("debug %s %u %u %s %s\n", snp->chrom, snp->chromStart, snp->chromEnd, snp->name, snp->strand ); } /* continue; */ /* add check for Heather for specific case locType=exact and size=1 */ if (sameString(snp->locType,"exact") && (snp->chromEnd == (snp->chromStart + 1))) { char *obs=cloneString(snp->observed); int n = chopString(obs, "/", NULL, 0); char **obsvd = needMem(n*sizeof(char*)); char *exactDna = checkAndFetchNib(chromSeq, snp, 0, 1); int i=0; boolean found = FALSE; if (sameString(snp->strand,"-")) { reverseComplement(exactDna,strlen(exactDna)); } chopString(obs, "/", obsvd, n); if (exactDna==NULL) { printf("exactDna=NULL for %s %s %u %u (%s) %s \n", snp->name, snp->chrom, snp->chromStart, snp->chromEnd, snp->observed, snp->locType ); } /* uglyf("%s: exactDna=%s obs=%s \n",snp->name,exactDna,snp->observed); */ for(i=0; i<n; i++) { if (strlen(obsvd[i]) > 1) { printf("%s: incorrect length of observed %s <> 1 \n", snp->name, obsvd[i] ); } if (sameWord(obsvd[i],exactDna)) { found = TRUE; } } if (found) { goodExact++; } else { badExact++; uglyf("id: %s exact %s not found in observed %s \n",snp->name,exactDna,snp->observed); writeExOut(snpExNotExact, snp); } freez(&obsvd); freez(&obs); freez(&exactDna); } while (cmp < 0) { while (cmp < 0) { next = flank; if (next == NULL) { cmp = 1; } else { break; } } if (cmp < 0) { id = flank->rsId; cmp=mystrcmp(id, snp->name); } if (cmp < 0) { flank = flank->next; } } if (cmp==0) { int strand=1; char *rc = NULL; int m = 0; int lf = 0; /* size of left flank context (lower case dna) */ int rf = 0; /* size of right flank context (lower case dna) */ int ls = 0; /* total size of assembly dna context plus actual region in dbSnpRs/affy */ char *origSeq = NULL; /* use to display the original dnSnpRs.assembly seq */ int flankSize = 0; lf=strlen(flank->leftFlank); rf=strlen(flank->rightFlank); /* if flanks exceed maxFlank, truncate */ if (lf>maxFlank) { char *temp=flank->leftFlank; flank->leftFlank=cloneString(temp+lf-maxFlank); lf = maxFlank; freez(&temp); } if (rf>maxFlank) { rf = maxFlank; flank->rightFlank[rf]=0; } /* at Daryl's request, try to make them same case */ toLowerN(flank->leftFlank , lf); toLowerN(flank->rightFlank, rf); flankSize = lf+1+rf; seq = needMem(flankSize+1); safef(seq,flankSize+1,"%s-%s",flank->leftFlank,flank->rightFlank); if (sameString(seq,"-")) { ++assemblyDash; if (Verbose) printf("(no assembly context) rsId=%s chrom=%s %u %u \n assembly=%s \n\n", id, snp->chrom, snp->chromStart, snp->chromEnd, seq ); continue; } origSeq = seq; seq = cloneString(origSeq); /* remove dashes indicating insert to simplify and correct processing of nib data */ stripDashes(seq); ls = lf + rf + (snp->chromEnd - snp->chromStart); nibDna = checkAndFetchNib(chromSeq, snp, lf, ls); if (nibDna==NULL) { ++noDna; printf("no dna for %s %s %u %u \n", snp->name, snp->chrom, snp->chromStart, snp->chromEnd ); continue; } if (allNs(nibDna)) { ++gapNib; ++mismatch; writeExOut(snpExMismatch, snp); if (Verbose) printf("(nib gap) rsId=%s chrom=%s %u %u \n assembly=%s \n snpMap=%s \n\n", id, snp->chrom, snp->chromStart, snp->chromEnd, seq, nibDna ); continue; } m = misses(seq,nibDna); if (m > 1) { int n = -1; rc = checkAndFetchNib(chromSeq, snp, rf, ls); if (rc==NULL) { ++noDna; printf("no dna for %s %s %u %u \n", snp->name, snp->chrom, snp->chromStart, snp->chromEnd ); continue; } reverseComplement(rc,strlen(rc)); n = misses(seq, rc); if (n < m) { strand=-1; m = n; } } if (m <= 1) { ++match; if (strand < 1) ++goodrc; } else { struct dnaSeq query, target; struct axt *axtAln = NULL; int bestScore = 0; ZeroVar(&query); query.dna = seq; query.size = strlen(query.dna); ZeroVar(&target); target.dna = nibDna; target.size = strlen(target.dna); axtAln = axtAffine(&query, &target, simpleDnaScheme); strand = 1; if (axtAln) { bestScore = axtAln->score / (lf+rf); } axtFree(&axtAln); if (bestScore < threshold) { ZeroVar(&target); target.dna = rc; target.size = strlen(target.dna); axtAln = axtAffine(&query, &target, simpleDnaScheme); if ((axtAln) && (bestScore < (axtAln->score / (lf+rf)))) { strand = -1; bestScore = axtAln->score / (lf+rf); } axtFree(&axtAln); } if (bestScore >= threshold) { ++match; if (strand < 1) ++goodrc; if (((strand == 1) && (!sameString(snp->strand,"+"))) || ((strand == -1) && (!sameString(snp->strand,"-")))) { strandMismatch++; /* reported strand differs from what we found */ writeExOut(snpExWrongStrand, snp); } } else { ++mismatch; writeExOut(snpExMismatch, snp); } if ((bestScore < threshold) || Verbose) { printf( "score=%d misses=%u strand=%d repstr=%s rsId=%s chrom=%s %u %u lf=%d ls=%d \n" " flanks=%s \n" " chr=%s \n" " rc chr=%s \n" "\n", bestScore, m, strand, /* best match found */ snp->strand, /* reported from snp table */ id, snp->chrom, snp->chromStart, snp->chromEnd, lf, ls, seq, nibDna, rc ); } } freez(&rc); freez(&seq); freez(&origSeq); } else /* if (cmp==0) */ { char snpLkup[10] = ""; /* this id is missing from flanks */ ++missing; safef(snpLkup,sizeof(snpLkup),"flnk%s",snp->chrom); if (Verbose) printf("snp.name=%s is missing from %s (now at %s) \n\n",snp->name,snpLkup,id); writeExOut(snpExNoFlanks, snp); } freez(&nibDna); } snpFreeList(&snps); flankFreeList(&flanks); dnaSeqFree(&chromSeq); printf("\n\n\n Total matches for chrom %s:\n ",cn->name); printf(" matches: %u \n ",match); printf(" mismatches: %u \n",mismatch); printf(" missing from flanks: %u \n",missing); printf(" rev compl matches: %u \n",goodrc); printf(" not rptd strand : %u \n",strandMismatch); printf(" assembly = -: %u \n",assemblyDash); printf(" nib in gap : %u \n",gapNib); printf("\n\n=========================================\n"); totalMatch += match; totalMismatch += mismatch; totalMissing += missing; totalGoodrc += goodrc; totalStrandMismatch += strandMismatch; totalAssemblyDash += assemblyDash; totalGapNib += gapNib; match = 0; mismatch = 0; missing = 0; goodrc = 0; strandMismatch = 0; assemblyDash = 0; gapNib = 0; printf("\n"); printf("\n"); } slFreeList(&cns); closeExOuts(); /* close Exception class output files */ axtScoreSchemeFree(&simpleDnaScheme); printf("\n\n\n Grand Totals: \n "); printf(" matches: %u \n",totalMatch); printf(" mismatches: %u \n",totalMismatch); printf(" missing from flanks: %u \n",totalMissing); printf(" rev compl matches: %u \n",totalGoodrc); printf(" not rptd strand : %u \n",totalStrandMismatch); printf(" assembly = -: %u \n",totalAssemblyDash); printf(" nib in gap : %u \n",totalGapNib); printf("\n Total rows in snp: %u \n ",snpRows); printf("\n # no dna found for : %u \n ",noDna); printf("\n Total goodExact: %u \n ",goodExact); printf("\n Total badExact: %u \n ",badExact); printf("\n\n=========================================\n"); } int main(int argc, char *argv[]) /* Process command line. */ { /* pushCarefulMemHandler(200000000); */ /* Set initial seed */ srand( (unsigned)time( NULL ) ); optionInit(&argc, argv, options); if (argc != 3) usage(); db = argv[1]; flankPath = argv[2]; chr = optionVal("chr" , chr ); threshold = optionInt("threshold", threshold); Verbose = optionInt("verbose" , Verbose ); maxFlank = optionInt("maxFlank" , maxFlank ); snpValid(); carefulCheckHeap(); return 0; }