extractRef.cpp

   1 #include<cstdio>
   2 #include<cstring>
   3 #include<cctype>
   4 #include<cstdlib>
   5 #include<fstream>
   6 #include<sstream>
   7 #include<map>
   8 #include<vector>
   9 #include<algorithm>
  10
  11 #include "utils.h"
  12 #include "GTFItem.h"
  13 #include "Transcript.h"
  14 #include "Transcripts.h"
  15
  16 using namespace std;
  17
  18 struct ChrInfo {
  19         string name;
  20         size_t len;
  21
  22         ChrInfo(const string& name, size_t len) {
  23                 this->name = name;
  24                 this->len = len;
  25         }
  26
  27         bool operator< (const ChrInfo& o) const {
  28                 return name < o.name;
  29         }
  30 };
  31
  32 int M;
  33
  34 vector<GTFItem> items;
  35 vector<string> seqs;
  36 vector<int> starts; // used to generate .grp
  37 map<string, vector<int> > sn2tr; // map from seqname to transcripts
  38 map<string, vector<int> >::iterator iter;
  39 vector<ChrInfo> chrvec;
  40
  41 Transcripts transcripts;
  42
  43 char groupF[STRLEN], tiF[STRLEN], refFastaF[STRLEN];
  44 char chromListF[STRLEN];
  45
  46 bool hasMappingFile;
  47 char mappingFile[STRLEN];
  48
  49 map<string, string> mi_table; // mapping info table
  50 map<string, string>::iterator mi_iter; //mapping info table's iterator
  51
  52 void loadMappingInfo(char* mappingF) {
  53   ifstream fin(mappingF);
  54   string line, key, value;
  55
  56   if (!fin.is_open()) {
  57           fprintf(stderr, "Cannot open %s! It may not exist.\n", mappingF);
  58           exit(-1);
  59   }
  60
  61   mi_table.clear();
  62   while (getline(fin, line)) {
  63     line = cleanStr(line);
  64     if (line[0] == '#') continue;
  65     istringstream strin(line);
  66     strin>>value>>key;
  67     mi_table[key] = value;
  68   }
  69
  70   fin.close();
  71 }
  72
  73 bool buildTranscript(int sp, int ep) {
  74         int cur_s, cur_e; // current_start, current_end
  75         vector<Interval> vec;
  76
  77         string transcript_id = items[sp].getTranscriptID();
  78         string gene_id = items[sp].getGeneID();
  79         char strand = items[sp].getStrand();
  80         string seqname = items[sp].getSeqName();
  81         string left = items[sp].getLeft();
  82
  83         vec.clear();
  84         cur_s = cur_e = -1;
  85         for (int i = sp; i <= ep; i++) {
  86                 int start = items[i].getStart();
  87                 int end = items[i].getEnd();
  88
  89                 assert(strand == items[i].getStrand());
  90                 assert(seqname == items[i].getSeqName());
  91
  92                 if (cur_e + 1 < start) {
  93                         if (cur_s > 0) vec.push_back(Interval(cur_s, cur_e));
  94                         cur_s = start;
  95                 }
  96                 cur_e = (cur_e < end ? end : cur_e);
  97         }
  98         if (cur_s > 0) vec.push_back(Interval(cur_s, cur_e));
  99
 100         transcripts.add(Transcript(transcript_id, gene_id, seqname, strand, vec, left));
 101
 102         return true;
 103 }
 104
 105 void parse_gtf_file(char* gtfF) {
 106         ifstream fin(gtfF);
 107         string line, curgid, tid, gid; //  curgid: current gene id;
 108         GTFItem item;
 109
 110         if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", gtfF); exit(-1); }
 111
 112         int cnt = 0;
 113
 114         items.clear();
 115         while (getline(fin, line)) {
 116                 if (line[0] == '#') continue; // if this line is comment, jump it
 117                 item.parse(line);
 118                 string feature = item.getFeature();
 119                 if (feature == "exon") {
 120                         if (item.getStart() > item.getEnd()) {
 121                                 fprintf(stderr, "Warning: exon's start position is larger than its end position! This exon is discarded.\n");
 122                                 fprintf(stderr, "\t%s\n\n", line.c_str());
 123                         }
 124                         else if (item.getStart() < 1) {
 125                                 fprintf(stderr, "Warning: exon's start position is less than 1! This exon is discarded.\n");
 126                                 fprintf(stderr, "\t%s\n\n", line.c_str());
 127                         }
 128                         else {
 129                                 if (hasMappingFile) {
 130                                         tid = item.getTranscriptID();
 131                                     mi_iter = mi_table.find(tid);
 132                                     if (mi_iter == mi_table.end()) {
 133                                         fprintf(stderr, "Mapping Info is not correct, cannot find %s's gene_id!\n", tid.c_str());
 134                                         exit(-1);
 135                                     }
 136                                     //assert(iter != table.end());
 137                                     gid = mi_iter->second;
 138                                     item.setGeneID(gid);
 139                                 }
 140                                 items.push_back(item);
 141                         }
 142                 }
 143
 144                 ++cnt;
 145                 if (verbose && cnt % 200000 == 0) { printf("Parsed %d lines\n", cnt); }
 146         }
 147         fin.close();
 148
 149         sort(items.begin(), items.end());
 150
 151         starts.clear();
 152         sn2tr.clear();
 153         curgid = "";
 154
 155         int sp = 0, ep; // start pointer, end pointer
 156         int nItems = items.size();
 157
 158         while (sp < nItems) {
 159                 tid = items[sp].getTranscriptID();
 160                 gid = items[sp].getGeneID();
 161
 162                 ep = sp + 1;
 163                 while (ep < nItems && items[ep].getTranscriptID() == tid) ep++;
 164                 ep--;
 165
 166                 buildTranscript(sp, ep);
 167
 168                 int sid = transcripts.getM();
 169                 const Transcript& transcript = transcripts.getTranscriptAt(sid);
 170
 171                 if (curgid != gid) {
 172                         starts.push_back(sid);
 173                         curgid = gid;
 174                 }
 175                 iter = sn2tr.find(transcript.getSeqName());
 176                 if (iter == sn2tr.end()) {
 177                         vector<int> vec(1, sid);
 178                         sn2tr[transcript.getSeqName()] = vec;
 179                 }
 180                 else {
 181                         iter->second.push_back(sid);
 182                 }
 183
 184                 sp = ep + 1;
 185         }
 186
 187         M = transcripts.getM();
 188         starts.push_back(M + 1);
 189         items.clear();
 190
 191         if (M < 1) {
 192                 fprintf(stderr, "Number of transcripts in the reference is less than 1!\n");
 193                 exit(-1);
 194         }
 195
 196         if (verbose) { printf("Parsing gtf File is done!\n"); }
 197 }
 198
 199 char check(char c) {
 200         if (!isalpha(c)) { fprintf(stderr, "Sequence contains unknown letter '%c'!\n", c); exit(-1); }
 201         //assert(isalpha(c));
 202         if (isupper(c) && c != 'A' && c != 'C' && c != 'G' && c != 'T') c = 'N';
 203         if (islower(c) && c != 'a' && c != 'c' && c != 'g' && c != 't') c = 'n';
 204         return c;
 205 }
 206
 207 void writeResults(char* refName) {
 208         int s;
 209         ofstream fout;
 210
 211         sprintf(groupF, "%s.grp", refName);
 212         sprintf(tiF, "%s.ti", refName);
 213         sprintf(refFastaF, "%s.transcripts.fa", refName);
 214         sprintf(chromListF, "%s.chrlist", refName);
 215
 216
 217         fout.open(groupF);
 218         s = starts.size();
 219         for (int i = 0; i < s; i++) fout<<starts[i]<<endl;
 220         fout.close();
 221         if (verbose) { printf("Group File is generated!\n"); }
 222
 223         transcripts.writeTo(tiF);
 224         if (verbose) { printf("Transcript Information File is generated!\n"); }
 225
 226         fout.open(chromListF);
 227         s = chrvec.size();
 228         for (int i = 0; i < s; i++) {
 229                 fout<<chrvec[i].name<<'\t'<<chrvec[i].len<<endl;
 230         }
 231         fout.close();
 232         if (verbose) { printf("Chromosome List File is generated!\n"); }
 233
 234         fout.open(refFastaF);
 235         for (int i = 1; i <= M; i++) {
 236                 fout<<">"<<transcripts.getTranscriptAt(i).getTranscriptID()<<endl;
 237                 fout<<seqs[i]<<endl;
 238         }
 239         fout.close();
 240         if (verbose) { printf("Extracted Sequences File is generated!\n"); }
 241 }
 242
 243 int main(int argc, char* argv[]) {
 244         if (argc < 6 || (hasMappingFile = atoi(argv[4])) && argc < 7) {
 245                 printf("Usage: rsem-extract-reference-transcripts refName quiet gtfF hasMappingFile [mappingFile] chromosome_file_1 [chromosome_file_2 ...]\n");
 246                 exit(-1);
 247         }
 248
 249         verbose = !atoi(argv[2]);
 250         if (hasMappingFile) {
 251                 loadMappingInfo(argv[5]);
 252         }
 253         parse_gtf_file(argv[3]);
 254
 255         ifstream fin;
 256         string line, gseq, seqname;
 257         void* pt;
 258
 259         chrvec.clear();
 260
 261         seqs.clear();
 262         seqs.resize(M + 1, "");
 263         int start = hasMappingFile ? 6 : 5;
 264         for (int i = start; i < argc; i++) {
 265                 fin.open(argv[i]);
 266                 if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", argv[i]); exit(-1); }
 267                 pt = getline(fin, line);
 268                 while (pt != 0 && line[0] == '>') {
 269                         istringstream strin(line.substr(1));
 270                         strin>>seqname;
 271
 272                         gseq = "";
 273                     while((pt = getline(fin, line)) && line[0] != '>') {
 274                         gseq += line;
 275                     }
 276
 277                     size_t len = gseq.length();
 278                     assert(len > 0);
 279                     for (size_t j = 0; j < len; j++) gseq[j] = check(gseq[j]);
 280
 281                     iter = sn2tr.find(seqname);
 282                     if (iter == sn2tr.end()) continue;
 283
 284                     chrvec.push_back(ChrInfo(seqname, len));
 285
 286                     vector<int>& vec = iter->second;
 287                     int s = vec.size();
 288                     for (int j = 0; j < s; j++) {
 289                         assert(vec[j] > 0 && vec[j] <= M);
 290                         transcripts.getTranscriptAt(vec[j]).extractSeq(gseq, seqs[vec[j]]);
 291                     }
 292                 }
 293                 fin.close();
 294
 295                 if (verbose) { printf("%s is processed!\n", argv[i]); }
 296         }
 297
 298         for (int i = 1; i <= M; i++) {
 299                 if (seqs[i] == "") {
 300                         fprintf(stderr, "%s's sequence is empty! You must provide all chromosome files of transcripts which are presented in the .gtf file!\n", transcripts.getTranscriptAt(i).getTranscriptID().c_str());
 301                         exit(-1);
 302                 }
 303         }
 304
 305         sort(chrvec.begin(), chrvec.end());
 306
 307         if (verbose) { printf("Extracting sequences is done!\n"); }
 308
 309         writeResults(argv[1]);
 310
 311         return 0;
 312 }