src/hg/instinct/bioInt2/bioGeneLevel.c 1.5
1.5 2009/04/27 06:15:48 jsanborn
updated lots of stuff, will break older implementation of database
Index: src/hg/instinct/bioInt2/bioGeneLevel.c
===================================================================
RCS file: /projects/compbio/cvsroot/kent/src/hg/instinct/bioInt2/bioGeneLevel.c,v
retrieving revision 1.4
retrieving revision 1.5
diff -b -B -U 1000000 -r1.4 -r1.5
--- src/hg/instinct/bioInt2/bioGeneLevel.c 4 Apr 2009 00:39:22 -0000 1.4
+++ src/hg/instinct/bioInt2/bioGeneLevel.c 27 Apr 2009 06:15:48 -0000 1.5
@@ -1,212 +1,175 @@
/* mapProbesToGenes - Will maps probes in BED format to overlapping gene(s). */
#include "common.h"
#include "linefile.h"
#include "hash.h"
#include "options.h"
#include "jksql.h"
#include "hPrint.h"
#include "hdb.h"
#include "dystring.h"
#include "bioIntDb.h"
#include "bioIntDriver.h"
#include "cprob.h"
#include "hgStatsLib.h"
#include "bioController.h"
/* Gene-level analysis functions */
-struct analysisResult *metaGene(struct biAnalysis *ba, void *data,
- char *sample, char *gene)
+struct analysisVals *metaGene(struct biAnalysis *ba, void *data,
+ int sample_id, int feature_id)
{
+if (!data)
+ return NULL;
+
+double total = 0.0, count = 0.0;
struct slPair *sp, *spList = data;
struct slDouble *sd, *sdList = NULL;
for (sp = spList; sp; sp = sp->next)
{
- struct slDouble *vals = sp->val;
- if (!vals)
+ struct analysisVals *av = sp->val;
+ if (!av)
continue;
- double val = slDoubleMedian(vals);
- sd = slDoubleNew(val);
+
+ total += av->val;
+ count += 1.0;
+ sd = slDoubleNew(av->conf);
slAddHead(&sdList, sd);
}
float chi2, metaP;
if (!fishersMetaSigned(sdList, &chi2, &metaP))
return NULL;
slFreeList(&sdList);
-struct analysisResult *ar;
-AllocVar(ar);
-ar->sample = cloneString(sample);
-ar->feature = cloneString(gene);
-ar->val = metaP;
-ar->conf = chi2;
-return ar;
+struct analysisVals *av;
+AllocVar(av);
+av->sample_id = sample_id;
+av->feature_id = feature_id;
+av->val = total / count;
+av->conf = metaP;
+return av;
}
/* Pipeline Stuff */
-void slPairDoubleFree(struct slPair **pEl)
-{
-struct slPair *el;
-if ((el = *pEl) == NULL) return;
-
-freeMem(el->name);
-struct slDouble *sdList = el->val;
-slFreeList(&sdList);
-freez(pEl);
-}
-
-void slPairDoubleFreeList(struct slPair **pList)
-{
-struct slPair *el, *next;
-
-for (el = *pList; el != NULL; el = next)
- {
- next = el->next;
- slPairDoubleFree(&el);
- }
-*pList = NULL;
-}
-
struct slName *getAvailableGenes(char *db, struct biResults *br)
{
-struct sqlConnection *biConn = hAllocConnProfile("localDb", db);
-char query[256];
-safef(query, sizeof(query),
- "select DISTINCT geneSymbol from kgXref join knownGene on kgXref.kgId = knownGene.name");
-
-struct slName *sl, *slList = sqlQuickList(biConn, query);
-
-struct slName *geneList = NULL;
-for (sl = slList; sl; sl = sl->next)
- {
- struct slName *probes = br->probesForGene(br, sl->name);
- int numProbes = slCount(probes);
- slNameFreeList(&probes);
- if (numProbes == 0)
- continue;
-
- slNameAddHead(&geneList, sl->name);
- }
-slReverse(&geneList);
-
-return geneList;
+return br->allFeatures(br);
}
struct slPair *geneLevelData(struct biResults *br, struct biData *bdList,
char *sample, char *gene)
{
struct biData *bd;
struct slPair *sp, *spList = NULL;
for (bd = bdList; bd; bd = bd->next)
{
- char *dataset = bd->name;
- struct slName *sl, *probes = br->probesForGeneInDataset(br, gene, dataset);
+ struct hashEl *el = hashLookup(bd->hash, gene);
+ if (!el)
+ continue;
AllocVar(sp);
sp->name = cloneString(bd->type);
sp->val = NULL;
- struct slDouble *sd, *sdList = NULL;
- for (sl = probes; sl; sl = sl->next)
- {
- struct hashEl *el = hashLookup(bd->hash, sl->name);
- if (!el)
- continue;
- struct slDouble *sdVal = el->val;
- sd = slDoubleNew(sdVal->val);
- slAddHead(&sdList, sd);
- }
- sp->val = sdList;
- slAddHead(&spList, sp);
- slNameFreeList(&probes);
+ struct analysisVals *newAv, *av = el->val;
+ newAv = cloneAnalysisVals(av);
+ sp->val = newAv;
+ slAddHead(&spList, sp);
}
return spList;
}
void geneLevelAnalysis(struct sqlConnection *biConn,
struct biAnalysis *ba, struct biResults *br,
struct slName *genes)
{
if (!ba->analyze)
return;
fprintf(stdout, "starting gene level analysis\n");
+struct hash *sampleHash = createIdHash(biConn, SA_TABLE, "name");
+struct hash *featureHash = createIdHash(biConn, AF_TABLE, "feature_name");
+
struct slName *gene, *sample, *samples = br->allSamplesInCommon(br);
-struct analysisResult *ar, *arList = NULL;
+int count = 0;
+int numSamples = slCount(samples);
+
+struct analysisVals *av, *avList = NULL;
for (sample = samples; sample; sample = sample->next)
{
+ int sample_id = hashIntValDefault(sampleHash, sample->name, -1);
+ if (sample_id == -1)
+ continue;
+
struct biData *bdList = br->dataForSample(br, sample->name);
for (gene = genes; gene; gene = gene->next)
{
+ int feature_id = hashIntValDefault(featureHash, gene->name, -1);
+ if (feature_id == -1)
+ continue;
+
struct slPair *spList = geneLevelData(br, bdList, sample->name, gene->name);
if (!spList)
continue;
- ar = ba->analyze(ba, spList, sample->name, gene->name);
- slPairDoubleFreeList(&spList);
- if (!ar)
+ av = ba->analyze(ba, spList, sample_id, feature_id); //e->name, gene->name);
+ slPairAnalysisValsFreeList(&spList);
+ if (!av)
continue;
- slAddHead(&arList, ar);
+ slAddHead(&avList, av);
}
- slReverse(&arList);
-
- storeAnalysisResultsInDb(biConn, ba, arList);
- analysisResultFreeList(&arList);
- arList = NULL;
-
- fprintf(stdout, ".");
+ count++;
+ fprintf(stdout, "%d of %d samples\n", count, numSamples);
fflush(stdout);
biDataFree(&bdList);
}
-fprintf(stdout, "\n");
+uglyTime(NULL);
+storeAnalysisValsInDb(biConn, ba->tableName, avList);
+uglyTime("stored ARs");
+
+hashFree(&sampleHash);
+hashFree(&featureHash);
}
struct biResults *retrieveRawData(char *db, struct slName *datasets, boolean toLogP)
{
struct slName *dataset;
struct biQuery *bq, *bqList = NULL;
for (dataset = datasets; dataset; dataset = dataset->next)
{
bq = biQueryNew(db, dataset->name); // Add Dataset
- bq->getAllProbes = TRUE; // Set flag to retrieve all probes
+ bq->getAllFeatures = TRUE; // Set flag to retrieve all data
biQueryAppend(&bqList, bq); // Append query to query list
}
/* Get results from all queries in list */
struct biResults *br = biQueryResults(bqList);
-/* Convert to log p */
-if (toLogP)
- br->toLogP(br);
-
return br;
}
void geneLevelPipeline(struct biAnalysis *ba)
{
/* Get raw gene/sample data for all overlapping samples in dataset list */
struct biResults *br = retrieveRawData(ba->db, ba->inputTables, TRUE);
if (!br)
return;
struct slName *genes = getAvailableGenes(ba->db, br);
struct sqlConnection *biConn = hAllocConnProfile("localDb", ba->db);
geneLevelAnalysis(biConn, ba, br, genes);
-//storeAnalysisResultsInDb(biConn, ba, arList);
hFreeConn(&biConn);
-//analysisResultFreeList(&arList);
biResultsFree(&br);
slNameFreeList(&genes);
}