Skip to content

Instantly share code, notes, and snippets.

@gatoravi
Last active September 30, 2019 23:33
Show Gist options
  • Save gatoravi/0e6bb97d2b1003470f2cf61e5f24c595 to your computer and use it in GitHub Desktop.
Save gatoravi/0e6bb97d2b1003470f2cf61e5f24c595 to your computer and use it in GitHub Desktop.
Takes a kmer and fasta as input and lists all positions //within the fasta for that kmer
/*
Copyright 2019 Avinash Ramu
Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/
//Takes a kmer and fasta as input and lists all positions
//within the fasta for that kmer.
//Example usage ./a.out hg38.fa AT
#include <iostream>
#include <fstream>
#include <sstream>
using namespace std;
int main(int argc, char* argv[]) {
if (argc != 3) {
cerr << "Example usage ./a.out hg38.fa AT";
return 1;
}
string kmer = argv[2];
// Read the FASTA
ifstream fasta_fh(argv[1]);
std::string line, contig;
long position;
char prev_line_last_char; // STore the last character of previous line
int line_number = 0;
while(std::getline(fasta_fh, line).good()) {
if (line[0] == '>') {
istringstream iss(line.substr(1));
iss >> contig; // Get the chromosome name
prev_line_last_char = '0';
position = 0;
continue;
}
for (std::string::size_type i = 0; i < line.size(); i++) {
position += 1;
if(i == 0 && prev_line_last_char != '0') { //CHeck concatenation with prev line
string first;
first += prev_line_last_char;
first += line[0];
if(first == kmer) {
cout << contig << " " << position << endl;
}
}
string twomer; //Compare consecutive characters
twomer += line[i];
twomer += line[i+1];
if(twomer == kmer) { // Check if the 2mer is what we want
cout << contig << " " << position << endl;
}
}
prev_line_last_char = line[line.size() - 1];
}
}
@gatoravi
Copy link
Author

To compile: g++ 2mer.cc
To run: ./a.out hs37d5.fa AT
Takes about 7 mins on the human genome, not extensively tested.
The positions are 1-based.

@gatoravi
Copy link
Author

One issue I can think of is some genome FASTAs have lower or upper cases depending on sequence complexity. This program currently only matches strings with the same case. So if the user provides lower-case 2mer only the matching lower-case positions will be printed out.

@arnavm
Copy link

arnavm commented Sep 30, 2019

The other thing it doesn't account for is strand. That's fine for 2-mers or palindromic k-mers, but not a general purpose solution.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment