| 1 |
/* |
| 2 |
* Author: Harold Pimentel |
| 3 |
* Contact: http://cs.berkeley.edu/~pimentel |
| 4 |
* Date: June 10, 2011 |
| 5 |
*/ |
| 6 |
|
| 7 |
#ifndef _MAP2GTF_H_ |
| 8 |
#define _MAP2GTF_H_ |
| 9 |
|
| 10 |
#ifdef HAVE_CONFIG_H |
| 11 |
#include <config.h> |
| 12 |
#endif |
| 13 |
|
| 14 |
#include <cassert> |
| 15 |
#include <cstdlib> |
| 16 |
#include <cstring> |
| 17 |
#include <fstream> |
| 18 |
#include <iostream> |
| 19 |
#include <string> |
| 20 |
#include <sstream> |
| 21 |
#include <vector> |
| 22 |
|
| 23 |
#include <bam/bam.h> |
| 24 |
#include <bam/sam.h> |
| 25 |
|
| 26 |
#include <seqan/sequence.h> |
| 27 |
|
| 28 |
#include <getopt.h> |
| 29 |
#include <unistd.h> |
| 30 |
|
| 31 |
#include "bwt_map.h" |
| 32 |
#include "common.h" |
| 33 |
#include "gff.h" |
| 34 |
|
| 35 |
#define MAX_READ_NAME_LEN 2048 |
| 36 |
|
| 37 |
/* |
| 38 |
* XXX: This class currently assumes someone used the script in TopHat to map |
| 39 |
* the reads already. It also depends on that same format. |
| 40 |
*/ |
| 41 |
class TranscriptomeHit; |
| 42 |
class Map2GTF |
| 43 |
{ |
| 44 |
public: |
| 45 |
Map2GTF(const std::string& gtf_fname, const std::string& sam_fname); |
| 46 |
~Map2GTF(); |
| 47 |
// Write out to a BAM file |
| 48 |
bool next_read_hits(std::vector<bam1_t*>& hits, size_t& num_hits); |
| 49 |
void convert_coords(const std::string& out_fname, const std::string& sam_header); |
| 50 |
void trans_to_genomic_coords(TranscriptomeHit& hit); |
| 51 |
|
| 52 |
private: |
| 53 |
GffReader gtfReader_; |
| 54 |
|
| 55 |
std::string gtf_fname_; |
| 56 |
std::string in_fname_; |
| 57 |
|
| 58 |
FILE* gtf_fhandle_; |
| 59 |
samfile_t* in_fhandle_; |
| 60 |
bam_header_t* in_sam_header_; |
| 61 |
|
| 62 |
map<string, int> ref_to_id_; |
| 63 |
bam_header_t* out_sam_header_; |
| 64 |
|
| 65 |
ReadTable readTable_; |
| 66 |
RefSequenceTable refSeqTable_; |
| 67 |
|
| 68 |
Map2GTF(); // Don't want anyone calling the constructor w/o options |
| 69 |
}; |
| 70 |
|
| 71 |
class TranscriptomeHit |
| 72 |
{ |
| 73 |
public: |
| 74 |
bam1_t* hit; |
| 75 |
GffObj* trans; |
| 76 |
TranscriptomeHit(bam1_t* h = NULL, GffObj* t = NULL) |
| 77 |
{ |
| 78 |
hit = h; |
| 79 |
trans = t; |
| 80 |
} |
| 81 |
bool operator==(const TranscriptomeHit& th) const |
| 82 |
{ |
| 83 |
if (hit->core.tid != th.hit->core.tid) |
| 84 |
return false; |
| 85 |
|
| 86 |
if (hit->core.pos != th.hit->core.pos) |
| 87 |
return false; |
| 88 |
|
| 89 |
if (hit->core.n_cigar != th.hit->core.n_cigar) |
| 90 |
return false; |
| 91 |
|
| 92 |
for (int i = 0; i < hit->core.n_cigar; ++i) |
| 93 |
{ |
| 94 |
if (bam1_cigar(hit)[i] != bam1_cigar(th.hit)[i]) |
| 95 |
return false; |
| 96 |
} |
| 97 |
|
| 98 |
return true; |
| 99 |
} |
| 100 |
bool operator<(const TranscriptomeHit& th) const |
| 101 |
{ |
| 102 |
if (hit->core.tid != th.hit->core.tid) |
| 103 |
return hit->core.tid < th.hit->core.tid; |
| 104 |
|
| 105 |
if (hit->core.pos != th.hit->core.pos) |
| 106 |
return hit->core.pos < th.hit->core.pos; |
| 107 |
|
| 108 |
if (hit->core.n_cigar != th.hit->core.n_cigar) |
| 109 |
return hit->core.n_cigar < th.hit->core.n_cigar; |
| 110 |
|
| 111 |
for (int i = 0; i < hit->core.n_cigar; ++i) |
| 112 |
{ |
| 113 |
if (bam1_cigar(hit)[i] != bam1_cigar(th.hit)[i]) |
| 114 |
return bam1_cigar(hit)[i] < bam1_cigar(th.hit)[i]; |
| 115 |
} |
| 116 |
|
| 117 |
return false; |
| 118 |
} |
| 119 |
}; |
| 120 |
|
| 121 |
bool get_read_start(GList<GffExon>* exon_list, size_t gtf_start, |
| 122 |
size_t& genome_start, int& exon_idx); |
| 123 |
|
| 124 |
void print_trans(GffObj* trans, const bam1_t* in, size_t rem_len, |
| 125 |
size_t match_len, size_t cur_pos, size_t start_pos); |
| 126 |
|
| 127 |
#endif /* _MAP2GTF_H_ */ |
