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