Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save gatoravi/454abf932c0b262058ab to your computer and use it in GitHub Desktop.
Save gatoravi/454abf932c0b262058ab to your computer and use it in GitHub Desktop.
Parse 'samtools depth' and obtain lines with sufficient coverage. Assumes two samples in the output. Same as, awk '$3 >= cutoff && $4 >= cutoff'
#include <stdint.h>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>
using namespace std;
void usage() {
cerr << "Usage: ./get_count_with_sufficient_coverage samtools_depth_file cutoff[INT]";
cerr << "Assumes 2 samples.";
cerr << endl;
return;
}
int count_number_satisfying_lines(string depth_file, int cutoff) {
ifstream fin;
fin.open(depth_file.c_str());
string line;
string chr;
uint32_t pos, s1_readcount, s2_readcount;
int lines_satisfying = 0;
while(getline(fin, line)) {
stringstream ss(line);
ss >> chr >> pos >> s1_readcount >> s2_readcount;
if(s1_readcount >= cutoff && s2_readcount >= cutoff) {
lines_satisfying++;
}
}
return lines_satisfying;
}
int main(int argc, char* argv[]) {
if (argc < 3) {
usage();
return 1;
}
string samtools_depth_file = string(argv[1]);
int cutoff = atoi(argv[2]);
cerr << "\nCutoff " << cutoff;
cerr << "\nSamtools file: " << samtools_depth_file;
cerr << "\nNumber of lines satisfying cutoff: ";
cout << count_number_satisfying_lines(samtools_depth_file, cutoff);
cout << endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment