Commands should output fastq rather than fasta (#60)

* Commands should output fastq rather than fasta

Author: kvg

Cargo.lock

@@ -42,4 +42,7 @@ plotters = "=0.3.7"
 
 [profile.release]
 opt-level = 2
-codegen-units = 16
+codegen-units = 16
+
+[dev-dependencies]
+tempfile = "3.8.1"

src/hidive/Cargo.toml

@@ -165,7 +165,7 @@ enum Commands {
     /// Stream loci from CRAM/BAM/FASTA files stored locally or in Google Cloud Storage.
     #[clap(arg_required_else_help = true)]
     Fetch {
-        /// Output path for FASTA file with reads spanning locus of interest.
+        /// Output path for FASTQ file with reads spanning locus of interest.
         #[clap(short, long, value_parser, default_value = "/dev/stdout")]
         output: PathBuf,
 
@@ -189,7 +189,7 @@ enum Commands {
     /// Find more sequences (both aligned and unaligned) overlapping previously fetched reads.
     #[clap(arg_required_else_help = true)]
     Rescue {
-        /// Output path for FASTA file with reads spanning locus of interest.
+        /// Output path for FASTQ file with reads spanning locus of interest.
         #[clap(short, long, value_parser, default_value = "/dev/stdout")]
         output: PathBuf,
 
@@ -217,7 +217,7 @@ enum Commands {
         #[clap(short, long, value_parser, required = true)]
         fasta_paths: Vec<PathBuf>,
 
-        /// Indexed WGS BAM, CRAM, or FASTA files from which to extract relevant sequences.
+        /// Indexed WGS BAM, CRAM, or FASTA files from which to find relevant sequences.
         #[clap(required = true, value_parser)]
         seq_paths: Vec<PathBuf>,
     },

src/hidive/src/main.rs

@@ -1,5 +1,5 @@
 use bio::data_structures::interval_tree::IntervalTree;
-use bio::io::fasta::Reader;
+use bio::io::fastq::{self, Reader};
 use bio::utils::Interval;
 use clap::ValueEnum;
 use flate2::write::GzEncoder;
@@ -10,7 +10,7 @@ use num_format::{Locale, ToFormattedString};
 use std::cmp::PartialEq;
 use std::collections::{BTreeMap, HashMap};
 use std::fs::File;
-use std::io::Write;
+use std::io::{BufWriter, Write};
 use std::sync::atomic::{AtomicUsize, Ordering};
 use std::sync::Arc;
 use std::{collections::HashSet, path::PathBuf};
@@ -42,10 +42,10 @@ pub fn start(
     search_option: SearchOption,
     ref_path: Option<PathBuf>,
     loci: Option<Vec<String>>,
-    fasta_paths: &Vec<PathBuf>,
+    fastq_paths: &Vec<PathBuf>,
     seq_paths: &Vec<PathBuf>,
 ) {
-    let fasta_urls = skydive::parse::parse_file_names(fasta_paths);
+    let fastq_urls = skydive::parse::parse_file_names(fastq_paths);
     let seq_urls = skydive::parse::parse_file_names(seq_paths);
 
     // Get the system's temporary directory path
@@ -57,8 +57,8 @@ pub fn start(
 
     // Read the FASTA files and prepare a hashmap of k-mers to search for in the other files.
     let mut kmer_set = HashSet::new();
-    for fasta_url in fasta_urls {
-        let reader = Reader::from_file(fasta_url.to_file_path().unwrap()).unwrap();
+    for fastq_url in fastq_urls {
+        let reader = Reader::from_file(fastq_url.to_file_path().unwrap()).unwrap();
 
         let mut reads = Vec::new();
         for record in reader.records().flatten() {
@@ -153,17 +153,35 @@ pub fn start(
             );
         }
 
+        /*
         let file = File::create(output).expect("Could not open file for writing.");
         let mut writer: Box<dyn Write> = if output.to_str().unwrap().ends_with(".gz") {
             Box::new(GzEncoder::new(file, Compression::default()))
         } else {
             Box::new(file)
         };
+        */
+
+        let mut buf_writer = BufWriter::new(File::create(output).unwrap());
+        let mut fastq_writer = fastq::Writer::new(&mut buf_writer);
 
         for record in all_records.iter() {
-            let fw_seq = record.seq().as_bytes();
-            let rc_seq = fw_seq.reverse_complement();
+            if record.is_reverse() {
+                let rv_seq = record.seq().as_bytes().reverse_complement();
+                let mut rv_qual = record.qual().iter().map(|&q| q + 33).collect::<Vec<u8>>();
+                rv_qual.reverse();
+                let rv_record = fastq::Record::with_attrs(&String::from_utf8_lossy(record.qname()), None, &rv_seq, &rv_qual);
+
+                fastq_writer.write_record(&rv_record).unwrap();
+            } else {
+                let fw_seq = record.seq().as_bytes();
+                let fw_qual = record.qual().iter().map(|&q| q + 33).collect::<Vec<u8>>();
+                let fw_record = fastq::Record::with_attrs(&String::from_utf8_lossy(record.qname()), None, &fw_seq, &fw_qual);
+
+                fastq_writer.write_record(&fw_record).unwrap();
+            }
 
+            /*
             writeln!(
                 writer,
                 ">read_{}_{}_{}_{}",
@@ -183,6 +201,7 @@ pub fn start(
                 }
             )
             .expect("Could not write to file");
+            */
         }
 
         if !all_records.is_empty() {

src/hidive/src/rescue.rs

@@ -0,0 +1,164 @@
+use std::path::PathBuf;
+use std::process::Command;
+use tempfile::TempDir;
+
+/// Validates that a FASTQ file exists and contains valid records
+fn validate_fastq_output(path: &PathBuf) {
+    let reader = bio::io::fastq::Reader::from_file(path).expect("Failed to open output file");
+    let mut records = Vec::new();
+
+    for result in reader.records() {
+        // Each record should be valid and parseable
+        let record = result.expect("Error reading FASTQ record");
+        records.push(record);
+    }
+
+    // Basic validation that we got some records
+    assert!(!records.is_empty(), "No FASTQ records found in output file");
+
+    // Validate each record has the expected components
+    for record in records {
+        // Check sequence name exists
+        assert!(!record.id().is_empty(), "FASTQ record missing sequence ID");
+        
+        // Check sequence exists and is not empty
+        assert!(!record.seq().is_empty(), "FASTQ record has empty sequence");
+        
+        // Check quality scores exist and match sequence length
+        assert_eq!(record.seq().len(), record.qual().len(), 
+            "Quality score length does not match sequence length");
+    }
+}
+
+#[test]
+fn test_fetch_basic() {
+    // Create a temporary directory for test files
+    let temp_dir = TempDir::new().unwrap();
+    let output_path = temp_dir.path().join("output.fastq");
+    
+    // Test BAM file
+    let test_bam = "gs://fc-8c3900db-633f-477f-96b3-fb31ae265c44/results/PBFlowcell/m84060_230907_210011_s2/reads/ccs/aligned/m84060_230907_210011_s2.bam";
+    
+    // Run the fetch command
+    let output = Command::new("cargo")
+        .args([
+            "run",
+            "--",
+            "fetch",
+            "--loci",
+            "chr15:23960193-23963918",
+            "--output",
+            output_path.to_str().unwrap(),
+            test_bam,
+        ])
+        .output()
+        .expect("Failed to execute command");
+
+    // Check if the command was successful
+    assert!(output.status.success(), "Command failed: {:?}", output);
+    
+    // Verify the output file exists
+    assert!(output_path.exists(), "Output file was not created");
+
+    // Validate the format of the output FASTQ file
+    validate_fastq_output(&output_path);
+}
+
+#[test]
+fn test_fetch_with_padding() {
+    let temp_dir = TempDir::new().unwrap();
+    let output_path = temp_dir.path().join("output.fastq");
+    let test_bam = "gs://fc-8c3900db-633f-477f-96b3-fb31ae265c44/results/PBFlowcell/m84060_230907_210011_s2/reads/ccs/aligned/m84060_230907_210011_s2.bam";
+    
+    let output = Command::new("cargo")
+        .args([
+            "run",
+            "--",
+            "fetch",
+            "--loci",
+            "chr15:23960193-23963918",
+            "--padding",
+            "100",
+            "--output",
+            output_path.to_str().unwrap(),
+            test_bam
+        ])
+        .output()
+        .expect("Failed to execute command");
+
+    assert!(output.status.success(), "Command failed: {:?}", output);
+    assert!(output_path.exists(), "Output file was not created");
+    validate_fastq_output(&output_path);
+}
+
+// #[test]
+// fn test_fetch_with_haplotype() {
+//     let temp_dir = TempDir::new().unwrap();
+//     let output_path = temp_dir.path().join("output.fastq");
+//     let test_bam = "gs://fc-8c3900db-633f-477f-96b3-fb31ae265c44/results/PBFlowcell/m84060_230907_210011_s2/reads/ccs/aligned/m84060_230907_210011_s2.bam";
+    
+//     let output = Command::new("cargo")
+//         .args([
+//             "run",
+//             "--",
+//             "fetch",
+//             "--loci",
+//             "chr15:23960193-23963918",
+//             "--haplotype",
+//             "1",
+//             "--output",
+//             output_path.to_str().unwrap(),
+//             test_bam,
+//         ])
+//         .output()
+//         .expect("Failed to execute command");
+
+//     assert!(output.status.success(), "Command failed: {:?}", output);
+//     assert!(output_path.exists(), "Output file was not created");
+//     validate_fastq_output(&output_path);
+// }
+
+#[test]
+fn test_fetch_invalid_locus() {
+    let temp_dir = TempDir::new().unwrap();
+    let output_path = temp_dir.path().join("output.fastq");
+    let test_bam = "gs://fc-8c3900db-633f-477f-96b3-fb31ae265c44/results/PBFlowcell/m84060_230907_210011_s2/reads/ccs/aligned/m84060_230907_210011_s2.bam";
+    
+    let output = Command::new("cargo")
+        .args([
+            "run",
+            "--",
+            "fetch",
+            "--loci",
+            "invalid_locus",
+            "--output",
+            output_path.to_str().unwrap(),
+            test_bam,
+        ])
+        .output()
+        .expect("Failed to execute command");
+
+    // This should fail because the locus format is invalid
+    assert!(!output.status.success(), "Command should have failed");
+}
+
+#[test]
+fn test_fetch_missing_required_args() {
+    let temp_dir = TempDir::new().unwrap();
+    let output_path = temp_dir.path().join("output.fastq");
+    
+    let output = Command::new("cargo")
+        .args([
+            "run",
+            "--",
+            "fetch",
+            "--output",
+            output_path.to_str().unwrap(),
+            // Missing required --loci argument and BAM file
+        ])
+        .output()
+        .expect("Failed to execute command");
+
+    // This should fail because required arguments are missing
+    assert!(!output.status.success(), "Command should have failed");
+}