Merge pull request #1 from mhaseeb123/revision

Changes for revision

.gitignore

@@ -16,7 +16,6 @@
 *.jekyll-metadata
 
 #IDE files
-/CMakeLists.txt.user*
 /.vscode
 /.vs
 /*.nja
@@ -99,7 +98,7 @@ __pycache__/
 # Distribution / packaging
 .Python
 env/
-build/
+build*/
 develop-eggs/
 dist/
 downloads/
@@ -183,3 +182,16 @@ ENV/
 # vscode
 .vscode/
 .vs/
+
+# cmake 
+CMakeLists.txt.user*
+CMakeCache.txt*
+CMakeFiles*
+CMakeScripts*
+Testing*
+Makefile*
+cmake_install.cmake*
+install_manifest.txt*
+compile_commands.json*
+CTestTestfile.cmake*
+_deps*

CMakeLists.txt

@@ -138,6 +138,7 @@ message(STATUS "Configuring...")
 
 option(TAILFIT "Use Tailfit method instead of Gumbelfit for e(x)" ON)
 option(PROGRESS "Enable Progress Marks" ON)
+option(MATCH_CHARGE "Ensure precursor charge - 1 in fragment-ion search" OFF)
 
 # Allowed maximum peptide sequence length
 if (NOT MAX_SEQ_LEN)

config.hpp.in

@@ -36,6 +36,10 @@
 // Print progress marks
 #cmakedefine PROGRESS
 
+// match fragment-ions to only the database ions with charge = precursor charge - 1
+// also do not match to the b1 ions that have low probability of generation.
+#cmakedefine MATCH_CHARGE
+
 // maximum peptide sequence length in the database: max 62 allowed
 #cmakedefine MAX_SEQ_LEN              @MAX_SEQ_LEN@
 

source/app/hicops.cpp

@@ -600,9 +600,9 @@ static status_t ParseParams(char_t* paramfile)
         getline(pfile, line);
         params.min_cpsm = std::atoi(line.c_str());
 
-        /* Base Intensity */
+        /* Base Intensity x 1000 */
         getline(pfile, line);
-        params.base_int = std::atoi(line.c_str());
+        params.base_int = std::atoi(line.c_str()) * 1000;
 
         /* Cutoff intensity ratio */
         getline(pfile, line);

source/core/dslim_fileout.cpp

@@ -65,10 +65,11 @@ status_t DFile_InitFiles()
                 string_t filename = common + "_" + std::to_string(f) + ".tsv";
                 tsvs[f].open(filename);
 
-                tsvs[f] << "file\t" << "scan_num\t" << "prec_mass\t" << "peptide\t"
-                        << "matched_ions\t" << "total_ions\t"
-                        << "calc_pep_mass\t" << "mass_diff\t" << "mod_info\t"
-                        << "hyperscore\t" << "expectscore\t" << "num_hits" << std::endl;
+                tsvs[f] << "file\t" << "scan_num\t" << "prec_mass\t" << "charge\t" 
+                        << "retention_time\t" << "peptide\t" << "matched_ions\t" 
+                        << "total_ions\t" << "calc_pep_mass\t" << "mass_diff\t" 
+                        << "mod_info\t" << "hyperscore\t" << "expectscore\t" 
+                        << "num_hits" << std::endl;
             }
         }
     }
@@ -133,12 +134,14 @@ status_t DFile_PrintScore(Index *index, uint_t specid, float_t pmass, hCell *psm
     tsvs[thno]         << queryfiles[psm->fileIndex];
     tsvs[thno] << '\t' << std::to_string(specid + 1);
     tsvs[thno] << '\t' << std::to_string(pmass);
+    tsvs[thno] << '\t' << std::to_string(psm->pchg);
+    tsvs[thno] << '\t' << std::to_string(psm->rtime);
     tsvs[thno] << '\t' << pep;
     tsvs[thno] << '\t' << std::to_string(psm->sharedions);
     tsvs[thno] << '\t' << std::to_string(psm->totalions);
     tsvs[thno] << '\t' << std::to_string(lclindex->pepEntries[pepid].Mass);
     tsvs[thno] << '\t' << std::to_string((pmass - lclindex->pepEntries[pepid].Mass));
-    tsvs[thno] << '\t'; // TODO: mod_info
+    tsvs[thno] << '\t'; // TODO: print (mod_info) here
     tsvs[thno] << '\t' << std::to_string(psm->hyperscore);
     tsvs[thno] << '\t' << std::to_string(e_x);
     tsvs[thno] << '\t' << std::to_string(npsms);

source/core/dslim_query.cpp

@@ -620,6 +620,8 @@ status_t DSLIM_QuerySpectrum(Queries *ss, Index *index, uint_t idxchunk)
             /* Pointer to each query spectrum */
             uint_t *QAPtr = ss->moz + ss->idx[queries];
             float_t pmass = ss->precurse[queries];
+            auto    pchg  = ss->charges[queries];
+            auto    rtime = ss->rtimes[queries];
             uint_t *iPtr = ss->intensity + ss->idx[queries];
             uint_t qspeclen = ss->idx[queries + 1] - ss->idx[queries];
             uint_t thno = omp_get_thread_num();
@@ -638,6 +640,9 @@ status_t DSLIM_QuerySpectrum(Queries *ss, Index *index, uint_t idxchunk)
             {
                 uint_t speclen = (index[ixx].pepIndex.peplen - 1) * maxz * iSERIES;
                 uint_t halfspeclen = speclen / 2;
+#ifdef MATCH_CHARGE
+                uint_t sixthspeclen = speclen / 6;
+#endif // MATCH_CHARGE
 
                 for (uint_t chno = 0; chno < index[ixx].nChunks; chno++)
                 {
@@ -696,8 +701,12 @@ status_t DSLIM_QuerySpectrum(Queries *ss, Index *index, uint_t idxchunk)
                                     int_t isY = residue / halfspeclen;
                                     int_t isB = 1 - isY;
 
-                                    while (isY < 0 || isY > 1);
+#ifdef MATCH_CHARGE
+                                    int_t ichg = residue / sixthspeclen;
 
+                                    isY *= residue && (ichg < (pchg + params.maxz));
+                                    isB *= residue && (ichg < pchg);
+#endif // MATCH_CHARGE
                                     /* Get the map element */
                                     BYC *elmnt = bycPtr + ppid;
 
@@ -728,13 +737,11 @@ status_t DSLIM_QuerySpectrum(Queries *ss, Index *index, uint_t idxchunk)
                             /* Create a heap cell */
                             hCell cell;
 
-                            ull_t pp = UTILS_Factorial(bcc) *
-                                    UTILS_Factorial(ycc);
+                            double_t pp = log10(UTILS_Factorial(bcc)) + 
+                                    log10(UTILS_Factorial(ycc));
 
                             /* Fill in the information */
-                            cell.hyperscore = 0.001 + pp * bycPtr[it].ibc * bycPtr[it].iyc;
-
-                            cell.hyperscore = log10(cell.hyperscore) - 6;
+                            cell.hyperscore = pp + log10(1 + (ull_t)bycPtr[it].ibc) + log10(1 + (ull_t)bycPtr[it].iyc) - 6;
 
                             /* hyperscore < 0 means either b- or y- ions were not matched */
                             if (cell.hyperscore > 0)
@@ -744,6 +751,8 @@ status_t DSLIM_QuerySpectrum(Queries *ss, Index *index, uint_t idxchunk)
                                 cell.sharedions = shpk;
                                 cell.totalions = speclen;
                                 cell.pmass = pmass;
+                                cell.pchg = pchg;
+                                cell.rtime = rtime;
                                 cell.fileIndex = ss->fileNum;
 
                                 /* Insert the cell in the heap dst */

source/core/include/msquery.h

@@ -31,6 +31,7 @@ typedef struct _Spectrum
     uint_t SpectrumSize;
     double_t prec_mz;
     uint_t Z;
+    double_t rtime;
 
     /* Overload the = operator - Required by MSQuery */
     _Spectrum &operator=(const _Spectrum &rhs)
@@ -40,6 +41,7 @@ typedef struct _Spectrum
         this->SpectrumSize = rhs.SpectrumSize;
         this->prec_mz = rhs.prec_mz;
         this->Z = rhs.Z;
+        this->rtime = rhs.rtime;
 
         return *this;
     }

source/core/include/slm_dsts.h

@@ -425,6 +425,8 @@ typedef struct _queries
     uint_t  *intensity; /* Stores the intensity values of the experimental spectra */
     uint_t        *idx; /* Row ptr. Starting index of each row */
     float_t  *precurse; /* Stores the precursor mass of each spectrum. */
+    int_t     *charges;
+    float_t    *rtimes;
     int_t     numPeaks;
     int_t     numSpecs; /* Number of theoretical spectra */
     int_t     batchNum;
@@ -444,6 +446,7 @@ typedef struct _queries
         this->idx       = NULL;
         this->precurse  = NULL;
         this->moz       = NULL;
+        this->charges   = NULL;
         this->intensity = NULL;
         numPeaks        = 0;
         numSpecs        = 0;
@@ -454,6 +457,8 @@ typedef struct _queries
     {
         this->idx       = new uint_t[QCHUNK + 1];
         this->precurse  = new float_t[QCHUNK];
+        this->charges   = new int_t[QCHUNK];
+        this->rtimes   = new float_t[QCHUNK];
         this->moz       = new uint_t[QCHUNK * QALEN];
         this->intensity = new uint_t[QCHUNK * QALEN];
         fileNum         = 0;
@@ -488,6 +493,18 @@ typedef struct _queries
             this->precurse = NULL;
         }
 
+        if (this->charges != NULL)
+        {
+            delete[] this->charges;
+            this->charges = NULL;
+        }
+
+        if (this->rtimes != NULL)
+        {
+            delete[] this->rtimes;
+            this->rtimes = NULL;
+        }
+
         if (this->idx != NULL)
         {
             delete[] this->idx;
@@ -520,6 +537,18 @@ typedef struct _queries
             this->precurse = NULL;
         }
 
+        if (this->charges != NULL)
+        {
+            delete[] this->charges;
+            this->charges = NULL;
+        }
+
+        if (this->rtimes != NULL)
+        {
+            delete[] this->rtimes;
+            this->rtimes = NULL;
+        }
+
         if (this->idx != NULL)
         {
             delete[] this->idx;
@@ -545,8 +574,9 @@ typedef struct _heapEntry
 
     /* Parent spectrum ID in the respective chunk of index */
     int_t        psid;
-
     float_t     pmass;
+    int_t        pchg;
+    float_t     rtime;
 
     /* Computed hyperscore */
     float_t hyperscore;
@@ -561,6 +591,8 @@ typedef struct _heapEntry
         sharedions = 0;
         totalions  = 0;
         pmass      = 0;
+        pchg       = 0;
+        rtime      = 0;
     }
 
     /* Copy constructor */
@@ -596,6 +628,8 @@ typedef struct _heapEntry
         this->sharedions = rhs;
         this->totalions = rhs;
         this->pmass = rhs;
+        this->pchg = rhs;
+        this->rtime = rhs;
 
         return *this;
     }