grumpy 1.1.3

Genetic analysis in Rust.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158

use std::time::SystemTime;

use clap::Parser;

pub mod common;
pub mod difference;
pub mod gene;
pub mod genome;
pub mod vcf;

use common::MinorType;
use difference::GeneDifference;
use difference::GenomeDifference;
use genome::mutate;
use genome::Genome;
use vcf::VCFFile;

#[derive(Parser, Debug)]
#[command(version, about, long_about = None)]
struct Args {
    /// Path to the reference genbank file
    #[arg(long)]
    reference: String,

    /// Path to the VCF file
    #[arg(long)]
    vcf: String,

    /// Target gene name. If given, pull out all variants and mutations within this gene.
    /// If omitted, pull out variants and mutations within all genes
    #[arg(long)]
    gene: Option<String>,

    /// Path to the output variable length FASTA file. If given, the sample genome will be written to this file
    #[arg(long)]
    fasta: Option<String>,
}

#[cfg(not(tarpaulin_include))]
fn main() {
    let args = Args::parse();

    let reference_path = args.reference;
    let vcf_path = args.vcf;

    let vcf_start = SystemTime::now();
    let vcf = VCFFile::new(vcf_path.to_string(), false, 3);
    let vcf_end = SystemTime::now();

    let reference_start = SystemTime::now();
    let mut reference = Genome::new(&reference_path);
    let reference_end = SystemTime::now();

    let sample_start = SystemTime::now();
    let mut sample = mutate(&reference, vcf);
    let sample_end = SystemTime::now();

    // If given a gene name, pull out the genome and gene level differences
    let target_genes;
    if let Some(target_gene) = args.gene {
        target_genes = vec![target_gene];
    } else {
        target_genes = reference.gene_names.clone();
    }

    let genome_start = SystemTime::now();
    let mut difference = GenomeDifference::new(reference.clone(), sample.clone(), MinorType::COV);
    let genome_end = SystemTime::now();

    let mut is_first = true;
    let mut gene_start = SystemTime::now();
    let mut gene_end = SystemTime::now();
    for target_gene in target_genes.iter() {
        for variant in difference.variants.iter_mut() {
            if variant.gene_name.clone().is_none()
                || (variant.gene_name.clone().is_some()
                    && variant.gene_name.clone().unwrap() != *target_gene)
            {
                continue;
            }
            println!(
                "{:?}@{:?} --> {:?}",
                variant.gene_name,
                variant.gene_position,
                variant.variant.clone()
            );
        }
        for variant in difference.minor_variants.iter_mut() {
            if variant.gene_name.clone().is_none()
                || (variant.gene_name.clone().is_some()
                    && variant.gene_name.clone().unwrap() != *target_gene)
            {
                continue;
            }
            println!(
                "{:?}@{:?} --> {:?}",
                variant.gene_name,
                variant.gene_position,
                variant.variant.clone()
            );
        }
        if is_first {
            gene_start = SystemTime::now();
        }
        for gene_name in sample.genes_with_mutations.clone().iter() {
            if gene_name != target_gene {
                continue;
            }
            println!("{}", gene_name);
            let gene_diff = GeneDifference::new(
                reference.get_gene(gene_name.clone()),
                sample.get_gene(gene_name.clone()),
                MinorType::COV,
            );
            println!(
                "{:?}",
                gene_diff
                    .mutations
                    .iter()
                    .map(|mutation| mutation.mutation.clone())
                    .collect::<Vec<String>>()
            );
            println!(
                "{:?}\n",
                gene_diff
                    .minor_mutations
                    .iter()
                    .map(|mutation| mutation.mutation.clone())
                    .collect::<Vec<String>>()
            );
        }
        gene_end = SystemTime::now();
        is_first = false;
    }
    println!("\n-----------------------------------\n");
    println!("VCF took {:?}", vcf_end.duration_since(vcf_start).unwrap());
    println!(
        "Reference took {:?}",
        reference_end.duration_since(reference_start).unwrap()
    );
    println!(
        "Sample took {:?}",
        sample_end.duration_since(sample_start).unwrap()
    );
    println!(
        "Genome diff took {:?}",
        genome_end.duration_since(genome_start).unwrap()
    );
    println!(
        "Gene diff took {:?}",
        gene_end.duration_since(gene_start).unwrap()
    );

    // If a FASTA path is provided, write the sample genome to it
    if let Some(fasta_path) = args.fasta {
        sample.write_fasta(&fasta_path);
    }
}