Picking up where we left off

We continue working on our DNA shotgun sequencing project

You can find the most up to date version in the project's solution/ folder

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
 CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
  CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
   CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
    CCCTCATCACAACTGAATGCTTGGGCTGAA

We have a match with an overlap of 5 bases!

• should we stop now?

Computing the overlap

Let's imagine that this is the current estimate of the sequence:

CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGCAGACG

and this is a candidate fragment:

CCCTCATCACAACTGAATGCTTGGGCTGAA

We could compute the overlap by comparing the overlapping region, starting from the smallest overlap:

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                 CCCTCATCACAACTGAATGCTTGGGCTGAA

No! If we carry on, we'll find a larger overlap

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                             CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                            CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                           CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                          CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                         CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                        CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                       CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                      CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                     CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                    CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                   CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                  CCCTCATCACAACTGAATGCTTGGGCTGAA

Computing the overlap

A better strategy is to start with the biggest overlap first, and gradually reduce it

• that way, if we do find a match, we can stop immediately

                             CTGAATGCTTGGGCTGAAAGGGCGCGAGACGTATTCCCCGGTTGC
                 CCCTCATCACAACTGAATGCTTGGGCTGAA

Now as soon as we find an overlap, it is guaranteed to be the largest

Exercises

• Add a function to compute the overlap between the current sequence and a candidate fragment

  • make sure it works for both ends of the string

• Use this function to identify the candidate fragment with the largest overlap with current sequence

• Add a function to merge this candidate fragment with the current sequence, given the computed overlap

• Use these functions to iteratively merge candidates fragments until no overlapping fragments remain

• Check that all unmerged fragments are already contained within the sequence

Computing the overlap

We create a header file overlap.h to hold our function declaration

#pragma once
#include <string>
int compute_overlap (const std::string& sequence, const std::string& fragment);

Computing the overlap

We create a header file overlap.h to hold our function declaration

#pragma once
#include <string>
int compute_overlap (const std::string& sequence, const std::string& fragment);

• we return an integer representing the number of bases in the overlap

Computing the overlap

We create a header file overlap.h to hold our function declaration

#pragma once
#include <string>
int compute_overlap (const std::string& sequence, const std::string& fragment);

• we return an integer representing the number of bases in the overlap
• we chose a suitable name for our function that clearly states its purpose

Computing the overlap

We create a header file overlap.h to hold our function declaration

#pragma once
#include <string>
int compute_overlap (const std::string& sequence, const std::string& fragment);

• we return an integer representing the number of bases in the overlap
• we chose a suitable name for our function that clearly states its purpose
• we pass both the sequence and candidate fragment as const references to strings

Computing the overlap

We create a header file overlap.h to hold our function declaration

#pragma once
#include <string>
int compute_overlap (const std::string& sequence, const std::string& fragment);

• we return an integer representing the number of bases in the overlap
• we chose a suitable name for our function that clearly states its purpose
• we pass both the sequence and candidate fragment as const references to strings
• we need to #include the <string> header since we are using std::string

Computing the overlap

We create a header file overlap.h to hold our function declaration

#pragma once
#include <string>
int compute_overlap (const std::string& sequence, const std::string& fragment);

• we return an integer representing the number of bases in the overlap
• we chose a suitable name for our function that clearly states its purpose
• we pass both the sequence and candidate fragment as const references to strings
• we need to #include the <string> header since we are using std::string
• we use the #pragma once directive to avoid double-inclusion

Computing the overlap

Next, we create the overlap.cpp file to hold our function definition

#include <iostream>
#include <stdexcept>
#include <string>
#include "overlap.h"
int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  // this is where our code will go
}

Computing the overlap

Next, we create the overlap.cpp file to hold our function definition

#include <iostream>
#include <stdexcept>
#include <string>
#include "overlap.h"
int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  // this is where our code will go
}

• we need to replicate the function declaration, making sure it matches exactly

Computing the overlap

Next, we create the overlap.cpp file to hold our function definition

#include <iostream>
#include <stdexcept>
#include <string>
#include "overlap.h"
int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  // this is where our code will go
}

• we need to replicate the function declaration, making sure it matches exactly
• we need to #include the matching header file

Computing the overlap

Next, we create the overlap.cpp file to hold our function definition

#include <iostream>
#include <stdexcept>
#include <string>
#include "overlap.h"
int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  // this is where our code will go
}

• we need to replicate the function declaration, making sure it matches exactly
• we need to #include the matching header file
• we need to #include the system headers for the functionality we will be using

Computing the overlap

Next, we create the overlap.cpp file to hold our function definition

#include <iostream>
#include <stdexcept>
#include <string>
#include "overlap.h"
int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  // this is where our code will go
}

• we need to replicate the function declaration, making sure it matches exactly
• we need to #include the matching header file
• we need to #include the system headers for the functionality we will be using
• now we can worry about what will go in the body of our function

Computing the overlap

Now we fill out the body of our compute_overlap() function

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  ...
  return 0;
}

Computing the overlap

Now we fill out the body of our compute_overlap() function

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  ...
  return 0;
}

• always check any assumptions and expectations upfront – this can save you a lot of frustration later...

Computing the overlap

Now we fill out the body of our compute_overlap() function

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  ...
  return 0;
}

• always check any assumptions and expectations upfront – this can save you a lot of frustration later...
• during the development phase, don't hesitate to display lots of information, it will often alert you to where things might be going wrong.

Computing the overlap

Now we fill out the body of our compute_overlap() function

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  ...
  return 0;
}

• always check any assumptions and expectations upfront – this can save you a lot of frustration later...
• during the development phase, don't hesitate to display lots of information, it will often alert you to where things might be going wrong.
• we need to return an integer, so let's return a value of zero (no overlap) by default

Computing the overlap

Now we fill out the body of our compute_overlap() function

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  ...
  return 0;
}

• always check any assumptions and expectations upfront – this can save you a lot of frustration later...
• during the development phase, don't hesitate to display lots of information, it will often alert you to where things might be going wrong.
• we need to return an integer, so let's return a value of zero (no overlap) by default

Computing the overlap

  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }

Computing the overlap

  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }

We iterate over the range of valid overlaps:

• starting from the largest possible overlap (the size of the candidate fragment)
• reducing the size of the overlap at each iteration
• until the overlap is zero

Computing the overlap

  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }

We iterate over the range of valid overlaps:

• starting from the largest possible overlap (the size of the candidate fragment)
• reducing the size of the overlap at each iteration
• until the overlap is zero

We extract the part of the current sequence that corresponds to that overlap

• to do this, we use the string substr() method

Computing the overlap

  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }

We iterate over the range of valid overlaps:

• starting from the largest possible overlap (the size of the candidate fragment)
• reducing the size of the overlap at each iteration
• until the overlap is zero

We extract the part of the current sequence that corresponds to that overlap

• to do this, we use the string substr() method
• and likewise for the candidate sequence – though it is now the end of the string

Computing the overlap

  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }

We iterate over the range of valid overlaps:

• starting from the largest possible overlap (the size of the candidate fragment)
• reducing the size of the overlap at each iteration
• until the overlap is zero

We extract the part of the current sequence that corresponds to that overlap

• to do this, we use the string substr() method
• and likewise for the candidate sequence – though it is now the end of the string

Finally, we check whether both substrings match

• If so, we can immediately return the size of the current overlap

Computing the overlap

  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }

We iterate over the range of valid overlaps:

• starting from the largest possible overlap (the size of the candidate fragment)
• reducing the size of the overlap at each iteration
• until the overlap is zero

We extract the part of the current sequence that corresponds to that overlap

• to do this, we use the string substr() method
• and likewise for the candidate sequence – though it is now the end of the string

Finally, we check whether both substrings match

• If so, we can immediately return the size of the current overlap

Computing the overlap

Here is our function definition in full:

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }
  return 0;
}

Exercises

• Add a function to compute the overlap between the current sequence and a candidate fragment

  • make sure it works for both ends of the string

• Use this function to identify the candidate fragment with the largest overlap with current sequence

• Add a function to merge this candidate fragment with the current sequence, given the computed overlap

• Use these functions to iteratively merge candidates fragments until no overlapping fragments remain

• Check that all unmerged fragments are already contained within the sequence

Computing the overlap

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }
  return 0;
}

Computing the overlap

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end)
      return overlap;
  }
  return 0;
}

Computing the overlap

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  int largest_overlap = 0;
  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  ...

Computing the overlap

int compute_overlap (const std::string& sequence, const std::string& fragment)
{
  if (fragment.size() > sequence.size())
    throw std::runtime_error ("fragment size larger than current sequence!");
  std::cerr << "sequence is \"" << sequence << "\"\n";
  std::cerr << "trying fragment \"" << fragment << "\"\n";
  int largest_overlap = 0;
  for (int overlap = fragment.size(); overlap > 0; --overlap) {
    const auto seq_start = sequence.substr(0, overlap);
    const auto frag_end = fragment.substr(fragment.size()-overlap);
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  ...

Computing the overlap

    ...
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  for (int overlap = fragment.size(); overlap > largest_overlap; --overlap) {
    const auto seq_end = sequence.substr(sequence.size() - overlap);
    const auto frag_start = fragment.substr(0, overlap);
    if (seq_end == frag_start) {
      largest_overlap = -overlap;
      break;
    }
  }
  return largest_overlap;
}

Computing the overlap

    ...
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  for (int overlap = fragment.size(); overlap > largest_overlap; --overlap) {
    const auto seq_end = sequence.substr(sequence.size() - overlap);
    const auto frag_start = fragment.substr(0, overlap);
    if (seq_end == frag_start) {
      largest_overlap = -overlap;
      break;
    }
  }
  return largest_overlap;
}

Computing the overlap

    ...
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  for (int overlap = fragment.size(); overlap > largest_overlap; --overlap) {
    const auto seq_end = sequence.substr(sequence.size() - overlap);
    const auto frag_start = fragment.substr(0, overlap);
    if (seq_end == frag_start) {
      largest_overlap = -overlap;
      break;
    }
  }
  return largest_overlap;
}

Computing the overlap

    ...
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  for (int overlap = fragment.size(); overlap > largest_overlap; --overlap) {
    const auto seq_end = sequence.substr(sequence.size() - overlap);
    const auto frag_start = fragment.substr(0, overlap);
    if (seq_end == frag_start) {
      largest_overlap = -overlap;
      break;
    }
  }
  return largest_overlap;
}

Computing the overlap

    ...
    if (seq_start == frag_end) {
      largest_overlap = overlap;
      break;
    }
  }
  for (int overlap = fragment.size(); overlap > largest_overlap; --overlap) {
    const auto seq_end = sequence.substr(sequence.size() - overlap);
    const auto frag_start = fragment.substr(0, overlap);
    if (seq_end == frag_start) {
      largest_overlap = -overlap;
      break;
    }
  }
  return largest_overlap;
}

Exercises

• Add a function to compute the overlap between the current sequence and a candidate fragment

  • make sure it works for both ends of the string

• Use this function to identify the candidate fragment with the largest overlap with current sequence

• Add a function to merge this candidate fragment with the current sequence, given the computed overlap

• Use these functions to iteratively merge candidates fragments until no overlapping fragments remain

• Check that all unmerged fragments are already contained within the sequence

Identifying fragment with the largest overlap

  ...
  int biggest_overlap = 0;
  int fragment_with_biggest_overlap = -1;
  for (int n = 0; n < fragments.size(); ++n) {
    const auto overlap = compute_overlap (sequence, fragments[n]);
    if (std::abs (biggest_overlap) < std::abs (overlap)) {
      biggest_overlap = overlap;
      fragment_with_biggest_overlap = n;
    }
  }
  if (fragment_with_biggest_overlap >= 0)
    std::cerr << "fragment with biggest overlap is at index " 
      << fragment_with_biggest_overlap
      << ", overlap = " << biggest_overlap << "\n";
  write_sequence (args[2], sequence);
}

Identifying fragment with the largest overlap

  ...
  int biggest_overlap = 0;
  int fragment_with_biggest_overlap = -1;
  for (int n = 0; n < fragments.size(); ++n) {
    const auto overlap = compute_overlap (sequence, fragments[n]);
    if (std::abs (biggest_overlap) < std::abs (overlap)) {
      biggest_overlap = overlap;
      fragment_with_biggest_overlap = n;
    }
  }
  if (fragment_with_biggest_overlap >= 0)
    std::cerr << "fragment with biggest overlap is at index " 
      << fragment_with_biggest_overlap
      << ", overlap = " << biggest_overlap << "\n";
  write_sequence (args[2], sequence);
}

Identifying fragment with the largest overlap

  ...
  int biggest_overlap = 0;
  int fragment_with_biggest_overlap = -1;
  for (int n = 0; n < fragments.size(); ++n) {
    const auto overlap = compute_overlap (sequence, fragments[n]);
    if (std::abs (biggest_overlap) < std::abs (overlap)) {
      biggest_overlap = overlap;
      fragment_with_biggest_overlap = n;
    }
  }
  if (fragment_with_biggest_overlap >= 0)
    std::cerr << "fragment with biggest overlap is at index " 
      << fragment_with_biggest_overlap
      << ", overlap = " << biggest_overlap << "\n";
  write_sequence (args[2], sequence);
}

Identifying fragment with the largest overlap

  ...
  int biggest_overlap = 0;
  int fragment_with_biggest_overlap = -1;
  for (int n = 0; n < fragments.size(); ++n) {
    const auto overlap = compute_overlap (sequence, fragments[n]);
    if (std::abs (biggest_overlap) < std::abs (overlap)) {
      biggest_overlap = overlap;
      fragment_with_biggest_overlap = n;
    }
  }
  if (fragment_with_biggest_overlap >= 0)
    std::cerr << "fragment with biggest overlap is at index " 
      << fragment_with_biggest_overlap
      << ", overlap = " << biggest_overlap << "\n";
  write_sequence (args[2], sequence);
}

Identifying fragment with the largest overlap

  ...
  int biggest_overlap = 0;
  int fragment_with_biggest_overlap = -1;
  for (int n = 0; n < fragments.size(); ++n) {
    const auto overlap = compute_overlap (sequence, fragments[n]);
    if (std::abs (biggest_overlap) < std::abs (overlap)) {
      biggest_overlap = overlap;
      fragment_with_biggest_overlap = n;
    }
  }
  if (fragment_with_biggest_overlap >= 0)
    std::cerr << "fragment with biggest overlap is at index " 
      << fragment_with_biggest_overlap
      << ", overlap = " << biggest_overlap << "\n";
  write_sequence (args[2], sequence);
}

Exercises

• Add a function to compute the overlap between the current sequence and a candidate fragment

  • make sure it works for both ends of the string

• Use this function to identify the candidate fragment with the largest overlap with current sequence

• Add a function to merge this candidate fragment with the current sequence, given the computed overlap

• Use these functions to iteratively merge candidates fragments until no overlapping fragments remain

• Check that all unmerged fragments are already contained within the sequence

Merge fragment with current sequence

Once we've identified the candidate fragment with the largest overlap, we need to merge it with the current sequence to produce a longer, updated sequence:

Merge fragment with current sequence

Next, add our function definition to overlap.cpp:

void merge (std::string& sequence, const std::string& fragment, const int overlap)
{
  if (overlap < 0) {
    sequence += fragment.substr (-overlap);
  }
  else if (overlap > 0) {
    sequence = fragment.substr (0, fragment.size()-overlap) + sequence;
  }
}

Merge fragment with current sequence

Use merge() function in main():

  ...
  if (fragment_with_biggest_overlap >= 0) {
    std::cerr << "fragment with biggest overlap is at index " 
      << fragment_with_biggest_overlap
      << ", overlap = " << biggest_overlap << "\n";
    merge (sequence, fragments[fragment_with_biggest_overlap], biggest_overlap);
  }
  std::cerr << "final sequence has length " << sequence.size() << "\n";
  write_sequence (args[2], sequence);
}

Exercises

• Add a function to compute the overlap between the current sequence and a candidate fragment

  • make sure it works for both ends of the string

• Use this function to identify the candidate fragment with the largest overlap with current sequence

• Add a function to merge this candidate fragment with the current sequence, given the computed overlap

• Use these functions to iteratively merge candidates fragments until no overlapping fragments remain

• Check that all unmerged fragments are already contained within the sequence

Iteratively merge fragments into sequence

  while (fragments.size()) {
    std::cerr << "---------------------------------------------------\n";
    std::cerr << fragments.size() << " fragments left\n";
    int biggest_overlap = 0;
    int fragment_with_biggest_overlap = -1;
    for (int n = 0; n < static_cast<int> (fragments.size()); ++n) {
      ...
    }
     if (fragment_with_biggest_overlap < 0)
       break;
     if (std::abs (biggest_overlap) < 10)
       break;
    std::cerr << "fragment with biggest overlap is at index " 
              << fragment_with_biggest_overlap
              << ", overlap = " << biggest_overlap << "\n";
    merge (sequence, fragments[fragment_with_biggest_overlap], biggest_overlap);
    fragments.erase (fragments.begin() + fragment_with_biggest_overlap);
  }

Exercises

• Add a function to compute the overlap between the current sequence and a candidate fragment

  • make sure it works for both ends of the string

• Use this function to identify the candidate fragment with the largest overlap with current sequence

• Add a function to merge this candidate fragment with the current sequence, given the computed overlap

• Use these functions to iteratively merge candidates fragments until no overlapping fragments remain

• Check that all unmerged fragments are already contained within the sequence

Check all unmerged fragments are already in sequence

    ...
    fragments.erase (fragments.begin() + fragment_with_biggest_overlap);
  }
  int num_unmatched = 0;
  for (const auto& frag : fragments) {
    if (sequence.find (frag) == std::string::npos)
      ++num_unmatched;
  }
  if (num_unmatched)
    std::cerr << "WARNING: " << num_unmatched << " fragments remain unmatched!\n";
  std::cerr << "final sequence has length " << sequence.size() << "\n";
  write_sequence (args[2], sequence);
}

Command-line options

We've already seen examples of command-line options in other commands:

$ ls -l
$ g++ -std=c++20 shotgun.cpp -o shotgun

Command-line options

There are many ways to handle command-line options

• indeed, there are many libraries dedicated to command-line parsing...

Detecting the -v option

Possible solution:

...
void run (std::vector<std::string>& args)
{
  bool verbose = std::erase (args, "-v");
  if (args.size() < 3)
    throw std::runtime_error ("expected 2 arguments: input_fragments output_sequence");
  ...
}

Global variables

One problem with the previous solution is that verbose is a local variable

• it is only accessible within the run() function

Global variables

A global variable is declared outside of any function and is accessible from anywhere in the code

• provided its declaration is encountered before the point of use within each translation unit

Global variables

We could declare our global variable like this:

...
bool verbose = false;
...
void run (std::vector<std::string>& args)
{
  verbose = std::erase (args, "-v");
  ...

But there are still issues with this:

• this is declared within our shotgun.cpp file – it won't be accessible in our other cpp files
  • they won't be in the same translation unit

Variable scope

The scope of a variable determines its lifetime and where it can be accessed from.

Global scope

• Variables are instantiated (created) before main() starts, and destroyed after main() returns
• they are accessible from any part of the code
  • provided their declaration occurs before their point of use in the current translation unit

Variable shadowing

Variable shadowing (also called name hiding) occurs when two variables of the same name co-exist in different scopes

int count = 0;   // <= global scope
...
void compute ()
{
  int count = 10;   // <= function scope
  ...
  for (int count = 0; count < 20; ++count) {    // <= block scope
    std::cout << count << "\n";
    ...
  }
  std::cout << count << "\n";
}

Variable shadowing

Variable shadowing (also called name hiding) occurs when two variables of the same name co-exist in different scopes

int count = 0;   // <= global scope
...
void compute ()
{
  int count = 10;   // <= function scope
  ...
  for (int count = 0; count < 20; ++count) {    // <= block scope
    std::cout << count << "\n";
    ...
  }
  std::cout << count << "\n";
}

This is perfectly legal in C++!

Variable shadowing

Variable shadowing (also called name hiding) occurs when two variables of the same name co-exist in different scopes

int count = 0;   // <= global scope
...
void compute ()
{
  int count = 10;   // <= function scope
  ...
  for (int count = 0; count < 20; ++count) {    // <= block scope
    std::cout << count << "\n";
    ...
  }
  std::cout << count << "\n";
}

This is perfectly legal in C++!

Variable shadowing

Variable shadowing (also called name hiding) occurs when two variables of the same name co-exist in different scopes

int count = 0;   // <= global scope
...
void compute ()
{
  int count = 10;   // <= function scope
  ...
  for (int count = 0; count < 20; ++count) {    // <= block scope
    std::cout << count << "\n";
    ...
  }
  std::cout << count << "\n";
}

This is perfectly legal in C++!

Variable shadowing

Variable shadowing (also called name hiding) occurs when two variables of the same name co-exist in different scopes

int count = 0;   // <= global scope
...
void compute ()
{
  int count = 10;   // <= function scope
  ...
  for (int count = 0; count < 20; ++count) {    // <= block scope
    std::cout << count << "\n";
    ...
  }
  std::cout << count << "\n";
}

This is perfectly legal in C++!

Variable shadowing

Variable shadowing allows a function to use local variables without worrying about whether some other global variable exists elsewhere with the same name

What is a namespace?

A namespace is essentially a named scope for your function and variable declarations

What is a namespace?

A namespace is essentially a named scope for your function and variable declarations

We have been using a namespace from the start: the std namespace

• this namespace is where all the functionality provided by the C++ standard is declared

std::vector

std::cout

• these are both declared within the std namespace

What is a namespace?

A namespace is essentially a named scope for your function and variable declarations

We have been using a namespace from the start: the std namespace

• this namespace is where all the functionality provided by the C++ standard is declared

std::vector

std::cout

• these are both declared within the std namespace
• accessing members of a namespace is done with the scope resolution operator (::)

Why would we want to use namespaces?

The main purpose of C++ namespace is:

• to organise your code into logical units
• to avoid name collisions

Declaring namespaces

Declaring variables or functions within a namespace is simply a matter of declaring them within the scope of our namespace.

Working with namespaces

Any code also declared within our debug namespace can already access our verbose variable directly.

Working with namespaces

We can also bring all identifiers declared in a namespace into the current scope with the using namespace directive:

using namespace debug;
...
if (verbose)
  ...

The problem with using namespace std;

In spite of the general recommendation to avoid blanket using namespace directives, you will find that many (if not most) C++ tutorials make use of this directive from the very beginning, starting at 'hello world':

#include <iostream>
using namespace std;
int main() {
  cout << "Hello World!\n";
  return 0;
}

The problem with using namespace std;

So why do we not follow this convention on this course?

• because its use is considered bad practice

Working with namespaces

For the reasons outlined in the previous slides, we recommend to always use the fully qualified name to access members of a namespace

Working with global variables

As previously mentioned, our debug::verbose variable must be declared prior to its point of use in any translation unit where it is needed

• this means it needs to be declared in a header file

Working with global variables

When trying to do this, you will most likely have encountered errors of this nature:

/usr/bin/ld: fragments.o:(.bss+0x0): multiple definition of `debug::verbose'; 
    shotgun.o:(.bss+0x0): first defined here
/usr/bin/ld: overlap.o:(.bss+0x0): multiple definition of `debug::verbose'; 
    shotgun.o:(.bss+0x0): first defined here
collect2: error: ld returned 1 exit status

Working with global variables

Second approach (much simpler):

Use the inline keyword in our debug.h header file to allow multiple definitions of the same variable to be provided across multiple translation units:

#pragma once
namespace debug {
  inline bool verbose = false;
}

We do not need an additional cpp file in this case, and the initial value can now be specified at the same time

Enabling the verbose option

We now have all the pieces required to implement our verbose option

• we add the debug.h header file as per the previous slide

Using the verbose option

We can now run our code as normal:

$ ./shotgun ../data/fragments-1.txt out
initial sequence has size 1000
final sequence has length 20108

which provides a clean output, showing just as much information as we might expect when everything works correctly

Using the verbose option

... or we can run it with the verbose option to see exactly what is going on:

$ ./shotgun ../data/fragments-1.txt out -v
reading fragments from file "../data/fragments-1.txt"...
read 190 fragments
190 fragments, mean fragment length: 529.158, range [ 51 1000 ]
initial sequence has size 1000
---------------------------------------------------
189 fragments left
fragment with biggest overlap is at index 41, overlap = -824
---------------------------------------------------
...
fragment with biggest overlap is at index 37, overlap = 51
---------------------------------------------------
104 fragments left
104 fragments remaining unmatched - checking whether already contained in sequence...
final sequence has length 20108
writing sequence to file "out"...

Using the verbose option

The availability of a verbose option allows you and your users to quickly narrow down problems

• the default output can be kept clean and minimal

• when issues occur, running in verbose mode provide useful information:

  • in the case of crashes, it provides a more precise indication of when the problem occurred (as indicated by the last message), and hence where to look in the code to figure out the problem
  • in the case of unexpected output, the information reported may point out issues with the input data, or provide clues as to where to look for errors in the code

• when there are performance issues, verbose mode can provide an indication of which stage is taking longer to complete than expected

• ...

Adding convenience functions for debugging

Now that we have the option to provide debugging output, we can wrap up this functionality more cleanly to avoid having lots of if (debug::verbose) statements throughout our code

Let's implement a debug::log() function, which will take a std::string as its argument, and print it to the standard error stream if verbose mode is enabled

• it can additionally format its output to more clearly label it as debugging information (for example, by prefixing each line with [DEBUG])

Adding convenience functions for debugging

In debug.h:

#pragma once
include <iostream>
include <string>
namespace debug {
  inline bool verbose = false;
  inline void log (const std::string& message) {
    if (verbose)
      std::cerr << "[DEBUG] " << message << "\n";
  }
}

Adding convenience functions for debugging

In debug.h:

#pragma once
include <iostream>
include <string>
namespace debug {
  inline bool verbose = false;
  inline void log (const std::string& message) {
    if (verbose)
      std::cerr << "[DEBUG] " << message << "\n";
  }
}

Note the use of the inline keyword here, this time to allow multiple definitions of the same function

A word about inline

Historically, the inline keyword was considered more as a performance optimisation feature

• it was used for functions only, to indicate that a given function was suitable for inlining
• this would suggest to the compiler that the required code could be substituted directly at the point of use, rather than calling that function explicitly
• There is always a small overhead to calling a function, and for functions that perform small operations and are called very often, this can impact performance
• declaring such functions as inline can help avoid the overhead of dedicated function calls, and so improve performance
  • this is only possible if the function definition is available, which is why the inline mechanism was introduced to allow definitions to be included in header files without causing the multiple definition problem
• For further details, refer to other online sources (e.g. GeeksForGeeks)

Adding convenience functions for debugging

We can now use our convenience function in our code, for example in fragments.cpp:

...
std::vector<std::string> load_fragments (const std::string& filename)
{
  if (debug::verbose)
    std::cerr << "reading fragments from file \"" << filename << "\"...\n";
  ...

becomes simply:

...
std::vector<std::string> load_fragments (const std::string& filename)
{
  debug::log ("reading fragments from file \"" + filename + "\"...");
  ...

Adding convenience functions for debugging

How do we handle this case?

  if (debug::verbose)
    std::cerr << "read " << fragments.size() << " fragments\n";

Exercise 1