|
|
@@ -181,108 +181,108 @@ pub fn create_bam_leave_two_out(input_path: &str) -> anyhow::Result<Vec<PathBuf>
|
|
|
Ok(output_paths)
|
|
|
}
|
|
|
|
|
|
-pub fn remap_bam(reference: &str, input_seq: &str, output_bam: &str) -> anyhow::Result<()> {
|
|
|
- duct::cmd!("bwa", "index", reference).run()?;
|
|
|
- let bwa = format!("bwa mem {reference} {input_seq}");
|
|
|
- let samtools = "samtools sort /dev/stdin";
|
|
|
- let pipe = format!("{bwa} | {samtools} > {output_bam}");
|
|
|
- duct::cmd!("bash", "-c", pipe).run()?;
|
|
|
- fs::remove_file(format!("{reference}.amb"))?;
|
|
|
- fs::remove_file(format!("{reference}.ann"))?;
|
|
|
- fs::remove_file(format!("{reference}.bwt"))?;
|
|
|
- fs::remove_file(format!("{reference}.pac"))?;
|
|
|
- fs::remove_file(format!("{reference}.sa"))?;
|
|
|
-
|
|
|
- Ok(())
|
|
|
-}
|
|
|
// pub fn remap_bam(reference: &str, input_seq: &str, output_bam: &str) -> anyhow::Result<()> {
|
|
|
-// info!("Remaping {input_seq} to {reference} into {output_bam}");
|
|
|
-// // Index the reference
|
|
|
-// let index_output = Command::new("bwa")
|
|
|
-// .args(["index", reference])
|
|
|
-// .output()
|
|
|
-// .context("Failed to run bwa index")?;
|
|
|
-//
|
|
|
-// if !index_output.status.success() {
|
|
|
-// anyhow::bail!(
|
|
|
-// "bwa index failed: {}",
|
|
|
-// String::from_utf8_lossy(&index_output.stderr)
|
|
|
-// );
|
|
|
-// }
|
|
|
-//
|
|
|
-// // Prepare the bwa mem command
|
|
|
-// let mut bwa_command = Command::new("bwa");
|
|
|
-// bwa_command.args(["mem", reference, input_seq]);
|
|
|
-// bwa_command.stdout(Stdio::piped());
|
|
|
-// bwa_command.stderr(Stdio::piped());
|
|
|
-//
|
|
|
-// let mut bwa_process = bwa_command
|
|
|
-// .spawn()
|
|
|
-// .context("Failed to spawn bwa mem command")?;
|
|
|
-//
|
|
|
-// // Prepare the samtools sort command
|
|
|
-// let mut samtools_command = Command::new("samtools");
|
|
|
-// samtools_command.args(["sort", "/dev/stdin"]);
|
|
|
-// samtools_command.stdin(Stdio::piped());
|
|
|
-// samtools_command.stdout(Stdio::from(
|
|
|
-// std::fs::File::create(output_bam).context("Failed to create output BAM file")?,
|
|
|
-// ));
|
|
|
-// samtools_command.stderr(Stdio::piped());
|
|
|
-//
|
|
|
-// let mut samtools_process = samtools_command
|
|
|
-// .spawn()
|
|
|
-// .context("Failed to spawn samtools command")?;
|
|
|
-//
|
|
|
-// // Connect bwa stdout to samtools stdin
|
|
|
-// if let (Some(mut bwa_stdout), Some(mut samtools_stdin)) =
|
|
|
-// (bwa_process.stdout.take(), samtools_process.stdin.take())
|
|
|
-// {
|
|
|
-// std::thread::spawn(move || std::io::copy(&mut bwa_stdout, &mut samtools_stdin));
|
|
|
-// }
|
|
|
-//
|
|
|
-// // Capture stderr from bwa
|
|
|
-// let mut bwa_stderr = String::new();
|
|
|
-// bwa_process
|
|
|
-// .stderr
|
|
|
-// .take()
|
|
|
-// .unwrap()
|
|
|
-// .read_to_string(&mut bwa_stderr)?;
|
|
|
-//
|
|
|
-// // Wait for bwa to finish
|
|
|
-// let bwa_status = bwa_process.wait().context("Failed to wait for bwa mem")?;
|
|
|
-// if !bwa_status.success() {
|
|
|
-// anyhow::bail!("bwa mem failed: {}", bwa_stderr);
|
|
|
-// }
|
|
|
-//
|
|
|
-// // Capture stderr from samtools
|
|
|
-// let mut samtools_stderr = String::new();
|
|
|
-// samtools_process
|
|
|
-// .stderr
|
|
|
-// .take()
|
|
|
-// .unwrap()
|
|
|
-// .read_to_string(&mut samtools_stderr)?;
|
|
|
-//
|
|
|
-// // Wait for samtools to finish
|
|
|
-// let samtools_status = samtools_process
|
|
|
-// .wait()
|
|
|
-// .context("Failed to wait for samtools")?;
|
|
|
-// if !samtools_status.success() {
|
|
|
-// anyhow::bail!("samtools sort failed: {}", samtools_stderr);
|
|
|
-// }
|
|
|
-//
|
|
|
+// duct::cmd!("bwa", "index", reference).run()?;
|
|
|
+// let bwa = format!("bwa mem {reference} {input_seq}");
|
|
|
+// let samtools = "samtools sort /dev/stdin";
|
|
|
+// let pipe = format!("{bwa} | {samtools} > {output_bam}");
|
|
|
+// duct::cmd!("bash", "-c", pipe).run()?;
|
|
|
// fs::remove_file(format!("{reference}.amb"))?;
|
|
|
// fs::remove_file(format!("{reference}.ann"))?;
|
|
|
// fs::remove_file(format!("{reference}.bwt"))?;
|
|
|
// fs::remove_file(format!("{reference}.pac"))?;
|
|
|
// fs::remove_file(format!("{reference}.sa"))?;
|
|
|
//
|
|
|
-// // Log the captured stderr
|
|
|
-// // println!("bwa mem stderr: {}", bwa_stderr);
|
|
|
-// // println!("samtools sort stderr: {}", samtools_stderr);
|
|
|
-//
|
|
|
// Ok(())
|
|
|
// }
|
|
|
-//
|
|
|
+pub fn remap_bam(reference: &str, input_seq: &str, output_bam: &str) -> anyhow::Result<()> {
|
|
|
+ info!("Remaping {input_seq} to {reference} into {output_bam}");
|
|
|
+ // Index the reference
|
|
|
+ // let index_output = Command::new("bwa")
|
|
|
+ // .args(["index", reference])
|
|
|
+ // .output()
|
|
|
+ // .context("Failed to run bwa index")?;
|
|
|
+ //
|
|
|
+ // if !index_output.status.success() {
|
|
|
+ // anyhow::bail!(
|
|
|
+ // "bwa index failed: {}",
|
|
|
+ // String::from_utf8_lossy(&index_output.stderr)
|
|
|
+ // );
|
|
|
+ // }
|
|
|
+
|
|
|
+ // Prepare the bwa mem command
|
|
|
+ let mut minimap2_command = Command::new("minimap2");
|
|
|
+ minimap2_command.args(["-a", reference, input_seq]);
|
|
|
+ minimap2_command.stdout(Stdio::piped());
|
|
|
+ minimap2_command.stderr(Stdio::piped());
|
|
|
+
|
|
|
+ let mut minimap2_process = minimap2_command
|
|
|
+ .spawn()
|
|
|
+ .context("Failed to spawn minimap2 command")?;
|
|
|
+
|
|
|
+ // Prepare the samtools sort command
|
|
|
+ let mut samtools_command = Command::new("samtools");
|
|
|
+ samtools_command.args(["sort", "/dev/stdin"]);
|
|
|
+ samtools_command.stdin(Stdio::piped());
|
|
|
+ samtools_command.stdout(Stdio::from(
|
|
|
+ std::fs::File::create(output_bam).context("Failed to create output BAM file")?,
|
|
|
+ ));
|
|
|
+ samtools_command.stderr(Stdio::piped());
|
|
|
+
|
|
|
+ let mut samtools_process = samtools_command
|
|
|
+ .spawn()
|
|
|
+ .context("Failed to spawn samtools command")?;
|
|
|
+
|
|
|
+ // Connect bwa stdout to samtools stdin
|
|
|
+ if let (Some(mut minimap2_stdout), Some(mut samtools_stdin)) =
|
|
|
+ (minimap2_process.stdout.take(), samtools_process.stdin.take())
|
|
|
+ {
|
|
|
+ std::thread::spawn(move || std::io::copy(&mut minimap2_stdout, &mut samtools_stdin));
|
|
|
+ }
|
|
|
+
|
|
|
+ // Capture stderr from bwa
|
|
|
+ let mut minimap2_stderr = String::new();
|
|
|
+ minimap2_process
|
|
|
+ .stderr
|
|
|
+ .take()
|
|
|
+ .unwrap()
|
|
|
+ .read_to_string(&mut minimap2_stderr)?;
|
|
|
+
|
|
|
+ // Wait for bwa to finish
|
|
|
+ let minimap2_status = minimap2_process.wait().context("Failed to wait for minimap2")?;
|
|
|
+ if !minimap2_status.success() {
|
|
|
+ anyhow::bail!("minimap2 failed: {}", minimap2_stderr);
|
|
|
+ }
|
|
|
+
|
|
|
+ // Capture stderr from samtools
|
|
|
+ let mut samtools_stderr = String::new();
|
|
|
+ samtools_process
|
|
|
+ .stderr
|
|
|
+ .take()
|
|
|
+ .unwrap()
|
|
|
+ .read_to_string(&mut samtools_stderr)?;
|
|
|
+
|
|
|
+ // Wait for samtools to finish
|
|
|
+ let samtools_status = samtools_process
|
|
|
+ .wait()
|
|
|
+ .context("Failed to wait for samtools")?;
|
|
|
+ if !samtools_status.success() {
|
|
|
+ anyhow::bail!("samtools sort failed: {}", samtools_stderr);
|
|
|
+ }
|
|
|
+
|
|
|
+ // fs::remove_file(format!("{reference}.amb"))?;
|
|
|
+ // fs::remove_file(format!("{reference}.ann"))?;
|
|
|
+ // fs::remove_file(format!("{reference}.bwt"))?;
|
|
|
+ // fs::remove_file(format!("{reference}.pac"))?;
|
|
|
+ // fs::remove_file(format!("{reference}.sa"))?;
|
|
|
+
|
|
|
+ // Log the captured stderr
|
|
|
+ // println!("bwa mem stderr: {}", bwa_stderr);
|
|
|
+ // println!("samtools sort stderr: {}", samtools_stderr);
|
|
|
+
|
|
|
+ Ok(())
|
|
|
+}
|
|
|
+
|
|
|
pub fn run_modkit(bam_path: &str, contig_fa: &str, pileup_path: &str) -> anyhow::Result<()> {
|
|
|
info!("Running Modkit pileup for {bam_path}");
|
|
|
let mut cmd = Command::new("modkit")
|