Fix HMM indexing errors

WGLab · Sep 16, 2024 · 87b3164 · 87b3164
1 parent 6d5bd3e
commit 87b3164
Show file tree

Hide file tree

Showing 5 changed files with 107 additions and 142 deletions.
diff --git a/__main__.py b/__main__.py
@@ -217,7 +217,7 @@ def main():
         if value is None:
             setattr(args, key, "")
 
-    # Set input parameters.
+    # Set input parameters
     input_data = contextsv.InputData()
     input_data.setVerbose(args.debug)
     input_data.setShortReadBam(args.short_read)
@@ -236,7 +236,7 @@ def main():
     input_data.saveCNVData(args.save_cnv)
     input_data.setWindowSize(args.window_size)
 
-    # Run the analysis.
+    # Run the analysis
     contextsv.run(input_data)
 
     # Determine the data paths for downstream analysis.

diff --git a/include/khmm.h b/include/khmm.h
@@ -48,7 +48,7 @@ CHMM ReadCHMM (const char *filename);
 std::pair<std::vector<int>, double> testVit_CHMM(CHMM hmm, int T, std::vector<double>& O1, std::vector<double>& O2, std::vector<double>& pfb);
 
 /// Viterbi algorithm
-std::pair<std::vector<int>, double> ViterbiLogNP_CHMM(CHMM phmm, int T, std::vector<double>& O1, std::vector<double>& O2, std::vector<double>& pfb, double **delta, int **psi, std::vector<double>& pprob);
+std::pair<std::vector<int>, double> ViterbiLogNP_CHMM(CHMM phmm, int T, std::vector<double>& O1, std::vector<double>& O2, std::vector<double>& pfb, double **delta, int **psi, double *pprob);
 
 /// O1 emission probability
 double b1iot (int state, double *mean, double *sd, double uf, double o);

diff --git a/python/cnv_plots.py b/python/cnv_plots.py
@@ -398,10 +398,6 @@ def run(cnv_data_file, output_html):
         # Append the figure to the output html file.
         output_html_file.write(fig.to_html(full_html=False, include_plotlyjs="cdn"))
         log.info("Plotted CNV %s %s:%d-%d.", 'SVType', chromosome, start_position, end_position)
-        log.info("Plot range: %d-%d", plot_start_position, plot_end_position)
-        log.info("Data before: %d", len(sv_data_before))
-        log.info("Data during: %d", len(sv_data))
-        log.info("Data after: %d", len(sv_data_after))
 
         # Increment the CNV count.
         # cnv_count += 1

diff --git a/src/cnv_caller.cpp b/src/cnv_caller.cpp
@@ -85,11 +85,11 @@ std::pair<SNPData, bool> CNVCaller::querySNPRegion(std::string chr, int64_t star
             // double baf_default = 0.5;
             double baf_default = -1.0;  // Use -1.0 to indicate no BAF data
             this->updateSNPData(snp_data, (window_start + window_end) / 2, pfb_default, baf_default, window_log2_ratio, false);
-            printMessage("No SNPs found in window " + chr + ":" + std::to_string((int)window_start) + "-" + std::to_string((int)window_end) + "...");
-            printMessage("Using position " + std::to_string((int)((window_start + window_end) / 2)) + " with default BAF and PFB values...");
+            // printMessage("No SNPs found in window " + chr + ":" + std::to_string((int)window_start) + "-" + std::to_string((int)window_end) + "...");
+            // printMessage("Using position " + std::to_string((int)((window_start + window_end) / 2)) + " with default BAF and PFB values...");
 
         } else {
-            printMessage("SNPs found in window " + chr + ":" + std::to_string((int)window_start) + "-" + std::to_string((int)window_end) + "...");
+            // printMessage("SNPs found in window " + chr + ":" + std::to_string((int)window_start) + "-" + std::to_string((int)window_end) + "...");
             snps_found = true;
 
             // Loop through the SNPs and calculate the log2 ratios
@@ -100,7 +100,7 @@ std::pair<SNPData, bool> CNVCaller::querySNPRegion(std::string chr, int64_t star
                 // Get the SNP position
                 int64_t snp_pos = snp_window_pos[j];
 
-                printMessage("SNP position: " + std::to_string((int)snp_pos));
+                // printMessage("SNP position: " + std::to_string((int)snp_pos));
 
                 // SNP bin starts at 1/2 the distance between the previous SNP
                 // and the current SNP, and ends at 1/2 the distance between
@@ -238,7 +238,7 @@ void CNVCaller::runCopyNumberPredictionChunk(std::string chr, std::map<SVCandida
 
         // Loop through the SV region, calculate the log2 ratios, and run the
         // Viterbi algorithm to predict the copy number states
-        printMessage("Querying SNPs for SV " + chr + ":" + std::to_string((int)start_pos) + "-" + std::to_string((int)end_pos) + "...");
+        // printMessage("Querying SNPs for SV " + chr + ":" + std::to_string((int)start_pos) + "-" + std::to_string((int)end_pos) + "...");
         std::pair<SNPData, bool> snp_call = this->querySNPRegion(chr, start_pos, end_pos, snp_info, pos_depth_map, mean_chr_cov);
         SNPData& sv_snps = snp_call.first;
         bool snps_found = snp_call.second;