sciencegenome

finished writing the code and had no internet connection and hence pushing it now. Final and binary release.

Author: Gaurav Sablok

README.md

@@ -6,29 +6,68 @@
  
  ```
  ```
- ➜  graph-clusters git:(main) ✗ ./target/debug/graph-kmers -h                                                                                       
- graphkmers                                                                                                                                         
-                                                                                                                                                   
- Usage: graph-kmers <COMMAND>                                                                                                                       
-                                                                                                                                                   
- Commands:                                                                                                                                          
-  sequence  profile the similarity based on the observed kmer                                                                                      
-  help      Print this message or the help of the given subcommand(s)                                                                              
-                                                                                                                                                   
- Options:                                                                                                                                           
-  -h, --help     Print help                                                                                                                        
-  -V, --version  Print version                                                                                                                     
- ➜  graph-clusters git:(main) ✗ ./target/debug/graph-kmers sequence -h                                                                              
- profile the similarity based on the observed kmer                                                                                                  
-                                                                                                                                                   
- Usage: graph-kmers sequence <SEQUENCEPATH> <SEQUENCEKMER>                                                                                          
-                                                                                                                                                   
- Arguments:                                                                                                                                         
-  <SEQUENCEPATH>  provide the path to sequence file                                                                                                
-  <SEQUENCEKMER>  provide the kmer to be profiled for the sequence similarity                                                                      
-                                                                                                                                                   
- Options:                                                                                                                                           
-  -h, --help  Print help 
+ 12:33:48 gauravsablok@genome graph-kmer main ? ./target/debug/graph-kmers -h
+ graphkmers
 
+ Usage: graph-kmers <COMMAND>
+
+ Commands:
+  sequence  identity kmer similarity index
+  filter    identity kmer filter
+  help      Print this message or the help of the given subcommand(s)
+
+ Options:
+  -h, --help     Print help
+  -V, --version  Print version
+ 12:33:53 gauravsablok@genome graph-kmer main ? ./target/debug/graph-kmers sequence -h
+ identity kmer similarity index
+
+ Usage: graph-kmers sequence <SEQUENCEPATH> <SEQUENCEKMER>
+
+ Arguments:
+  <SEQUENCEPATH>  provide the path to sequence file
+  <SEQUENCEKMER>  provide the kmer to be profiled for the sequence similarity
+
+ Options:
+  -h, --help  Print help
+ 12:33:57 gauravsablok@genome graph-kmer main ? ./target/debug/graph-kmers filter -h
+ identity kmer filter
+
+ Usage: graph-kmers filter <SEQUENCE> <KMER> <THRESHOLD>
+
+ Arguments:
+  <SEQUENCE>   provide the path to the sequence file
+  <KMER>       sequence kmer for the identity kmer
+  <THRESHOLD>  provide the threshold
+
+ Options:
+  -h, --help  Print help
+ 
+ ```
+ - to run the compiled library 
+ ```
+ 12:53:04 gauravsablok@genome graph-kmer main ? \ 
+           ./target/debug/graph-kmers filter \
+                           ./samplefile/sample.fasta 4 10
+ 
+ 12:56:27 gauravsablok@genome samplefile main ? 
+           ./target/debug/graph-kmers sequence 
+                          ./samplefile/sample.fasta 4
+
+
+ 12:55:34 gauravsablok@genome samplefile main ? head -n 4 frequencies-t*
+ ==> frequencies-tab.txt <== filtering without the threshold
+ "chr11:66478458"        "chr11:66478459"        22      44      50.0
+ "chr11:66478458"        "chr11:66478460"        12      44      27.272728
+ "chr11:66478458"        "chr11:66478461"        8       42      19.047619
+ "chr11:66478458"        "chr11:66478462"        9       41      21.95122
+
+ ==> frequencies-threshold.txt <== filtering after the threshold
+ "chr11:66478458"        "chr11:66478459"        22      44      50.0
+ "chr11:66478458"        "chr11:66478460"        12      44      27.272728
+ "chr11:66478458"        "chr11:66478461"        8       42      19.047619
+ "chr11:66478458"        "chr11:66478462"        9       41      21.95122
  ```
 
+ Gaurav Sablok
+

graph-kmers

@@ -0,0 +1,10 @@
+"chr11:66478458"	"chr11:66478459"	22	44	50.0
+"chr11:66478458"	"chr11:66478460"	12	44	27.272728
+"chr11:66478458"	"chr11:66478461"	8	42	19.047619
+"chr11:66478458"	"chr11:66478462"	9	41	21.95122
+"chr11:66478459"	"chr11:66478460"	12	44	27.272728
+"chr11:66478459"	"chr11:66478461"	8	42	19.047619
+"chr11:66478459"	"chr11:66478462"	9	41	21.95122
+"chr11:66478460"	"chr11:66478461"	8	42	19.047619
+"chr11:66478460"	"chr11:66478462"	9	41	21.95122
+"chr11:66478461"	"chr11:66478462"	11	39	28.20513

graph-kmers.tar

@@ -0,0 +1,10 @@
+"chr11:66478458"	"chr11:66478459"	22	44	50.0
+"chr11:66478458"	"chr11:66478460"	12	44	27.272728
+"chr11:66478458"	"chr11:66478461"	8	42	19.047619
+"chr11:66478458"	"chr11:66478462"	9	41	21.95122
+"chr11:66478459"	"chr11:66478460"	12	44	27.272728
+"chr11:66478459"	"chr11:66478461"	8	42	19.047619
+"chr11:66478459"	"chr11:66478462"	9	41	21.95122
+"chr11:66478460"	"chr11:66478461"	8	42	19.047619
+"chr11:66478460"	"chr11:66478462"	9	41	21.95122
+"chr11:66478461"	"chr11:66478462"	11	39	28.20513

samplefile/frequencies-tab.txt

@@ -0,0 +1,10 @@
+>chr11:66478458
+ATGGCGGACACCCAGTACATCCTGC
+>chr11:66478459
+ATGGCGGACACCCAGTACATCCTGC
+>chr11:66478460
+ATGGCGGACACCCTATAGTACTAGA
+>chr11:66478461
+ATGGCGGACACTGACTGAATGCATA
+>chr11:66478462
+ATGGCGGACATAGATGCATGCATCA

samplefile/frequencies-threshold.txt

@@ -0,0 +1,10 @@
+"chr11:66478458"	"chr11:66478459"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACCCAGTACATCCTGC"	22	44	50.0
+"chr11:66478458"	"chr11:66478460"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACCCTATAGTACTAGA"	12	44	27.272728
+"chr11:66478458"	"chr11:66478461"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACTGACTGAATGCATA"	8	42	19.047619
+"chr11:66478458"	"chr11:66478462"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACATAGATGCATGCATCA"	9	41	21.95122
+"chr11:66478459"	"chr11:66478460"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACCCTATAGTACTAGA"	12	44	27.272728
+"chr11:66478459"	"chr11:66478461"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACTGACTGAATGCATA"	8	42	19.047619
+"chr11:66478459"	"chr11:66478462"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACATAGATGCATGCATCA"	9	41	21.95122
+"chr11:66478460"	"chr11:66478461"	"ATGGCGGACACCCTATAGTACTAGA"	"ATGGCGGACACTGACTGAATGCATA"	8	42	19.047619
+"chr11:66478460"	"chr11:66478462"	"ATGGCGGACACCCTATAGTACTAGA"	"ATGGCGGACATAGATGCATGCATCA"	9	41	21.95122
+"chr11:66478461"	"chr11:66478462"	"ATGGCGGACACTGACTGAATGCATA"	"ATGGCGGACATAGATGCATGCATCA"	11	39	28.20513

samplefile/sample.fasta

@@ -0,0 +1,10 @@
+"chr11:66478458"	"chr11:66478459"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACCCAGTACATCCTGC"	22	44	50.0
+"chr11:66478458"	"chr11:66478460"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACCCTATAGTACTAGA"	12	44	27.272728
+"chr11:66478458"	"chr11:66478461"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACTGACTGAATGCATA"	8	42	19.047619
+"chr11:66478458"	"chr11:66478462"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACATAGATGCATGCATCA"	9	41	21.95122
+"chr11:66478459"	"chr11:66478460"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACCCTATAGTACTAGA"	12	44	27.272728
+"chr11:66478459"	"chr11:66478461"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACACTGACTGAATGCATA"	8	42	19.047619
+"chr11:66478459"	"chr11:66478462"	"ATGGCGGACACCCAGTACATCCTGC"	"ATGGCGGACATAGATGCATGCATCA"	9	41	21.95122
+"chr11:66478460"	"chr11:66478461"	"ATGGCGGACACCCTATAGTACTAGA"	"ATGGCGGACACTGACTGAATGCATA"	8	42	19.047619
+"chr11:66478460"	"chr11:66478462"	"ATGGCGGACACCCTATAGTACTAGA"	"ATGGCGGACATAGATGCATGCATCA"	9	41	21.95122
+"chr11:66478461"	"chr11:66478462"	"ATGGCGGACACTGACTGAATGCATA"	"ATGGCGGACATAGATGCATGCATCA"	11	39	28.20513

samplefile/sequence-clusters.fasta

@@ -9,11 +9,20 @@ pub struct CommandParse {
 
 #[derive(Subcommand, Debug)]
 pub enum Commands {
-    /// profile the similarity based on the observed kmer
+    /// identity kmer similarity index
     Sequence {
         /// provide the path to sequence file
         sequencepath: String,
         /// provide the kmer to be profiled for the sequence similarity
         sequencekmer: String,
     },
+    /// identity kmer filter
+    Filter {
+        /// provide the path to the sequence file
+        sequence: String,
+        /// sequence kmer for the identity kmer
+        kmer: String,
+        /// provide the threshold
+        threshold: String,
+    },
 }

samplefile/sequence-threshold.fasta

@@ -7,18 +7,11 @@ pub struct Genomeiter {
 
 #[derive(Debug, Clone, PartialOrd, PartialEq)]
 
-pub struct Collector {
+pub struct CollectIter {
     pub name: String,
+    pub namenext: String,
     pub id: String,
-    pub count: usize,
-}
-
-#[derive(Debug, Clone, PartialOrd, PartialEq)]
-
-pub struct ProfileKmer {
-    pub name: String,
-    pub sequence: Vec<String>,
+    pub idnext: String,
     pub count: usize,
     pub shared: usize,
-    pub ratio: usize,
 }

src/args.rs

@@ -1,10 +1,11 @@
 mod args;
 mod filestruct;
+mod simfilter;
 mod similarity;
 use crate::args::CommandParse;
 use crate::args::Commands;
+use crate::simfilter::simfilterarg;
 use crate::similarity::profilesimilarity;
-
 use clap::Parser;
 
 /*
@@ -26,5 +27,13 @@ fn main() {
                 command
             );
         }
+        Commands::Filter {
+            sequence,
+            kmer,
+            threshold,
+        } => {
+            let command = simfilterarg(sequence, kmer, threshold).unwrap();
+            println!("The filtered files have been written: {:?}", command);
+        }
     }
 }

src/filestruct.rs

@@ -0,0 +1,134 @@
+use crate::filestruct::CollectIter;
+use crate::filestruct::Genomeiter;
+use std::collections::HashSet;
+use std::error::Error;
+use std::fs::File;
+use std::io::{BufRead, BufReader, Write};
+/*
+
+ Author Gaurav Sablok
+ SLB Potsdam
+ Date: 2025-2-7
+
+  profiling the kmer on the sequences and then based on the
+  sequence kmers implementing a similarity ratio and that says
+
+  kmerratio = observer unique kmer /
+          number of the total unique kmer * 100
+*/
+
+pub fn simfilterarg(path: &str, kmer: &str, threshold: &str) -> Result<String, Box<dyn Error>> {
+    let fileopen = File::open(path).expect("file not found");
+    let fileread = BufReader::new(fileopen);
+    let mut sequencevector: Vec<String> = Vec::new();
+    let mut headervector: Vec<String> = Vec::new();
+    let thresh: f32 = threshold.parse::<f32>().unwrap();
+    for i in fileread.lines() {
+        let line = i.expect("file not found");
+        if line.starts_with(">") {
+            headervector.push(line.replace(">", ""));
+        } else if !line.starts_with(">") {
+            sequencevector.push(line);
+        }
+    }
+    let mut combinedinfo: Vec<Genomeiter> = Vec::new();
+    for i in 0..headervector.len() {
+        combinedinfo.push(Genomeiter {
+            header: headervector[i].clone(),
+            sequence: sequencevector[i].clone(),
+        })
+    }
+
+    let mut seqbtreemap: Vec<(String, (String, Vec<String>))> = Vec::new();
+    for i in combinedinfo.iter() {
+        let windowkmer: Vec<_> = i
+            .sequence
+            .chars()
+            .map(String::from)
+            .collect::<Vec<_>>()
+            .windows(kmer.parse::<usize>().unwrap())
+            .map(|x| x.join("").to_string())
+            .collect::<Vec<_>>();
+        let sequencehash: Vec<String> = windowkmer
+            .into_iter()
+            .collect::<HashSet<String>>()
+            .into_iter()
+            .collect::<Vec<_>>();
+        seqbtreemap.push((i.header.clone(), (i.sequence.clone(), sequencehash)));
+    }
+
+    let mut newbase: Vec<Vec<CollectIter>> = Vec::new();
+    for i in 0..seqbtreemap.len() - 1 {
+        let mut shared: Vec<CollectIter> = Vec::new();
+        let vec = seqbtreemap[i].clone();
+        let restvec: Vec<_> = seqbtreemap[i + 1usize..seqbtreemap.len()]
+            .iter()
+            .collect::<Vec<_>>();
+        for restvectiter in 0..restvec.len() {
+            let mut countkmer: usize = 0usize;
+            for itercount in 0..vec.1 .1.len() {
+                if restvec[restvectiter]
+                    .1
+                     .1
+                    .contains(&vec.1 .1[itercount].to_string())
+                {
+                    countkmer += 1usize;
+                }
+            }
+            shared.push(CollectIter {
+                name: vec.0.clone(),
+                namenext: restvec[restvectiter].0.clone(),
+                id: vec.1 .0.clone(),
+                idnext: restvec[restvectiter].1 .0.clone(),
+                count: countkmer,
+                shared: vec.1 .1.len() + restvec[restvectiter].1 .1.len(),
+            });
+        }
+        newbase.push(shared);
+    }
+
+    let mut filewrite = File::create("sequence-threshold.fasta").expect("file not found");
+    for i in newbase.iter() {
+        for j in i.iter() {
+            let calibrationpoint: f32 = j.count as f32 / j.shared as f32 * 100.0;
+            if calibrationpoint >= thresh {
+                writeln!(
+                    filewrite,
+                    "{:?}\t{:?}\t{:?}\t{:?}\t{:?}\t{:?}\t{:?}",
+                    j.name,
+                    j.namenext,
+                    j.id,
+                    j.idnext,
+                    j.count,
+                    j.shared,
+                    j.count as f32 / j.shared as f32 * 100.0
+                )
+                .expect("file not found");
+            }
+        }
+    }
+
+    let mut filesecondwrite = File::create("frequencies-threshold.txt").expect("file not present");
+    for i in newbase.iter() {
+        for j in i.iter() {
+            let calibrationpoint: f32 = j.count as f32 / j.shared as f32 * 100.0;
+            if calibrationpoint >= thresh {
+                writeln!(
+                    filesecondwrite,
+                    "{:?}\t{:?}\t{:?}\t{:?}\t{:?}",
+                    j.name,
+                    j.namenext,
+                    j.count,
+                    j.shared,
+                    j.count as f32 / j.shared as f32 * 100.0
+                )
+                .expect("file not found");
+            }
+        }
+    }
+
+    Ok(
+        "The sequence similarity scores and the cluster of the sequences have been written"
+            .to_string(),
+    )
+}