5
\$\begingroup\$

This program takes a simple nucleotide sequence and finds the most common "k-mers" in the sequence, as determined by the supplied dataset (see below). The goal of the program is to find the origin of replication of a DNA sequence, as it is presumably the most common sequence, or in other words, the substring with the most occurrences.

#include <iostream>
#include <string>
#include <fstream>
#include <stdlib.h>
#include <map>

//-------------------------------------------
//FUNCTION PROTYPES
//Dataset files include k in them
//void chooseK();

//Open file and read sequence
void enterFileToOpen();
void openFile();
void readFile();

//Generate substrings and save them to file
void enterOFileName();
void openOFile();
void printSubstrToFile();

//Count substrings and determine modal substrings
//Step 1: Open substring output file
//Step 2: Read and aggregate substrings
void setOutPutReadFiletoOutputFile();
void openOutputFile();
void countSubStrings();


//---------------------------------
//VARIABLES
std::string sequence;
int k;

int kmerCount;
int i=0;
std::string subStrPrint;

std::ifstream inFile;
std::string fileName;

std::ofstream outFile;
std::string ofileName;

std::ifstream outputReadFile;
std::string outputReadFileName;

std::map<std::string, int> subSequenceCountMap;
std::string stringMapElement;

//------------------------------------------

int main()
{
    enterFileToOpen();
    openFile();
    readFile();

    enterOFileName();
    openOFile();
    printSubstrToFile();

    setOutPutReadFiletoOutputFile();
    openOutputFile();
    countSubStrings();

    inFile.close();
    outFile.close();

    return 0;
}

void enterFileToOpen()
{
    std::cout << "Please enter file to open: (remember to include full filepath if file is not located within this program's directory)" << std::endl;
    std::cin >> fileName;
    std::cout << "File: " << fileName << std::endl;
}


void openFile()
{
    inFile.open(fileName.c_str());

    //Check for Error during file open process
    if(inFile.fail())
    {
        std::cerr << "Error during Open File Process" << std::endl;
        exit(1);
    }
}


void readFile()
{
    inFile >> sequence >> k;

    std::cout << "\nNucleotide sequence: " << sequence << "\n" << std::endl;
    std::cout << "k = " << k << std::endl;

}

void enterOFileName()
{
    std::cout << "Please enter name for nucleotide subsequence output file: ";
    std::cin >> ofileName;
    std::cout << "Nucleotide subsequence output file: " << ofileName << std::endl;
}

void openOFile()
{
    outFile.open(ofileName.c_str());
}

void printSubstrToFile()
{
    //variables declared globally: kmerCount, i, subStrPrint
    kmerCount = sequence.length()-k;

    i = 0;

    for(i=0; i<kmerCount; i++)
    {
        subStrPrint = sequence.substr(i,k);
        outFile << subStrPrint << std::endl;
        std::cout << "Iteration " << i+1 << ": sub-sequence " << subStrPrint << std::endl;
    }

    outFile.close();
    inFile.close();
}

void setOutPutReadFiletoOutputFile()
{
    //No need to specify file path for read file, since we just created it
    outputReadFileName = ofileName;
}

void openOutputFile()
{
    inFile.open(outputReadFileName.c_str());

    //Check for Errors
    if(inFile.fail())
    {
        std::cerr << "Error during File Open" << std::endl;
        exit(1);
    }
    else
    {
        std::cout << "\nFile opened successfully\n" << std::endl;
    }
}

void countSubStrings()
{
    i = 0;

    for(i = 0; i < kmerCount; i++)
    {
        inFile >> stringMapElement;

        if(subSequenceCountMap.find(stringMapElement)!=subSequenceCountMap.end())
        {
            subSequenceCountMap[stringMapElement]++;
        }
        else
        {
            subSequenceCountMap[stringMapElement]=1;
        }
    }

    for (std::map<std::string, int>::iterator itr = subSequenceCountMap.begin(); itr!=subSequenceCountMap.end(); ++itr) //changed const_iterator to iterator, cbegin to begin, cend to end
    {
        if(itr->second >1)
        {
            std::cout << itr->first << " occurs " << itr->second << " times" << std::endl;
        }
        else
        {
            std::cout << itr->first << " occurs " << itr->second << " time" << std::endl;
        }
    }

}

Sample dataset:

TGAATAGGAATCCCAACGTCTTGTGAATAGGAATCCCAAGAATCCCAAGAATCCCAAGAATCCCAAGAATCCCAACGTCTTGTGCAAGGTGCAAGGGAATCCCAAGAATCCCAACGTCTTGTGCAAGGGAATCCCAACGTCTTGTGGCCTAGAATCCCAACGTCTTGGAATCCCAACGTCTTGGAATCCCAATGCAAGGGAATCCCAACGTCTTGTGCAAGGTGCAAGGTGAATAGCGTCTTGGAATCCCAATGCAAGGTGGCCTATGCAAGGTGAATAGTGAATAGTGGCCTATGGCCTATGCAAGGTGCAAGGTGAATAGTGCAAGGGAATCCCAAGAATCCCAATGAATAGGAATCCCAAGAATCCCAATGGCCTATGGCCTACGTCTTGTGCAAGGTGGCCTACGTCTTGGAATCCCAATGCAAGGGAATCCCAATGAATAGTGCAAGGGAATCCCAAGAATCCCAATGCAAGGTGGCCTATGCAAGGCGTCTTGCGTCTTGTGGCCTATGCAAGGGAATCCCAATGCAAGGGAATCCCAATGCAAGGTGCAAGGTGCAAGGGAATCCCAACGTCTTGGAATCCCAACGTCTTGCGTCTTGGAATCCCAAGAATCCCAATGGCCTACGTCTTGGAATCCCAATGCAAGGTGAATAGTGAATAGGAATCCCAATGGCCTACGTCTTGTGGCCTATGAATAGTGAATAGTGGCCTAGAATCCCAATGAATAGGAATCCCAATGGCCTACGTCTTGTGCAAGGTGCAAGGTGGCCTACGTCTTGCGTCTTGGAATCCCAAGAATCCCAATGCAAGGTGGCCTATGGCCTATGAATAGTGCAAGGGAATCCCAACGTCTTGCGTCTTGCGTCTTGGAATCCCAA
12

When this file is entered into the program, the program will generate every possible substring of length 12 and calculate the occurrence of each distinct one.

Sample data set output:

AACGTCTTGCGT occurs 2 times
AACGTCTTGGAA occurs 3 times
AACGTCTTGTGA occurs 1 time
AACGTCTTGTGC occurs 3 times
AACGTCTTGTGG occurs 1 time
AAGAATCCCAAC occurs 2 times
AAGAATCCCAAG occurs 3 times
AAGAATCCCAAT occurs 5 times
AAGGCGTCTTGC occurs 1 time
AAGGGAATCCCA occurs 10 times
AAGGTGAATAGC occurs 1 time
AAGGTGAATAGT occurs 3 times
AAGGTGCAAGGG occurs 2 times
AAGGTGCAAGGT occurs 4 times
AAGGTGGCCTAC occurs 2 times
AAGGTGGCCTAT occurs 3 times
AATAGCGTCTTG occurs 1 time
AATAGGAATCCC occurs 5 times
AATAGTGAATAG occurs 3 times
AATAGTGCAAGG occurs 3 times
AATAGTGGCCTA occurs 2 times
AATCCCAACGTC occurs 10 times
AATCCCAAGAAT occurs 10 times
AATCCCAATGAA occurs 3 times
AATCCCAATGCA occurs 8 times
AATCCCAATGGC occurs 4 times
AATGAATAGGAA occurs 2 times
AATGAATAGTGC occurs 1 time
AATGCAAGGGAA occurs 3 times
AATGCAAGGTGA occurs 1 time
AATGCAAGGTGC occurs 1 time
AATGCAAGGTGG occurs 3 times
AATGGCCTACGT occurs 3 times
AATGGCCTATGG occurs 1 time
ACGTCTTGCGTC occurs 3 times
ACGTCTTGGAAT occurs 5 times
ACGTCTTGTGAA occurs 1 time
ACGTCTTGTGCA occurs 5 times
ACGTCTTGTGGC occurs 2 times
AGAATCCCAACG occurs 3 times
AGAATCCCAAGA occurs 3 times
AGAATCCCAATG occurs 6 times
AGCGTCTTGGAA occurs 1 time
AGGAATCCCAAC occurs 1 time
AGGAATCCCAAG occurs 2 times
AGGAATCCCAAT occurs 2 times
AGGCGTCTTGCG occurs 1 time
AGGGAATCCCAA occurs 10 times
AGGTGAATAGCG occurs 1 time
AGGTGAATAGTG occurs 3 times
AGGTGCAAGGGA occurs 2 times
AGGTGCAAGGTG occurs 4 times
AGGTGGCCTACG occurs 2 times
AGGTGGCCTATG occurs 3 times
AGTGAATAGGAA occurs 1 time
AGTGAATAGTGG occurs 2 times
AGTGCAAGGGAA occurs 3 times
AGTGGCCTAGAA occurs 1 time
AGTGGCCTATGG occurs 1 time
ATAGCGTCTTGG occurs 1 time
ATAGGAATCCCA occurs 5 times
ATAGTGAATAGG occurs 1 time
ATAGTGAATAGT occurs 2 times
ATAGTGCAAGGG occurs 3 times
ATAGTGGCCTAG occurs 1 time
ATAGTGGCCTAT occurs 1 time
ATCCCAACGTCT occurs 10 times
ATCCCAAGAATC occurs 10 times
ATCCCAATGAAT occurs 3 times
ATCCCAATGCAA occurs 8 times
ATCCCAATGGCC occurs 4 times
ATGAATAGGAAT occurs 2 times
ATGAATAGTGAA occurs 1 time
ATGAATAGTGCA occurs 2 times
ATGCAAGGCGTC occurs 1 time
ATGCAAGGGAAT occurs 4 times
ATGCAAGGTGAA occurs 2 times
ATGCAAGGTGCA occurs 2 times
ATGCAAGGTGGC occurs 3 times
ATGGCCTACGTC occurs 4 times
ATGGCCTATGAA occurs 1 time
ATGGCCTATGCA occurs 1 time
ATGGCCTATGGC occurs 1 time
CAACGTCTTGCG occurs 2 times
CAACGTCTTGGA occurs 3 times
CAACGTCTTGTG occurs 5 times
CAAGAATCCCAA occurs 10 times
CAAGGCGTCTTG occurs 1 time
CAAGGGAATCCC occurs 10 times
CAAGGTGAATAG occurs 4 times
CAAGGTGCAAGG occurs 6 times
CAAGGTGGCCTA occurs 5 times
CAATGAATAGGA occurs 2 times
CAATGAATAGTG occurs 1 time
CAATGCAAGGGA occurs 3 times
CAATGCAAGGTG occurs 5 times
CAATGGCCTACG occurs 3 times
CAATGGCCTATG occurs 1 time
CCAACGTCTTGC occurs 2 times
CCAACGTCTTGG occurs 3 times
CCAACGTCTTGT occurs 5 times
CCAAGAATCCCA occurs 10 times
CCAATGAATAGG occurs 2 times
CCAATGAATAGT occurs 1 time
CCAATGCAAGGG occurs 3 times
CCAATGCAAGGT occurs 5 times
CCAATGGCCTAC occurs 3 times
CCAATGGCCTAT occurs 1 time
CCCAACGTCTTG occurs 10 times
CCCAAGAATCCC occurs 10 times
CCCAATGAATAG occurs 3 times
CCCAATGCAAGG occurs 8 times
CCCAATGGCCTA occurs 4 times
CCTACGTCTTGC occurs 1 time
CCTACGTCTTGG occurs 2 times
CCTACGTCTTGT occurs 3 times
CCTAGAATCCCA occurs 2 times
CCTATGAATAGT occurs 2 times
CCTATGCAAGGC occurs 1 time
CCTATGCAAGGG occurs 1 time
CCTATGCAAGGT occurs 2 times
CCTATGGCCTAC occurs 1 time
CCTATGGCCTAT occurs 2 times
CGTCTTGCGTCT occurs 5 times
CGTCTTGGAATC occurs 9 times
CGTCTTGTGAAT occurs 1 time
CGTCTTGTGCAA occurs 5 times
CGTCTTGTGGCC occurs 3 times
CTACGTCTTGCG occurs 1 time
CTACGTCTTGGA occurs 2 times
CTACGTCTTGTG occurs 3 times
CTAGAATCCCAA occurs 2 times
CTATGAATAGTG occurs 2 times
CTATGCAAGGCG occurs 1 time
CTATGCAAGGGA occurs 1 time
CTATGCAAGGTG occurs 2 times
CTATGGCCTACG occurs 1 time
CTATGGCCTATG occurs 2 times
CTTGCGTCTTGC occurs 1 time
CTTGCGTCTTGG occurs 3 times
CTTGCGTCTTGT occurs 1 time
CTTGGAATCCCA occurs 9 times
CTTGTGAATAGG occurs 1 time
CTTGTGCAAGGG occurs 1 time
CTTGTGCAAGGT occurs 4 times
CTTGTGGCCTAG occurs 1 time
CTTGTGGCCTAT occurs 2 times
GAATAGCGTCTT occurs 1 time
GAATAGGAATCC occurs 5 times
GAATAGTGAATA occurs 3 times
GAATAGTGCAAG occurs 3 times
GAATAGTGGCCT occurs 2 times
GAATCCCAACGT occurs 10 times
GAATCCCAAGAA occurs 10 times
GAATCCCAATGA occurs 3 times
GAATCCCAATGC occurs 8 times
GAATCCCAATGG occurs 4 times
GCAAGGCGTCTT occurs 1 time
GCAAGGGAATCC occurs 10 times
GCAAGGTGAATA occurs 4 times
GCAAGGTGCAAG occurs 6 times
GCAAGGTGGCCT occurs 5 times
GCCTACGTCTTG occurs 6 times
GCCTAGAATCCC occurs 2 times
GCCTATGAATAG occurs 2 times
GCCTATGCAAGG occurs 4 times
GCCTATGGCCTA occurs 3 times
GCGTCTTGCGTC occurs 2 times
GCGTCTTGGAAT occurs 4 times
GCGTCTTGTGGC occurs 1 time
GGAATCCCAACG occurs 7 times
GGAATCCCAAGA occurs 7 times
GGAATCCCAATG occurs 9 times
GGCCTACGTCTT occurs 6 times
GGCCTAGAATCC occurs 2 times
GGCCTATGAATA occurs 2 times
GGCCTATGCAAG occurs 4 times
GGCCTATGGCCT occurs 3 times
GGCGTCTTGCGT occurs 1 time
GGGAATCCCAAC occurs 4 times
GGGAATCCCAAG occurs 3 times
GGGAATCCCAAT occurs 3 times
GGTGAATAGCGT occurs 1 time
GGTGAATAGTGA occurs 2 times
GGTGAATAGTGC occurs 1 time
GGTGCAAGGGAA occurs 2 times
GGTGCAAGGTGA occurs 2 times
GGTGCAAGGTGC occurs 1 time
GGTGCAAGGTGG occurs 1 time
GGTGGCCTACGT occurs 2 times
GGTGGCCTATGC occurs 2 times
GGTGGCCTATGG occurs 1 time
GTCTTGCGTCTT occurs 5 times
GTCTTGGAATCC occurs 9 times
GTCTTGTGAATA occurs 1 time
GTCTTGTGCAAG occurs 5 times
GTCTTGTGGCCT occurs 3 times
GTGAATAGCGTC occurs 1 time
GTGAATAGGAAT occurs 2 times
GTGAATAGTGAA occurs 2 times
GTGAATAGTGCA occurs 1 time
GTGAATAGTGGC occurs 2 times
GTGCAAGGGAAT occurs 6 times
GTGCAAGGTGAA occurs 2 times
GTGCAAGGTGCA occurs 4 times
GTGCAAGGTGGC occurs 2 times
GTGGCCTACGTC occurs 2 times
GTGGCCTAGAAT occurs 2 times
GTGGCCTATGAA occurs 1 time
GTGGCCTATGCA occurs 3 times
GTGGCCTATGGC occurs 2 times
TACGTCTTGCGT occurs 1 time
TACGTCTTGGAA occurs 2 times
TACGTCTTGTGC occurs 2 times
TACGTCTTGTGG occurs 1 time
TAGAATCCCAAC occurs 1 time
TAGAATCCCAAT occurs 1 time
TAGCGTCTTGGA occurs 1 time
TAGGAATCCCAA occurs 5 times
TAGTGAATAGGA occurs 1 time
TAGTGAATAGTG occurs 2 times
TAGTGCAAGGGA occurs 3 times
TAGTGGCCTAGA occurs 1 time
TAGTGGCCTATG occurs 1 time
TATGAATAGTGA occurs 1 time
TATGAATAGTGC occurs 1 time
TATGCAAGGCGT occurs 1 time
TATGCAAGGGAA occurs 1 time
TATGCAAGGTGA occurs 1 time
TATGCAAGGTGC occurs 1 time
TATGGCCTACGT occurs 1 time
TATGGCCTATGA occurs 1 time
TATGGCCTATGC occurs 1 time
TCCCAACGTCTT occurs 10 times
TCCCAAGAATCC occurs 10 times
TCCCAATGAATA occurs 3 times
TCCCAATGCAAG occurs 8 times
TCCCAATGGCCT occurs 4 times
TCTTGCGTCTTG occurs 5 times
TCTTGGAATCCC occurs 9 times
TCTTGTGAATAG occurs 1 time
TCTTGTGCAAGG occurs 5 times
TCTTGTGGCCTA occurs 3 times
TGAATAGCGTCT occurs 1 time
TGAATAGGAATC occurs 5 times
TGAATAGTGAAT occurs 3 times
TGAATAGTGCAA occurs 3 times
TGAATAGTGGCC occurs 2 times
TGCAAGGCGTCT occurs 1 time
TGCAAGGGAATC occurs 10 times
TGCAAGGTGAAT occurs 4 times
TGCAAGGTGCAA occurs 6 times
TGCAAGGTGGCC occurs 5 times
TGCGTCTTGCGT occurs 1 time
TGCGTCTTGGAA occurs 3 times
TGCGTCTTGTGG occurs 1 time
TGGAATCCCAAC occurs 2 times
TGGAATCCCAAG occurs 2 times
TGGAATCCCAAT occurs 4 times
TGGCCTACGTCT occurs 6 times
TGGCCTAGAATC occurs 2 times
TGGCCTATGAAT occurs 2 times
TGGCCTATGCAA occurs 4 times
TGGCCTATGGCC occurs 3 times
TGTGAATAGGAA occurs 1 time
TGTGCAAGGGAA occurs 1 time
TGTGCAAGGTGC occurs 3 times
TGTGCAAGGTGG occurs 1 time
TGTGGCCTAGAA occurs 1 time
TGTGGCCTATGA occurs 1 time
TGTGGCCTATGC occurs 1 time
TTGCGTCTTGCG occurs 1 time
TTGCGTCTTGGA occurs 3 times
TTGCGTCTTGTG occurs 1 time
TTGGAATCCCAA occurs 8 times
TTGTGAATAGGA occurs 1 time
TTGTGCAAGGGA occurs 1 time
TTGTGCAAGGTG occurs 4 times
TTGTGGCCTAGA occurs 1 time
TTGTGGCCTATG occurs 2 times
Improve this question
\$\endgroup\$

1 Answer 1

Reset to default
6
\$\begingroup\$

Code Review

Replace:

    if(subSequenceCountMap.find(stringMapElement)!=subSequenceCountMap.end())
    {
        subSequenceCountMap[stringMapElement]++;
    }
    else
    {
        subSequenceCountMap[stringMapElement]=1;
    }

With:

    subSequenceCountMap[stringMapElement]++;

This is because operator[] will automatically insert the element 0 if it does not already exist before returning a reference. Thus if you just increment the value you get the correct result.

I don't see the point of using an intermediate file to store all the sub strings. Just put them dirextly into the map.

Replace:

for(i=0; i<kmerCount; i++)
{
    subStrPrint = sequence.substr(i,k);
    outFile << subStrPrint << std::endl;
    std::cout << "Iteration " << i+1 << ": sub-sequence " << subStrPrint << std::endl;
}

With:

for(i=0; i<kmerCount; i++)
{
    subSequenceCountMap[sequence.substr(i,k)]++;
}

Not sure why you need this test:

    if(itr->second >1)
    {
        std::cout << itr->first << " occurs " << itr->second << " times" << std::endl;
    }
    else
    {
        std::cout << itr->first << " occurs " << itr->second << " time" << std::endl;
    }

Both sides of the if statement are exactly the same. Remove all the ones that have only happened once by removing the whole else block!

Ohh. I see it now. The extra s on the end of time. Yea. If this were a million dollar project maybe I would go to the extra effort of pluralization. But I would do it simpler.

    std::cout << itr->first << " occurs "
              << itr->second << " time" << ((itr->second > 1)?"s":"")
              << "\n";

Design:

Don't use global variables.

//---------------------------------
//VARIABLES
std::string sequence;
int k;

int kmerCount;
int i=0;
std::string subStrPrint;

std::ifstream inFile;
std::string fileName;

std::ofstream outFile;
std::string ofileName;

std::ifstream outputReadFile;
std::string outputReadFileName;

std::map<std::string, int> subSequenceCountMap;
std::string stringMapElement;

Have local variables and pass them as parameters to functions and get results as a return value from the function.

Loading the whole gene into memory seems like it could be very expensive. Especially since you only need the k characters at a time. Why not load the first k characters (that's your first string). Then move everything down once then read one more character to the end.

Error

You have an off by one error here:

 kmerCount = sequence.length()-k;

Your code reports:

 TTGGAATCCCAA occurrs 8 times

But if you manually count it. You will find that it happens 9 times.

Writting in C++ style

I would have used an object to encapsulate the gene fragment.

#include <iostream>
#include <fstream>
#include <map>
#include <string>

class Gene
{
    std::istream&   in;
    std::string     section;
    public:
        Gene(std::istream& in, std::size_t size)
            : in(in)
        {
            section.resize(size);
            in.read(&section[0], size);
        }

        bool more()
        {
            return in.good();
        }

        void next()
        {
            std::move(std::begin(section) + 1, std::end(section), std::begin(section));
            in.read(&section.back(), 1);
        }
        operator std::string const&()
        {
            return section;
        }
};

int main()
{
    std::string fileName;
    std::size_t size;

    std::cout << "Enter File Name with Gene and size\n";
    if (!(std::cin >> fileName >> size))
    {
        std::cerr << "Failed to read info\n";
        return 1;
    }

    std::ifstream   file(fileName);
    if (!file)
    {
        std::cerr << "Failed to open file\n";
        return 1;
    }

    std::map<std::string, int>  count;
    for(Gene gene(file, size);gene.more();gene.next())
    {
        ++count[gene];
    }

    for(auto const& item: count)
    {
        if (item.second > 1)
        {
            std::cout << item.first << " " << item.second << "\n";
        }
    }
}
Improve this answer
\$\endgroup\$
5
  • \$\begingroup\$ Thank you for taking the time to review my program. I think it's pretty funny that I never stopped using global variables since I started doing it while I was learning Visual Basic in '05-'06. I don't think I ever would have realized it was wrong if you hadn't pointed it out, although I must admit it felt weird doing it. Thank you again, that was truly insightful. You recommended encapsulating the gene fragment using an object, but I must admit I'm a bit at a loss as to how to do that right now, but having your template for the main function for now helped me out a lot. \$\endgroup\$
    – Jose Lopez
    Nov 30, 2016 at 17:55
  • \$\begingroup\$ I especially liked your mention of not loading the entire gene into memory. I hadn't considered the performance repercussions, about which you're absolutely right. The sample dataset I used for this project wasn't even a tenth of the size of a real sequence, so thank you for that. Thank you again for the review, it was remarkably enlightening. Happy holidays, Jose \$\endgroup\$
    – Jose Lopez
    Nov 30, 2016 at 17:55
  • 1
    \$\begingroup\$ @JoseLopez Added the Gene class for you. \$\endgroup\$
    – Martin York
    Nov 30, 2016 at 18:57
  • 1
    \$\begingroup\$ Probably want to write errors to std:cerr, and return 1; instead of exit(1); (exit doesn't guarantee flushing output streams, among other things). \$\endgroup\$
    – Jerry Coffin
    Dec 1, 2016 at 0:29
  • \$\begingroup\$ @JerryCoffin That makes sense, I didn't know that, thank you \$\endgroup\$
    – Jose Lopez
    Dec 2, 2016 at 18:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.