| 1 |
/* |
| 2 |
* GFaIdx.cpp |
| 3 |
* |
| 4 |
* Created on: Aug 25, 2010 |
| 5 |
* Author: gpertea |
| 6 |
*/ |
| 7 |
|
| 8 |
#include "GFastaIndex.h" |
| 9 |
#define ERR_FAIDXLINE "Error parsing fasta index line: \n%s\n" |
| 10 |
#define ERR_FALINELEN "Error: sequence lines in a FASTA record must have the same length!\n" |
| 11 |
void GFastaIndex::addRecord(const char* seqname, uint seqlen, off_t foffs, int llen, int llen_full) { |
| 12 |
GFastaRec* farec=records.Find(seqname); |
| 13 |
if (farec!=NULL) { |
| 14 |
GMessage("Warning: duplicate sequence ID (%s) added to the fasta index! Only last entry data will be kept.\n"); |
| 15 |
farec->seqlen=seqlen; |
| 16 |
farec->fpos=foffs; |
| 17 |
farec->line_len=llen; |
| 18 |
farec->line_blen=llen_full; |
| 19 |
} |
| 20 |
else { |
| 21 |
farec=new GFastaRec(seqlen,foffs,llen,llen_full); |
| 22 |
records.Add(seqname,farec); |
| 23 |
farec->seqname=records.getLastKey(); |
| 24 |
} |
| 25 |
} |
| 26 |
|
| 27 |
int GFastaIndex::loadIndex(const char* finame) { //load record info from existing fasta index |
| 28 |
if (finame==NULL) finame=fai_name; |
| 29 |
if (finame!=fai_name) { |
| 30 |
fai_name=Gstrdup(finame); |
| 31 |
} |
| 32 |
if (fai_name==NULL) GError("Error: GFastaIndex::loadIndex() called with no file name!\n"); |
| 33 |
records.Clear(); |
| 34 |
haveFai=false; |
| 35 |
FILE* fi=fopen(fai_name,"rb"); |
| 36 |
if (fi==NULL) { |
| 37 |
GMessage("Warning: cannot open fasta index file: %s!\n",fai_name); |
| 38 |
return 0; |
| 39 |
} |
| 40 |
GLineReader fl(fi); |
| 41 |
char* s=NULL; |
| 42 |
while ((s=fl.nextLine())!=NULL) { |
| 43 |
if (*s=='#') continue; |
| 44 |
char* p=strchrs(s,"\t "); |
| 45 |
if (p==NULL) GError(ERR_FAIDXLINE,s); |
| 46 |
*p=0; //s now holds the genomic sequence name |
| 47 |
p++; |
| 48 |
uint len=0; |
| 49 |
int line_len=0, line_blen=0; |
| 50 |
#ifdef __WIN32__ |
| 51 |
long offset=-1; |
| 52 |
sscanf(p, "%d%ld%d%d", &len, &offset, &line_len, &line_blen); |
| 53 |
#else |
| 54 |
long long offset=-1; |
| 55 |
sscanf(p, "%d%lld%d%d", &len, &offset, &line_len, &line_blen); |
| 56 |
#endif |
| 57 |
if (len==0 || line_len==0 || line_blen==0 || line_blen<line_len) |
| 58 |
GError(ERR_FAIDXLINE,p); |
| 59 |
addRecord(s,len,offset,line_len, line_blen); |
| 60 |
} |
| 61 |
fclose(fi); |
| 62 |
haveFai=(records.Count()>0); |
| 63 |
return records.Count(); |
| 64 |
} |
| 65 |
|
| 66 |
int GFastaIndex::buildIndex() { |
| 67 |
//this parses the whole fasta file, so it could be slow |
| 68 |
if (fa_name==NULL) |
| 69 |
GError("Error: GFastaIndex::buildIndex() called with no fasta file!\n"); |
| 70 |
FILE* fa=fopen(fa_name,"rb"); |
| 71 |
if (fa==NULL) { |
| 72 |
GMessage("Warning: cannot open fasta index file: %s!\n",fa_name); |
| 73 |
return 0; |
| 74 |
} |
| 75 |
records.Clear(); |
| 76 |
GLineReader fl(fa); |
| 77 |
char* s=NULL; |
| 78 |
uint seqlen=0; |
| 79 |
int line_len=0,line_blen=0; |
| 80 |
bool newSeq=false; //set to true after defline |
| 81 |
off_t newSeqOffset=0; |
| 82 |
int prevOffset=0; |
| 83 |
char* seqname=NULL; |
| 84 |
int last_len=0; |
| 85 |
bool mustbeLastLine=false; //true if the line length decreases |
| 86 |
while ((s=fl.nextLine())!=NULL) { |
| 87 |
if (s[0]=='>') { |
| 88 |
if (seqname!=NULL) { |
| 89 |
if (seqlen==0) |
| 90 |
GError("Warning: empty FASTA record skipped (%s)!\n",seqname); |
| 91 |
else { //seqlen!=0 |
| 92 |
addRecord(seqname, seqlen,newSeqOffset, line_len, line_blen); |
| 93 |
} |
| 94 |
} |
| 95 |
char *p=s; |
| 96 |
while (*p > 32) p++; |
| 97 |
*p=0; |
| 98 |
GFREE(seqname); |
| 99 |
seqname=Gstrdup(&s[1]); |
| 100 |
newSeq=true; |
| 101 |
newSeqOffset=fl.getfpos(); |
| 102 |
last_len=0; |
| 103 |
line_len=0; |
| 104 |
line_blen=0; |
| 105 |
seqlen=0; |
| 106 |
mustbeLastLine=false; |
| 107 |
} //defline parsing |
| 108 |
else { //sequence line |
| 109 |
int llen=fl.length(); |
| 110 |
int lblen=fl.getFpos()-prevOffset; |
| 111 |
if (newSeq) { //first sequence line after defline |
| 112 |
line_len=llen; |
| 113 |
line_blen=lblen; |
| 114 |
} |
| 115 |
else {//next seq lines after first |
| 116 |
if (mustbeLastLine || llen>last_len) |
| 117 |
GError(ERR_FALINELEN); |
| 118 |
if (llen<last_len) mustbeLastLine=true; |
| 119 |
} |
| 120 |
seqlen+=llen; |
| 121 |
last_len=llen; |
| 122 |
newSeq=false; |
| 123 |
} //sequence line |
| 124 |
prevOffset=fl.getfpos(); |
| 125 |
}//for each line of the fasta file |
| 126 |
if (seqlen>0) |
| 127 |
addRecord(seqname, seqlen, newSeqOffset, line_len, line_blen); |
| 128 |
GFREE(seqname); |
| 129 |
fclose(fa); |
| 130 |
return records.Count(); |
| 131 |
} |
| 132 |
|
| 133 |
|
| 134 |
int GFastaIndex::storeIndex(const char* finame) { //write the hash to a file |
| 135 |
if (records.Count()==0) |
| 136 |
GError("Error at GFastaIndex:storeIndex(): no records found!\n"); |
| 137 |
FILE* fai=fopen(finame, "w"); |
| 138 |
if (fai==NULL) GError("Error creating fasta index file: %s\n",finame); |
| 139 |
int rcount=storeIndex(fai); |
| 140 |
GFREE(fai_name); |
| 141 |
fai_name=Gstrdup(finame); |
| 142 |
return rcount; |
| 143 |
} |
| 144 |
|
| 145 |
int GFastaIndex::storeIndex(FILE* fai) { |
| 146 |
int rcount=0; |
| 147 |
GList<GFastaRec> reclist(true,false,true); //sorted, don't free members, unique |
| 148 |
records.startIterate(); |
| 149 |
GFastaRec* rec=NULL; |
| 150 |
while ((rec=records.NextData())!=NULL) { |
| 151 |
reclist.Add(rec); |
| 152 |
} |
| 153 |
//reclist has records sorted by file offset |
| 154 |
for (int i=0;i<reclist.Count();i++) { |
| 155 |
#ifdef __WIN32__ |
| 156 |
int written=fprintf(fai, "%s\t%d\t%ld\t%d\t%d\n", |
| 157 |
reclist[i]->seqname,reclist[i]->seqlen,(long)reclist[i]->fpos, |
| 158 |
reclist[i]->line_len, reclist[i]->line_blen); |
| 159 |
#else |
| 160 |
int written=fprintf(fai, "%s\t%d\t%lld\t%d\t%d\n", |
| 161 |
reclist[i]->seqname, reclist[i]->seqlen, (long long)(reclist[i]->fpos), |
| 162 |
reclist[i]->line_len, reclist[i]->line_blen); |
| 163 |
#endif |
| 164 |
if (written>0) rcount++; |
| 165 |
else break; //couldn't write anymore |
| 166 |
} |
| 167 |
fclose(fai); |
| 168 |
haveFai=(rcount>0); |
| 169 |
return rcount; |
| 170 |
} |
