Commit 9881decf by Eric Coissac

--no commit message

parent a127f573
......@@ -44,6 +44,8 @@ buffer_t *newBuffer(size_t maxsize,size_t readSize)
rep->index1 = MALLOC(sizeof(uint32_t) * (1 << 16));
rep->index2 = MALLOC(sizeof(uint32_t) * (1 << 16));
rep->complement = NULL;
return rep;
#undef MASK
......@@ -58,6 +60,9 @@ void freeBuffer(buffer_t *buffer)
FREE(buffer->arena);
FREE(buffer->index1);
FREE(buffer->index2);
if (buffer->complement)
FREE(buffer->complement);
FREE(buffer);
}
};
......
......@@ -18,7 +18,7 @@ buffer_t *buildIndex(const char *forwardFileName, const char* reverseFileName,
char* indexFileName;
uint32_t indexcount;
// uint32_t indexcount;
uint32_t i;
buffer_t *reads;
......@@ -93,28 +93,9 @@ buffer_t *buildIndex(const char *forwardFileName, const char* reverseFileName,
fclose(index);
free(indexFileName);
indexcount = reads->readSize - minword;
reads->minword =minword;
reads->suffixidx=indexcount;
// for (i=1; i <= indexcount; i++)
// {
// sortSuffix(reads,i);
// swapOrder(reads);
//
// asprintf(&indexFileName,"%s.o%02d",indexname,i);
// if (indexFileName == NULL)
// FATALERROR("Cannot allocate memory for file name %s",indexname);
//
// index = fopen(indexFileName,"w");
// if (index == NULL)
// FATALERROR("Cannot open file %s",indexFileName);
//
// writeOrderData(reads,index);
//
// fclose(index);
//
// }
// indexcount = reads->readSize - minword;
// reads->minword =minword;
// reads->suffixidx=indexcount;
asprintf(&indexFileName,"%s.ogx",indexname);
if (indexFileName == NULL)
......
......@@ -133,6 +133,8 @@ char *getRead(buffer_t *buffer, int32_t recordid, uint32_t begin, int32_t length
{
recordid=-recordid;
recordid--;
recordid=buffer->complement[recordid];
return decodeComp(buffer, recordid, begin, length, dest);
}
......
......@@ -7,6 +7,18 @@
#include "orgasm.h"
void indexReverseComplement(buffer_t *reads)
{
int32_t i;
if (reads->complement)
FREE(reads->complement);
reads->complement=(int32_t*) MALLOC(sizeof(int32_t) * reads->readCount);
for (i=0;i < reads->readCount; i++)
reads->complement[reads->order1[i]]=i;
}
buffer_t *loadIndexedReads(const char *indexname)
{
buffer_t tmp;
......@@ -15,6 +27,10 @@ buffer_t *loadIndexedReads(const char *indexname)
FILE *index;
size_t maxsize;
//
// Read the general data about the sequence index
//
asprintf(&indexFileName,"%s.ogx",indexname);
if (indexFileName == NULL)
FATALERROR("Cannot allocate memory for file name %s",indexname);
......@@ -35,6 +51,10 @@ buffer_t *loadIndexedReads(const char *indexname)
reads->readCount = tmp.readCount;
//
// Read the forward sequence data
//
asprintf(&indexFileName,"%s.ofx",indexname);
if (indexFileName == NULL)
FATALERROR("Cannot allocate memory for file name %s",indexname);
......@@ -48,6 +68,10 @@ buffer_t *loadIndexedReads(const char *indexname)
fclose(index);
free(indexFileName);
//
// Read the pair-end data --> into order2
//
asprintf(&indexFileName,"%s.opx",indexname);
if (indexFileName == NULL)
FATALERROR("Cannot allocate memory for file name %s",indexname);
......@@ -61,6 +85,10 @@ buffer_t *loadIndexedReads(const char *indexname)
fclose(index);
free(indexFileName);
//
// Read index for the reverse complemented sequences --> into order1
//
asprintf(&indexFileName,"%s.orx",indexname);
if (indexFileName == NULL)
FATALERROR("Cannot allocate memory for file name %s",indexname);
......@@ -76,8 +104,10 @@ buffer_t *loadIndexedReads(const char *indexname)
fclose(index);
free(indexFileName);
reads->minword = tmp.minword;
reads->suffixidx= tmp.suffixidx;
reads->complement=NULL; // Provisoire...
fprintf(stderr,"\nIndexing reverse complement sequences ...\n\n");
indexReverseComplement(reads);
fprintf(stderr,"\nFast indexing forward reads...\n\n");
indexForward(reads);
......
......@@ -453,13 +453,15 @@ int32_t lookForReads(buffer_t *buffer, pnuc key, size_t length, int32_t* endofli
{
int32_t pos = lookForForward(buffer, key, length, endoflist);
if (*endoflist==0)
pos++;
else
{
pos = lookForReverse(buffer, key, length, endoflist);
if (*endoflist==0)
pos = -pos - 1;
pos = - buffer->order1[pos] - 1;
}
return pos;
......@@ -467,9 +469,12 @@ int32_t lookForReads(buffer_t *buffer, pnuc key, size_t length, int32_t* endofli
int32_t nextRead(buffer_t *buffer, int32_t current, size_t length, int32_t* endoflist)
{
int32_t next;
int32_t next=0;
pnuc from;
if (current == 0)
*endoflist=1;
if (current > 0)
{
current--;
......@@ -484,23 +489,20 @@ int32_t nextRead(buffer_t *buffer, int32_t current, size_t length, int32_t* endo
next = lookForReverse(buffer, from, length, endoflist);
if (*endoflist==0)
next = -next -1;
next = - buffer->order1[next] -1;
}
}
if (current < 0)
{
current = -current -1;
current = buffer->complement[-current -1];
next = nextReverse(buffer, current, length, endoflist);
if (*endoflist==0)
next = -next -1;
next = - buffer->order1[next] -1;
}
if (current == 0)
*endoflist=0;
return next;
}
......
......@@ -70,6 +70,7 @@ typedef struct {
// ex: gives the num of the pair in records pair(n) => records[order2[n]]
uint32_t* index1;
uint32_t* index2;
uint32_t* complement; // convert index to reverse complement sequence
uint32_t minword; // minimal size of word to consider
uint32_t suffixidx; // count of suffix indices stored on disk
} buffer_t;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment