Create another check on updating variant based on long read coverage.

seqan · Nov 29, 2019 · b30c062 · b30c062
1 parent d7652ea
commit b30c062
Show file tree

Hide file tree

Showing 4 changed files with 35 additions and 5 deletions.
diff --git a/src/alignment.cpp b/src/alignment.cpp
@@ -80,8 +80,8 @@ bool AlignmentManager::load()
             this->pairedEndRead->prepAfterHeaderParsing(this->bamHeader, this->bamFileIn);
         if ( optionManager->doClippedReadAnalysis() )
             this->clippedRead->prepAfterHeaderParsing(this->bamHeader, this->bamFileIn);
-        if ( optionManager->doReadDepthAnalysis(this->isLongRead) )
-            this->readDepth->prepAfterHeaderParsing(this->bamHeader, this->bamFileIn);
+        if ( optionManager->doReadDepthAnalysis() )
+            this->readDepth->prepAfterHeaderParsing(this->bamHeader, this->bamFileIn, this->isLongRead);
 
         seqan::CharString qNameWithPairInfo;
         seqan::BamAlignmentRecord record;

diff --git a/src/readdepth.cpp b/src/readdepth.cpp
@@ -63,6 +63,35 @@ void ReadDepth::prepAfterHeaderParsing(seqan::BamHeader& header, seqan::BamFileI
             std::fill(readDepthInfo[rID].begin(), readDepthInfo[rID].end(), 0);
         }
     }
+    if (isLongRead)
+    {
+        std::vector<double> backgroundDepth{};
+        this->setRandomSeed(0);
+        for (unsigned i = 0; i < this->getOptionManager()->getSamplingNum(); i++)
+        {
+            TPosition p;
+            TTemplateID t;
+            double kl = 0.0, depth = 0.0;
+            double leftDepthAvg, rightDepthAvg;
+
+            // get kl and depth at random position
+            while (depth <= 0.0)
+            {
+                this->getRandomPos(t, p);
+                getAvgReadDepth(leftDepthAvg, rightDepthAvg, t, p, this->getOptionManager()->getReadDepthWindowSize(),
+                                this->getOptionManager()->getReadDepthWindowSize(), BreakpointEvidence::SIDE::LEFT);
+
+                depth = leftDepthAvg + rightDepthAvg / 2.0;
+            }
+            backgroundDepth.push_back(depth);
+        }
+
+        std::sort(backgroundDepth.begin(), backgroundDepth.end(), [](auto &left, auto &right) {return left < right;} );
+        double medianDepth = MID_ELEMENT(backgroundDepth);
+        this->setMedianDepth( medianDepth );
+
+        printMessage("Depth median for long reads: " + std::to_string(medianDepth));
+    }
 }
 
 void ReadDepth::setRandomSeed(int seed)

diff --git a/src/sv.hpp b/src/sv.hpp
@@ -122,9 +122,9 @@ class VcfRecordEnhanced : public seqan::VcfRecord
             return tmp;
         }
 
-        void update(sviper::Variant const & tmp)
+        void update(sviper::Variant const & tmp, int lr_depth)
         {
-            if (this->se + this->ce + this->pe < 8)
+            if (lr_depth >= 15 || (this->se + this->ce + this->pe < 8))
             {
                 this->status = (tmp.filter == "FAIL1" || tmp.filter == "FAIL2" || tmp.filter == "FAIL5") ? STATUS::FILTERED : STATUS::PASS;
                 this->filter = (tmp.filter == "FAIL1" || tmp.filter == "FAIL2" || tmp.filter == "FAIL5") ? tmp.filter : STATUS_PASS();

diff --git a/src/vaquita.cpp b/src/vaquita.cpp
@@ -159,6 +159,7 @@ int callMain(int argc, char const ** argv)
         options.mean_insert_size_of_short_reads = alnMgr.getInsMedian();
         options.stdev_insert_size_of_short_reads = alnMgr.getInsSD();
         sviper::input_output_information info{options};
+        int lr_depth = static_cast<int>(std::round(bpMgrLR.getReadDepth()->getMedianDepth()));
 
         // Check files
         if (!sviper::open_file_success(info.long_read_bai, (info.cmd_options.long_read_file_name + ".bai").c_str()) ||
@@ -186,7 +187,7 @@ int callMain(int argc, char const ** argv)
                 sviper::Variant tmp{deletions_sviper[vidx]};
                 tmp.ref_chrom = std::string{seqan::toCString(alnMgr.getRefName(std::stoi(tmp.ref_chrom)))};
                 sviper::polish_variant(tmp, info);
-                deletions_sviper[vidx].update(tmp);
+                deletions_sviper[vidx].update(tmp, lr_depth);
             }
         } // parallel for loop
         endTimeMessage("SViper: POLISHING VARIANTS");