6d437c28a5ff82a2af59ca516559053c4bf7c3f9 max Thu Jun 18 01:18:34 2026 -0700 hgTracks: amino-acid name in codon + ruler mouseovers; exon-length label; ruler "Complement bases" config toggle Three genome browser display changes around AA display: 1) amino acids shown on genePred codons and on rule codons, and making the "complement option" easier to find. 1) The genePred/bigGenePred codon mouseover (zoomed to the codon level) now shows the codon's amino acid on its own "Amino acid:" line, as the three-letter abbreviation plus full name (e.g. "Ala (alanine)"). Rather than reverse-decode the amino acid out of the codon's packed grayIx, the one-letter code is now stored on the codon when it is translated: struct simpleFeature gains a codonAa field, codonToGrayIx() reports the letter through a new out-param, and the codon mouseover reads simpleFeature.codonAa directly. The drawn codon letter and the mouseover then derive from the same translation and cannot drift apart, without the mouseover having to understand the grayIx encoding (the drawing path, colorAndCodonFromGrayIx(), still decodes its own grayIx inline, unchanged). Adds a one-letter->full-name aaToName() in lib/dnautil.c (using the previously unused name field of aminoAcidTable). Stop codons show "Ter (termination)". Also relabels the exon "Length" field as "Exon Length" in the codon and zoomed-out exon mouseovers; introns keep "Length". 2) The base-position ruler's three-reading-frame translation (hgt.baseShowCodons) now gives each codon box a mouse-over with the same three-letter abbreviation and full name, reading the stored codonAa via aaToName()/aaToAbbr() (baseColorDrawRulerCodons in cds.c). 3) Adds a "Complement the bases" checkbox to the base position (ruler) track configuration page (hgTrackUi rulerUi). It toggles the same per-assembly COMPLEMENT_BASES_VAR cart variable that was previously reachable only by clicking the "Click to complement" arrow next to the ruler. Also fixes the indentation of the adjacent drawComplementArrow() text assignments. refs #37779 Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com> diff --git src/lib/dnautil.c src/lib/dnautil.c index 56c46142460..5e27251f5a2 100644 --- src/lib/dnautil.c +++ src/lib/dnautil.c @@ -1,1297 +1,1311 @@ /* Some stuff that you'll likely need in any program that works with * DNA. Includes stuff for amino acids as well. * * Assumes that DNA is stored as a character. * The DNA it generates will include the bases * as lowercase tcag. It will generally accept * uppercase as well, and also 'n' or 'N' or '-' * for unknown bases. * * Amino acids are stored as single character upper case. * * This file is copyright 2002 Jim Kent, but license is hereby * granted for all use - public, private or commercial. */ #include "common.h" #include "dnautil.h" struct codonTable /* The dread codon table. */ { DNA *codon; /* Lower case. */ AA protCode; /* Upper case. The "Standard" code */ AA mitoCode; /* Upper case. Vertebrate Mitochondrial translations */ AA uniqCode; /* unique code for each codon */ }; struct codonTable codonTable[] = /* The master codon/protein table. */ { {"ttt", 'F', 'F', 'a'}, {"ttc", 'F', 'F', 'b'}, {"tta", 'L', 'L', 'c'}, {"ttg", 'L', 'L', 'd'}, {"tct", 'S', 'S', 'e'}, {"tcc", 'S', 'S', 'f'}, {"tca", 'S', 'S', 'g'}, {"tcg", 'S', 'S', 'h'}, {"tat", 'Y', 'Y', 'i'}, {"tac", 'Y', 'Y', 'j'}, {"taa", 0, 0, 'k'}, {"tag", 0, 0, 'l'}, {"tgt", 'C', 'C', 'm'}, {"tgc", 'C', 'C', 'n'}, {"tga", 0, 'W', 'o'}, {"tgg", 'W', 'W', 'p'}, {"ctt", 'L', 'L', 'q'}, {"ctc", 'L', 'L', 'r'}, {"cta", 'L', 'L', 's'}, {"ctg", 'L', 'L', 't'}, {"cct", 'P', 'P', 'u'}, {"ccc", 'P', 'P', 'v'}, {"cca", 'P', 'P', 'w'}, {"ccg", 'P', 'P', 'x'}, {"cat", 'H', 'H', 'y'}, {"cac", 'H', 'H', 'z'}, {"caa", 'Q', 'Q', 'A'}, {"cag", 'Q', 'Q', 'B'}, {"cgt", 'R', 'R', 'C'}, {"cgc", 'R', 'R', 'D'}, {"cga", 'R', 'R', 'E'}, {"cgg", 'R', 'R', 'F'}, {"att", 'I', 'I', 'G'}, {"atc", 'I', 'I', 'H'}, {"ata", 'I', 'M', 'I'}, {"atg", 'M', 'M', 'J'}, {"act", 'T', 'T', 'K'}, {"acc", 'T', 'T', 'L'}, {"aca", 'T', 'T', 'M'}, {"acg", 'T', 'T', 'N'}, {"aat", 'N', 'N', 'O'}, {"aac", 'N', 'N', 'P'}, {"aaa", 'K', 'K', 'Q'}, {"aag", 'K', 'K', 'R'}, {"agt", 'S', 'S', 'S'}, {"agc", 'S', 'S', 'T'}, {"aga", 'R', 0, 'U'}, {"agg", 'R', 0, 'V'}, {"gtt", 'V', 'V', 'W'}, {"gtc", 'V', 'V', 'X'}, {"gta", 'V', 'V', 'Y'}, {"gtg", 'V', 'V', 'Z'}, {"gct", 'A', 'A', '1'}, {"gcc", 'A', 'A', '2'}, {"gca", 'A', 'A', '3'}, {"gcg", 'A', 'A', '4'}, {"gat", 'D', 'D', '5'}, {"gac", 'D', 'D', '6'}, {"gaa", 'E', 'E', '7'}, {"gag", 'E', 'E', '8'}, {"ggt", 'G', 'G', '9'}, {"ggc", 'G', 'G', '0'}, {"gga", 'G', 'G', '@'}, {"ggg", 'G', 'G', '$'}, }; /* A table that gives values 0 for t 1 for c 2 for a 3 for g * (which is order aa's are in biochemistry codon tables) * and gives -1 for all others. */ int ntVal[256]; int ntValLower[256]; /* NT values only for lower case. */ int ntValUpper[256]; /* NT values only for upper case. */ int ntVal5[256]; int ntValNoN[256]; /* Like ntVal, but with T_BASE_VAL in place of -1 for nonexistent ones. */ DNA valToNt[(N_BASE_VAL|MASKED_BASE_BIT)+1]; /* convert tables for bit-4 indicating masked */ int ntValMasked[256]; DNA valToNtMasked[256]; static boolean inittedNtVal = FALSE; static void initNtVal() { if (!inittedNtVal) { int i; for (i=0; i<ArraySize(ntVal); i++) { ntValUpper[i] = ntValLower[i] = ntVal[i] = -1; ntValNoN[i] = T_BASE_VAL; if (isspace(i) || isdigit(i)) ntVal5[i] = ntValMasked[i] = -1; else { ntVal5[i] = N_BASE_VAL; ntValMasked[i] = (islower(i) ? (N_BASE_VAL|MASKED_BASE_BIT) : N_BASE_VAL); } } ntVal5['t'] = ntVal5['T'] = ntValNoN['t'] = ntValNoN['T'] = ntVal['t'] = ntVal['T'] = ntValLower['t'] = ntValUpper['T'] = T_BASE_VAL; ntVal5['u'] = ntVal5['U'] = ntValNoN['u'] = ntValNoN['U'] = ntVal['u'] = ntVal['U'] = ntValLower['u'] = ntValUpper['U'] = U_BASE_VAL; ntVal5['c'] = ntVal5['C'] = ntValNoN['c'] = ntValNoN['C'] = ntVal['c'] = ntVal['C'] = ntValLower['c'] = ntValUpper['C'] = C_BASE_VAL; ntVal5['a'] = ntVal5['A'] = ntValNoN['a'] = ntValNoN['A'] = ntVal['a'] = ntVal['A'] = ntValLower['a'] = ntValUpper['A'] = A_BASE_VAL; ntVal5['g'] = ntVal5['G'] = ntValNoN['g'] = ntValNoN['G'] = ntVal['g'] = ntVal['G'] = ntValLower['g'] = ntValUpper['G'] = G_BASE_VAL; valToNt[T_BASE_VAL] = valToNt[T_BASE_VAL|MASKED_BASE_BIT] = 't'; valToNt[C_BASE_VAL] = valToNt[C_BASE_VAL|MASKED_BASE_BIT] = 'c'; valToNt[A_BASE_VAL] = valToNt[A_BASE_VAL|MASKED_BASE_BIT] = 'a'; valToNt[G_BASE_VAL] = valToNt[G_BASE_VAL|MASKED_BASE_BIT] = 'g'; valToNt[N_BASE_VAL] = valToNt[N_BASE_VAL|MASKED_BASE_BIT] = 'n'; /* masked values */ ntValMasked['T'] = T_BASE_VAL; ntValMasked['U'] = U_BASE_VAL; ntValMasked['C'] = C_BASE_VAL; ntValMasked['A'] = A_BASE_VAL; ntValMasked['G'] = G_BASE_VAL; ntValMasked['t'] = T_BASE_VAL|MASKED_BASE_BIT; ntValMasked['u'] = U_BASE_VAL|MASKED_BASE_BIT; ntValMasked['c'] = C_BASE_VAL|MASKED_BASE_BIT; ntValMasked['a'] = A_BASE_VAL|MASKED_BASE_BIT; ntValMasked['g'] = G_BASE_VAL|MASKED_BASE_BIT; valToNtMasked[T_BASE_VAL] = 'T'; valToNtMasked[C_BASE_VAL] = 'C'; valToNtMasked[A_BASE_VAL] = 'A'; valToNtMasked[G_BASE_VAL] = 'G'; valToNtMasked[N_BASE_VAL] = 'N'; valToNtMasked[T_BASE_VAL|MASKED_BASE_BIT] = 't'; valToNtMasked[C_BASE_VAL|MASKED_BASE_BIT] = 'c'; valToNtMasked[A_BASE_VAL|MASKED_BASE_BIT] = 'a'; valToNtMasked[G_BASE_VAL|MASKED_BASE_BIT] = 'g'; valToNtMasked[N_BASE_VAL|MASKED_BASE_BIT] = 'n'; inittedNtVal = TRUE; } } /* Returns one letter code for protein, * 0 for stop codon or X for bad input, * The "Standard" Code */ AA lookupCodon(DNA *dna) { int ix; int i; char c; if (!inittedNtVal) initNtVal(); ix = 0; for (i=0; i<3; ++i) { int bv = ntVal[(int)dna[i]]; if (bv<0) return 'X'; ix = (ix<<2) + bv; } c = codonTable[ix].protCode; return c; } boolean isStopCodon(DNA *dna) /* Return TRUE if it's a stop codon. */ { return lookupCodon(dna) == 0; } boolean isKozak(char *dna, int dnaSize, int pos) /* Return TRUE if it's a Kozak compatible start. */ { if (lookupCodon(dna+pos) != 'M') { return FALSE; } if (pos + 3 < dnaSize) { if (ntVal[(int)dna[pos+3]] == G_BASE_VAL) return TRUE; } if (pos >= 3) { int c = ntVal[(int)dna[pos-3]]; if (c == A_BASE_VAL || c == G_BASE_VAL) return TRUE; } return FALSE; } boolean isReallyStopCodon(char *dna, boolean selenocysteine) /* Return TRUE if it's really a stop codon, even considering * possibilility of selenocysteine. */ { if (selenocysteine) { /* Luckily the mitochondria *also* replaces TGA with * something else, even though it isn't selenocysteine */ return lookupMitoCodon(dna) == 0; } else { return lookupCodon(dna) == 0; } } /* Returns one letter code for protein, * 0 for stop codon or X for bad input, * Vertebrate Mitochondrial Code */ AA lookupMitoCodon(DNA *dna) { int ix; int i; char c; if (!inittedNtVal) initNtVal(); ix = 0; for (i=0; i<3; ++i) { int bv = ntVal[(int)dna[i]]; if (bv<0) return 'X'; ix = (ix<<2) + bv; } c = codonTable[ix].mitoCode; c = toupper(c); return c; } AA lookupUniqCodon(DNA *dna) { int ix; int i; char c; if (!inittedNtVal) initNtVal(); ix = 0; for (i=0; i<3; ++i) { int bv = ntVal[(int)dna[i]]; if (bv<0) return 'X'; ix = (ix<<2) + bv; } c = codonTable[ix].uniqCode; c = toupper(c); return c; } Codon codonVal(DNA *start) /* Return value from 0-63 of codon starting at start. * Returns -1 if not a codon. */ { int v1,v2,v3; if ((v1 = ntVal[(int)start[0]]) < 0) return -1; if ((v2 = ntVal[(int)start[1]]) < 0) return -1; if ((v3 = ntVal[(int)start[2]]) < 0) return -1; return ((v1<<4) + (v2<<2) + v3); } DNA *valToCodon(int val) /* Return codon corresponding to val (0-63) */ { assert(val >= 0 && val < 64); return codonTable[val].codon; } void dnaTranslateSome(DNA *dna, char *out, int outSize) /* Translate DNA upto a stop codon or until outSize-1 amino acids, * whichever comes first. Output will be zero terminated. */ { int i; int dnaSize; int protSize = 0; outSize -= 1; /* Room for terminal zero */ dnaSize = strlen(dna); for (i=0; i<dnaSize-2; i+=3) { if (protSize >= outSize) break; if ((out[protSize++] = lookupCodon(dna+i)) == 0) break; } out[protSize] = 0; } /* A little array to help us decide if a character is a * nucleotide, and if so convert it to lower case. */ char ntChars[256]; static void initNtChars() { static boolean initted = FALSE; if (!initted) { zeroBytes(ntChars, sizeof(ntChars)); ntChars['a'] = ntChars['A'] = 'a'; ntChars['c'] = ntChars['C'] = 'c'; ntChars['g'] = ntChars['G'] = 'g'; ntChars['t'] = ntChars['T'] = 't'; ntChars['n'] = ntChars['N'] = 'n'; ntChars['u'] = ntChars['U'] = 'u'; ntChars['-'] = 'n'; initted = TRUE; } } char ntMixedCaseChars[256]; static void initNtMixedCaseChars() { static boolean initted = FALSE; if (!initted) { zeroBytes(ntMixedCaseChars, sizeof(ntMixedCaseChars)); ntMixedCaseChars['a'] = 'a'; ntMixedCaseChars['A'] = 'A'; ntMixedCaseChars['c'] = 'c'; ntMixedCaseChars['C'] = 'C'; ntMixedCaseChars['g'] = 'g'; ntMixedCaseChars['G'] = 'G'; ntMixedCaseChars['t'] = 't'; ntMixedCaseChars['T'] = 'T'; ntMixedCaseChars['n'] = 'n'; ntMixedCaseChars['N'] = 'N'; ntMixedCaseChars['u'] = 'u'; ntMixedCaseChars['U'] = 'U'; ntMixedCaseChars['-'] = 'n'; initted = TRUE; } } /* Another array to help us do complement of DNA */ DNA ntCompTable[256]; static boolean inittedCompTable = FALSE; static void initNtCompTable() { zeroBytes(ntCompTable, sizeof(ntCompTable)); ntCompTable[' '] = ' '; ntCompTable['-'] = '-'; ntCompTable['='] = '='; ntCompTable['a'] = 't'; ntCompTable['c'] = 'g'; ntCompTable['g'] = 'c'; ntCompTable['t'] = 'a'; ntCompTable['u'] = 'a'; ntCompTable['n'] = 'n'; ntCompTable['-'] = '-'; ntCompTable['.'] = '.'; ntCompTable['A'] = 'T'; ntCompTable['C'] = 'G'; ntCompTable['G'] = 'C'; ntCompTable['T'] = 'A'; ntCompTable['U'] = 'A'; ntCompTable['N'] = 'N'; ntCompTable['R'] = 'Y'; ntCompTable['Y'] = 'R'; ntCompTable['M'] = 'K'; ntCompTable['K'] = 'M'; ntCompTable['S'] = 'S'; ntCompTable['W'] = 'W'; ntCompTable['V'] = 'B'; ntCompTable['H'] = 'D'; ntCompTable['D'] = 'H'; ntCompTable['B'] = 'V'; ntCompTable['X'] = 'N'; ntCompTable['r'] = 'y'; ntCompTable['y'] = 'r'; ntCompTable['s'] = 's'; ntCompTable['w'] = 'w'; ntCompTable['m'] = 'k'; ntCompTable['k'] = 'm'; ntCompTable['v'] = 'b'; ntCompTable['h'] = 'd'; ntCompTable['d'] = 'h'; ntCompTable['b'] = 'v'; ntCompTable['x'] = 'n'; ntCompTable['('] = ')'; ntCompTable[')'] = '('; inittedCompTable = TRUE; } /* Complement DNA (not reverse). */ void complement(DNA *dna, long length) { int i; if (!inittedCompTable) initNtCompTable(); for (i=0; i<length; ++i) { *dna = ntCompTable[(int)*dna]; ++dna; } } /* Reverse complement DNA. */ void reverseComplement(DNA *dna, long length) { reverseBytes(dna, length); complement(dna, length); } /* Reverse offset - return what will be offset (0 based) to * same member of array after array is reversed. */ long reverseOffset(long offset, long arraySize) { return arraySize-1 - offset; } /* Switch start/end (zero based half open) coordinates * to opposite strand. */ void reverseIntRange(int *pStart, int *pEnd, int size) { int temp; temp = *pStart; *pStart = size - *pEnd; *pEnd = size - temp; } /* Switch start/end (zero based half open) coordinates * to opposite strand. */ void reverseUnsignedRange(unsigned *pStart, unsigned *pEnd, int size) { unsigned temp; temp = *pStart; *pStart = size - *pEnd; *pEnd = size - temp; } char *reverseComplementSlashSeparated(char *alleleStr) /* Given a slash-separated series of sequences (a common representation of variant alleles), * returns a slash-sep series with the reverse complement of each sequence (if it is a * nucleotide sequence). * Special behavior to support dbSNP's variant allele conventions: * 1. Reverse the order of sequences (to maintain alphabetical ordering). * 2. If alleleStr begins with "-/", then after reversing, move "-/" back to the beginning. */ { int len = strlen(alleleStr); char choppyCopy[len+1]; safecpy(choppyCopy, sizeof(choppyCopy), alleleStr); char *alleles[len]; int alCount = chopByChar(choppyCopy, '/', alleles, ArraySize(alleles)); char *outStr = needMem(len+1); int i; for (i = alCount-1; i >= 0; i--) { char *allele = alleles[i]; int alLen = strlen(allele); if (isAllNt(allele, alLen)) reverseComplement(allele, alLen); if (i != alCount-1) safecat(outStr, len+1, "/"); safecat(outStr, len+1, allele); } if (startsWith("-/", alleleStr)) { // Keep "-/" at the beginning: memmove(outStr+2, outStr, len-2); outStr[0] = '-'; outStr[1] = '/'; } return outStr; } int cmpDnaStrings(DNA *a, DNA *b) /* Compare using screwy non-alphabetical DNA order TCGA */ { for (;;) { DNA aa = *a++; DNA bb = *b++; if (aa != bb) return ntVal[(int)aa] - ntVal[(int)bb]; if (aa == 0) break; } return 0; } /* Convert U's to T's */ void toDna(DNA *rna) { DNA c; for (;;) { c = *rna; if (c == 'u') *rna = 't'; else if (c == 'U') *rna = 'T'; else if (c == 0) break; ++rna; } } /* Convert T's to U's */ void toRna(DNA *dna) { DNA c; for (;;) { c = *dna; if (c == 't') *dna = 'u'; else if (c == 'T') *dna = 'U'; else if (c == 0) break; ++dna; } } char *skipIgnoringDash(char *a, int size, bool skipTrailingDash) /* Count size number of characters, and any * dash characters. */ { while (size > 0) { if (*a++ != '-') --size; } if (skipTrailingDash) while (*a == '-') ++a; return a; } int countNonDash(char *a, int size) /* Count number of non-dash characters. */ { int count = 0; int i; for (i=0; i<size; ++i) if (a[i] != '-') ++count; return count; } int nextPowerOfFour(long x) /* Return next power of four that would be greater or equal to x. * For instance if x < 4, return 1, if x < 16 return 2.... * (From biological point of view how many bases are needed to * code this number.) */ { int count = 1; while (x > 4) { count += 1; x >>= 2; } return count; } long dnaOrAaFilteredSize(char *raw, char filter[256]) /* Return how long DNA will be after non-DNA is filtered out. */ { char c; long count = 0; dnaUtilOpen(); while ((c = *raw++) != 0) { if (filter[(int)c]) ++count; } return count; } void dnaOrAaFilter(char *in, char *out, char filter[256]) /* Run chars through filter. */ { char c; dnaUtilOpen(); while ((c = *in++) != 0) { if ((c = filter[(int)c]) != 0) *out++ = c; } *out++ = 0; } long dnaFilteredSize(char *rawDna) /* Return how long DNA will be after non-DNA is filtered out. */ { return dnaOrAaFilteredSize(rawDna, ntChars); } void dnaFilter(char *in, DNA *out) /* Filter out non-DNA characters and change to lower case. */ { dnaOrAaFilter(in, out, ntChars); } void dnaFilterToN(char *in, DNA *out) /* Change all non-DNA characters to N. */ { DNA c; initNtChars(); while ((c = *in++) != 0) { if ((c = ntChars[(int)c]) != 0) *out++ = c; else *out++ = 'n'; } *out++ = 0; } void dnaMixedCaseFilter(char *in, DNA *out) /* Filter out non-DNA characters but leave case intact. */ { dnaOrAaFilter(in, out, ntMixedCaseChars); } long aaFilteredSize(char *raw) /* Return how long aa will be after non-aa chars is filtered out. */ { return dnaOrAaFilteredSize(raw, aaChars); } void aaFilter(char *in, DNA *out) /* Filter out non-aa characters and change to upper case. */ { dnaOrAaFilter(in, out, aaChars); } void upperToN(char *s, int size) /* Turn upper case letters to N's. */ { char c; int i; for (i=0; i<size; ++i) { c = s[i]; if (isupper(c)) s[i] = 'n'; } } void lowerToN(char *s, int size) /* Turn lower case letters to N's. */ { char c; int i; for (i=0; i<size; ++i) { c = s[i]; if (islower(c)) s[i] = 'N'; } } void dnaBaseHistogram(DNA *dna, int dnaSize, int histogram[4]) /* Count up frequency of occurance of each base and store * results in histogram. */ { int val; zeroBytes(histogram, 4*sizeof(int)); while (--dnaSize >= 0) { if ((val = ntVal[(int)*dna++]) >= 0) ++histogram[val]; } } bits64 basesToBits64(char *dna, int size) /* Convert dna of given size (up to 32) to binary representation */ { if (size > 32) errAbort("basesToBits64 called on %d bases, can only go up to 32", size); bits64 result = 0; int i; for (i=0; i<size; ++i) { result <<= 2; result += ntValNoN[(int)dna[i]]; } return result; } bits32 packDna16(DNA *in) /* pack 16 bases into a word */ { bits32 out = 0; int count = 16; int bVal; while (--count >= 0) { bVal = ntValNoN[(int)*in++]; out <<= 2; out += bVal; } return out; } bits16 packDna8(DNA *in) /* Pack 8 bases into a short word */ { bits16 out = 0; int count = 8; int bVal; while (--count >= 0) { bVal = ntValNoN[(int)*in++]; out <<= 2; out += bVal; } return out; } UBYTE packDna4(DNA *in) /* Pack 4 bases into a UBYTE */ { UBYTE out = 0; int count = 4; int bVal; while (--count >= 0) { bVal = ntValNoN[(int)*in++]; out <<= 2; out += bVal; } return out; } void unpackDna(bits32 *tiles, int tileCount, DNA *out) /* Unpack DNA. Expands to 16x tileCount in output. */ { int i, j; bits32 tile; for (i=0; i<tileCount; ++i) { tile = tiles[i]; for (j=15; j>=0; --j) { out[j] = valToNt[tile & 0x3]; tile >>= 2; } out += 16; } } void unpackDna4(UBYTE *tiles, int byteCount, DNA *out) /* Unpack DNA. Expands to 4x byteCount in output. */ { int i, j; UBYTE tile; for (i=0; i<byteCount; ++i) { tile = tiles[i]; for (j=3; j>=0; --j) { out[j] = valToNt[tile & 0x3]; tile >>= 2; } out += 4; } } static void checkSizeTypes() /* Make sure that some of our predefined types are the right size. */ { assert(sizeof(UBYTE) == 1); assert(sizeof(WORD) == 2); assert(sizeof(bits32) == 4); assert(sizeof(bits16) == 2); } int intronOrientationMinSize(DNA *iStart, DNA *iEnd, int minIntronSize) /* Given a gap in genome from iStart to iEnd, return * Return 1 for GT/AG intron between left and right, -1 for CT/AC, 0 for no * intron. Assumes DNA is lower cased. */ { if (iEnd - iStart < minIntronSize) return 0; if (iStart[0] == 'g' && iStart[1] == 't' && iEnd[-2] == 'a' && iEnd[-1] == 'g') { return 1; } else if (iStart[0] == 'c' && iStart[1] == 't' && iEnd[-2] == 'a' && iEnd[-1] == 'c') { return -1; } else return 0; } int intronOrientation(DNA *iStart, DNA *iEnd) /* Given a gap in genome from iStart to iEnd, return * Return 1 for GT/AG intron between left and right, -1 for CT/AC, 0 for no * intron. Assumes DNA is lower cased. */ { return intronOrientationMinSize(iStart, iEnd, 32); } int dnaScore2(DNA a, DNA b) /* Score match between two bases (relatively crudely). */ { if (a == 'n' || b == 'n') return 0; if (a == b) return 1; else return -1; } int dnaOrAaScoreMatch(char *a, char *b, int size, int matchScore, int mismatchScore, char ignore) /* Compare two sequences (without inserts or deletions) and score. */ { int i; int score = 0; for (i=0; i<size; ++i) { char aa = a[i]; char bb = b[i]; if (aa == ignore || bb == ignore) continue; if (aa == bb) score += matchScore; else score += mismatchScore; } return score; } int dnaScoreMatch(DNA *a, DNA *b, int size) /* Compare two pieces of DNA base by base. Total mismatches are * subtracted from total matches and returned as score. 'N's * neither hurt nor help score. */ { return dnaOrAaScoreMatch(a, b, size, 1, -1, 'n'); } int aaScore2(AA a, AA b) /* Score match between two amino acids (relatively crudely). */ { if (a == 'X' || b == 'X') return 0; if (a == b) return 2; else return -1; } int aaScoreMatch(AA *a, AA *b, int size) /* Compare two peptides aa by aa. */ { return dnaOrAaScoreMatch(a, b, size, 2, -1, 'X'); } void writeSeqWithBreaks(FILE *f, char *letters, int letterCount, int maxPerLine) /* Write out letters with newlines every maxLine. */ { int lettersLeft = letterCount; int lineSize; while (lettersLeft > 0) { lineSize = lettersLeft; if (lineSize > maxPerLine) lineSize = maxPerLine; mustWrite(f, letters, lineSize); fputc('\n', f); letters += lineSize; lettersLeft -= lineSize; } } static int findTailPolyAMaybeMask(DNA *dna, int size, boolean doMask, boolean loose) /* Identify PolyA at end; mask to 'n' if specified. This allows a few * non-A's as noise to be trimmed too. Returns number of bases trimmed. * Leaves first two bases of PolyA in case there's a taa stop codon. */ { int i; int score = 10; int bestScore = 10; int bestPos = -1; int trimSize = 0; for (i=size-1; i>=0; --i) { DNA b = dna[i]; if (b == 'n' || b == 'N') continue; if (score > 20) score = 20; if (b == 'a' || b == 'A') { score += 1; if (score >= bestScore) { bestScore = score; bestPos = i; } else if (loose && score >= (bestScore - 8)) { /* If loose, keep extending even if score isn't back up to best. */ bestPos = i; } } else { score -= 10; } if (score < 0) { break; } } if (bestPos >= 0) { trimSize = size - bestPos - 2; // Leave two for aa in taa stop codon if (trimSize > 0) { if (doMask) for (i=size - trimSize; i<size; ++i) dna[i] = 'n'; } else trimSize = 0; } return trimSize; } int tailPolyASizeLoose(DNA *dna, int size) /* Return size of PolyA at end (if present). This allows a few non-A's as * noise to be trimmed too, but skips first two aa for taa stop codon. * It is less conservative in extending the polyA region than maskTailPolyA. */ { return findTailPolyAMaybeMask(dna, size, FALSE, TRUE); } int maskTailPolyA(DNA *dna, int size) /* Convert PolyA at end to n. This allows a few non-A's as noise to be * trimmed too. Returns number of bases trimmed. Leaves very last a. */ { return findTailPolyAMaybeMask(dna, size, TRUE, FALSE); } static int findHeadPolyTMaybeMask(DNA *dna, int size, boolean doMask, boolean loose) /* Return size of PolyT at start (if present); mask to 'n' if specified. * This allows a few non-T's as noise to be trimmed too, but skips last * two tt for revcomp'd taa stop codon. */ { int i; int score = 10; int bestScore = 10; int bestPos = -1; int trimSize = 0; for (i=0; i<size; ++i) { DNA b = dna[i]; if (b == 'n' || b == 'N') continue; if (score > 20) score = 20; if (b == 't' || b == 'T') { score += 1; if (score >= bestScore) { bestScore = score; bestPos = i; } else if (loose && score >= (bestScore - 8)) { /* If loose, keep extending even if score isn't back up to best. */ bestPos = i; } } else { score -= 10; } if (score < 0) { break; } } if (bestPos >= 0) { trimSize = bestPos+1 - 2; // Leave two for aa in taa stop codon if (trimSize > 0) { if (doMask) memset(dna, 'n', trimSize); } else trimSize = 0; } return trimSize; } int headPolyTSizeLoose(DNA *dna, int size) /* Return size of PolyT at start (if present). This allows a few non-T's as * noise to be trimmed too, but skips last two tt for revcomp'd taa stop * codon. * It is less conservative in extending the polyA region than maskHeadPolyT. */ { return findHeadPolyTMaybeMask(dna, size, FALSE, TRUE); } int maskHeadPolyT(DNA *dna, int size) /* Convert PolyT at start. This allows a few non-T's as noise to be * trimmed too. Returns number of bases trimmed. */ { return findHeadPolyTMaybeMask(dna, size, TRUE, FALSE); } boolean isDna(char *poly, int size) /* Return TRUE if letters in poly are at least 90% ACGTNU- */ { int i; int dnaCount = 0; dnaUtilOpen(); for (i=0; i<size; ++i) { if (ntChars[(int)poly[i]]) dnaCount += 1; } return (dnaCount >= round(0.9 * size)); } boolean isAllNt(char *seq, int size) /* Return TRUE if all letters in seq are ACGTNU-. */ { int i; dnaUtilOpen(); for (i = 0; i < size; ++i) { if (ntChars[(int)seq[i]] == 0) return FALSE; } return TRUE; } /* Tables to convert from 0-20 to ascii single letter representation * of proteins. */ int aaVal[256]; AA valToAa[21]; AA aaChars[256]; /* 0 except for value aa characters. Converts to upper case rest. */ struct aminoAcidTable /* A little info about each amino acid. */ { int ix; char letter; char abbreviation[3]; char *name; }; struct aminoAcidTable aminoAcidTable[] = { {0, 'A', "ala", "alanine"}, {1, 'C', "cys", "cysteine"}, {2, 'D', "asp", "aspartic acid"}, {3, 'E', "glu", "glutamic acid"}, {4, 'F', "phe", "phenylalanine"}, {5, 'G', "gly", "glycine"}, {6, 'H', "his", "histidine"}, {7, 'I', "ile", "isoleucine"}, {8, 'K', "lys", "lysine"}, {9, 'L', "leu", "leucine"}, {10, 'M', "met", "methionine"}, {11, 'N', "asn", "asparagine"}, {12, 'P', "pro", "proline"}, {13, 'Q', "gln", "glutamine"}, {14, 'R', "arg", "arginine"}, {15, 'S', "ser", "serine"}, {16, 'T', "thr", "threonine"}, {17, 'V', "val", "valine"}, {18, 'W', "trp", "tryptophan"}, {19, 'Y', "tyr", "tyrosine"}, {20, 'X', "ter", "termination"}, }; char *aaAbbr(int i) /* return pointer to AA abbrevation */ { return(aminoAcidTable[i].abbreviation); } char aaLetter(int i) /* return AA letter */ { return(aminoAcidTable[i].letter); } static void initAaVal() /* Initialize aaVal and valToAa tables. */ { int i; char c, lowc; for (i=0; i<ArraySize(aaVal); ++i) aaVal[i] = -1; for (i=0; i<ArraySize(aminoAcidTable); ++i) { c = aminoAcidTable[i].letter; lowc = tolower(c); aaVal[(int)c] = aaVal[(int)lowc] = i; aaChars[(int)c] = aaChars[(int)lowc] = c; valToAa[i] = c; } aaChars['x'] = aaChars['X'] = 'X'; } void dnaUtilOpen() /* Initialize stuff herein. */ { static boolean opened = FALSE; if (!opened) { checkSizeTypes(); initNtVal(); initAaVal(); initNtChars(); initNtMixedCaseChars(); initNtCompTable(); opened = TRUE; } } char aaAbbrToLetter(char *abbr) /* Convert an AA abbreviation such as "Ala", "Asp" etc., to its single letter code * such as "A", "D" etc. Return the null char '\0' if abbr is not found. */ { // Lowercase for comparison. char abbrLC[4]; safencpy(abbrLC, sizeof(abbrLC), abbr, 3); toLowerN(abbrLC, 3); int ix; for (ix = 0; ix < ArraySize(aminoAcidTable); ix++) { if (sameStringN(abbrLC, aminoAcidTable[ix].abbreviation, 3)) return aminoAcidTable[ix].letter; } return '\0'; } void aaToAbbr(char aa, char *abbrBuf, size_t abbrBufSize) /* Convert an AA single letter such as "A", "D" etc. to its abbreviation such as "Ala", "Asp" etc. * abbrBufSize must be at least 4. If aa is not found, "?%c?",aa is written into abbrBuf. */ { // Uppercase for comparison. char aaUC = toupper(aa); int ix; for (ix = 0; ix < ArraySize(aminoAcidTable); ix++) { if (aaUC == aminoAcidTable[ix].letter) { // safencpy(...3) is required here because aminoAcidTable.abbreviation is char[3] not [4] safencpy(abbrBuf, abbrBufSize, aminoAcidTable[ix].abbreviation, 3); abbrBuf[0] = toupper(abbrBuf[0]); return; } } safef(abbrBuf, abbrBufSize, "?%c?", aa); } +char *aaToName(char aa) +/* Convert an AA single letter such as 'A', 'D' etc. to its full name such as "alanine", + * "aspartic acid" etc. Returns NULL if aa is not found. */ +{ +char aaUC = toupper(aa); +int ix; +for (ix = 0; ix < ArraySize(aminoAcidTable); ix++) + { + if (aaUC == aminoAcidTable[ix].letter) + return aminoAcidTable[ix].name; + } +return NULL; +} + void trimRefAltDir(char *ref, char *alt, uint *pStart, uint *pEnd, int *pRefLen, int *pAltLen, boolean leftJustify) /* If ref and alt have identical bases at beginning and/or end, trim those & update all params. */ { int trimStart = 0, trimEnd = 0; int refLen = strlen(ref); int altLen = strlen(alt); if (leftJustify) { // first trim end, then start while (refLen > 0 && altLen > 0 && ref[refLen-1] == alt[altLen-1]) { (*pEnd)--; refLen--; altLen--; ref[refLen] = 0; alt[altLen] = 0; } while (ref[trimStart] != '\0' && alt[trimStart] != '\0' && ref[trimStart] == alt[trimStart]) trimStart++; } else { // first trim start, then end while (ref[trimStart] != '\0' && alt[trimStart] != '\0' && ref[trimStart] == alt[trimStart]) trimStart++; int iR = refLen - 1, iA = altLen - 1; while (iR >= trimStart && iA >= trimStart && ref[iR] == alt[iA]) { iR--; iA--; trimEnd++; } } if (trimEnd) { *pEnd -= trimEnd; refLen -= trimEnd; altLen -= trimEnd; ref[refLen] = '\0'; alt[altLen] = '\0'; } if (trimStart) { *pStart += trimStart; refLen -= trimStart; altLen -= trimStart; memmove(ref, ref+trimStart, refLen+1); memmove(alt, alt+trimStart, altLen+1); } *pRefLen = refLen; *pAltLen = altLen; } void trimRefAlt(char *ref, char *alt, uint *pStart, uint *pEnd, int *pRefLen, int *pAltLen) /* If ref and alt have identical bases at beginning and/or end, trim those & update all params. */ { trimRefAltDir(ref, alt, pStart, pEnd, pRefLen, pAltLen, FALSE); } void trimRefAltLeft(char *ref, char *alt, uint *pStart, uint *pEnd, int *pRefLen, int *pAltLen) /* If ref and alt have identical bases at beginning and/or end, trim those starting on the right * so we get the leftmost representation & update all params. */ { trimRefAltDir(ref, alt, pStart, pEnd, pRefLen, pAltLen, TRUE); }