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                 item.parse(line);
 117                 string feature = item.getFeature();
 118                 if (feature == "exon") {
 119                         if (item.getStart() > item.getEnd()) {
 120                                 fprintf(stderr, "Warning: exon's start position is larger than its end position! This exon is discarded.\n");
 121                                 fprintf(stderr, "\t%s\n\n", line.c_str());
 122                         }
 123                         else if (item.getStart() < 1) {
 124                                 fprintf(stderr, "Warning: exon's start position is less than 1! This exon is discarded.\n");
 125                                 fprintf(stderr, "\t%s\n\n", line.c_str());
 126                         }
 127                         else {
 128                                 if (hasMappingFile) {
 129                                         tid = item.getTranscriptID();
 130                                     mi_iter = mi_table.find(tid);
 131                                     if (mi_iter == mi_table.end()) {
 132                                         fprintf(stderr, "Mapping Info is not correct, cannot find %s's gene_id!\n", tid.c_str());
 133                                         exit(-1);
 134                                     }
 135                                     //assert(iter != table.end());
 136                                     gid = mi_iter->second;
 137                                     item.setGeneID(gid);
 138                                 }
 139                                 items.push_back(item);
 140                         }
 141                 }
 142
 143                 ++cnt;
 144                 if (verbose && cnt % 200000 == 0) { printf("Parsed %d lines\n", cnt); }
 145         }
 146         fin.close();
 147
 148         sort(items.begin(), items.end());
 149
 150         starts.clear();
 151         sn2tr.clear();
 152         curgid = "";
 153
 154         int sp = 0, ep; // start pointer, end pointer
 155         int nItems = items.size();
 156
 157         while (sp < nItems) {
 158                 tid = items[sp].getTranscriptID();
 159                 gid = items[sp].getGeneID();
 160
 161                 ep = sp + 1;
 162                 while (ep < nItems && items[ep].getTranscriptID() == tid) ep++;
 163                 ep--;
 164
 165                 buildTranscript(sp, ep);
 166
 167                 int sid = transcripts.getM();
 168                 const Transcript& transcript = transcripts.getTranscriptAt(sid);
 169
 170                 if (curgid != gid) {
 171                         starts.push_back(sid);
 172                         curgid = gid;
 173                 }
 174                 iter = sn2tr.find(transcript.getSeqName());
 175                 if (iter == sn2tr.end()) {
 176                         vector<int> vec(1, sid);
 177                         sn2tr[transcript.getSeqName()] = vec;
 178                 }
 179                 else {
 180                         iter->second.push_back(sid);
 181                 }
 182
 183                 sp = ep + 1;
 184         }
 185
 186         M = transcripts.getM();
 187         starts.push_back(M + 1);
 188         items.clear();
 189
 190         if (M < 1) {
 191                 fprintf(stderr, "Number of transcripts in the reference is less than 1!\n");
 192                 exit(-1);
 193         }
 194
 195         if (verbose) { printf("Parsing gtf File is done!\n"); }
 196 }
 197
 198 char check(char c) {
 199         if (!isalpha(c)) { fprintf(stderr, "Sequence contains unknown letter '%c'!\n", c); exit(-1); }
 200         //assert(isalpha(c));
 201         if (isupper(c) && c != 'A' && c != 'C' && c != 'G' && c != 'T') c = 'N';
 202         if (islower(c) && c != 'a' && c != 'c' && c != 'g' && c != 't') c = 'n';
 203         return c;
 204 }
 205
 206 void writeResults(char* refName) {
 207         int s;
 208         ofstream fout;
 209
 210         sprintf(groupF, "%s.grp", refName);
 211         sprintf(tiF, "%s.ti", refName);
 212         sprintf(refFastaF, "%s.transcripts.fa", refName);
 213         sprintf(chromListF, "%s.chrlist", refName);
 214
 215
 216         fout.open(groupF);
 217         s = starts.size();
 218         for (int i = 0; i < s; i++) fout<<starts[i]<<endl;
 219         fout.close();
 220         if (verbose) { printf("Group File is generated!\n"); }
 221
 222         transcripts.writeTo(tiF);
 223         if (verbose) { printf("Transcript Information File is generated!\n"); }
 224
 225         fout.open(chromListF);
 226         s = chrvec.size();
 227         for (int i = 0; i < s; i++) {
 228                 fout<<chrvec[i].name<<'\t'<<chrvec[i].len<<endl;
 229         }
 230         fout.close();
 231         if (verbose) { printf("Chromosome List File is generated!\n"); }
 232
 233         fout.open(refFastaF);
 234         for (int i = 1; i <= M; i++) {
 235                 fout<<">"<<transcripts.getTranscriptAt(i).getTranscriptID()<<endl;
 236                 fout<<seqs[i]<<endl;
 237         }
 238         fout.close();
 239         if (verbose) { printf("Extracted Sequences File is generated!\n"); }
 240 }
 241
 242 int main(int argc, char* argv[]) {
 243         if (argc < 6 || (hasMappingFile = atoi(argv[4])) && argc < 7) {
 244                 printf("Usage: rsem-extract-reference-transcripts refName quiet gtfF hasMappingFile [mappingFile] chromosome_file_1 [chromosome_file_2 ...]\n");
 245                 exit(-1);
 246         }
 247
 248         verbose = !atoi(argv[2]);
 249         if (hasMappingFile) {
 250                 loadMappingInfo(argv[5]);
 251         }
 252         parse_gtf_file(argv[3]);
 253
 254         ifstream fin;
 255         string line, gseq, seqname;
 256         void* pt;
 257
 258         chrvec.clear();
 259
 260         seqs.clear();
 261         seqs.resize(M + 1, "");
 262         int start = hasMappingFile ? 6 : 5;
 263         for (int i = start; i < argc; i++) {
 264                 fin.open(argv[i]);
 265                 if (!fin.is_open()) { fprintf(stderr, "Cannot open %s! It may not exist.\n", argv[i]); exit(-1); }
 266                 pt = getline(fin, line);
 267                 while (pt != 0 && line[0] == '>') {
 268                         istringstream strin(line.substr(1));
 269                         strin>>seqname;
 270
 271                         gseq = "";
 272                     while((pt = getline(fin, line)) && line[0] != '>') {
 273                         gseq += line;
 274                     }
 275
 276                     size_t len = gseq.length();
 277                     assert(len > 0);
 278                     for (size_t j = 0; j < len; j++) gseq[j] = check(gseq[j]);
 279
 280                     iter = sn2tr.find(seqname);
 281                     if (iter == sn2tr.end()) continue;
 282
 283                     chrvec.push_back(ChrInfo(seqname, len));
 284
 285                     vector<int>& vec = iter->second;
 286                     int s = vec.size();
 287                     for (int j = 0; j < s; j++) {
 288                         assert(vec[j] > 0 && vec[j] <= M);
 289                         transcripts.getTranscriptAt(vec[j]).extractSeq(gseq, seqs[vec[j]]);
 290                     }
 291                 }
 292                 fin.close();
 293
 294                 if (verbose) { printf("%s is processed!\n", argv[i]); }
 295         }
 296
 297         for (int i = 1; i <= M; i++) {
 298                 if (seqs[i] == "") {
 299                         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());
 300                         exit(-1);
 301                 }
 302         }
 303
 304         sort(chrvec.begin(), chrvec.end());
 305
 306         if (verbose) { printf("Extracting sequences is done!\n"); }
 307
 308         writeResults(argv[1]);
 309
 310         return 0;
 311 }