File: SeqPlot.cpp

package info (click to toggle)
libseqlib 1.1.2+dfsg-1~bpo9+1
  • links: PTS, VCS
  • area: main
  • in suites: stretch-backports
  • size: 1,460 kB
  • sloc: cpp: 7,176; sh: 805; makefile: 60
file content (85 lines) | stat: -rw-r--r-- 2,211 bytes parent folder | download | duplicates (3)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
#include "SeqLib/SeqPlot.h"

namespace SeqLib {

std::string SeqPlot::PlotAlignmentRecords(const BamRecordVector& brv) const {

  PlottedReadVector plot_vec;

  for (BamRecordVector::const_iterator i = brv.begin(); i != brv.end(); ++i) {
    
    // get the position in the view window
    if (i->ChrID() != m_view.chr)
      continue;

    int pos = i->Position() - m_view.pos1;
    if (pos < 0)
      continue;

    if (i->PositionEnd() > m_view.pos2)
      continue;

    // plot with gaps
    std::string tseq = i->Sequence();
    std::string gapped_seq;

    size_t p = i->AlignmentPosition(); // move along on sequence, starting at first non-clipped base
    Cigar cc = i->GetCigar();
    for (Cigar::const_iterator c = cc.begin(); c != cc.end(); ++c) {
      if (c->Type() == 'M') { // 
	assert(p + c->Length() <= tseq.length());
	gapped_seq += tseq.substr(p, c->Length());
      } else if (c->Type() == 'D') {
	gapped_seq += std::string(c->Length(), '-');
      }

      if (c->Type() == 'I' || c->Type() == 'M')
	p += c->Length();
    }

    std::stringstream msg;
    msg << i->Qname() << ">>>" << (i->ChrID() + 1) << ":" << i->Position();
      
    // add to the read plot
    plot_vec.push_back(PlottedRead(pos, gapped_seq, msg.str()));
    
  }

  // sort them
  std::sort(plot_vec.begin(), plot_vec.end());

  // make a list of lines
  PlottedReadLineVector line_vec;

  // plot the reads from the ReadPlot vector
  for (PlottedReadVector::iterator i = plot_vec.begin(); i != plot_vec.end(); ++i) {
    bool found = false;
    for (PlottedReadLineVector::iterator j = line_vec.begin(); j != line_vec.end(); ++j) {
      if (j->readFits(*i)) { // it fits here
	j->addRead(&(*i));
	found = true;
	break;
      }
    }
    if (!found) { // didn't fit anywhere, so make a new line
      PlottedReadLine prl;
      prl.pad = m_pad;
      prl.contig_len = m_view.Width(); //ac.getSequence().length();
      prl.addRead(&(*i));
      line_vec.push_back(prl);
    }
  }
  
  std::stringstream ss;
  
  // plot the lines. Add contig identifier to each
  for (PlottedReadLineVector::const_iterator i = line_vec.begin(); i != line_vec.end(); ++i) 
    ss << (*i) << std::endl;

  return ss.str();

  
}


}