gclib/tophat_cpp/insertions.cpp

/*
 *  insertions.cpp
 *  TopHat
 *
 *  Created by Ryan Kelley on 11/04/2010.
 *
 */

#ifdef HAVE_CONFIG_H
#include <config.h>
#else
#define PACKAGE_VERSION "INTERNAL"
#define SVN_REVISION "XXX"
#endif


#include <cassert>
#include <cstdio>
#include <cstring>
#include <vector>
#include <string>
#include <map>
#include <algorithm>
#include <set>
#include <stdexcept>
#include <iostream>
#include <fstream>
#include <seqan/sequence.h>
#include <seqan/find.h>
#include <seqan/file.h>
#include <getopt.h>

#include "common.h"
#include "bwt_map.h"
#include "junctions.h"
#include "insertions.h"
#include "fragments.h"
#include "wiggles.h"
#include "tokenize.h"
#include "reads.h"

#include "inserts.h"

/**
 * Add insertions from an alignment to an InsertionSet.
 * This will look for insertion in the alignment specified by bh. If the 
 * insertion is already in insertions, it will updated the count. Otherwise,
 * it will add the insertion to the set and initialize the count to 1.
 * @param bh The bowtie hit to be used to specify alignment infromation.
 * @param insertions The InsertionSet that will be updated with the insertion information from teh alignment.
 */
void insertions_from_alignment(const BowtieHit& bh, InsertionSet& insertions) {
  vector<pair<Insertion, InsertionStats> > new_insertions;
  insertions_from_spliced_hit(bh, new_insertions);

  for(size_t i = 0; i < new_insertions.size(); ++i) {
    const pair<Insertion, InsertionStats>& insertion = new_insertions[i];
    InsertionSet::iterator itr = insertions.find(insertion.first);
    if (itr != insertions.end()) {
      itr->second.supporting_hits += 1;
      itr->second.left_extent = max(itr->second.left_extent, insertion.second.left_extent);
      itr->second.right_extent = max(itr->second.right_extent, insertion.second.right_extent);
    }
    else {
      assert(insertion.first.refid != VMAXINT32);
      insertions[insertion.first] = insertion.second;
    }
  }
  return;
}

/**
 * Print insertions in BED format.
 * Note, as per the BED-standard (http://genome.ucsc.edu/FAQ/FAQformat)
 *   -The coordinates should be 0-based
 *   -The chromEnd field should not include the actual feature
 *   -The name will be the inserted sequence
 *   -The score will be the number of supporing counts, which is capped at 1,000
 * By (my) convention, the chromStart will be the last genome postion
 * before hte insertio.
 * 
 * <chrom>\t<left>\t<left>\t<inserted sequence>\t<read count>\n
 * @param insertions_out The output file
 * @param insertions Maps from insertions to number of supporting reads
 * @param ref_sequences The table of reference sequences 
 */
void print_insertions(FILE* insertions_out, const InsertionSet& insertions, RefSequenceTable& ref_sequences) {
        fprintf(insertions_out, "track name=insertions description=\"TopHat insertions\"\n");
        for(InsertionSet::const_iterator i = insertions.begin(); i != insertions.end(); ++i){
          int counts = i->second.supporting_hits;
          if(counts > 1000){
            counts = 1000;
          }
          fprintf(insertions_out, "%s\t%d\t%d\t%s\t%d\n",
                  ref_sequences.get_name(i->first.refid),
                  i->first.left,
                  i->first.left,
                  (i->first.sequence).c_str(),
                  counts);
        }
}

/**
 * Extract a list of insertions from a bowtie hit.
 * Given a bowtie hit, extract a vector of insertions.  
 * @param bh The bowtie hit to use for alignment information.
 * @param insertions Used to store the resultant vector of insertions.
 */
void insertions_from_spliced_hit(const BowtieHit& bh, vector<pair<Insertion, InsertionStats> >& insertions){
        const vector<CigarOp>& cigar = bh.cigar();
        unsigned int positionInGenome = bh.left();
        unsigned int positionInRead = 0;

        bool bSawFusion = false;
        for(size_t c = 0; c < cigar.size(); ++c){
          switch(cigar[c].opcode){
          case REF_SKIP:
            positionInGenome += cigar[c].length;
            break;
          case rEF_SKIP:
            positionInGenome -= cigar[c].length;
            break;
          case MATCH:
          case mATCH:
            if (cigar[c].opcode == MATCH)
              positionInGenome += cigar[c].length;
            else
              positionInGenome -= cigar[c].length;
            positionInRead += cigar[c].length;
            break;
          case DEL:
            positionInGenome += cigar[c].length;
            break;
          case dEL:
            positionInGenome -= cigar[c].length;
            break;
          case INS:
          case iNS:
            {
              Insertion insertion;
              InsertionStats stats;
              
              /*
               * Note that the reported position in the genome from the SAM
               * alignment is 1-based, since the insertion object is expecting
               * a 0-based co-ordinate, we need to subtract 1
               */
              if (bSawFusion)
                insertion.refid = bh.ref_id2();
              else
                insertion.refid = bh.ref_id();
              
              if (cigar[c].opcode == INS)
                insertion.left = positionInGenome - 1;
              else
                insertion.left = positionInGenome + 1;
              
              insertion.sequence = bh.seq().substr(positionInRead, cigar[c].length);
              stats.supporting_hits = 1;
              if (c > 0)
                stats.left_extent = cigar[c-1].length;
              if (c + 1 < cigar.size())
                stats.right_extent = cigar[c+1].length;
              
              insertions.push_back(make_pair(insertion, stats));
              positionInRead += cigar[c].length;
            }
            break;
          case FUSION_FF:
          case FUSION_FR:
          case FUSION_RF:
            bSawFusion = true;
            positionInGenome = cigar[c].length;
            break;
          default:
            break;
                }       
        }       
        return;
} 

void merge_with(InsertionSet& insertions, const InsertionSet& other)
{
  for (InsertionSet::const_iterator insertion = other.begin(); insertion != other.end(); ++insertion)
    {
      InsertionSet::iterator itr = insertions.find(insertion->first);
      if (itr != insertions.end())
        {
          itr->second.merge_with(insertion->second);
        }
      else
        {
          insertions[insertion->first] = insertion->second;
        }
    }
}
Revision:	239
Committed:	Mon May 14 21:01:08 2012 UTC (12 years, 3 months ago) by gpertea
File size:	5798 byte(s)
Log Message:	wip prep_reads.cpp, tophat
Line	File contents
1	/*
2	* insertions.cpp
3	* TopHat
4	*
5	* Created by Ryan Kelley on 11/04/2010.
6	*
7	*/
8
9	#ifdef HAVE_CONFIG_H
10	#include <config.h>
11	#else
12	#define PACKAGE_VERSION "INTERNAL"
13	#define SVN_REVISION "XXX"
14	#endif
15
16
17	#include <cassert>
18	#include <cstdio>
19	#include <cstring>
20	#include <vector>
21	#include <string>
22	#include <map>
23	#include <algorithm>
24	#include <set>
25	#include <stdexcept>
26	#include <iostream>
27	#include <fstream>
28	#include <seqan/sequence.h>
29	#include <seqan/find.h>
30	#include <seqan/file.h>
31	#include <getopt.h>
32
33	#include "common.h"
34	#include "bwt_map.h"
35	#include "junctions.h"
36	#include "insertions.h"
37	#include "fragments.h"
38	#include "wiggles.h"
39	#include "tokenize.h"
40	#include "reads.h"
41
42	#include "inserts.h"
43
44	/**
45	* Add insertions from an alignment to an InsertionSet.
46	* This will look for insertion in the alignment specified by bh. If the
47	* insertion is already in insertions, it will updated the count. Otherwise,
48	* it will add the insertion to the set and initialize the count to 1.
49	* @param bh The bowtie hit to be used to specify alignment infromation.
50	* @param insertions The InsertionSet that will be updated with the insertion information from teh alignment.
51	*/
52	void insertions_from_alignment(const BowtieHit& bh, InsertionSet& insertions) {
53	vector<pair<Insertion, InsertionStats> > new_insertions;
54	insertions_from_spliced_hit(bh, new_insertions);
55
56	for(size_t i = 0; i < new_insertions.size(); ++i) {
57	const pair<Insertion, InsertionStats>& insertion = new_insertions[i];
58	InsertionSet::iterator itr = insertions.find(insertion.first);
59	if (itr != insertions.end()) {
60	itr->second.supporting_hits += 1;
61	itr->second.left_extent = max(itr->second.left_extent, insertion.second.left_extent);
62	itr->second.right_extent = max(itr->second.right_extent, insertion.second.right_extent);
63	}
64	else {
65	assert(insertion.first.refid != VMAXINT32);
66	insertions[insertion.first] = insertion.second;
67	}
68	}
69	return;
70	}
71
72	/**
73	* Print insertions in BED format.
74	* Note, as per the BED-standard (http://genome.ucsc.edu/FAQ/FAQformat)
75	* -The coordinates should be 0-based
76	* -The chromEnd field should not include the actual feature
77	* -The name will be the inserted sequence
78	* -The score will be the number of supporing counts, which is capped at 1,000
79	* By (my) convention, the chromStart will be the last genome postion
80	* before hte insertio.
81	*
82	* <chrom>\t<left>\t<left>\t<inserted sequence>\t<read count>\n
83	* @param insertions_out The output file
84	* @param insertions Maps from insertions to number of supporting reads
85	* @param ref_sequences The table of reference sequences
86	*/
87	void print_insertions(FILE* insertions_out, const InsertionSet& insertions, RefSequenceTable& ref_sequences) {
88	fprintf(insertions_out, "track name=insertions description=\"TopHat insertions\"\n");
89	for(InsertionSet::const_iterator i = insertions.begin(); i != insertions.end(); ++i){
90	int counts = i->second.supporting_hits;
91	if(counts > 1000){
92	counts = 1000;
93	}
94	fprintf(insertions_out, "%s\t%d\t%d\t%s\t%d\n",
95	ref_sequences.get_name(i->first.refid),
96	i->first.left,
97	i->first.left,
98	(i->first.sequence).c_str(),
99	counts);
100	}
101	}
102
103	/**
104	* Extract a list of insertions from a bowtie hit.
105	* Given a bowtie hit, extract a vector of insertions.
106	* @param bh The bowtie hit to use for alignment information.
107	* @param insertions Used to store the resultant vector of insertions.
108	*/
109	void insertions_from_spliced_hit(const BowtieHit& bh, vector<pair<Insertion, InsertionStats> >& insertions){
110	const vector<CigarOp>& cigar = bh.cigar();
111	unsigned int positionInGenome = bh.left();
112	unsigned int positionInRead = 0;
113
114	bool bSawFusion = false;
115	for(size_t c = 0; c < cigar.size(); ++c){
116	switch(cigar[c].opcode){
117	case REF_SKIP:
118	positionInGenome += cigar[c].length;
119	break;
120	case rEF_SKIP:
121	positionInGenome -= cigar[c].length;
122	break;
123	case MATCH:
124	case mATCH:
125	if (cigar[c].opcode == MATCH)
126	positionInGenome += cigar[c].length;
127	else
128	positionInGenome -= cigar[c].length;
129	positionInRead += cigar[c].length;
130	break;
131	case DEL:
132	positionInGenome += cigar[c].length;
133	break;
134	case dEL:
135	positionInGenome -= cigar[c].length;
136	break;
137	case INS:
138	case iNS:
139	{
140	Insertion insertion;
141	InsertionStats stats;
142
143	/*
144	* Note that the reported position in the genome from the SAM
145	* alignment is 1-based, since the insertion object is expecting
146	* a 0-based co-ordinate, we need to subtract 1
147	*/
148	if (bSawFusion)
149	insertion.refid = bh.ref_id2();
150	else
151	insertion.refid = bh.ref_id();
152
153	if (cigar[c].opcode == INS)
154	insertion.left = positionInGenome - 1;
155	else
156	insertion.left = positionInGenome + 1;
157
158	insertion.sequence = bh.seq().substr(positionInRead, cigar[c].length);
159	stats.supporting_hits = 1;
160	if (c > 0)
161	stats.left_extent = cigar[c-1].length;
162	if (c + 1 < cigar.size())
163	stats.right_extent = cigar[c+1].length;
164
165	insertions.push_back(make_pair(insertion, stats));
166	positionInRead += cigar[c].length;
167	}
168	break;
169	case FUSION_FF:
170	case FUSION_FR:
171	case FUSION_RF:
172	bSawFusion = true;
173	positionInGenome = cigar[c].length;
174	break;
175	default:
176	break;
177	}
178	}
179	return;
180	}
181
182	void merge_with(InsertionSet& insertions, const InsertionSet& other)
183	{
184	for (InsertionSet::const_iterator insertion = other.begin(); insertion != other.end(); ++insertion)
185	{
186	InsertionSet::iterator itr = insertions.find(insertion->first);
187	if (itr != insertions.end())
188	{
189	itr->second.merge_with(insertion->second);
190	}
191	else
192	{
193	insertions[insertion->first] = insertion->second;
194	}
195	}
196	}