diff --git a/Cargo.toml b/Cargo.toml
index fd9f9b9..90dccdd 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -12,7 +12,9 @@ bio = "1.6.0"
 clap = "4.5.4"
 multimap = "0.9.1"
 rayon = "1.5.0"
+path-absolutize = "3.0.14"
 
 [dev-dependencies]
 env_logger = "0.11.3"
-tempfile = "3.2"
\ No newline at end of file
+tempfile = "3.2"
+assert_cmd = "2.0"
diff --git a/assets/metajunction_m6A.png b/assets/metajunction_m6A.png
index 846a35c..a45b2c4 100644
Binary files a/assets/metajunction_m6A.png and b/assets/metajunction_m6A.png differ
diff --git a/assets/metatranscript_m6A.png b/assets/metatranscript_m6A.png
index da211db..8b43631 100644
Binary files a/assets/metatranscript_m6A.png and b/assets/metatranscript_m6A.png differ
diff --git a/scripts/R2_plotMetaCodon.R b/scripts/R2_plotMetaCodon.R
index 8fdb3c2..1e8a1fb 100644
--- a/scripts/R2_plotMetaCodon.R
+++ b/scripts/R2_plotMetaCodon.R
@@ -5,6 +5,9 @@ suppressMessages({
   library(binom)
 })
 
+cat("Arguments received:\n")
+cat(paste(commandArgs(trailingOnly = TRUE), collapse = " "), "\n")
+
 help_message <- function() {
   cat("\nUsage: Rscript script.R '/path/to/annotated.bed' '/path/to/output.png' '<probability field>' '<cutoff>' '<upper/lower>' [-c 'loess'/'binom'] [-o '/path/to/output.tsv'] [-l] (-s | -e)\n")
   cat("Options:\n")
@@ -17,55 +20,65 @@ help_message <- function() {
 }
 
 args <- commandArgs(trailingOnly = TRUE)
-options <- list(ci_method = "loess", output_path = NULL, add_labels = FALSE, plot_from = NULL)
+required_args <- args[1:5]
+optional_args <- args[6:length(args)]
 
-# help message 
-if ("-h" %in% args || length(args) == 0) {
+if (length(required_args) != 5) {
+  cat("Error: Incorrect number of required arguments.\n")
+  cat("Received arguments:", paste(args, collapse = " "), "\n")
   help_message()
-  quit(status = 0)
+  quit(status = 1)
 }
 
-remaining_args <- list()
-plot_from_specified <- FALSE
+input_file <- required_args[1]
+output_file <- required_args[2]
+col_name <- required_args[3]
+cutoff <- as.numeric(required_args[4])
+direction <- required_args[5]
+
+options <- list(ci_method = "loess", output_path = NULL, add_labels = FALSE, plot_from = NULL)
 
-for (arg in args) {
-  if (arg %in% c("-s", "-e")) {
-    if (!is.null(options$plot_from)) {
+# parse optional arguments
+i <- 1
+while (i <= length(optional_args)) {
+  flag <- optional_args[i]
+  
+  if (flag == "-c" && i + 1 <= length(optional_args)) {
+    options$ci_method <- optional_args[i + 1]
+    i <- i + 2
+  } else if (flag == "-o" && i + 1 <= length(optional_args)) {
+    options$output_path <- optional_args[i + 1]
+    i <- i + 2
+  } else if (flag == "-l") {
+    options$add_labels <- TRUE
+    i <- i + 1
+  } else if (flag == "-s") {
+    if (is.null(options$plot_from)) {
+      options$plot_from <- "abs_cds_start"
+    } else {
       stop("Error: Both -s and -e flags cannot be used simultaneously.")
     }
-    options$plot_from <- ifelse(arg == "-s", "abs_cds_start", "abs_cds_end")
-    plot_from_specified <- TRUE
-  } else if (arg %in% c("-c", "-o")) {
-    next_index <- which(args == arg) + 1
-    if (next_index > length(args)) {
-      stop(paste("Error: The", arg, "flag requires a following argument."))
+    i <- i + 1
+  } else if (flag == "-e") {
+    if (is.null(options$plot_from)) {
+      options$plot_from <- "abs_cds_end"
+    } else {
+      stop("Error: Both -s and -e flags cannot be used simultaneously.")
     }
-    options[[sub("-", "", arg)]] <- args[next_index]
-  } else if (arg == "-l") {
-    options$add_labels <- TRUE
+    i <- i + 1
   } else {
-    remaining_args <- c(remaining_args, arg)
+    warning(paste("Unrecognized flag:", flag))
+    i <- i + 1
   }
 }
 
-# check for codon flag 
-if (!plot_from_specified) {
+# check for start or end flag 
+if (is.null(options$plot_from)) {
+  cat("Error: Either -s or -e flag must be specified.\n")
   help_message()
-  stop("Error: Either -s or -e flag must be specified.")
+  quit(status = 1)
 }
 
-# check remaining positional arguments
-if (length(remaining_args) != 5) {
-  help_message()
-  stop("Error: Incorrect number of positional arguments.")
-}
-
-input_file <- as.character(remaining_args[1])
-output_file <- as.character(remaining_args[2])
-col_name <- as.character(remaining_args[3])
-cutoff <- as.numeric(remaining_args[4])
-direction <- as.character(remaining_args[5])
-
 if (!file.exists(input_file)) {
   stop("Input file does not exist")
 }
@@ -93,11 +106,11 @@ filter_calls <- function(file, col, cutoff, direction) {
 
 # calculate ratio of significant sites 
 compute_ratio <- function(calls, interval) {
-  # Filter calls to include only those within the specified range
+
+
   calls <- calls %>%
     filter(.data[[interval]] >= -100 & .data[[interval]] <= 100)
 
-  # Calculate the count of 'sig' and 'ns' at each position
   out_ratio <- calls %>%
     group_by(.data[[interval]], filter) %>%
     summarise(n = n(), .groups = 'drop') %>%
@@ -107,7 +120,7 @@ compute_ratio <- function(calls, interval) {
            ratio = sig / (sig + ns + 1e-9)) %>%
     rename(interval = 1)
 
-  # Calculate confidence intervals if the binom method is chosen
+
   if (options$ci_method == "binom") {
     conf_int <- binom.confint(out_ratio$sig, out_ratio$sig + out_ratio$ns, methods = "wilson")
     out_ratio <- mutate(out_ratio, lower = conf_int$lower, upper = conf_int$upper)
@@ -157,6 +170,8 @@ calls <- filter_calls(input_file, col_name, cutoff, direction)
 out_ratio <- compute_ratio(calls, options$plot_from)
 
 # optionally, write the data shown in the plot to a file 
+print(options$output_path)
+
 if (!is.null(options$output_path)) {
   write_tsv(out_ratio, options$output_path)
 }
diff --git a/scripts/manuscript_demo/process_data_for_manuscript.sh b/scripts/manuscript_demo/process_data_for_manuscript.sh
new file mode 100644
index 0000000..eed0e8a
--- /dev/null
+++ b/scripts/manuscript_demo/process_data_for_manuscript.sh
@@ -0,0 +1,78 @@
+#!/bin/bash 
+
+# using minimap v2.1 and samtools v1.19
+export PATH="$PATH:/g/data/lf10/as7425/apps/minimap2/"
+module load samtools
+
+# map hela DRS reads to GRCh38 cDNA transcriptome
+export wd="/g/data/lf10/as7425/R2DTool_demo"; mkdir -p ${wd} 2>/dev/null
+export reads="/g/data/lf10/as7425/2023_mrna-biogenesis-maps/analysis/2024-03-28_ASR012_HeLa-total-RNA004-rep1/ASR012_HeLa-total-RNA004-rep1_dorado_polyA_m6A_all_reads.fastq"
+export genome="/g/data/lf10/as7425/genomes/human_genome/ensembl_release_110/Homo_sapiens.GRCh38.cdna.all.fa"
+
+minimap2 -ax map-ont -y -k 14 -N 20 ${genome} ${reads} -t 104 | samtools view -bh -@ 104 > ${wd}/aligned_reads.bam
+
+# filter for plus strand primary alignments with nonzero mapq
+samtools view -b -F 2308 -@ 104 -q 1 ${wd}/aligned_reads.bam | samtools sort > ${wd}/primary.bam
+samtools view -@ 104 -c ${wd}/primary.bam
+samtools index ${wd}/primary.bam
+
+# pileup the modification calls using modkit 
+export PATH="$PATH:/g/data/lf10/as7425/apps/dist"
+modkit pileup -t 104 ${wd}/primary.bam "${wd}/pileup.bed" --log-filepath "${wd}/pileup.bed.log"
+
+# convert to table with header for R2Dtool 
+# filter for coverage >=10
+printf "transcript\tstart\tend\tmod\tcov\tstrand\tstart\tend\tcolor\tX1\tX2\tcov\tfrac_mod\tn_mod\tn_canonical\tn_other_mod\tn_delete\tn_fail\tn_diff\tn_no_call\n" > ${wd}/R2D_input.bed
+awk 'BEGIN {FS=OFS="\t"} {gsub(",", OFS, $0); $9="."; print}' "${wd}/pileup.bed" | tr " " "\t" | awk '($12 > 9)' >> ${wd}/R2D_input.bed
+
+# if not filtering for coverage 
+# printf "transcript\tstart\tend\tmod\tcov\tstrand\tstart\tend\tcolor\tX1\tX2\tcov\tfrac_mod\tn_mod\tn_canonical\tn_other_mod\tn_delete\tn_fail\tn_diff\tn_no_call\n" > ${wd}/R2D_input.bed
+# awk 'BEGIN {FS=OFS="\t"} {gsub(",", OFS, $0); $9="."; print}' "${wd}/pileup.bed" | tr " " "\t" >> ${wd}/R2D_input.bed
+
+
+# verify that all sites mapped to 'A' nucleotide as expected
+export wd="/g/data/lf10/as7425/R2DTool_demo"; mkdir -p ${wd} 2>/dev/null
+export genome="/g/data/lf10/as7425/genomes/human_genome/ensembl_release_110/Homo_sapiens.GRCh38.cdna.all.fa"
+tail -n +2 ${wd}/R2D_input.bed > ${wd}/R2D_input_noheader.bed
+bedtools getfasta -s -fi ${genome} -bed ${wd}/R2D_input_noheader.bed | tail -n +2 | awk 'NR%2==1' | sort | uniq -c | sort -nr > ${wd}/input_sequence_context.txt
+cat ${wd}/input_sequence_context.txt
+
+# 336160 A
+#     782 T
+#     333 G
+#     300 C
+
+
+##################################################
+
+# run R2Dtool 
+
+export wd="/g/data/lf10/as7425/R2DTool_demo";
+export sites="${wd}/R2D_input.bed"
+export anno="/g/data/lf10/as7425/genomes/human_genome/ensembl_release_110/Homo_sapiens.GRCh38.110.chr_patch_hapl_scaff.gtf"
+
+# R2Dtool rust implementation 
+cd ~/R2Dtool && rm -rf target; cargo build --release 
+export PATH="${PATH}:/home/150/as7425/R2Dtool/target/release/"
+
+time r2d annotate -i "${sites}" -g ${anno} -H > "${wd}/methylation_calls_annotated.bed"
+time r2d liftover -i "${sites}" -g ${anno} -H > "${wd}/methylation_calls_annotated_liftover.bed"
+
+##################################################
+
+# check sequence context in the liftover 
+
+export genome="/g/data/lf10/as7425/genomes/human_genome/ensembl_release_110/Homo_sapiens.GRCh38.dna_sm.toplevel.fa"
+tail -n +2 "${wd}/methylation_calls_annotated_liftover.bed" | cut -f1-6 | awk -F "\t" '($6 == "-")' > "${wd}/methylation_calls_annotated_liftover_noheader.bed" 
+bedtools getfasta -s -fi ${genome} -bed "${wd}/methylation_calls_annotated_liftover_noheader.bed" | tail -n +2 | awk 'NR%2==1' | sort | uniq -c | sort -nr > ${wd}/liftover_sequence_context.txt
+cat ${wd}/liftover_sequence_context.txt
+
+##################################################
+
+# make R2Dtool plots 
+
+time Rscript ~/R2Dtool/scripts/R2_plotMetaTranscript.R "${wd}/methylation_calls_annotated.bed" "${wd}/metatranscript_out_Rust.svg" "frac_mod" "9" "upper" -c loess -l -o "${wd}/transcript_data_Rust.tsv"
+time Rscript ~/R2Dtool/scripts/R2_plotMetaJunction.R "${wd}/methylation_calls_annotated.bed" "${wd}/metajunction_out_binom_rust.svg" "frac_mod" "9" "upper" -c loess -o "$wd/table"
+
+
+##################################################
\ No newline at end of file
diff --git a/src/main.rs b/src/main.rs
index 861e071..0957453 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -2,16 +2,14 @@
 use clap::{Arg,Command};
 use std::process::Command as ProcessCommand;
 use std::env;
+use std::path::{Path,PathBuf};
+use path_absolutize::Absolutize;
 
 // modules
 pub mod parse_annotation;
 pub mod annotate;
 pub mod liftover;
 pub mod parse_gtf;
-use std::path::Path;
-
-#[cfg(test)]
-mod tests;
 
 const RELATIVE_SCRIPT_PATH: &str = "../../scripts/";
 
@@ -376,7 +374,7 @@ fn main() {
         println!("Arguments being passed to R2_plotMetaTranscript.R:");
         println!("Rscript R2_plotMetaTranscript.R {}", args.join(" "));
 
-        match generate_r_plots("R2_plotMetaTranscript.R", &args) {
+        match generate_r_plots("R2_plotMetaTranscript.R", &args, None) {
             Ok(_) => println!("PlotMetaTranscript generated successfully."),
             Err(e) => {
                 eprintln!("Error: {}", e);
@@ -406,7 +404,7 @@ fn main() {
         println!("Arguments being passed to R2_plotMetaJunction.R:");
         println!("Rscript R2_plotMetaJunction.R {}", args.join(" "));
 
-        match generate_r_plots("R2_plotMetaJunction.R", &args) {
+        match generate_r_plots("R2_plotMetaJunction.R", &args, None) {
             Ok(_) => println!("PlotMetaJunction generated successfully."),
             Err(e) => {
                 eprintln!("Error: {}", e);
@@ -416,37 +414,49 @@ fn main() {
         }
     }
     
-// plotMetaCodon 
+// plotMetaCodon
 if let Some(matches) = matches.subcommand_matches("plotMetaCodon") {
-        let codon_flag = if matches.get_flag("start_codon") { "-s" } else { "-e" };
+    let codon_flag = if matches.get_flag("start_codon") { "-s" } else { "-e" };
 
-        let mut args: Vec<String> = vec![
-            matches.get_one::<String>("input_file").unwrap().clone(),
-            matches.get_one::<String>("output_file").unwrap().clone(),
-            matches.get_one::<String>("field_name").unwrap().clone(),
-            matches.get_one::<String>("cutoff_value").unwrap().clone(),
-            matches.get_one::<String>("cutoff_type").unwrap().clone(),
-            codon_flag.to_string(),
-        ];
+    let current_dir = std::env::current_dir().expect("Failed to get current directory");
+    let input_file = current_dir.join(matches.get_one::<String>("input_file").unwrap())
+        .canonicalize()
+        .expect("Failed to canonicalize input file path");
+    let output_file_path = PathBuf::from(matches.get_one::<String>("output_file").unwrap());
+    let output_file = output_file_path.absolutize()
+        .expect("Failed to absolutize output file path");
 
-        if let Some(method) = matches.get_one::<String>("confidence_method") {
-            args.extend_from_slice(&["-c".to_string(), method.clone()]);
-        }
+    let mut args: Vec<String> = vec![
+        input_file.to_string_lossy().into_owned(),
+        output_file.to_string_lossy().into_owned(),
+        matches.get_one::<String>("field_name").unwrap().clone(),
+        matches.get_one::<String>("cutoff_value").unwrap().clone(),
+        matches.get_one::<String>("cutoff_type").unwrap().clone(),
+    ];
 
-        if let Some(table) = matches.get_one::<String>("save_table") {
-            args.extend_from_slice(&["-o".to_string(), table.clone()]);
-        }
+    args.push(codon_flag.to_string());
 
-        println!("Arguments being passed to R2_plotMetaCodon.R:");
-        println!("Rscript R2_plotMetaCodon.R {}", args.join(" "));
-    
-        match generate_r_plots("R2_plotMetaCodon.R", &args) {
-            Ok(_) => println!("PlotMetaCodon generated successfully."),
-            Err(e) => {
-                eprintln!("Error: {}", e);
-                eprintln!("PlotMetaCodon generation failed. Please check the error messages above.");
-                std::process::exit(1);
-            }
+    if let Some(method) = matches.get_one::<String>("confidence_method") {
+        args.extend_from_slice(&["-c".to_string(), method.clone()]);
+    }
+
+    if let Some(table) = matches.get_one::<String>("save_table") {
+        let save_table_path_buf = PathBuf::from(table);
+        let save_table_path = save_table_path_buf.absolutize()
+            .expect("Failed to absolutize save table path");    
+        args.extend_from_slice(&["-o".to_string(), save_table_path.to_string_lossy().into_owned()]);
+    }
+
+    println!("Arguments being passed to R2_plotMetaCodon.R:");
+    println!("Rscript R2_plotMetaCodon.R {}", args.join(" "));
+
+    match generate_r_plots("R2_plotMetaCodon.R", &args, None) {
+        Ok(_) => println!("PlotMetaCodon generated successfully."),
+        Err(e) => {
+            eprintln!("Error: {}", e);
+            eprintln!("PlotMetaCodon generation failed. Please check the error messages above.");
+            std::process::exit(1);
         }
     }
 }
+}
diff --git a/src/tests/mod.rs b/src/tests/mod.rs
new file mode 100644
index 0000000..28f7bf0
--- /dev/null
+++ b/src/tests/mod.rs
@@ -0,0 +1 @@
+mod plot_integration_tests;
diff --git a/src/tests/plot_integration_tests.rs b/src/tests/plot_integration_tests.rs
index a0a3918..501a4c7 100644
--- a/src/tests/plot_integration_tests.rs
+++ b/src/tests/plot_integration_tests.rs
@@ -2,11 +2,30 @@ use std::path::{PathBuf, Path};
 use tempfile::NamedTempFile;
 use crate::{annotate, generate_r_plots};
 use std::env;
+use std::fs;
 
 fn get_project_root() -> PathBuf {
     Path::new(&env!("CARGO_MANIFEST_DIR")).to_path_buf()
 }
 
+fn create_output_directory() -> std::io::Result<PathBuf> {
+    let output_dir = get_project_root().join("test_outputs");
+    fs::create_dir_all(&output_dir)?;
+    Ok(output_dir)
+}
+
+fn generate_output_filename(plot_type: &str, cutoff: &str, direction: &str, labels: bool, ci_method: &str, codon_type: Option<&str>) -> String {
+    let mut filename = format!("{}_cutoff{}_{}_{}", plot_type, cutoff, direction, ci_method);
+    if labels {
+        filename.push_str("_with_labels");
+    }
+    if let Some(codon) = codon_type {
+        filename.push_str(&format!("_{}_codon", codon));
+    }
+    filename.push_str(".png");
+    filename
+}
+
 fn run_annotate(input: &Path, gtf: &Path, output: &Path) -> Result<(), Box<dyn std::error::Error>> {
     let matches = clap::Command::new("test")
         .arg(clap::Arg::new("gtf").short('g').long("gtf").required(true))
@@ -22,13 +41,11 @@ fn run_annotate(input: &Path, gtf: &Path, output: &Path) -> Result<(), Box<dyn s
     annotate::run_annotate(&matches, true, false)
 }
 
-#[test]
-fn test_plotting_functions() -> Result<(), Box<dyn std::error::Error>> {
+fn setup_test() -> Result<(PathBuf, PathBuf, NamedTempFile), Box<dyn std::error::Error>> {
     let project_root = get_project_root();
     let script_dir = project_root.join("scripts");
     let test_dir = project_root.join("test");
 
-    // Run annotate and save to a temp file
     let annotate_output = NamedTempFile::new()?;
     run_annotate(
         &test_dir.join("m6A_isoform_sites_GRCh38_subset.bed"),
@@ -36,69 +53,229 @@ fn test_plotting_functions() -> Result<(), Box<dyn std::error::Error>> {
         annotate_output.path()
     )?;
 
-    // Test plotMetaTranscript with different configurations
-    let meta_transcript_plots = [
-        ("metatranscript_m6A_upper.png", "upper", true),
-        ("metatranscript_m6A_upper_no_labels.png", "upper", false),
-        ("metatranscript_m6A_lower.png", "lower", true),
+    Ok((script_dir, test_dir, annotate_output))
+}
+
+fn run_plot_test(script_name: &str, args: Vec<String>, script_dir: &Path) -> Result<(), Box<dyn std::error::Error>> {
+    generate_r_plots(script_name, &args, Some(script_dir))?;
+    let output_path = Path::new(&args[1]);
+    assert!(output_path.exists());
+    Ok(())
+}
+
+// metatranscript tests 
+#[test]
+fn test_plot_meta_transcript_cutoff10_upper_loess_with_labels() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metatranscript", "10", "upper", true, "loess", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "10".to_string(),
+        "upper".to_string(),
+        "-l".to_string(),
     ];
+    run_plot_test("R2_plotMetaTranscript.R", args, &script_dir)
+}
 
-    for (_, cutoff_type, show_labels) in &meta_transcript_plots {
-        let output_path = NamedTempFile::new()?;
-        let mut args = vec![
-            annotate_output.path().to_str().unwrap().to_string(),
-            output_path.path().to_str().unwrap().to_string(),
-            "fraction_modified".to_string(),
-            "10".to_string(),
-            cutoff_type.to_string(),
-        ];
-        if *show_labels {
-            args.push("-l".to_string());
-        }
-        generate_r_plots("R2_plotMetaTranscript.R", &args, Some(&script_dir))?;
-        assert!(output_path.path().exists());
-    }
+#[test]
+fn test_plot_meta_transcript_cutoff0_lower_binom_without_labels() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metatranscript", "0", "lower", false, "binom", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "0".to_string(),
+        "lower".to_string(),
+        "-c".to_string(),
+        "binom".to_string(),
+    ];
+    run_plot_test("R2_plotMetaTranscript.R", args, &script_dir)
+}
 
-    // Test plotMetaJunction with different configurations
-    let meta_junction_plots = [
-        ("metajunction_m6A_upper.png", "upper"),
-        ("metajunction_m6A_lower.png", "lower"),
+#[test]
+fn test_plot_meta_transcript_cutoff100_upper_loess_without_labels() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metatranscript", "100", "upper", false, "loess", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "100".to_string(),
+        "upper".to_string(),
     ];
+    run_plot_test("R2_plotMetaTranscript.R", args, &script_dir)
+}
 
-    for (_, cutoff_type) in &meta_junction_plots {
-        let output_path = NamedTempFile::new()?;
-        let args = vec![
-            annotate_output.path().to_str().unwrap().to_string(),
-            output_path.path().to_str().unwrap().to_string(),
-            "fraction_modified".to_string(),
-            "10".to_string(),
-            cutoff_type.to_string(),
-        ];
-        generate_r_plots("R2_plotMetaJunction.R", &args, Some(&script_dir))?;
-        assert!(output_path.path().exists());
-    }
+// metajunction tests 
+#[test]
+fn test_plot_meta_junction_cutoff10_upper_loess() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metajunction", "10", "upper", false, "loess", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "10".to_string(),
+        "upper".to_string(),
+    ];
+    run_plot_test("R2_plotMetaJunction.R", args, &script_dir)
+}
 
-    // Test plotMetaCodon for start and stop codons
-    let meta_codon_plots = [
-        ("metacodon_start_upper.png", "upper", "-s"),
-        ("metacodon_start_lower.png", "lower", "-s"),
-        ("metacodon_stop_upper.png", "upper", "-e"),
-        ("metacodon_stop_lower.png", "lower", "-e"),
+#[test]
+fn test_plot_meta_junction_cutoff0_lower_binom() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metajunction", "0", "lower", false, "binom", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "0".to_string(),
+        "lower".to_string(),
+        "-c".to_string(),
+        "binom".to_string(),
     ];
+    run_plot_test("R2_plotMetaJunction.R", args, &script_dir)
+}
 
-    for (_, cutoff_type, codon_type) in &meta_codon_plots {
-        let output_path = NamedTempFile::new()?;
-        let args = vec![
-            annotate_output.path().to_str().unwrap().to_string(),
-            output_path.path().to_str().unwrap().to_string(),
-            "fraction_modified".to_string(),
-            "10".to_string(),
-            cutoff_type.to_string(),
-            codon_type.to_string(),
-        ];
-        generate_r_plots("R2_plotMetaCodon.R", &args, Some(&script_dir))?;
-        assert!(output_path.path().exists());
-    }
+// metacodon tests
+#[test]
+fn test_plot_meta_codon_start_cutoff10_upper_loess() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metacodon", "10", "upper", false, "loess", Some("start"));
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "10".to_string(),
+        "upper".to_string(),
+        "-s".to_string(),
+    ];
+    run_plot_test("R2_plotMetaCodon.R", args, &script_dir)
+}
+
+#[test]
+fn test_plot_meta_codon_stop_cutoff0_lower_binom() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metacodon", "0", "lower", false, "binom", Some("stop"));
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "0".to_string(),
+        "lower".to_string(),
+        "-e".to_string(),
+        "-c".to_string(),
+        "binom".to_string(),
+    ];
+    run_plot_test("R2_plotMetaCodon.R", args, &script_dir)
+}
+
+#[test]
+fn test_plot_meta_codon_start_cutoff100_upper_loess() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metacodon", "100", "upper", false, "loess", Some("start"));
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "100".to_string(),
+        "upper".to_string(),
+        "-s".to_string(),
+    ];
+    run_plot_test("R2_plotMetaCodon.R", args, &script_dir)
+}
+
+#[test]
+fn test_plot_meta_codon_start_cutoff100_lower_loess() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metacodon", "100", "lower", false, "loess", Some("start"));
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "100".to_string(),
+        "lower".to_string(),
+        "-s".to_string(),
+    ];
+    run_plot_test("R2_plotMetaCodon.R", args, &script_dir)
 
+// helper function for saving plot tables 
+fn run_plot_test_with_table(script_name: &str, args: Vec<String>, script_dir: &Path) -> Result<(), Box<dyn std::error::Error>> {
+    let mut table_args = args.clone();
+    let table_path = Path::new(&args[1]).with_extension("tsv");
+    table_args.extend_from_slice(&["-o".to_string(), table_path.to_str().unwrap().to_string()]);
+    generate_r_plots(script_name, &table_args, Some(script_dir))?;
+    assert!(Path::new(&args[1]).exists());
+    assert!(table_path.exists());
     Ok(())
+}
+
+#[test]
+fn test_plot_meta_transcript_with_table() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metatranscript", "10", "upper", false, "loess", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "10".to_string(),
+        "upper".to_string(),
+    ];
+    run_plot_test_with_table("R2_plotMetaTranscript.R", args, &script_dir)
+}
+
+#[test]
+fn test_plot_meta_junction_with_table() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metajunction", "10", "upper", false, "loess", None);
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "10".to_string(),
+        "upper".to_string(),
+    ];
+    run_plot_test_with_table("R2_plotMetaJunction.R", args, &script_dir)
+}
+
+#[test]
+fn test_plot_meta_codon_with_table() -> Result<(), Box<dyn std::error::Error>> {
+    let (script_dir, _, annotate_output) = setup_test()?;
+    let output_dir = create_output_directory()?;
+    let filename = generate_output_filename("metacodon", "10", "upper", false, "loess", Some("start"));
+    let output_path = output_dir.join(filename);
+    let args = vec![
+        annotate_output.path().to_str().unwrap().to_string(),
+        output_path.to_str().unwrap().to_string(),
+        "fraction_modified".to_string(),
+        "10".to_string(),
+        "upper".to_string(),
+        "-s".to_string(),
+    ];
+    run_plot_test_with_table("R2_plotMetaCodon.R", args, &script_dir)
 }
\ No newline at end of file