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/mouseStuff/chainInfo/chainInfo.c src/hg/mouseStuff/chainInfo/chainInfo.c
index 7133b09..4ca70a7 100644
--- src/hg/mouseStuff/chainInfo/chainInfo.c
+++ src/hg/mouseStuff/chainInfo/chainInfo.c
@@ -1,595 +1,595 @@
/* Copyright (C) 2005 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 "hCommon.h"
#include "linefile.h"
#include "hash.h"
#include "options.h"
//#include "jksql.h"
#include "hdb.h"
#include "dnautil.h"
#include "chain.h"
#include "chainNet.h"
#include "axt.h"
#include "bed.h"
#include "axtInfo.h"
#include "rmskOut.h"
#include "binRange.h"
#include "obscure.h"
#include "genePred.h"
int minGap = 60; /* minimum intron size */
float minOverlap = 0.60;
int bigChain = 1000000; /* max chain size to consider */
static struct hash *chainHash = NULL; /* hash of all chains indexed by chainId*/
struct intronChain {
/* scoring of aligned introns in a chain */
struct intronChain *next; /* Next in list. */
struct chain *chain; /* scored chain */
int qIntron ; /* count of introns on query side */
int tIntron ; /* count of introns on target side */
int tGap ; /* count of bases in gaps on target */
int qGap ; /* count of bases in gaps on query */
int tReps; /* total number of repeats in target */
int qReps; /* total number of repeats in target */
int total; /* total size of chain */
float score1; /* score2 * ratio of bases in target/query gaps */
float score2; /* ratio of introns target/query , negative if query has more */
int overlap1; /* bases that overlap with gene */
int overlap2;
char strand; /* strand where putative pseudogene is */
};
void usage()
/* Explain usage and exit. */
{
errAbort(
"chainInfo - parse chains to determine number of introns on query and target\n"
"usage:\n"
" chainInfo query_db target_chrom chainFile syntenicNet.bed target.lst query.lst target_repeatMasker.out query_repeatDir output refGene.tab genePred2.tab genePred3.tab\n"
"options:\n"
" -minGap=N minimum size to be considered a gap (intron). Default is 60.\n"
" -minOverlap=0.N threshold for repeat overlap with gap to be considered a repeat instead of intron. Default is 0.60.\n"
"db is query database to match genes to pseudogene\n"
"chrom is target chromosome to match genes to pseudogene\n"
"syntenicNet.bed is a filter net file using netFilter -syn converted to bed format\n"
"query_repeatDir is the dirctory containing repeatMasker out files chrXX.fa.out\n"
"target.lst, query.lst is a list of chromosome names followed by size of each chromosome\n"
"refGene.tab, genePred2,3 are bed12 files of genes in query database\n"
);
}
int qStart(struct chain *chain)
/* convert coordinates */
{
return (chain->qStrand == '+') ? chain->qStart : chain->qSize-chain->qEnd;
}
int qEnd(struct chain *chain)
/* convert coordinates */
{
return (chain->qStrand == '+') ? chain->qEnd : chain->qSize-chain->qStart ;
}
int qStartBlk(struct chain *chain, int start, int end)
/* convert coordinates */
{
return (chain->qStrand == '+') ? start : chain->qSize-end;
}
int qEndBlk(struct chain *chain, int start, int end)
/* convert coordinates */
{
return (chain->qStrand == '+') ? end : (chain->qSize)-start ;
}
struct hash *readSizes(char *fileName)
/* Read tab-separated file into hash with
* name key size value. */
{
struct lineFile *lf = lineFileOpen(fileName, TRUE);
struct hash *hash = newHash(0);
char *row[2];
while (lineFileRow(lf, row))
{
char *name = row[0];
int size = lineFileNeedNum(lf, row, 1);
if (hashLookup(hash, name) != NULL)
warn("Duplicate %s, ignoring all but first\n", name);
else
hashAdd(hash, name, intToPt(size));
}
lineFileClose(&lf);
return hash;
}
struct intronChain *scoreNextAlignedIntron(struct chain *chain, FILE *f, struct binKeeper *bk, struct hash *bkHash)
/* get next chain from sorted list. Return NULL at end of list. attempt to chain chains together if possible
* count introns on query and target side , blocks with < minGap bases between will be combined.
* Note that chain block scores are not filled in by this .
* repeats are filtered out and are passed in the bk structure
* score 2 = ratio of intron target/query; score1 = score2 * ratio bases in target/query gaps */
{
char *row[13];
int wordCount;
int q = qStart(chain),t = chain->tStart;
static int id = 0;
struct binElement *el;
struct rmskOut *rmsk ;
struct intronChain *iChain;
struct cBlock *b, *nextB;
int tGapStart = 0;
int qGapStart = 0;
int tGapEnd = 0;
int qGapEnd = 0, tg,qg;
int tReps, qReps;
float percentOverlap, percentQOverlap;
struct binKeeper *bkQ;
AllocVar(iChain);
iChain->chain = chain;
iChain->qIntron = 0; iChain->tIntron = 0;
iChain->tGap = 0; iChain->qGap = 0;
iChain->tReps = 0; iChain->qReps = 0; iChain->total = 0;
iChain->strand = ' ';
for (b = chain->blockList; b != NULL; b = nextB)
{
nextB = b->next;
if (nextB == NULL)
break;
tGapStart = b->tEnd;
tGapEnd = nextB->tStart;
if (chain->qStrand == '+')
{
qGapStart = b->qEnd;
qGapEnd = nextB->qStart;
}
else
{
qGapStart = chain->qSize-nextB->qStart;
qGapEnd = chain->qSize-b->qEnd;
}
tg = tGapEnd - tGapStart;
qg = qGapEnd - qGapStart;
t += tg;
q += qg;
/* remove repeasts so we don't call them 'introns' on target*/
tReps = 0;
for (el = binKeeperFindSorted(bk, tGapStart, tGapEnd); el != NULL; el = el->next)
{
rmsk = el->val;
tReps += positiveRangeIntersection(tGapStart, tGapEnd, rmsk->genoStart, rmsk->genoEnd);
}
percentOverlap = (float)tReps / (float)(tGapEnd - tGapStart);
if (((tGapEnd - tGapStart) > 0) && percentOverlap > minOverlap)
{
if (tReps > tg) tReps=tg;
tg -= tReps;
if (tg < 0)
assert(tg >=0);
iChain->tReps += tReps;
}
/* remove repeasts so we don't call them 'introns' on query */
/* get bk for query */
qReps = 0;
if (bkHash != NULL)
{
bkQ = hashFindVal(bkHash, chain->qName);
for (el = binKeeperFindSorted(bkQ, qGapStart, qGapEnd); el != NULL; el = el->next)
{
rmsk = el->val;
qReps += positiveRangeIntersection(qGapStart, qGapEnd, rmsk->genoStart, rmsk->genoEnd);
}
}
percentQOverlap = (float)qReps / (float)(qGapEnd - qGapStart);
if (((qGapEnd - qGapStart) > 0) && percentQOverlap > minOverlap)
{
if (qReps > qg) qReps=qg;
qg -= qReps;
if (qg < 0)
assert(qg >= 0);
iChain->qReps += qReps;
}
/* if (chain->qStrand == '+')
* printf("\n---%d %d %3.2f %d gap=%s:%d-%d q=%s:%d-%d %d",tg, tReps , percentOverlap, tGapEnd - tGapStart, chain->tName,
* tGapStart, tGapEnd, chain->qName, b->qStart, b->qEnd , chain->id);
* else
* printf("\n---%d %d %3.2f %d gap=%s:%d-%d q=%s:%d-%d %d",tg, tReps , percentOverlap, tGapEnd - tGapStart, chain->tName,
* tGapStart, tGapEnd, chain->qName, chain->qSize-b->qEnd, chain->qSize-b->qStart , chain->id);
*/
iChain->tGap += tg;
if(iChain->tGap <0)
assert(iChain->tGap >= 0);
iChain->qGap += qg;
if(iChain->qGap <0)
assert(iChain->qGap >= 0);
if (tg >= minGap)
iChain->tIntron++;
if (qg >= minGap)
iChain->qIntron++;
}
/*
if (q != chain->qEnd)
errAbort("q end mismatch %d vs %d in chain %d \n",
q, chain->qEnd, chain->id);
if (t != chain->tEnd)
errAbort("t end mismatch %d vs %d in chain %d \n",
t, chain->tEnd,chain->id );
*/
//slReverse(&chain->blockList);
assert(iChain->tIntron>=0);
assert(iChain->qIntron>=0);
iChain->score2 = iChain->qIntron > iChain->tIntron ? (iChain->qIntron+1.0)/(iChain->tIntron+1.0) : -(iChain->tIntron+1.0)/(iChain->qIntron+1.0) ;
iChain->score1 = (iChain->qIntron > iChain->tIntron) ? ((iChain->qGap+1)/(iChain->tGap+1)) : -(iChain->tGap+1)/(iChain->qGap+1);
return iChain;
}
#ifdef DUPE
struct binKeeper *readRepeats(char *chrom, char *rmskFileName, struct hash *tSizeHash)
/* read all repeats for a chromosome of size size, returns results in binKeeper structure for fast query*/
{
boolean rmskRet;
struct lineFile *rmskF = NULL;
struct rmskOut *rmsk;
struct binKeeper *bk;
int size;
size = hashIntVal(tSizeHash, chrom);
bk = binKeeperNew(0, size);
assert(size > 1);
rmskOutOpenVerify(rmskFileName ,&rmskF , &rmskRet);
while ((rmsk = rmskOutReadNext(rmskF)) != NULL)
{
binKeeperAdd(bk, rmsk->genoStart, rmsk->genoEnd, rmsk);
}
lineFileClose(&rmskF);
return bk;
}
struct hash *oldreadRepeatsAll(char *chrom, char *sizeFileName, char *rmskDir)
/* read all repeats for a all chromosome of size size, returns results in hash of binKeeper structure for fast query*/
{
boolean rmskRet;
struct binKeeper *bk;
int size;
struct lineFile *rmskF = NULL;
struct rmskOut *rmsk;
struct lineFile *lf = lineFileOpen(sizeFileName, TRUE);
struct hash *hash = newHash(0);
char *row[2];
char rmskFileName[256];
while (lineFileRow(lf, row))
{
char *name = row[0];
int size = lineFileNeedNum(lf, row, 1);
// char *repSubDir = cloneString(skipChr(name));
// chopSuffixAt(repSubDir,'_');
if (hashLookup(hash, name) != NULL)
warn("Duplicate %s, ignoring all but first\n", name);
else
{
bk = binKeeperNew(0, size);
assert(size > 1);
sprintf(rmskFileName, "%s/%s.fa.out",rmskDir,name);
rmskOutOpenVerify(rmskFileName ,&rmskF , &rmskRet);
while ((rmsk = rmskOutReadNext(rmskF)) != NULL)
{
binKeeperAdd(bk, rmsk->genoStart, rmsk->genoEnd, rmsk);
}
lineFileClose(&rmskF);
hashAdd(hash, name, bk);
}
}
lineFileClose(&lf);
return hash;
}
#endif
struct chain *readSortChains(char *fileName)
{
struct lineFile *lf = lineFileOpen(fileName, TRUE);
struct chain *chainList = NULL, *chain;
char chainId[128];
chainHash = newHash(0);
while ((chain = chainRead(lf)) != NULL)
{
slAddHead(&chainList, chain);
safef(chainId, sizeof(chainId), "%d", chain->id);
hashAddUnique(chainHash, chainId, chain);
}
lineFileClose(&lf);
/* Sort. */
printf("Sorting Chains\n");
slSort(&chainList, chainCmpQuery);
return chainList;
}
struct chain *chainFromId(int id)
/** Return a chain given the id. */
{
char key[128];
struct chain *chain = NULL;
if(chainHash == NULL)
errAbort("chainFromId() - no chains found");
safef(key, sizeof(key), "%d", id);
chain = hashFindVal(chainHash, key);
if(chain == NULL)
warn("Chain not found for id: %d", id);
return chain;
}
boolean checkChain(struct chain *chain, int start, int end)
/* Return true if chain covers start, end, false otherwise. */
{
struct chain *subChain=NULL, *toFree=NULL;
boolean good = FALSE;
if (chain == NULL) return good;
chainSubsetOnT(chain, start, end, &subChain, &toFree);
if(subChain != NULL)
good = TRUE;
chainFree(&toFree);
return good;
}
struct chain *lookForChain(struct cnFill *list, int start, int end)
/* Recursively look for a chain in this list containing the coordinates
desired. */
{
struct cnFill *fill=NULL;
struct cnFill *gap=NULL;
struct chain *chain = NULL;
for(fill = list; fill != NULL; fill = fill->next)
{
if (positiveRangeIntersection(fill->tStart, fill->tStart+fill->tSize, start, end))
{
chain = chainFromId(fill->chainId);
if(checkChain(chain, start,end))
return chain;
for(gap = fill->children; gap != NULL; gap = gap->next)
{
chain = chainFromId(fill->chainId);
if(checkChain(chain, start,end))
return chain;
if(gap->children)
{
chain = lookForChain(gap->children, start, end);
if(checkChain(chain, start,end))
return chain;
}
}
}
}
return chain;
}
struct cnFill *getNetList(struct chain *chain, struct chainNet *netList)
/* return all cnFills from a netList that overlap a particular chain by at least 100 bases*/
{
struct cnFill *fillList = NULL, *fill, *fillClone;
struct chainNet *net;
for (net = netList; net != NULL; net = net->next)
{
for (fill = net->fillList; fill != NULL; fill = fill->next)
{
if (positiveRangeIntersection(fill->tStart,fill->tStart+fill->tSize, chain->tStart, chain->tEnd) > 10)
{
fillClone = cnFillNew();
fillClone->next = NULL ;
fillClone->children = fill->children ;
fillClone->tStart = fill->tStart ;
fillClone->tSize = fill->tSize ;
fillClone->qName = cloneString(fill->qName) ;
fillClone->type = cloneString(fill->type) ;
fillClone->qStrand = fill->qStrand ;
fillClone->qStart = fill->qStart ;
fillClone->qSize = fill->qSize ;
fillClone->chainId = fill->chainId ;
fillClone->ali = fill->ali = fill->ali ;
fillClone->tN = fill->tN = fill->tN ;
fillClone->qN = fill->qN = fill->qN ;
fillClone->tR = fill->tR ;
fillClone->qR = fill->qR ;
fillClone->tNewR = fill->tNewR ;
fillClone->qNewR = fill->qNewR ;
fillClone->tOldR = fill->tOldR ;
fillClone->qOldR = fill->qOldR ;
fillClone->tTrf = fill->tTrf ;
fillClone->qTrf = fill->qTrf ;
fillClone->qOver = fill->qOver ;
fillClone->qFar = fill->qFar ;
fillClone->qDup = fill->qDup ;
slAddHead(&fillList, fillClone);
printf(" %s %d-%d ", fill->type, fill->tStart, fill->tStart+fill->tSize);
}
//else printf(" %d-%d ", fill->tStart, fill->tStart+fill->tSize);
}
}
//slReverse(&fillList);
if (fillList == NULL)
printf("return null for %d %s %s:%d-%d\n",chain->id, chain->qName, chain->tName, chain->tStart, chain->tEnd);
return fillList;
}
void scoreChains(struct chain *chainList, struct bed *bedList, char *db, FILE *f, struct binKeeper *bk, struct hash *bkHash, struct genePred *gpList1, struct genePred *gpList2, struct genePred *gpList3, struct genePred *gpList4)
/* score all chains in a chain file based on ration of introns and gaps , mask repeats first */
/* reduce score for chains that are on the main syntenic net */
/* return best overlapping gene in the target sequence */
{
struct chain *chain, *retSubChain,*retChainToFree, *netChain;
struct cnFill *fillList, *fill;
struct intronChain *iChain=NULL, *iChainList = NULL;
struct genePred *gp = NULL, *vg;
struct bed *bed;
int overlap1, overlapVega;
boolean onDiagonal;
int overlapDiagonal;
for (chain = chainList; chain != NULL; chain = chain->next)
{
/* skip long chains - they are not pseudogenes */
if ((chain->tEnd-chain->tStart)> bigChain)
continue;
overlapDiagonal = 0; overlap1 = 0; overlapVega = 0;
for (bed = bedList ; bed != NULL ; bed = bed->next)
{
if( positiveRangeIntersection(bed->chromStart, bed->chromEnd, chain->tStart, chain->tEnd) > overlapDiagonal )
overlapDiagonal = positiveRangeIntersection(bed->chromStart, bed->chromEnd, chain->tStart, chain->tEnd);
}
/* look for gene on query side */
gp = getOverlappingGene(&gpList1, "refGene",chain->qName, qStart(chain), qEnd(chain) , &overlap1);
if (gp == NULL)
{
gp = getOverlappingGene(&gpList2, "softberryGene",chain->qName, qStart(chain), qEnd(chain) , &overlap1);
if (gp == NULL)
{
gp = getOverlappingGene(&gpList4, "all_mrna",chain->qName, qStart(chain), qEnd(chain) , &overlap1);
}
}
/* check for annontated pseudogene */
vg = getOverlappingGene(&gpList3, "vegaPseudoGene",chain->tName, chain->tStart, chain->tEnd, &overlapVega );
/* if we found something , score it and output record */
// if (gp != NULL || vg != NULL)
{
assert (chain->tEnd > 0);
retSubChain = NULL;
if (gp != NULL)
if (chain->qStrand =='+')
chainSubsetOnQ(chain, gp->txStart, gp->txEnd, &retSubChain, &retChainToFree);
else
chainSubsetOnQ(chain, chain->qSize - gp->txEnd, chain->qSize - gp->txStart, &retSubChain, &retChainToFree);
if (retSubChain != NULL)
{
iChain = scoreNextAlignedIntron(retSubChain,f, bk, bkHash);
iChain->overlap1 = overlap1;
iChain->overlap1 = overlapVega;
if (gp != NULL)
iChain->strand = gp->strand[0] == retSubChain->qStrand ? '+' : '-';
else
iChain->strand = '0';
slAddHead(&iChainList, iChain);
if (gp != NULL)
fprintf(f," %7.1f %7.1f %6d %2d %2d %5d %5d %4d %4d %s %d %d %c %s %d %d %c %d %d" ,
iChain->score1, iChain->score2, retSubChain->tEnd-retSubChain->tStart, iChain->tIntron, iChain->qIntron,
iChain->tGap, iChain->qGap, iChain->tReps, iChain->qReps, chain->tName,
retSubChain->tStart, retSubChain->tEnd, iChain->strand ,
retSubChain->qName, qStart(retSubChain), qEnd(retSubChain), retSubChain->qStrand,
chain->id, overlapDiagonal);
else
fprintf(f," %7.1f %7.1f %6d %2d %2d %5d %5d %4d %4d %s %d %d %c %s %d %d %c %d %d",
iChain->score1, iChain->score2, retSubChain->tEnd-retSubChain->tStart, iChain->tIntron, iChain->qIntron,
iChain->tGap, iChain->qGap, iChain->tReps, iChain->qReps, chain->tName,
retSubChain->tStart, retSubChain->tEnd, iChain->strand ,
retSubChain->qName, qStart(retSubChain), qEnd(retSubChain), retSubChain->qStrand,
chain->id, overlapDiagonal);
if (gp != NULL && genePredBases(gp) != 0)
fprintf(f," %s %s %d %d %d %2.2f %c",
gp->name, gp->chrom, gp->txStart, gp->txEnd, iChain->overlap1,
(float)iChain->overlap1/(float)genePredBases(gp), gp->strand[0]);
else
fprintf(f," NO gene 0 0 0 0.0 0");
if (vg != NULL && genePredBases(vg) != 0 && gp != NULL)
fprintf(f," %s %s %d %d %d %2.2f %c",
vg->name, vg->chrom, vg->txStart, vg->txEnd, iChain->overlap2,
(float)iChain->overlap2/(float)genePredBases(vg), vg->strand[0]);
fprintf(f,"\n");
chainFree(&retChainToFree);
}
else
{
fprintf(f," 0 0 0 0 0 0 0 0 0 %s %d %d + %s %d %d %c %d %d NO gene 0 0 0 0.0 0\n",
chain->tName, chain->tStart, chain->tEnd,
chain->qName, qStart(chain), qEnd(chain), chain->qStrand,
chain->id, overlapDiagonal);
}
}
}
}
int main(int argc, char *argv[])
/* Process command line. */
{
struct hash *chainHash;
FILE *outFile, *f;
struct binKeeper *tBk = NULL;
struct hash *bkHash = NULL;
struct hash *tSizeHash;
struct chain *chainList;
struct genePred *genePred1;
struct genePred *genePred2;
struct genePred *genePred3;
struct genePred *genePred4 = NULL, *gp;
struct lineFile *lf = NULL;
struct chainNet *netList = NULL, *net;
struct bed *bedList = NULL, *bed;
char *row[3];
int overlap;
optionHash(&argc, argv);
if (argc < 13)
usage();
outFile = mustOpen(argv[9], "w");
/*hSetDb(argv[1]);
hSetDb2(argv[2]);
db = hGetDb();
db2 = hGetDb2();*/
minGap = optionInt("minGap", minGap);
minOverlap = optionFloat("minOverlap", minOverlap);
tSizeHash = readSizes(argv[5]);
hSetDb(argv[1]);
printf("Loading Syntenic Bed\n");
lf = lineFileOpen(argv[4], TRUE);
while (lineFileRow(lf,row))
{
bed = bedLoad3(row);
slAddHead(&bedList, bed);
}
printf("Loading %s\n",argv[10]);
genePred1 = genePredLoadAll(argv[10]);
printf("Loading %s\n",argv[11]);
genePred2 = genePredLoadAll(argv[11]);
printf("Loading %s\n",argv[12]);
genePred3 = genePredLoadAll(argv[12]);
if (argc > 13)
{
printf("Loading %s\n",argv[13]);
genePred4 = genePredLoadAll(argv[13]);
}
printf("Loading Chains\n");
chainList = readSortChains(argv[3]);
printf("Loading Target Repeats\n");
tBk = readRepeats(argv[2], argv[7], tSizeHash);
printf("Loading Query Repeats\n");
bkHash = readRepeatsAll(argv[6], argv[8]);
printf("Scoring Chains\n");
if (genePred4 == NULL)
{
gp = getOverlappingGene(&genePred4, "all_mrna","chr1", 1, 100 , &overlap);
f = mustOpen("mrna.tab", "w");
for (gp = genePred4; gp != NULL; gp = gp->next)
{
genePredTabOut(gp, f);
}
fclose(f);
}
scoreChains(chainList, bedList, argv[1], outFile, tBk, bkHash, genePred1, genePred2, genePred3, genePred4);
return 0;
}