bpp-seq-omics  2.4.1
ConcatenateMafIterator.cpp
Go to the documentation of this file.
1 //
2 // File: ConcatenateMafIterator.cpp
3 // Authors: Julien Dutheil
4 // Created: Tue Sep 07 2010
5 //
6 
7 /*
8 Copyright or © or Copr. Bio++ Development Team, (2010)
9 
10 This software is a computer program whose purpose is to provide classes
11 for sequences analysis.
12 
13 This software is governed by the CeCILL license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's author, the holder of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL license and that you accept its terms.
38 */
39 
40 #include "ConcatenateMafIterator.h"
41 
42 using namespace bpp;
43 
44 //From the STL:
45 #include <string>
46 #include <numeric>
47 
48 //From bpp-core:
50 
51 using namespace std;
52 
54 {
55  if (!incomingBlock_) return 0;
56  currentBlock_ = incomingBlock_;
57  incomingBlock_ = iterator_->nextBlock();
58  size_t count = 1;
59  if (verbose_)
60  ApplicationTools::displayMessage("Concatenating new block...");
61  while (incomingBlock_ &&
62  (refSpecies_ == "" ||
63  (incomingBlock_->hasSequenceForSpecies(refSpecies_) &&
64  currentBlock_->hasSequenceForSpecies(refSpecies_) &&
65  incomingBlock_->getSequenceForSpecies(refSpecies_).getChromosome() ==
66  currentBlock_->getSequenceForSpecies(refSpecies_).getChromosome()
67  )
68  )
69  )
70  {
71  if (currentBlock_->getNumberOfSites() >= minimumSize_) {
72  return currentBlock_;
73  }
74  if (verbose_) {
75  ApplicationTools::displayUnlimitedGauge(count++, "Concatenating...");
76  }
77 
78  //We merge the two blocks:
79  vector<string> sp1 = currentBlock_->getSpeciesList();
80  vector<string> sp2 = incomingBlock_->getSpeciesList();
81  vector<string> allSp = VectorTools::unique(VectorTools::vectorUnion(sp1, sp2));
82  //We need to create a new MafBlock:
83  MafBlock* mergedBlock = new MafBlock();
84  //We average the score and pass values:
85  unsigned int p1 = currentBlock_->getPass();
86  unsigned int p2 = incomingBlock_->getPass();
87  if (p1 == p2) mergedBlock->setPass(p1);
88  double s1 = currentBlock_->getScore();
89  double n1 = static_cast<double>(currentBlock_->getNumberOfSites());
90  double s2 = incomingBlock_->getScore();
91  double n2 = static_cast<double>(incomingBlock_->getNumberOfSites());
92  mergedBlock->setScore((s1 * n1 + s2 * n2) / (n1 + n2));
93 
94  //Now fill the new block:
95  for (size_t i = 0; i < allSp.size(); ++i) {
96  unique_ptr<MafSequence> seq;
97  try {
98  seq.reset(new MafSequence(currentBlock_->getSequenceForSpecies(allSp[i])));
99 
100  //Check is there is a second sequence:
101  try {
102  unique_ptr<MafSequence> tmp(new MafSequence(incomingBlock_->getSequenceForSpecies(allSp[i])));
103  string ref1 = seq->getDescription(), ref2 = tmp->getDescription();
104  if (seq->getChromosome() != tmp->getChromosome()) {
105  seq->setChromosome("fus");
106  seq->removeCoordinates();
107  }
108  if (seq->getStrand() != tmp->getStrand()) {
109  seq->setStrand('?');
110  seq->removeCoordinates();
111  }
112  if (seq->getName() != tmp->getName())
113  tmp->setName(seq->getName()); //force name conversion to prevent exception in 'merge'.
114  seq->merge(*tmp);
115  if (logstream_) {
116  (*logstream_ << "BLOCK CONCATENATE: merging " << ref1 << " with " << ref2 << " into " << seq->getDescription()).endLine();
117  }
118  } catch (SequenceNotFoundException& snfe2) {
119  //There was a first sequence, we just extend it:
120  string ref1 = seq->getDescription();
121  seq->setToSizeR(seq->size() + incomingBlock_->getNumberOfSites());
122  if (logstream_) {
123  (*logstream_ << "BLOCK CONCATENATE: extending " << ref1 << " with " << incomingBlock_->getNumberOfSites() << " gaps on the right.").endLine();
124  }
125  }
126  } catch (SequenceNotFoundException& snfe1) {
127  //There must be a second sequence then:
128  seq.reset(new MafSequence(incomingBlock_->getSequenceForSpecies(allSp[i])));
129  string ref2 = seq->getDescription();
130  seq->setToSizeL(seq->size() + currentBlock_->getNumberOfSites());
131  if (logstream_) {
132  (*logstream_ << "BLOCK CONCATENATE: adding " << ref2 << " and extend it with " << currentBlock_->getNumberOfSites() << " gaps on the left.").endLine();
133  }
134  }
135  mergedBlock->addSequence(*seq);
136  }
137  //Cleaning stuff:
138  delete currentBlock_;
139  delete incomingBlock_;
140  currentBlock_ = mergedBlock;
141  //We check if we can also merge the next block:
142  incomingBlock_ = iterator_->nextBlock();
143  }
144  return currentBlock_;
145 }
146 
static void displayMessage(const std::string &text)
static void displayUnlimitedGauge(size_t iter, const std::string &mes="")
A synteny block data structure, the basic unit of a MAF alignement file.
Definition: MafBlock.h:57
void setScore(double score)
Definition: MafBlock.h:102
void setPass(unsigned int pass)
Definition: MafBlock.h:103
void addSequence(const MafSequence &sequence)
Definition: MafBlock.h:115
A sequence class which is used to store data from MAF files.
Definition: MafSequence.h:64
static std::vector< T > vectorUnion(const std::vector< T > &vec1, const std::vector< T > &vec2)
static std::vector< T > unique(const std::vector< T > &v)
std::size_t count(const std::string &s, const std::string &pattern)