[ViewVC] Diff of: gclib/tophat_cpp/bwt

+        return false;
-+
+bool SplicedSAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,
-+
+                                               BowtieHit& bh,
-+
+                                               bool strip_slash,
-+
+                                               char* name_out,
-+
+                                               char* name_tags,
-+
+                                               char* seq,
-+
+                                               char* qual)
-+
+{
-+
+        if (!orig_bwt_buf || !*orig_bwt_buf)
-+
+                return false;
-+
-+
+        char bwt_buf[2048];
-+
+        strcpy(bwt_buf, orig_bwt_buf);
-+
+        // Are we still in the header region?
-+
+        if (bwt_buf[0] == '@')
-+
+                return false;
-+
-+
+        //char orientation;
-+
+        /*
-+
+        char mismatches[buf_size];
-+
-+
+        //char* orientation_str = get_token((char**)&buf,"\t");
-+
-+
+        mismatches[0] = 0;
-+
+        char* mismatches_str = get_token((char**)&buf, "\t");
-+
+        if (mismatches_str)
-+
+                strcpy(mismatches, mismatches_str);
-+
-+
+        orientation = orientation_str[0];
-+
+        text_offset = atoi(text_offset_str);
-+
+    */
-+
-+
+        char* buf = bwt_buf;
-+
+        char* name = get_token((char**)&buf,"\t");
-+
+        char* sam_flag_str = get_token((char**)&buf,"\t");
-+
+        char* text_name = get_token((char**)&buf,"\t");
-+
+        char* text_offset_str = get_token((char**)&buf,"\t");
-+
+        const char* map_qual_str = get_token((char**)&buf,"\t");
-+
+        char* cigar_str = get_token((char**)&buf,"\t");
-+
+        const char* mate_ref_str =  get_token((char**)&buf,"\t");
-+
+        const char* mate_pos_str =  get_token((char**)&buf,"\t");
-+
+        const char* inferred_insert_sz_str =  get_token((char**)&buf,"\t");
-+
-+
+        const char* seq_str =  get_token((char**)&buf,"\t");
-+
+        if (seq)
-+
+          strcpy(seq, seq_str);
-+
-+
+        const char* qual_str =  get_token((char**)&buf,"\t");
-+
+        if (qual)
-+
+          strcpy(qual, qual_str);
-+
-+
+        if (!name ||
-+
+                !sam_flag_str ||
-+
+                !text_name ||
-+
+                !text_offset_str ||
-+
+                !map_qual_str ||
-+
+                !cigar_str ||
-+
+                !mate_ref_str ||
-+
+                !mate_pos_str ||
-+
+                !inferred_insert_sz_str ||
-+
+                !seq_str ||
-+
+                !qual_str)
-+
+        {
-+
+                // truncated or malformed SAM record
-+
+                return false;
-+
+        }
-+
-+
-+
+        int sam_flag = atoi(sam_flag_str);
-+
+        int text_offset = atoi(text_offset_str);
-+
-+
+        //##############################################
-+
+        bool end = true;
-+
+        unsigned int seg_offset = 0;
-+
+        unsigned int seg_num = 0;
-+
+        unsigned int num_segs = 0;
-+
-+
+        // Copy the tag out of the name field before we might wipe it out
-+
+        char* pipe = strrchr(name, '|');
-+
+        if (pipe)
-+
+        {
-+
+                if (name_tags)
-+
+                  strcpy(name_tags, pipe);
-+
-+
+                char* tag_buf = pipe + 1;
-+
+                if (strchr(tag_buf, ':'))
-+
+                  {
-+
+                    sscanf(tag_buf, "%u:%u:%u", &seg_offset, &seg_num, &num_segs);
-+
+                    if (seg_num + 1 == num_segs)
-+
+                      end = true;
-+
+                    else
-+
+                      end = false;
-+
+                  }
-+
-+
+                *pipe = 0;
-+
+        }
-+
+        // Stripping the slash and number following it gives the insert name
-+
+        char* slash = strrchr(name, '/');
-+
+        if (strip_slash && slash)
-+
+                *slash = 0;
-+
-+
+        //int read_len = strlen(sequence);
-+
-+
+        // Add this alignment to the table of hits for this half of the
-+
+        // Bowtie map
-+
-+
+        // Parse the text_name field to recover the splice coords
-+
+        vector<string> toks;
-+
-+
+        tokenize_strict(text_name, "|", toks);
-+
-+
+        int num_extra_toks = (int)toks.size() - 6;
-+
-+
+        if (num_extra_toks >= 0)
-+
+        {
-+
+                static const uint8_t left_window_edge_field = 1;
-+
+                static const uint8_t splice_field = 2;
-+
+                //static const uint8_t right_window_edge_field = 3;
-+
+                static const uint8_t junction_type_field = 4;
-+
+                static const uint8_t strand_field = 5;
-+
-+
+                string contig = toks[0];
-+
+                for (int t = 1; t <= num_extra_toks; ++t)
-+
+                {
-+
+                        contig += "|";
-+
+                        contig += toks[t];
-+
+                }
-+
-+
+                vector<string> splice_toks;
-+
+                tokenize(toks[num_extra_toks + splice_field], "-", splice_toks);
-+
+                if (splice_toks.size() != 2)
-+
+                {
-+
+                        fprintf(stderr, "Warning: found malformed splice record, skipping:\n");
-+
+                        //fprintf(stderr, "%s (token: %s)\n", text_name,
-+
+                        //        toks[num_extra_toks + splice_field].c_str());
-+
+                        return false;
-+
+                }
-+
-+
+                //
-+
+                // check for an insertion hit
-+
+                //
-+
+                if(toks[num_extra_toks + junction_type_field] == "ins")
-+
+                  {
-+
+                        int8_t spliced_read_len = strlen(seq_str);
-+
+                        /*
-+
+                         * The 0-based position of the left edge of the alignment. Note that this
-+
+                         * value may need to be futher corrected to account for the presence of
-+
+                         * of the insertion.
-+
+                        */
-+
+                        uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
-+
+                        uint32_t right = left + spliced_read_len - 1;
-+
-+
+                        /*
-+
+                         * The 0-based position of the last genomic sequence before the insertion
-+
+                         */
-+
+                        uint32_t left_splice_pos = atoi(splice_toks[0].c_str());
-+
-+
+                        string insertedSequence = splice_toks[1];
-+
+                        /*
-+
+                         * The 0-based position of the first genomic sequence after teh insertion
-+
+                         */
-+
+                        uint32_t right_splice_pos = left_splice_pos + 1;
-+
+                        if(left > left_splice_pos){
-+
+                                /*
-+
+                                 * The genomic position of the left edge of the alignment needs to be corrected
-+
+                                 * If the alignment does not extend into the insertion, simply subtract the length
-+
+                                 * of the inserted sequence, otherwise, just set it equal to the right edge
-+
+                                 */
-+
+                                left = left - insertedSequence.length();
-+
+                                if(left < right_splice_pos){
-+
+                                        left = right_splice_pos;
-+
+                                }
-+
+                        }
-+
+                        if(right > left_splice_pos){
-+
+                                right = right - insertedSequence.length();
-+
+                                if(right < left_splice_pos){
-+
+                                        right = left_splice_pos;
-+
+                                }
-+
+                        }
-+
+                        /*
-+
+                         * Now, right and left should be properly transformed into genomic coordinates
-+
+                         * We should be able to deduce how much the alignment matches the insertion
-+
+                         * simply based on the length of the read
-+
+                         */
-+
+                        int left_match_length = 0;
-+
+                        if(left <= left_splice_pos){
-+
+                                left_match_length = left_splice_pos - left + 1;
-+
+                        }
-+
+                        int right_match_length = 0;
-+
+                        if(right >= right_splice_pos){
-+
+                                right_match_length = right - right_splice_pos + 1;
-+
+                        }
-+
+                        int insertion_match_length = spliced_read_len - left_match_length - right_match_length;
-+
-+
+                        if(left_match_length <= 0 || right_match_length <= 0 || insertion_match_length <= 0)
-+
+                          return false;
-+
-+
+                        string junction_strand = toks[num_extra_toks + strand_field];
-+
+                        if(junction_strand != "rev" && junction_strand != "fwd"){
-+
+                                fprintf(stderr, "Malformed insertion record\n");
-+
+                                return false;
-+
+                        }
-+
-+
+                        char* pch = strtok( mismatches, ",");
-+
+                        unsigned char num_mismatches = 0;
-+
-+
+                        /*
-+
+                         * remember that text_offset holds the left end of the
-+
+                         * alignment, relative to the start of the contig
-+
+                         */
-+
-+
+                        /*
-+
+                         * The 0-based relative position of the left-most character
-+
+                         * before the insertion in the contig
-+
+                         */
-+
+                        int relative_splice_pos = left_splice_pos - atoi(toks[num_extra_toks + left_window_edge_field].c_str());
-+
+                        while (pch != NULL)
-+
+                        {
-+
+                                char* colon = strchr(pch, ':');
-+
+                                if (colon)
-+
+                                {
-+
+                                        *colon = 0;
-+
+                                        int mismatch_pos = atoi(pch);
-+
-+
+                                        /*
-+
+                                         * for reversely mapped reads,
-+
+                                         * find the correct mismatched position.
-+
+                                         */
-+
+                                        if (sam_flag & 0x0010){
-+
+                                          mismatch_pos = spliced_read_len - mismatch_pos - 1;
-+
+                                        }
-+
-+
+                                        /*
-+
+                                         * Only count mismatches outside of the insertion region
-+
+                                         * If there is a mismatch within the insertion,
-+
+                                         * disallow this hit
-+
+                                         */
-+
+                                        if(mismatch_pos + text_offset <= relative_splice_pos || mismatch_pos + text_offset > relative_splice_pos + insertedSequence.length()){
-+
+                                          num_mismatches++;
-+
+                                        }else{
-+
+                                          return false;
-+
+                                        }
-+
+                                }
-+
+                                pch = strtok (NULL, ",");
-+
+                        }
-+
-+
-+
+                        vector<CigarOp> cigar;
-+
+                        cigar.push_back(CigarOp(MATCH, left_match_length));
-+
+                        cigar.push_back(CigarOp(INS, insertion_match_length));
-+
+                        cigar.push_back(CigarOp(MATCH, right_match_length));
-+
-+
+                        /*
-+
+                         * For now, disallow hits that don't span
-+
+                         * the insertion. If we allow these types of hits,
-+
+                         * then long_spanning.cpp needs to be updated
-+
+                         * in order to intelligently merge these kinds
-+
+                         * of reads back together
-+
+                         *
-+
+                         * Following code has been changed to allow segment that end
-+
+                         * in an insertion
-+
+                         */
-+
+                        bh = create_hit(name,
-+
+                                        contig,
-+
+                                        left,
-+
+                                        cigar,
-+
+                                        sam_flag & 0x0010,  //#######
-+
+                                        junction_strand == "rev",
-+
+                                        num_mismatches,
-+
+,
-+
+                                        end);
-+
+                        return true;
-+
+                }
-+
-+
+                else
-+
+                  {
-+
+                    uint32_t left = atoi(toks[num_extra_toks + left_window_edge_field].c_str()) + text_offset;
-+
+                    int spliced_read_len = strlen(seq_str);
-+
+                    int8_t left_splice_pos = atoi(splice_toks[0].c_str()) - left + 1;
-+
+                    if(left_splice_pos > spliced_read_len) left_splice_pos = spliced_read_len;
-+
+                    int8_t right_splice_pos = spliced_read_len - left_splice_pos;
-+
-+
+                    if (right_splice_pos <= 0 || left_splice_pos <= 0)
-+
+                      return false;
-+
-+
+                    if ((sam_flag & 0x0010) == 0) //######
-+
+                      {
-+
+                        if (left_splice_pos + seg_offset < _anchor_length){
-+
+                          return false;
-+
+                        }
-+
+                      }
-+
+                    else
-+
+                      {
-+
+                        if (right_splice_pos + seg_offset < _anchor_length)
-+
+                          return false;
-+
+                      }
-+
+                    //uint32_t right = atoi(splice_toks[1].c_str()) + right_splice_pos;
-+
+                    //atoi(toks[right_window_edge_field].c_str());
-+
+                    int gap_len = atoi(splice_toks[1].c_str()) - atoi(splice_toks[0].c_str()) - 1;
-+
-+
+                    string junction_strand = toks[num_extra_toks + strand_field];
-+
+                    if (!(junction_strand == "rev" || junction_strand == "fwd")||
-+
+                        !(orientation == '-' || orientation == '+'))
-+
+                      {
-+
+                        fprintf(stderr, "Warning: found malformed splice record, skipping\n");
-+
+                        //fprintf(stderr, "junction_strand=%s, orientation='%c'\n",
-+
+                        //           junction_strand.c_str(), orientation);
-+
+                        return false;
-+
+                      }
-+
-+
+                    //vector<string> mismatch_toks;
-+
+                    char* pch = strtok (mismatches,",");
-+
+                    int mismatches_in_anchor = 0;
-+
+                    unsigned char num_mismatches = 0;
-+
+                    while (pch != NULL)
-+
+                      {
-+
+                        char* colon = strchr(pch, ':');
-+
+                        if (colon)
-+
+                          {
-+
+                            *colon = 0;
-+
+                            num_mismatches++;
-+
+                            int mismatch_pos = atoi(pch);
-+
+                            if ((orientation == '+' && abs(mismatch_pos - left_splice_pos) < (int)min_anchor_len) ||
-+
+                                (orientation == '-' && abs(((int)spliced_read_len - left_splice_pos + 1) - mismatch_pos)) < (int)min_anchor_len)
-+
+                              mismatches_in_anchor++;
-+
+                          }
-+
+                        //mismatch_toks.push_back(pch);
-+
+                        pch = strtok (NULL, ",");
-+
+                      }
-+
-+
+                    // FIXME: we probably should exclude these hits somewhere, but this
-+
+                    // isn't the right place
-+
+                    vector<CigarOp> cigar;
-+
+                    cigar.push_back(CigarOp(MATCH, left_splice_pos));
-+
+                    if(toks[num_extra_toks + junction_type_field] == "del"){
-+
+                      cigar.push_back(CigarOp(DEL, gap_len));
-+
+                    }else{
-+
+                      cigar.push_back(CigarOp(REF_SKIP, gap_len));
-+
+                    }
-+
+                    cigar.push_back(CigarOp(MATCH, right_splice_pos));
-+
-+
+                    bh = create_hit(name,
-+
+                                    contig,
-+
+                                    left,
-+
+                                    cigar,
-+
+                                    orientation == '-',
-+
+                                    junction_strand == "rev",
-+
+                                    num_mismatches,
-+
+                                    mismatches_in_anchor,
-+
+                                    end);
-+
+                    return true;
-+
+                  }
-+
+        }
-+
+        else
-+
+        {
-+
+          fprintf(stderr, "Warning: found malformed splice record, skipping\n");
-+
+          //fprintf(stderr, "%s\n", orig_bwt_buf);
-+
+          //                    continue;
-+
+                return false;
-+
+        }
-+
-+
+        return false;
-+
+}
+bool BAMHitFactory::get_hit_from_buf(const char* orig_bwt_buf,

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines