v1.1

Author: Yichen-HNU

README.md

@@ -4,6 +4,7 @@ Repeat and haplotype aware error correction in nanopore sequencing reads with De
 
 Error correction is the canonical first step in long-read sequencing data analysis. Nanopore R10 reads have error rates below 2\%. we introduce DeChat, a novel approach specifically designed for Nanopore R10 reads.DeChat enables repeat- and haplotype-aware error correction, leveraging the strengths of both de Bruijn graphs and variant-aware multiple sequence alignment to create a synergistic approach. This approach avoids read overcorrection, ensuring that variants in repeats and haplotypes are preserved while sequencing errors are accurately corrected.
 
+DeChat can use HIFi or NGS to correct ONT now
 
 Dechat is implemented with C++.
 
@@ -48,7 +49,7 @@ Usage: dechat [options] -o <output> -t <thread>  -i <reads> <...>
 Options:
   Input/Output:
        -o STR       prefix of output files [(null)]
-                    The output for the stage 1 of correction is "recorrected.fa", 
+                   The output for the stage 1 of correction is "recorrected.fa", 
                     The final corrected file is "file name".ec.fa;
        -t INT       number of threads [1]
        -h           show help information
@@ -57,9 +58,10 @@ Options:
        -k INT       k-mer length (must be <64) [21]
   Error correction stage 1 (dBG):
        -r1           set the maximal abundance threshold for a k-mer in dBG [2]
+       -d           input reads file for building dBG (Default use input ONT reads) 
   Error correction stage 2 (MSA):
        -r            round of correction in alignment [3]
-       -e            maximum allowed error rate used for filtering overlaps [0.04]       
+       -e            maximum allowed error rate used for filtering overlaps [0.04]     
 ```
 
 ## Examples
@@ -69,6 +71,10 @@ The example folder contains test data, including the 10X depth sim-ont10.4 data
 cd example
 dechat -i reads.fa.gz -o reads -t 8
 ```
-
+### Using HIFi or NGS to correct ONT
+```
+cd example
+dechat -i reads.fa.gz -o reads -t 8 -d HiFi-reads.fq.gz/NGS-reads.fq.g
+```
 
 

aligner-correct/correct_round1.h

@@ -9,11 +9,13 @@
 #include "ketopt.h"
 #include <iostream>
 
-void correct_round2(chat_opt_t *chat_opt,hifiasm_opt_t* asm_opt);
+void correct_round2(chat_opt_t *chat_opt, hifiasm_opt_t *asm_opt);
 void correct_round1(chat_opt_t *chat_opt)
 {
-    std::cout << "correct_round1" << chat_opt->thread_num << std::endl;
-
+    std::cout << "correct_round1 thread:" << chat_opt->thread_num << std::endl;
+    // PRINT_LINE_FUNC();
+    // if (chat_opt->dBGFile != NULL)
+    //     std::cout << "dBG use " << chat_opt->dBGFile << std::endl;
     std::string pacbioFile = chat_opt->read_file_names;
     BankFasta bsize(pacbioFile);
     BankFasta::Iterator itSeqSize(bsize);
@@ -29,12 +31,22 @@ void correct_round1(chat_opt_t *chat_opt)
         }
         nbSeq++;
     }
+    // PRINT_LINE_FUNC();
     max_read_len = max_read_len * 1.25;
     std::string ONTGraph;
     int comaPosition = std::string::npos;
-    PRINT_LINE_FUNC();
-
-    ONTGraph = chat_opt->read_file_names;
+    // PRINT_LINE_FUNC();
+    if (chat_opt->dBGFile == NULL)
+    {
+        std::cout << "dBG use ONT reads"<<chat_opt->read_file_names<< std::endl;
+        ONTGraph = chat_opt->read_file_names;
+        std::cerr << "creating the graph from file(s): " << chat_opt->read_file_names << std::endl;
+    }else{
+        std::cout << "dBG use " << chat_opt->dBGFile << std::endl;
+        std::cerr << "creating the graph from file(s): " << chat_opt->dBGFile << std::endl;
+        ONTGraph = chat_opt->dBGFile;
+    }
+    
     comaPosition = ONTGraph.find(",");
 
     if (comaPosition != std::string::npos)
@@ -62,21 +74,32 @@ void correct_round1(chat_opt_t *chat_opt)
     }
 
     Graph graph;
-    PRINT_LINE_FUNC();
+    // PRINT_LINE_FUNC();
     if (DEBUG_l)
     {
-        std::cerr << "creating the graph from file(s): " << chat_opt->read_file_names << std::endl;
+        
     }
     try
     {
         // v106: open IBank from 1/ list of filenames 2/ a file of filenames
-        IBank *b = Bank::open(chat_opt->read_file_names);
+        IBank *b;
+        if (chat_opt->dBGFile == NULL)
+        {
+            std::cout << "dBG use ONT reads" << std::endl;
+            b = Bank::open(chat_opt->read_file_names);
+        }
+        else
+        {
+            std::cout << "dBG use " << chat_opt->dBGFile << std::endl;
+            b = Bank::open(chat_opt->dBGFile);
+        }
+
         std::string outTmpPath = "";
-        //outTmpPath = outTmpPath + dirname(std::string(chat_opt->read_file_names)) + "recorrected.fa";
-        PRINT_LINE_FUNC();
+        // outTmpPath = outTmpPath + dirname(std::string(chat_opt->read_file_names)) + "recorrected.fa";
+        // PRINT_LINE_FUNC();
         std::cout << outTmpPath.c_str() << "   ONTGraph:" << ONTGraph << std::endl;
-        graph = Graph::create(b, (const char *)"%s -out %s -kmer-size %d -abundance-min %d -bloom cache -debloom original -debloom-impl basic -nb-cores %d -abundance-max 2147483647", outTmpPath.c_str(), ONTGraph.c_str(), chat_opt->k_mer_length,chat_opt->abundance_min, chat_opt->thread_num);
-        PRINT_LINE_FUNC();
+        graph = Graph::create(b, (const char *)"%s -out %s -kmer-size %d -abundance-min %d -bloom cache -debloom original -debloom-impl basic -nb-cores %d -abundance-max 2147483647", outTmpPath.c_str(), ONTGraph.c_str(), chat_opt->k_mer_length, chat_opt->abundance_min, chat_opt->thread_num);
+        // PRINT_LINE_FUNC();
         if (is_readable(ONTGraph))
         {
             std::cerr << "!!! file present : " << ONTGraph << std::endl;
@@ -201,7 +224,7 @@ void correct_round1(chat_opt_t *chat_opt)
     }
 }
 
-void correct_round2(chat_opt_t *chat_opt,hifiasm_opt_t* asm_opt)
+void correct_round2(chat_opt_t *chat_opt, hifiasm_opt_t *asm_opt)
 {
     int ret;
     yak_reset_realtime();
@@ -210,25 +233,25 @@ void correct_round2(chat_opt_t *chat_opt,hifiasm_opt_t* asm_opt)
     std::cout << "chenggong" << std::endl;
     //.................传入参数..................
     asm_opt->num_reads = 1;
-    asm_opt->read_file_names = new char*[asm_opt->num_reads];
+    asm_opt->read_file_names = new char *[asm_opt->num_reads];
     asm_opt->thread_num = chat_opt->thread_num;
     asm_opt->read_file_names[0] = strdup(chat_opt->output_dir_ec.c_str());
 
     //
     std::cout << asm_opt->read_file_names[0] << std::endl;
     std::string outputname = chat_opt->outReadFile;
-    std::cout<<"chat_opt->outReadFile:"<<chat_opt->outReadFile<<std::endl;
+    std::cout << "chat_opt->outReadFile:" << chat_opt->outReadFile << std::endl;
     asm_opt->output_file_name = chat_opt->outReadFile;
-    std::cout<<"asm_opt->output_file_name:"<<asm_opt->output_file_name<<std::endl;
+    std::cout << "asm_opt->output_file_name:" << asm_opt->output_file_name << std::endl;
 
     asm_opt->max_ov_diff_ec = chat_opt->max_ov_diff_ec;
     asm_opt->number_of_round = chat_opt->second_number_of_round;
-    std::cout<<"asm_opt->max_ov_diff_ec:"<<asm_opt->max_ov_diff_ec<<std::endl;
-    std::cout<<"asm_opt->number_of_round:"<<asm_opt->number_of_round<<std::endl;
+    std::cout << "asm_opt->max_ov_diff_ec:" << asm_opt->max_ov_diff_ec << std::endl;
+    std::cout << "asm_opt->number_of_round:" << asm_opt->number_of_round << std::endl;
     check_option1(asm_opt);
 
     ret = ha_assemble();
-    
+
     return;
 }
 

bin/dechat

@@ -86,6 +86,7 @@ void Print_help(chat_opt_t *chat_opt)
     fprintf(stderr, "       -k INT       k-mer length (must be <64) [%d]\n", chat_opt->k_mer_length);
     fprintf(stderr, "  Error correction stage 1 (dBG):\n");
     fprintf(stderr, "       -r1           set the maximal abundance threshold for a k-mer in dBG [%d]\n",chat_opt->abundance_min);
+    fprintf(stderr, "       -d           input reads file for building dBG (Default use input ONT reads) \n");
     fprintf(stderr, "  Error correction stage 2 (MSA):\n");
     fprintf(stderr, "       -r            round of correction in alignment [%d]\n",chat_opt->second_number_of_round);
     fprintf(stderr, "       -e            maximum allowed error rate used for filtering overlaps [%.2f]\n",chat_opt->max_ov_diff_ec);
@@ -97,6 +98,7 @@ void init_opt(chat_opt_t *chat_opt)
     /// chat_opt->flag = 0;
     // chat_opt->flag = HA_F_PARTITION;
     chat_opt->abundance_min = 2;
+    chat_opt->dBGFile = NULL;
     chat_opt->coverage = -1;
     chat_opt->num_reads = 0;
     chat_opt->fast = 0;
@@ -226,7 +228,7 @@ int check_option(chat_opt_t *chat_opt)
         fprintf(stderr, "[ERROR] the number of threads must be > 0 (-t)\n");
         return 0;
     }
-    std::cout << chat_opt->second_number_of_round << std::endl;
+    //std::cout << chat_opt->second_number_of_round << std::endl;
     if (chat_opt->second_number_of_round < 1)
     {
         fprintf(stderr, "[ERROR] the number of rounds for correction must be > 0 (-r)\n");
@@ -242,8 +244,9 @@ int Dechat_command(int argc, char *argv[], chat_opt_t *chat_opt, hifiasm_opt_t *
     int c;
     int option_index = 0;
     // PRINT_LINE_FUNC();
-    while ((c = ketopt(&opt, argc, argv, 1, "hvt:i:k:o:r:r1:e:hifi", long_options1)) >= 0)
+    while ((c = ketopt(&opt, argc, argv, 1, "hvt:i:d:k:o:r:r1:e:hifi", long_options1)) >= 0)
     {
+        // std::cout<<c<<std::endl;
         if (c == 'h')
         {
             Print_help(chat_opt);
@@ -260,6 +263,11 @@ int Dechat_command(int argc, char *argv[], chat_opt_t *chat_opt, hifiasm_opt_t *
         }
         else if (c == 't')
             chat_opt->thread_num = atoi(opt.arg);
+        else if (c == 'd'){
+            chat_opt->dBGFile = opt.arg;
+            // std::cout<<"chat_opt->dBGFile:"<<chat_opt->dBGFile<<std::endl;
+        }
+            
         else if (c == 'k')
             chat_opt->k_mer_length = atoi(opt.arg);
         else if (c == 'i')

src/CommandLines.cpp

@@ -13,6 +13,7 @@
 typedef struct {
     //...................................Test.......................................
     char* read_file_names;
+    char* dBGFile;
     int thread_num;
     int k_mer_length;
     char* k_mer_length_str;