Method calculateFlankLength

src/BAFpileup.cpp:62–156

}

float BAFpileup::calculateFlankLength(std::string const& mateFileName, std::string const& inputFormat_str, std::string const& matesOrientation_str, std::string pathToSamtools_, std::string pathToSambamba, std::string SambambaThreads)
{
        if (matesOrientation_str=="0") return 0; // do not add anything in case of single end data
        if (getInputFormat(inputFormat_str)!=SAM_INPUT_FORMAT)  return 0;
        std::ifstream fileMates (mateFileName.c_str());
        vector<float> insertSizeVector;
        string line;
        if (!fileMates.is_open())        {
            cerr << "Error: unable to open "+mateFileName+"\n" ;
            exit(-1);
        }
        fileMates.close();

        #ifdef PROFILE_TRACE
            time_t t0 = time(NULL);
        #endif

       // MateOrientation matesOrientation = getMateOrientation(matesOrientation_str); // will not consider read orientation to increase speed
        char* line_buffer;
        FILE *stream;
        char buffer[MAX_BUFFER];
        bool zgOrbam=0;

        int numberOfReadsToCheck=3000000;
        int j = 0;
        int fragmentLength=0;
        if(mateFileName.substr(mateFileName.size()-3,3).compare("bam")==0 || mateFileName.substr(mateFileName.size()-3,3).compare(".gz")==0) {
            string command;
            if (mateFileName.substr(mateFileName.size()-3,3).compare("bam")==0) {
                if (pathToSambamba != "")     {
                     command = pathToSambamba + " view -t " + SambambaThreads + " " + mateFileName;
                }    else    {
                     command = pathToSamtools_ + " view "+mateFileName;
                }
            }
            if (mateFileName.substr(mateFileName.size()-3,3).compare(".gz")==0) {
                          command = "gzip -c -d "+mateFileName;
            }
            stream =
            #if defined(_WIN32)
                _popen(command.c_str(), "r");
            #else
                popen(command.c_str(), "r");
            #endif
            zgOrbam=1;

            while (((line_buffer = getLine(buffer, MAX_BUFFER, stream, line)) != NULL) && j<numberOfReadsToCheck) {
                if (line_buffer[0] == '@')
                    continue;

                char* strs[32];
                unsigned int strs_cnt = split((char*)line_buffer, '\t', strs);

                if (strs_cnt > 8 && strs[6][0]=='=') {
                    int frLen=atoi(strs[8]);
                    if(frLen>0 && frLen< 10000) {
                        fragmentLength+=frLen;
                        j++;