|
@@ -75,8 +75,9 @@ pub fn cp_mod_tags(input_records: &Vec<Record>, output_path: &str) -> anyhow::Re
|
|
|
let tags = record.tags_mut();
|
|
let tags = record.tags_mut();
|
|
|
|
|
|
|
|
if let Some(from) = tags_hm.get(&record_name) {
|
|
if let Some(from) = tags_hm.get(&record_name) {
|
|
|
- from.iter().filter(|(name, _)| !omit_tags.contains(&name)).for_each(|(name, tag)| {
|
|
|
|
|
- match tag {
|
|
|
|
|
|
|
+ from.iter()
|
|
|
|
|
+ .filter(|(name, _)| !omit_tags.contains(&name))
|
|
|
|
|
+ .for_each(|(name, tag)| match tag {
|
|
|
TagValue::Char(v) => tags.push_char(&name, v),
|
|
TagValue::Char(v) => tags.push_char(&name, v),
|
|
|
TagValue::Int(v, iv) => match iv {
|
|
TagValue::Int(v, iv) => match iv {
|
|
|
bam::record::tags::IntegerType::I8 => tags.push_num(&name, v as i8),
|
|
bam::record::tags::IntegerType::I8 => tags.push_num(&name, v as i8),
|
|
@@ -93,8 +94,7 @@ pub fn cp_mod_tags(input_records: &Vec<Record>, output_path: &str) -> anyhow::Re
|
|
|
},
|
|
},
|
|
|
TagValue::IntArray(arr) => tags.push_array(&name, arr.raw()),
|
|
TagValue::IntArray(arr) => tags.push_array(&name, arr.raw()),
|
|
|
TagValue::FloatArray(arr) => tags.push_array(&name, arr.raw()),
|
|
TagValue::FloatArray(arr) => tags.push_array(&name, arr.raw()),
|
|
|
- }
|
|
|
|
|
- })
|
|
|
|
|
|
|
+ })
|
|
|
}
|
|
}
|
|
|
// if let Some(mm) = mm_hm.get(&record_name) {
|
|
// if let Some(mm) = mm_hm.get(&record_name) {
|
|
|
// tags.push_string(b"MM", mm)
|
|
// tags.push_string(b"MM", mm)
|
|
@@ -181,102 +181,108 @@ pub fn create_bam_leave_two_out(input_path: &str) -> anyhow::Result<Vec<PathBuf>
|
|
|
Ok(output_paths)
|
|
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()?;
|
|
|
|
|
-// Ok(())
|
|
|
|
|
-// }
|
|
|
|
|
pub fn remap_bam(reference: &str, input_seq: &str, output_bam: &str) -> anyhow::Result<()> {
|
|
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}.amb"))?;
|
|
|
fs::remove_file(format!("{reference}.ann"))?;
|
|
fs::remove_file(format!("{reference}.ann"))?;
|
|
|
fs::remove_file(format!("{reference}.bwt"))?;
|
|
fs::remove_file(format!("{reference}.bwt"))?;
|
|
|
fs::remove_file(format!("{reference}.pac"))?;
|
|
fs::remove_file(format!("{reference}.pac"))?;
|
|
|
fs::remove_file(format!("{reference}.sa"))?;
|
|
fs::remove_file(format!("{reference}.sa"))?;
|
|
|
|
|
|
|
|
- // Log the captured stderr
|
|
|
|
|
- // println!("bwa mem stderr: {}", bwa_stderr);
|
|
|
|
|
- // println!("samtools sort stderr: {}", samtools_stderr);
|
|
|
|
|
-
|
|
|
|
|
Ok(())
|
|
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);
|
|
|
|
|
+// }
|
|
|
|
|
+//
|
|
|
|
|
+// 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<()> {
|
|
pub fn run_modkit(bam_path: &str, contig_fa: &str, pileup_path: &str) -> anyhow::Result<()> {
|
|
|
info!("Running Modkit pileup for {bam_path}");
|
|
info!("Running Modkit pileup for {bam_path}");
|
|
|
let mut cmd = Command::new("modkit")
|
|
let mut cmd = Command::new("modkit")
|