add subcommnad join for merge overlap pe reads

Author: sharkLoc

src/cli/join.rs

@@ -0,0 +1,111 @@
+use crate::utils::*;
+use anyhow::{Ok, Result};
+use bio::io::fastq;
+use log::*;
+
+fn reverse_complement(seq: &[u8]) -> Vec<u8> {
+    seq.iter()
+        .rev()
+        .map(|b| match b {
+            b'A' => b'T',
+            b'T' => b'A',
+            b'C' => b'G',
+            b'G' => b'C',
+            _ => *b, // unknow base like N
+        })
+        .collect()
+}
+
+pub fn join_overlap(
+    read1: &str,
+    read2: &str,
+    max_mismatch_rate: f64,
+    min_overlap_len: usize,
+    overlap_merge: Option<&String>,
+    nonoverlap_pe: Option<&String>,
+    compression_level: u32,
+    stdout_type: char,
+) -> Result<()> {
+
+    let reader1 = fastq::Reader::new(file_reader(Some(read1))?);
+    let reader2 = fastq::Reader::new(file_reader(Some(read2))?);
+    
+    let mut count_join = 0usize;
+    let mut count_total = 0usize;
+    let mut writer_single = fastq::Writer::new(file_writer(overlap_merge, compression_level, stdout_type)?);
+    let mut nuoverlap_writer = fastq::Writer::new(file_writer(nonoverlap_pe, compression_level, stdout_type)?);
+
+    for (rec1,rec2) in reader1.records().map_while(Result::ok).zip(reader2.records().map_while(Result::ok)) {
+        count_total += 1;
+        let max_overlap_len = rec1.seq().len().min(rec2.seq().len());
+        let rec2_seq_rev = reverse_complement(rec2.seq());
+    
+        let mut fine_overlap_len = 0;
+        if min_overlap_len <= max_overlap_len {
+            for overlap_len in (min_overlap_len..=max_overlap_len).rev() {
+                assert!(rec1.seq().len() >= overlap_len && rec2_seq_rev.len() >= overlap_len);
+                // overlap region in PE reads
+                let over1 = &rec1.seq()[rec1.seq().len()-overlap_len..];
+                let over2 = &rec2_seq_rev[..overlap_len];
+
+                let mismatch = over1.iter().zip(over2.iter()).filter(|(x,y)| x != y).count();
+                let mismatch_rate = mismatch as f64 / overlap_len as f64;
+
+                if max_mismatch_rate < mismatch_rate { continue; } // mismatch count too much
+                if overlap_len < min_overlap_len { continue; } // overlap length is too short
+                // pe reads overlaped
+                fine_overlap_len = overlap_len;
+                break;
+            }
+        }
+
+        // build longer single read
+        if fine_overlap_len > 0{
+            count_join += 1; 
+            let mut single_seq = vec![];
+            let mut single_qual = vec![];
+            single_seq.extend_from_slice(&rec1.seq()[..rec1.seq().len()-fine_overlap_len]);
+            single_qual.extend_from_slice(&rec1.qual()[..rec1.qual().len()-fine_overlap_len]);
+        
+            let overlap_r1_qual = &rec1.qual()[rec1.qual().len()-fine_overlap_len..];
+            let overlap_r1_seq = &rec1.seq()[rec1.seq().len()-fine_overlap_len..];
+
+            let rec2_qual_rev = rec2.qual().iter().copied().rev().collect::<Vec<u8>>();
+            let overlap_r2_qual = &rec2_qual_rev[..fine_overlap_len];
+            let overlap_r2_seq = &rec2_seq_rev[..fine_overlap_len];
+
+            for i in 0..fine_overlap_len {
+                // Prefer R1 if quality is equal
+                if overlap_r1_qual[i] >= overlap_r2_qual[i] {
+                    single_seq.push(overlap_r1_seq[i]);
+                    single_qual.push(overlap_r1_qual[i]);
+                } else {
+                    single_seq.push(overlap_r2_seq[i]); 
+                    single_qual.push(overlap_r2_qual[i]);
+                };
+            }
+            single_seq.extend_from_slice(&rec2_seq_rev[fine_overlap_len..]);
+            single_qual.extend_from_slice(&rec2_qual_rev[fine_overlap_len..]);
+
+            let read_long = fastq::Record::with_attrs(rec1.id(), rec1.desc(), &single_seq, &single_qual);
+            writer_single.write_record(&read_long)?;
+        } else {
+            // no overlap
+            if nonoverlap_pe.is_some() {
+                nuoverlap_writer.write_record(&rec1)?;
+                nuoverlap_writer.write_record(&rec2)?;
+            }
+
+        }
+    }
+    writer_single.flush()?;
+    nuoverlap_writer.flush()?;
+    
+    let rate = count_join as f64 / count_total as f64;
+    info!("total pe reads overlaped and joined number: {}", count_join);
+    info!("total pe reads number: {}", count_total);
+    info!("pe reads overlap rate: {:.2}%", rate*100.0);
+
+    Ok(())
+}
+

src/cli/mod.rs

@@ -9,6 +9,7 @@ pub mod fq2sam;
 pub mod fqscore;
 pub mod gcplot;
 pub mod grep;
+pub mod join;
 pub mod kmer;
 pub mod length;
 pub mod mask;

src/command.rs

@@ -227,6 +227,36 @@ pub enum Subcli {
         #[arg(short = 'r', long = "out2", value_name = "FILE")]
         out2: String,
     },
+    /// join paired end reads that are overlapping into a single longer read
+    #[command(before_help=r"Note:
+    1. if overlapping regions are detected, low quality bases are corrected by high quality paired bases. 
+    2. if a base is corrected, the quality value is also corrected.
+    3. only paired end reads such as the following will be detected and merged, overlap mode:
+       r1: GCAAGCGTTAATCGGAATTTATGGGCGTAAAGCGCACGCAGGA
+                            |\|||||||||||||||||||||||\ 
+                            TCTATGGGCGTAAAGCGCACGCAGGCATGCTGGGCGTAAAGCGCACGCAGGC  r2: reverse complement                       
+    ")]
+    join {
+        /// input read1 fastq file
+        #[arg(short = '1', long = "read1", value_name = "FILE")]
+        read1: String,
+        /// input read2 fastq file
+        #[arg(short = '2', long = "read2", value_name = "FILE")]
+        read2: String,
+        /// minimum overlap length in PE reads
+        #[arg(short = 'l', long = "length", default_value_t=30)]
+        length: usize,
+        /// maximum mismatch rate count in overlap region
+        #[arg(short = 'm', long = "miss", default_value_t=0.1)]
+        miss: f64,
+        /// output joinde long fastq file name or write to stdout, file ending in .gz/.bz2/.xz will be compressed automatically
+        #[arg(short = 'o', long = "output")]
+        output: Option<String>,
+        /// output interleaved fastq file name for non-overlap pe reads
+        #[arg(short = 'n', long = "non-overlap", value_name = "FILE")]
+        non: Option<String>,
+
+    },
     /// print fastq records in a range
     range {
         /// input fastq file, or read from stdin

src/main.rs

@@ -11,11 +11,11 @@ mod utils;
 use command::*;
 mod cli;
 use cli::{
-    barcode::*, check::*, concat::*, cutadapter::cut_adapter, filter::*, flatten::*, fq2fa::*,
+    barcode::*, check::*, concat::*, cutadapter::*, filter::*, flatten::*, fq2fa::*,
     fq2sam::*, fqscore::*, gcplot::*, grep::*, kmer::*, length::*, mask::*, merge::*, plot::*,
     range::*, remove::*, rename::*, reverse::*, search::*, select::*, shuffle::*, size::*,
     slide::*, sort::*, split::*, split2::*, stats::*, subfq::*, tail::*, top::*, trimfq::*,
-    view::*,
+    view::*, join::*,
 };
 
 fn main() -> Result<(), Error> {
@@ -557,6 +557,18 @@ fn main() -> Result<(), Error> {
                 arg.stdout_type,
             )?;
         }
+        Subcli::join { read1, read2, length, miss, output, non } => {
+            join_overlap(
+                &read1, 
+                &read2,
+                miss,
+                length,
+                output.as_ref(),
+                non.as_ref(),
+                arg.compression_level,
+                arg.stdout_type
+            )?;
+        }
         Subcli::kmer {
             input,
             size,