aavar computing for flu. make dais_vars.json and aavars.csv. started data_processing.rs.

Author: mandysulli

src/processes/prepare_mira_reports.rs

@@ -1,5 +1,5 @@
 #![allow(dead_code, unused_imports)]
-use crate::utils::{data_ingest::*, writing_outputs::*};
+use crate::utils::{data_ingest::*, data_processing::*, writing_outputs::*};
 use clap::Parser;
 use csv::ReaderBuilder;
 use either::Either;
@@ -131,10 +131,10 @@ pub fn prepare_mira_reports_process(args: ReportsArgs) -> Result<(), Box<dyn Err
         }));
 
     //Read in DAIS-ribosome data
-    let dais_ins_data = dais_insertion_data_collection(&args.irma_path);
-    let dais_del_data = dias_deletion_data_collection(&args.irma_path);
-    let dais_seq_data = dias_sequence_data_collection(&args.irma_path);
-    let dais_ref_data = dais_ref_seq_data_collection(&args.workdir_path, "flu");
+    let dais_ins_data = dais_insertion_data_collection(&args.irma_path)?;
+    let dais_del_data = dias_deletion_data_collection(&args.irma_path)?;
+    let dais_seq_data = dias_sequence_data_collection(&args.irma_path)?;
+    let dais_ref_data = dais_ref_seq_data_collection(&args.workdir_path, "flu")?;
 
     //TODO: remove these at end
     //println!("{samplesheet:?}");
@@ -149,6 +149,12 @@ pub fn prepare_mira_reports_process(args: ReportsArgs) -> Result<(), Box<dyn Err
     //println!("dais ref data: {dais_ref_data:#?}");
     //println!("ref length data: {ref_lengths:#?}");
 
+    //Calculate AA variants for aavars.csv and dais_vars.json
+    //TODO: circle back for the rsv and sc2 situations
+    let dais_vars_data = compute_dais_variants(&dais_ref_data, &dais_seq_data)?;
+
+    println!("{dais_vars_data:?}");
+
     /////////////// Write the structs to JSON files and CSV files ///////////////
     // Writing out coverage data
     let coverage_struct_values = vec![
@@ -227,7 +233,7 @@ pub fn prepare_mira_reports_process(args: ReportsArgs) -> Result<(), Box<dyn Err
         &vtype_columns,
     )?;
 
-    // Writing out allele csv and josn file
+    // Writing out allele csv and json file
     let allele_struct_values = vec![
         "Sample",
         "Upstream_Position",
@@ -301,6 +307,7 @@ pub fn prepare_mira_reports_process(args: ReportsArgs) -> Result<(), Box<dyn Err
         &indels_struct_values,
     )?;
 
+    // Write out ref_data.json
     write_ref_data_json(
         "/home/xpa3/mira-oxide/test/ref_data.json",
         &ref_lengths,
@@ -309,6 +316,29 @@ pub fn prepare_mira_reports_process(args: ReportsArgs) -> Result<(), Box<dyn Err
         &segcolor,
     )?;
 
+    // write out the dais_vars.json and the {runid}_aavars.csv
+    let aavars_columns = vec![
+        "sample_id",
+        "reference_id",
+        "protein",
+        "aa_variant_count",
+        "aa_variants",
+    ];
+    write_structs_to_split_json_file(
+        "/home/xpa3/mira-oxide/test/dais_vars.json",
+        &dais_vars_data,
+        &aavars_columns,
+        &aavars_columns,
+    )?;
+
+    write_structs_to_csv_file(
+        &format!("/home/xpa3/mira-oxide/test/{}_aavars.csv", &args.runid),
+        &dais_vars_data,
+        &aavars_columns,
+        &aavars_columns,
+    )?;
+
+    /////////////// Write the structs to parquet files if flag invoked ///////////////
     // Write fields to parq if flag given
     write_reads_to_parquet(&read_data, "/home/xpa3/mira-oxide/test/read_data.parquet")?;
 

src/utils/data_ingest.rs

@@ -197,7 +197,7 @@ pub struct DeletionsData {
     pub del_cds_length: Option<String>,
 }
 
-/// Sequence Data
+/// Dais Sequence Data
 #[derive(Serialize, Deserialize, Debug)]
 pub struct DaisSeqData {
     #[serde(rename = "ID")]

src/utils/data_processing.rs

@@ -0,0 +1,90 @@
+use crate::utils::data_ingest;
+use either::Either;
+use glob::glob;
+use serde::{self, Deserialize, Serialize, de::DeserializeOwned};
+use std::{collections::HashSet, error::Error};
+
+use super::data_ingest::DaisSeqData;
+
+/// Dais Variants
+#[derive(Serialize, Deserialize, Debug)]
+pub struct DaisVarsData {
+    pub sample_id: Option<String>,
+    pub reference_id: String,
+    pub protein: String,
+    pub aa_variant_count: i32,
+    pub aa_variants: String,
+}
+
+// Function to append a new string with a comma
+pub fn append_with_comma(base: &mut String, new_entry: &str) {
+    if !base.is_empty() {
+        base.push(',');
+        base.push_str(new_entry);
+    } else {
+        base.push_str(new_entry);
+    }
+}
+
+pub fn compute_dais_variants(
+    ref_seqs_data: &Vec<DaisSeqData>,
+    sample_seqs_data: &Vec<DaisSeqData>,
+) -> Result<Vec<DaisVarsData>, Box<dyn Error>> {
+    let mut dais_vars_data: Vec<DaisVarsData> = Vec::new();
+
+    // Compute aa variants
+    for sample_entry in sample_seqs_data {
+        for ref_entry in ref_seqs_data {
+            if sample_entry.reference == ref_entry.reference
+                && sample_entry.protein == ref_entry.protein
+            {
+                let sample_aa_seq = sample_entry.aa_aln.clone();
+                let ref_aa_seq = ref_entry.aa_aln.clone();
+                let mut var_aa_count = 0;
+
+                // Use .chars() to iterate over aa
+                let mut aa_vars = String::new();
+                for (index, (sample_aa, ref_aa)) in
+                    sample_aa_seq.chars().zip(ref_aa_seq.chars()).enumerate()
+                {
+                    if sample_aa != ref_aa {
+                        let pos = index + 1;
+                        var_aa_count += 1;
+                        let variant = format!("{ref_aa}{pos}{sample_aa}");
+                        append_with_comma(&mut aa_vars, &variant);
+                    }
+                }
+
+                let dais_vars_entry = DaisVarsData {
+                    sample_id: sample_entry.sample_id.clone(),
+                    reference_id: sample_entry.reference.clone(),
+                    protein: sample_entry.protein.clone(),
+                    aa_variant_count: var_aa_count,
+                    aa_variants: aa_vars,
+                };
+                dais_vars_data.push(dais_vars_entry);
+            }
+        }
+    }
+
+    // Sort by protein, sample_id, and aa_variant_count
+    dais_vars_data.sort_by(|a, b| {
+        a.protein
+            .cmp(&b.protein)
+            .then(a.sample_id.cmp(&b.sample_id))
+            .then(a.aa_variant_count.cmp(&b.aa_variant_count))
+    });
+
+    // Remove duplicates based on sample_id and protein, keeping the first entry
+    let mut unique_entries = Vec::new();
+    let mut seen = HashSet::new();
+
+    for entry in dais_vars_data {
+        let key = (entry.sample_id.clone(), entry.protein.clone());
+        if seen.insert(key) {
+            unique_entries.push(entry);
+        }
+    }
+
+    Ok(unique_entries)
+}

src/utils/mod.rs

@@ -1,2 +1,3 @@
 pub mod data_ingest;
+pub mod data_processing;
 pub mod writing_outputs;

src/utils/writing_outputs.rs

@@ -11,8 +11,6 @@ use std::io::Write;
 use std::sync::Arc;
 use std::{error::Error, fs::File};
 
-use super::data_ingest::CoverageData;
-
 /////////////// Functions to write to json and csv files ///////////////
 /// Function to serialize a vector of structs into split-oriented JSON with precision and indexing
 pub fn write_structs_to_split_json_file<T: Serialize>(