");
+ printf("");
printf(" ACTIONS QUERY SCORE START END QSIZE IDENTITY CHRO STRAND START END SPAN\n");
printf("---------------------------------------------------------------------------------------------------\n");
for (psl = pslList; psl != NULL; psl = psl->next)
{
if (customText)
printf("",
browserUrl, psl->tName, psl->tStart + 1, psl->tEnd, database,
customText, uiState, unhideTrack);
else
printf("",
browserUrl, psl->tName, psl->tStart + 1, psl->tEnd, database,
pslName, faName, uiState, unhideTrack);
printf("browser ");
printf("",
hgcUrl, psl->tStart, pslName, cgiEncode(faName), psl->qName, psl->tName,
psl->tStart, psl->tEnd, database, uiState);
printf("details ");
printf("%-14s %5d %5d %5d %5d %5.1f%% %4s %2s %9d %9d %6d\n",
psl->qName, pslScore(psl), psl->qStart+1, psl->qEnd, psl->qSize,
100.0 - pslCalcMilliBad(psl, TRUE) * 0.1,
skipChr(psl->tName), psl->strand, psl->tStart+1, psl->tEnd,
psl->tEnd - psl->tStart);
}
printf("\n");
puts("Missing a match?
\n");
puts("\n");
}
pslFreeList(&pslList);
}
void trimUniq(bioSeq *seqList)
/* Check that all seq's in list have a unique name. Try and
* abbreviate longer sequence names. */
{
struct hash *hash = newHash(0);
bioSeq *seq;
for (seq = seqList; seq != NULL; seq = seq->next)
{
char *saferString = needMem(strlen(seq->name)+1);
char *c, *s;
/* Some chars are safe to allow through, other chars cause
* problems. It isn't necessarily a URL safe string that is
* being calculated here. The original problem was a user had
* the fasta header line of:
* chr8|59823648:59825047|+
* The plus sign was being taken as the query name and this
* created problems as that name was passed on to hgc via
* the ss cart variable. The + sign became part of a URL
* eventually. This loop allows only isalnum and =_/.:;_|
* to get through as part of the header name. These characters
* all proved to be safe as single character names, or all
* together.
*/
s = saferString;
for (c = seq->name; *c != '\0'; ++c)
{
if (c && (*c != '\0'))
{
if ( isalnum(*c) || (*c == '=') || (*c == '-') || (*c == '/') ||
(*c == '.') || (*c == ':') || (*c == ';') || (*c == '_') ||
(*c == '|') )
*s++ = *c;
}
}
*s = '\0';
freeMem(seq->name);
if (*saferString == '\0')
{
freeMem(saferString);
saferString = cloneString("YourSeq");
}
seq->name = saferString;
if (strlen(seq->name) > 14) /* Try and get rid of long NCBI .fa cruft. */
{
char *nameClone = NULL;
char *abbrv = NULL;
char *words[32];
int wordCount;
boolean isEns = (stringIn("ENSEMBL:", seq->name) != NULL);
nameClone = cloneString(seq->name);
wordCount = chopString(nameClone, "|", words, ArraySize(words));
if (wordCount > 1) /* Looks like it's an Ensembl/NCBI
* long name alright. */
{
if (isEns)
{
abbrv = words[0];
if (abbrv[0] == 0) abbrv = words[1];
}
else if (sameString(words[1], "dbSNP"))
{
if (wordCount > 2)
abbrv = words[2];
else
abbrv = nameClone;
}
else
{
abbrv = words[wordCount-1];
if (abbrv[0] == 0) abbrv = words[wordCount-2];
}
if (hashLookup(hash, abbrv) == NULL)
{
freeMem(seq->name);
seq->name = cloneString(abbrv);
}
freez(&nameClone);
}
}
hashAddUnique(hash, seq->name, hash);
}
freeHash(&hash);
}
int realSeqSize(bioSeq *seq, boolean isDna)
/* Return size of sequence without N's or (for proteins)
* X's. */
{
char unknown = (isDna ? 'n' : 'X');
int i, size = seq->size, count = 0;
char *s = seq->dna;
for (i=0; inext)
subChar(seq->dna, 'u', 't');
}
static struct slName *namesInPsl(struct psl *psl)
/* Find all the unique names in a list of psls. */
{
struct hash *hash = newHash(3);
struct slName *nameList = NULL;
struct slName *name;
for(; psl; psl = psl->next)
{
if (hashLookup(hash, psl->qName) == NULL)
{
name = slNameNew(psl->qName);
slAddHead(&nameList, name);
hashStore(hash, psl->qName);
}
}
slReverse(&nameList);
return nameList;
}
static char *makeNameUnique(char *name, char *database, struct cart *cart)
/* Make sure track name will create a unique custom track. */
{
struct slName *browserLines = NULL;
struct customTrack *ctList = customTracksParseCart(database, cart, &browserLines, NULL);
struct customTrack *ct;
int count = 0;
char buffer[4096];
safef(buffer, sizeof buffer, "%s", name);
for(;;count++)
{
char *customName = customTrackTableFromLabel(buffer);
for (ct=ctList;
ct != NULL;
ct=ct->next)
{
if (startsWith(customName, ct->tdb->track))
// Found a track with this name.
break;
}
if (ct == NULL)
break;
safef(buffer, sizeof buffer, "%s (%d)", name, count + 1);
}
return cloneString(buffer);
}
static void getCustomName(char *database, struct cart *cart, struct psl *psl, char **pName, char **pDescription)
// Find a track name that isn't currently a custom track. Also fill in description.
{
struct slName *names = namesInPsl(psl);
char shortName[4096];
char description[4096];
unsigned count = slCount(names);
if (count == 1)
{
safef(shortName, sizeof shortName, "blat %s", names->name);
safef(description, sizeof description, "blat on %s", names->name);
}
else if (count == 2)
{
safef(shortName, sizeof shortName, "blat %s+%d", names->name, count - 1);
safef(description, sizeof description, "blat on %d queries (%s, %s)", count, names->name, names->next->name);
}
else
{
safef(shortName, sizeof shortName, "blat %s+%d", names->name, count - 1);
safef(description, sizeof description, "blat on %d queries (%s, %s, ...)", count, names->name, names->next->name);
}
*pName = makeNameUnique(shortName, database, cart);
*pDescription = cloneString(description);
}
void blatSeq(char *userSeq, char *organism, char *database)
/* Blat sequence user pasted in. */
{
FILE *f;
struct dnaSeq *seqList = NULL, *seq;
struct tempName pslTn, faTn;
int maxSingleSize, maxTotalSize, maxSeqCount;
int minSingleSize = minMatchShown;
char *genome, *db;
char *type = cgiString("type");
char *seqLetters = cloneString(userSeq);
struct serverTable *serve;
int conn;
int oneSize, totalSize = 0, seqCount = 0;
boolean isTx = FALSE;
boolean isTxTx = FALSE;
boolean txTxBoth = FALSE;
struct gfOutput *gvo;
boolean qIsProt = FALSE;
enum gfType qType, tType;
struct hash *tFileCache = gfFileCacheNew();
boolean feelingLucky = cgiBoolean("Lucky");
getDbAndGenome(cart, &db, &genome, oldVars);
if(!feelingLucky)
cartWebStart(cart, db, "%s BLAT Results", trackHubSkipHubName(organism));
/* Load user sequence and figure out if it is DNA or protein. */
if (sameWord(type, "DNA"))
{
seqList = faSeqListFromMemText(seqLetters, TRUE);
uToT(seqList);
isTx = FALSE;
}
else if (sameWord(type, "translated RNA") || sameWord(type, "translated DNA"))
{
seqList = faSeqListFromMemText(seqLetters, TRUE);
uToT(seqList);
isTx = TRUE;
isTxTx = TRUE;
txTxBoth = sameWord(type, "translated DNA");
}
else if (sameWord(type, "protein"))
{
seqList = faSeqListFromMemText(seqLetters, FALSE);
isTx = TRUE;
qIsProt = TRUE;
}
else
{
seqList = faSeqListFromMemTextRaw(seqLetters);
isTx = !seqIsDna(seqList);
if (!isTx)
{
for (seq = seqList; seq != NULL; seq = seq->next)
{
seq->size = dnaFilteredSize(seq->dna);
dnaFilter(seq->dna, seq->dna);
toLowerN(seq->dna, seq->size);
subChar(seq->dna, 'u', 't');
}
}
else
{
for (seq = seqList; seq != NULL; seq = seq->next)
{
seq->size = aaFilteredSize(seq->dna);
aaFilter(seq->dna, seq->dna);
toUpperN(seq->dna, seq->size);
}
qIsProt = TRUE;
}
}
if (seqList != NULL && seqList->name[0] == 0)
{
freeMem(seqList->name);
seqList->name = cloneString("YourSeq");
}
trimUniq(seqList);
/* If feeling lucky only do the first one. */
if(feelingLucky && seqList != NULL)
{
seqList->next = NULL;
}
/* Figure out size allowed. */
maxSingleSize = (isTx ? 10000 : 75000);
maxTotalSize = maxSingleSize * 2.5;
#ifdef LOWELAB
maxSeqCount = 200;
#else
maxSeqCount = 25;
#endif
/* Create temporary file to store sequence. */
trashDirFile(&faTn, "hgSs", "hgSs", ".fa");
faWriteAll(faTn.forCgi, seqList);
/* Create a temporary .psl file with the alignments against genome. */
trashDirFile(&pslTn, "hgSs", "hgSs", ".pslx");
f = mustOpen(pslTn.forCgi, "w");
gvo = gfOutputPsl(0, qIsProt, FALSE, f, FALSE, TRUE);
serve = findServer(db, isTx);
/* Write header for extended (possibly protein) psl file. */
if (isTx)
{
if (isTxTx)
{
qType = gftDnaX;
tType = gftDnaX;
}
else
{
qType = gftProt;
tType = gftDnaX;
}
}
else
{
qType = gftDna;
tType = gftDna;
}
pslxWriteHead(f, qType, tType);
if (qType == gftProt)
{
minSingleSize = 14;
}
else if (qType == gftDnaX)
{
minSingleSize = 36;
}
/* Loop through each sequence. */
for (seq = seqList; seq != NULL; seq = seq->next)
{
printf(" "); fflush(stdout); /* prevent apache cgi timeout by outputting something */
oneSize = realSeqSize(seq, !isTx);
if ((seqCount&1) == 0) // Call bot delay every 2nd time starting with first time
hgBotDelay();
if (++seqCount > maxSeqCount)
{
warn("More than 25 input sequences, stopping at %s.",
seq->name);
break;
}
if (oneSize > maxSingleSize)
{
warn("Sequence %s is %d letters long (max is %d), skipping",
seq->name, oneSize, maxSingleSize);
continue;
}
if (oneSize < minSingleSize)
{
warn("Warning: Sequence %s is only %d letters long (%d is the recommended minimum)",
seq->name, oneSize, minSingleSize);
// we could use "continue;" here to actually enforce skipping,
// but let's give the short sequence a chance, it might work.
// minimum possible length = tileSize+stepSize, so mpl=16 for dna stepSize=5, mpl=10 for protein.
if (qIsProt && oneSize < 1) // protein does not tolerate oneSize==0
continue;
}
totalSize += oneSize;
if (totalSize > maxTotalSize)
{
warn("Sequence %s would take us over the %d letter limit, stopping here.",
seq->name, maxTotalSize);
break;
}
conn = gfConnect(serve->host, serve->port);
if (isTx)
{
gvo->reportTargetStrand = TRUE;
if (isTxTx)
{
gfAlignTransTrans(&conn, serve->nibDir, seq, FALSE, 5,
tFileCache, gvo, !txTxBoth);
if (txTxBoth)
{
reverseComplement(seq->dna, seq->size);
conn = gfConnect(serve->host, serve->port);
gfAlignTransTrans(&conn, serve->nibDir, seq, TRUE, 5,
tFileCache, gvo, FALSE);
}
}
else
{
gfAlignTrans(&conn, serve->nibDir, seq, 5, tFileCache, gvo);
}
}
else
{
gfAlignStrand(&conn, serve->nibDir, seq, FALSE, minMatchShown, tFileCache, gvo);
reverseComplement(seq->dna, seq->size);
conn = gfConnect(serve->host, serve->port);
gfAlignStrand(&conn, serve->nibDir, seq, TRUE, minMatchShown, tFileCache, gvo);
}
gfOutputQuery(gvo, f);
}
carefulClose(&f);
showAliPlaces(pslTn.forCgi, faTn.forCgi, NULL, serve->db, qType, tType,
organism, feelingLucky);
if(!feelingLucky)
cartWebEnd();
gfFileCacheFree(&tFileCache);
}
void askForSeq(char *organism, char *db)
/* Put up a little form that asks for sequence.
* Call self.... */
{
/* ignore struct serverTable* return, but can error out if not found */
findServer(db, FALSE);
/* JavaScript to update form when org changes */
char *onChangeText = ""
"document.mainForm.changeInfo.value='orgChange';"
"document.mainForm.submit();";
char *userSeq = NULL;
printf(
"\n");
webNewSection("About BLAT");
printf(
"BLAT on DNA is designed to\n"
"quickly find sequences of 95%% and greater similarity of length 25 bases or\n"
"more. It may miss more divergent or shorter sequence alignments. It will find\n"
"perfect sequence matches of 20 bases.\n"
"BLAT on proteins finds sequences of 80%% and greater similarity of length 20 amino\n"
"acids or more. In practice DNA BLAT works well on primates, and protein\n"
"blat on land vertebrates."
);
printf("%s",
"\n
BLAT is not BLAST. DNA BLAT works by keeping an index of the entire genome\n"
"in memory. The index consists of all overlapping 11-mers stepping by 5 except for\n"
"those heavily involved in repeats. The index takes up about\n"
"2 gigabytes of RAM. RAM can be further reduced to less than 1 GB by increasing step size to 11.\n"
"The genome itself is not kept in memory, allowing\n"
"BLAT to deliver high performance on a reasonably priced Linux box.\n"
"The index is used to find areas of probable homology, which are then\n"
"loaded into memory for a detailed alignment. Protein BLAT works in a similar\n"
"manner, except with 4-mers rather than 11-mers. The protein index takes a little\n"
"more than 2 gigabytes.
\n"
"BLAT was written by Jim Kent.\n"
"Like most of Jim's software, interactive use on this web server is free to all.\n"
"Sources and executables to run batch jobs on your own server are available free\n"
"for academic, personal, and non-profit purposes. Non-exclusive commercial\n"
"licenses are also available. See the \n"
"Kent Informatics\n"
"website for details.
\n"
"\n"
"For more information on the graphical version of BLAT, click the Help \n"
"button on the top menu bar");
if (hIsGsidServer())
printf(".
\n");
else
printf(" or see the Genome Browser FAQ. \n");
}
void doMiddle(struct cart *theCart)
/* Write header and body of html page. */
{
char *userSeq;
char *db, *organism;
boolean clearUserSeq = cgiBoolean("Clear");
cart = theCart;
dnaUtilOpen();
orgChange = sameOk(cgiOptionalString("changeInfo"),"orgChange");
if (orgChange)
{
cgiVarSet("db", hDefaultDbForGenome(cgiOptionalString("org")));
}
getDbAndGenome(cart, &db, &organism, oldVars);
char *oldDb = cloneString(db);
findClosestServer(&db, &organism);
/* Get sequence - from userSeq variable, or if
* that is empty from a file. */
if (clearUserSeq)
{
cartSetString(cart, "userSeq", "");
cartSetString(cart, "seqFile", "");
}
userSeq = cartUsualString(cart, "userSeq", "");
if (isEmpty(userSeq))
{
userSeq = cartOptionalString(cart, "seqFile");
}
if (isEmpty(userSeq) || orgChange)
{
cartWebStart(theCart, db, "%s BLAT Search", trackHubSkipHubName(organism));
if (differentString(oldDb, db))
printf("
Note: BLAT search is not available for %s %s; "
"defaulting to %s %s
\n",
hGenome(oldDb), hFreezeDate(oldDb), organism, hFreezeDate(db));
askForSeq(organism,db);
cartWebEnd();
}
else
{
blatSeq(skipLeadingSpaces(userSeq), organism, db);
}
}
/* Null terminated list of CGI Variables we don't want to save
* permanently. */
char *excludeVars[] = {"Submit", "submit", "Clear", "Lucky", "type", "userSeq", "seqFile", "showPage", "changeInfo", NULL};
int main(int argc, char *argv[])
/* Process command line. */
{
long enteredMainTime = clock1000();
oldVars = hashNew(10);
cgiSpoof(&argc, argv);
/* org has precedence over db when changeInfo='orgChange' */
cartEmptyShell(doMiddle, hUserCookie(), excludeVars, oldVars);
cgiExitTime("hgBlat", enteredMainTime);
return 0;
}