Method addInfoFromAPileUp

src/SNPinGenome.cpp:476–549

}

float SNPinGenome::addInfoFromAPileUp (int totalLetterCount, int minimalTotalLetterCountPerPosition,char whatToLook, int index,int &positionCount, int &sNPpositionToProceed, const char * pileup,int minimalQualityPerPosition, const char * quality){
    float value=NA;
    float status = NA;
    string pileupShort = "";
    if (totalLetterCount>=minimalTotalLetterCountPerPosition) {
        //check that there is not deletions/insertions at this position:
        if (strccnt(pileup, int('+'))+strccnt(pileup, int('-'))==0) {


            if (minimalQualityPerPosition>0) {

                pileupShort=pileup;
          //  cout << quality << "\t"<< pileupShort <<"\t";
                string qualityS = quality;
                chomp(qualityS);
                deleteChar(pileupShort, '^', 1);
                strkeepOnly(pileupShort, ".,ACGTNacgtn*");
         //   cout << pileupShort <<"\n";
                if (pileupShort.length()==qualityS.length()) {
                    filterWithQualities(pileupShort,qualityS,minimalQualityPerPosition);
                    deleteChar(pileupShort, '*');
                    pileup = pileupShort.c_str();
                    totalLetterCount = pileupShort.length();
                } else {
                    cout << "position to skip: " << pileupShort << " " << qualityS << "\n";
                    //skip this position..
                    totalLetterCount = -1;
                }
            }

            if (totalLetterCount>=minimalTotalLetterCountPerPosition) {

                //we will recalculate totalLetterCount = refCount+countForOverLetter
                int refCount =strccnt(pileup, '.')+strccnt(pileup, ',');
                int countForOverLetter=strccnt(pileup, whatToLook)+strccnt(pileup, tolower(whatToLook));

                //value = countForOverLetter*1./totalLetterCount;


                if (refCount+countForOverLetter>=minimalTotalLetterCountPerPosition) {
                    if (countForOverLetter+refCount>totalLetterCount) {
                            countForOverLetter=totalLetterCount-refCount;
                            if (countForOverLetter<0)
                                countForOverLetter=0;
                    }
                    value = countForOverLetter*1./(refCount+countForOverLetter);

                    //value = fabs(countForOverLetter*1./totalLetterCount-0.5);
                  //  cout << "pileup: "<< pileup << " ; whatToLook: "<< whatToLook<< "\n";
                   // cout << "countForOverLetter-1: "<< countForOverLetter-1 << "; totalLetterCount: "<< totalLetterCount << "; ERROR_PER_POS: "<< ERROR_PER_POS << "\n";
                    double prob = binomialcdistribution(countForOverLetter-1, totalLetterCount, ERROR_PER_POS); //probability that there > k Beta Allel bases if the SNP is HomoZ.
                    status = 0;
                    if (prob<0.05) { //0.01

                        if (binomialcdistribution(totalLetterCount-countForOverLetter-1, totalLetterCount, ERROR_PER_POS)<0.01)
                            //decide that SNP is not homoZ.
                            status = 0.5; //at this point, to indicate possible AB
                    }