【生物信息学】使用genome 作为ref时，由bam格式 or pileup格式 计算 depth 的 cpp程序

主要是首次测试 map引用指针，构建多重map。
main.cpp

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <map>
#include <set>
#include <cmath>
#include <inttypes.h>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>
#include "depth.h"

using namespace std;

void usage()
{
cout << "./depth_pileup Sample.pileup <genome.fa.fai> <output> " << endl;
cout << " -d depth for count; sep by comma " << endl;
cout << " -h get help information" << endl;
cout << "Usage: "<< endl;
cout << " samtools -O -f <genome.fa> Sample.sort.bam >Sample.pileup "<< endl;
cout << " ./depth_pileup -d 0,4,9 Sample.pileup <genome.fa.fai> <output> " << endl;

exit (0);
}

int main(int argc, char *argv[])
{
int c;
vector<int> Depth;
char *p;
while ( (c=getopt(argc,argv,"d:h")) != -1 )
{
switch(c)
{
case 'd' :
p = optarg;
break;
case 'h' :
usage();break;
default : usage();
}
}
if (argc < 3 ) usage();

stringstream ss;
ss << p;
string s = ss.str();

vector<string> Indepth;
boost::split(Indepth,s,boost::is_any_of(","),boost::token_compress_on);

for ( int i=0; i<Indepth.size(); ++i){
Depth.push_back(boost::lexical_cast<int>(Indepth[i]));
}
sort(Depth.begin(),Depth.end());

string sam_file = argv[optind++];
string fai_file = argv[optind++];

map<string,int> ChrSize;
map<string,int> ChrDepth;
map<string,map<int,int> > Chr_Dep_Cnt;

load_fai(fai_file,ChrSize);

ifstream infile;
infile.open(sam_file.c_str());
if (! infile )
{
cerr << "fail to open input file " << sam_file << endl;
}
string lineStr;
int total_line=0;

while ( getline(infile,lineStr,'\n' ))
{
vector<string> lineVec;
total_line++;
boost::split(lineVec,lineStr,boost::is_any_of("\t"),boost::token_compress_on);
ChrDepth[lineVec[0]] += boost::lexical_cast<int>(lineVec[3]);
for ( int i=0; i<Depth.size(); ++i){
if (Depth[i] <= boost::lexical_cast<int>(lineVec[3]) ){
Chr_Dep_Cnt[lineVec[0]][Depth[i]] += 1;
}
else {
Chr_Dep_Cnt[lineVec[0]][Depth[i]] += 0;
}
}
}
infile.close();

cout << "#Chr\tTotal_length\tCovered-sites\tTotal-covered-bases\tmean-Depth";
for ( int i=0;i<Depth.size();++i){
cout << "\tabove_" << Depth[i] << "x";
}
cout << endl;
map<string,map<int,int> >::iterator multi_tr;
map<int,int>::iterator inter_tr;
for ( multi_tr=Chr_Dep_Cnt.begin(); multi_tr!=Chr_Dep_Cnt.end(); multi_tr++ ){
double mean = boost::lexical_cast<double>(ChrDepth[multi_tr->first]) / boost::lexical_cast<double>(total_line);

cout << multi_tr->first << "\t" << ChrSize[multi_tr->first] << "\t" << total_line << "\t" << ChrDepth[multi_tr->first] << "\t" << mean << "\t";
for ( inter_tr=multi_tr->second.begin(); inter_tr!=multi_tr->second.end();inter_tr++ ){
cout << "\t" << inter_tr->second;
}
cout << endl;
}
}

depth.h

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <algorithm>
#include <map>
#include <set>
#include <cmath>
#include <inttypes.h>
#include <boost/lexical_cast.hpp>
#include <boost/algorithm/string.hpp>

void load_fai(std::string &matrix_file,std::map<std::string,int> &GenSize);
depth.cpp

#include "depth.h"
using namespace std;
void load_fai(string &matrix_file, map<string,int> &GenSize){
ifstream infile;
infile.open(matrix_file.c_str());
if (! infile ) {
cerr << "fail to open input file " << matrix_file << endl;
exit(0);
}
string lineStr;

while (getline(infile,lineStr,'\n')){
if (lineStr[0] == ' ' || lineStr[0] == '\n'){
continue;
}

vector<string> lineVec;
boost::split(lineVec,lineStr, boost::is_any_of(":, \t\n"), boost::token_compress_on);
GenSize[lineVec[0]] = boost::lexical_cast<int>(lineVec[1]);
}
infile.close();
}

Makefile

cc=g++
exe=depth_pileup
obj=main.o depth.o

$(exe):$(obj)
$(cc) -o $(exe) $(obj)
main.o:main.cpp depth.h
$(cc) -c main.cpp
depth.o:depth.cpp depth.h
$(cc) -c depth.cpp
clean:
rm -rf *.o $(exe)