IBCHgenomics

codebiogenomics · codebiogenomics · commit e26b3a271193 · 2025-06-10T13:06:01.000+02:00
final release.
diff --git a/README.md b/README.md
@@ -29,6 +29,42 @@ Options:
   -V, --version  Print version
 ```
 
+- Annotate only the ref variant with A
+```
+./target/debug/varlinker variant-trefanno ./sample-files/sample.vcf A
+```
+
+- Annotate only the alt variant with A
+```
+./target/debug/varlinker variant-taltanno ./sample-files/sample.vcf T
+```
+- Annotate all the variants in the vcf
+```
+./target/debug/varlinker variant-linker ./sample-files/sample.vcf
+```
+
+- it will produce three output files classifying variant annotation to gene, exon and transcript level.
+
+```
+==> annotationfile-exon.txt <==
+19      111     .       A       C       exon    ENSG00000284900.2       KBTBD4
+19      111     .       A       C       exon    ENSG00000284900.2       KBTBD4
+19      111     .       A       C       exon    ENSG00000284900.2       KBTBD4
+19      111     .       A       C       exon    ENSG00000284900.2       KBTBD4
+
+==> annotationfile-gene.txt <==
+19      111     .       A       C       gene    ENSG00000284900.2       KBTBD4
+19      111     .       A       C       gene    ENSG00000278550.4       SLC43A2
+19      111     .       A       C       gene    ENSG00000264450.1       ENSG00000264450
+19      111     .       A       C       gene    ENSG00000273929.1       U6
+
+==> annotationfile-transcript.txt <==
+19      111     .       A       C       transcript      ENSG00000284900.2       KBTBD4
+19      111     .       A       C       transcript      ENSG00000284900.2       KBTBD4
+19      111     .       A       C       transcript      ENSG00000284900.2       KBTBD4
+19      111     .       A       C       transcript      ENSG00000284900.2       KBTBD4
+```
+
 - Acknowledgements: MOSAIC platform, developed as part of the ECBiG-MOSAIC project (POIR.04.02.00-00-D017/20), co-financed by the European Regional Development Fund (ERDF) under the Smart Growth Operational Programme 2014-2020, Measure 4.2 for the development of modern research infrastructure in the science sector.
 - Project PI and Informal queries: Prof. Luiza Handschuh: luizahan@ibch.poznan.pl.
 - Code related queries: Dr. Gaurav Sablok: gsablok@ibch.poznan.pl.
diff --git a/sample-files/sample.vcf b/sample-files/sample.vcf
@@ -28,5 +28,3 @@
 20	1230237	.	T	.	47	PASS	NS=3;DP=13;AA=T	GT:GQ:DP:HQ	0|0:54:.:56,60	0|0:48:4:51,51	0/0:61:2:.,.
 20	1234567	microsat1	G	GA,GAC	50	PASS	NS=3;DP=9;AA=G;AN=6;AC=3,1	GT:GQ:DP	0/1:.:4	0/2:17:2	1/1:40:3
 20	1235237	.	T	.	.	.	.	GT	0/0	0|0	./.
-X	10	rsTest	AC	A,ATG	10	PASS	.	GT	0	0/1	0|2
-
diff --git a/src/store.rs b/src/store.rs
@@ -10,21 +10,43 @@ pub struct VCF {
 }
 
 #[derive(Debug, Clone, PartialOrd, PartialEq)]
+pub struct GenCodeGene {
+    pub chrom: String,
+    pub typeannotate: String,
+    pub start: usize,
+    pub stop: usize,
+    pub geneid: String,
+    pub genename: String,
+}
+
+#[derive(Debug, Clone, PartialOrd, PartialEq)]
+pub struct GenCodeExon {
+    pub chrom: String,
+    pub typeannotate: String,
+    pub start: usize,
+    pub stop: usize,
+    pub geneid: String,
+    pub genename: String,
+}
 
-pub struct GenCode {
+#[derive(Debug, Clone, PartialOrd, PartialEq)]
+pub struct GenCodeTranscript {
     pub chrom: String,
     pub typeannotate: String,
     pub start: usize,
     pub stop: usize,
+    pub geneid: String,
     pub genename: String,
 }
 
+#[derive(Debug, Clone, PartialOrd, PartialEq)]
 pub struct OUTPUT {
     pub chrom: String,
     pub pos: String,
     pub id: String,
     pub refnuc: String,
     pub altnuc: String,
     pub typeannotate: String,
+    pub geneid: String,
     pub genename: String,
 }
diff --git a/src/varaltannot.rs b/src/varaltannot.rs
@@ -1,5 +1,6 @@
+use crate::store::OUTPUT;
 use crate::store::VCF;
-use crate::store::{GenCode, OUTPUT};
+use crate::store::{GenCodeExon, GenCodeGene, GenCodeTranscript};
 use async_std::task;
 use std::error::Error;
 use std::fs::File;
@@ -16,7 +17,7 @@ use std::process::Command;
 */
 
 pub async fn varaltanno(pathfile: &str, variant: &str) -> Result<String, Box<dyn Error>> {
-    let _ = Command::new("wegt").
+    let _ = Command::new("wget").
         arg("https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_48/gencode.v48.chr_patch_hapl_scaff.annotation.gtf.gz")
         .output()
         .expect("command failed");
@@ -26,82 +27,240 @@ pub async fn varaltanno(pathfile: &str, variant: &str) -> Result<String, Box<dyn
         .expect("command failed");
     let fileopen = File::open(pathfile).expect("file not present");
     let fileread = BufReader::new(fileopen);
-    let gtfresults: Vec<GenCode> =
-        task::block_on(gtfread("gencode.v48.chr_patch_hapl_scaff.annotation.gtf")).unwrap();
+    let gtfresults_gene: Vec<GenCodeGene> = task::block_on(gtfread_gene_annotation(
+        "gencode.v48.chr_patch_hapl_scaff.annotation.gtf",
+    ))
+    .unwrap();
+    let gtfresults_exon: Vec<GenCodeExon> = task::block_on(gtfread_exon_annotation(
+        "gencode.v48.chr_patch_hapl_scaff.annotation.gtf",
+    ))
+    .unwrap();
+    let gtfresults_transcript: Vec<GenCodeTranscript> = task::block_on(
+        gtfread_transcript_annotation("gencode.v48.chr_patch_hapl_scaff.annotation.gtf"),
+    )
+    .unwrap();
     let mut vcstring_file: Vec<VCF> = Vec::new();
     for i in fileread.lines() {
         let linevcf = i.expect("file not present");
-        let linevec: Vec<String> = linevcf
-            .split("\t")
-            .map(|x| x.to_string())
-            .collect::<Vec<_>>();
-        vcstring_file.push(VCF {
-            chrom: linevec[0].to_string(),
-            pos: linevec[1].parse::<usize>().unwrap(),
-            id: linevec[2].to_string(),
-            refnuc: linevec[3].to_string(),
-            altnuc: linevec[4].to_string(),
-            qual: linevec[5].to_string(),
-        });
+        if !linevcf.starts_with("#") {
+            let linevec: Vec<String> = linevcf
+                .split("\t")
+                .map(|x| x.to_string())
+                .collect::<Vec<_>>();
+            vcstring_file.push(VCF {
+                chrom: linevec[0].to_string(),
+                pos: linevec[1].parse::<usize>().unwrap(),
+                id: linevec[2].to_string(),
+                refnuc: linevec[3].to_string(),
+                altnuc: linevec[4].to_string(),
+                qual: linevec[5].to_string(),
+            });
+        }
     }
 
-    let mut output: Vec<OUTPUT> = Vec::new();
+    let mut output_gene: Vec<OUTPUT> = Vec::new();
+    let mut output_exon: Vec<OUTPUT> = Vec::new();
+    let mut output_transcript: Vec<OUTPUT> = Vec::new();
+
     for i in vcstring_file.iter() {
-        for j in gtfresults.iter() {
-            if i.pos > j.start && i.pos <= j.stop && j.typeannotate == variant {
-                output.push(OUTPUT {
+        for j in gtfresults_gene.iter() {
+            if i.pos > j.start && i.pos <= j.stop && i.altnuc == variant {
+                output_gene.push(OUTPUT {
                     chrom: i.chrom.clone(),
                     pos: i.pos.clone().to_string(),
                     id: i.id.clone(),
-                    refnuc: i.id.clone(),
+                    refnuc: i.refnuc.clone(),
                     altnuc: i.altnuc.clone(),
                     typeannotate: j.typeannotate.clone(),
+                    geneid: j.geneid.clone(),
                     genename: j.genename.clone(),
                 });
             }
         }
     }
 
-    let mut mutwrite = File::create("annotationfile.txt").expect("file not present");
-    writeln!(
-        mutwrite,
-        "{}\t{}\t{}\t{}\t{}\t{}\t{}",
-        "chrom", "pos", "id", "refnuc", "altnuc", "typeannotate", "genename"
-    )
-    .expect("file not found");
-    for i in output.iter() {
+    for i in vcstring_file.iter() {
+        for j in gtfresults_exon.iter() {
+            if i.pos > j.start && i.pos <= j.stop && i.altnuc == variant {
+                output_exon.push(OUTPUT {
+                    chrom: i.chrom.clone(),
+                    pos: i.pos.clone().to_string(),
+                    id: i.id.clone(),
+                    refnuc: i.refnuc.clone(),
+                    altnuc: i.altnuc.clone(),
+                    typeannotate: j.typeannotate.clone(),
+                    geneid: j.geneid.clone(),
+                    genename: j.genename.clone(),
+                });
+            }
+        }
+    }
+
+    for i in vcstring_file.iter() {
+        for j in gtfresults_transcript.iter() {
+            if i.pos > j.start && i.pos <= j.stop && i.altnuc == variant {
+                output_transcript.push(OUTPUT {
+                    chrom: i.chrom.clone(),
+                    pos: i.pos.clone().to_string(),
+                    id: i.id.clone(),
+                    refnuc: i.refnuc.clone(),
+                    altnuc: i.altnuc.clone(),
+                    typeannotate: j.typeannotate.clone(),
+                    geneid: j.geneid.clone(),
+                    genename: j.genename.clone(),
+                });
+            }
+        }
+    }
+
+    let mut mutwrite_gene = File::create("annotationfile-gene.txt").expect("file not present");
+    let mut mutwrite_exon = File::create("annotationfile-exon.txt").expect("file not present");
+    let mut mutwrite_transcript =
+        File::create("annotationfile-transcript.txt").expect("file not present");
+
+    for i in output_gene.iter() {
         writeln!(
-            mutwrite,
-            "{}\t{}\t{}\t{}\t{}\t{}\t{}",
-            i.chrom, i.pos, i.id, i.refnuc, i.altnuc, i.typeannotate, i.genename
+            mutwrite_gene,
+            "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
+            i.chrom, i.pos, i.id, i.refnuc, i.altnuc, i.typeannotate, i.geneid, i.genename
         )
         .expect("line not found");
     }
+
+    for i in output_exon.iter() {
+        writeln!(
+            mutwrite_exon,
+            "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
+            i.chrom, i.pos, i.id, i.refnuc, i.altnuc, i.typeannotate, i.geneid, i.genename
+        )
+        .expect("line not found");
+    }
+
+    for i in output_transcript.iter() {
+        writeln!(
+            mutwrite_transcript,
+            "{}\t{}\t{}\t{}\t{}\t{}\t{}\t{}",
+            i.chrom, i.pos, i.id, i.refnuc, i.altnuc, i.typeannotate, i.geneid, i.genename
+        )
+        .expect("line not found");
+    }
+
     Ok("The regions have been annotated".to_string())
 }
 
-pub async fn gtfread(gtffile: &str) -> Result<Vec<GenCode>, Box<dyn Error>> {
+pub async fn gtfread_gene_annotation(gtffile: &str) -> Result<Vec<GenCodeGene>, Box<dyn Error>> {
     let fileopen = File::open(gtffile).expect("file not found");
     let fileread = BufReader::new(fileopen);
-    let mut gtf_vector: Vec<GenCode> = Vec::new();
+    let mut gtf_vector: Vec<GenCodeGene> = Vec::new();
     for i in fileread.lines() {
         let linegtf = i.expect("line not found");
-        let linevec: Vec<String> = linegtf
-            .split("\t")
-            .map(|x| x.to_string())
-            .collect::<Vec<String>>();
-        let linecollect: String = linevec[9].split(";").collect::<Vec<_>>()[2]
-            .replace(" ", "-")
-            .split("-")
-            .collect::<Vec<_>>()[1]
-            .to_string();
-        gtf_vector.push(GenCode {
-            chrom: linevec[0].clone(),
-            typeannotate: linevec[2].clone(),
-            start: linevec[3].parse::<usize>().unwrap(),
-            stop: linevec[4].parse::<usize>().unwrap(),
-            genename: linecollect,
-        })
+        if !linegtf.starts_with("#") {
+            let linevec: Vec<String> = linegtf
+                .split("\t")
+                .map(|x| x.to_string())
+                .collect::<Vec<String>>();
+            if linevec[2] == "gene" {
+                let linecollect: String = linevec[8].split(";").collect::<Vec<_>>()[2]
+                    .replace(" ", "-")
+                    .split("-")
+                    .collect::<Vec<_>>()[2]
+                    .to_string()
+                    .replace("\"", "");
+                let genecollect: String = linevec[8].split(";").collect::<Vec<_>>()[0]
+                    .replace(" ", "-")
+                    .split("-")
+                    .collect::<Vec<_>>()[1]
+                    .to_string()
+                    .replace("\"", "");
+                gtf_vector.push(GenCodeGene {
+                    chrom: linevec[0].clone(),
+                    typeannotate: linevec[2].clone(),
+                    start: linevec[3].parse::<usize>().unwrap(),
+                    stop: linevec[4].parse::<usize>().unwrap(),
+                    geneid: genecollect,
+                    genename: linecollect,
+                })
+            }
+        }
+    }
+    Ok(gtf_vector)
+}
+
+pub async fn gtfread_exon_annotation(gtffile: &str) -> Result<Vec<GenCodeExon>, Box<dyn Error>> {
+    let fileopen = File::open(gtffile).expect("file not found");
+    let fileread = BufReader::new(fileopen);
+    let mut gtf_vector: Vec<GenCodeExon> = Vec::new();
+    for i in fileread.lines() {
+        let linegtf = i.expect("line not found");
+        if !linegtf.starts_with("#") {
+            let linevec: Vec<String> = linegtf
+                .split("\t")
+                .map(|x| x.to_string())
+                .collect::<Vec<String>>();
+            if linevec[2] == "exon" {
+                let linecollect: String = linevec[8].split(";").collect::<Vec<_>>()[3]
+                    .replace(" ", "-")
+                    .split("-")
+                    .collect::<Vec<_>>()[2]
+                    .to_string()
+                    .replace("\"", "");
+                let genecollect: String = linevec[8].split(";").collect::<Vec<_>>()[0]
+                    .replace(" ", "-")
+                    .split("-")
+                    .collect::<Vec<_>>()[1]
+                    .to_string()
+                    .replace("\"", "");
+
+                gtf_vector.push(GenCodeExon {
+                    chrom: linevec[0].clone(),
+                    typeannotate: linevec[2].clone(),
+                    start: linevec[3].parse::<usize>().unwrap(),
+                    stop: linevec[4].parse::<usize>().unwrap(),
+                    geneid: genecollect,
+                    genename: linecollect,
+                })
+            }
+        }
+    }
+    Ok(gtf_vector)
+}
+
+pub async fn gtfread_transcript_annotation(
+    gtffile: &str,
+) -> Result<Vec<GenCodeTranscript>, Box<dyn Error>> {
+    let fileopen = File::open(gtffile).expect("file not found");
+    let fileread = BufReader::new(fileopen);
+    let mut gtf_vector: Vec<GenCodeTranscript> = Vec::new();
+    for i in fileread.lines() {
+        let linegtf = i.expect("line not found");
+        if !linegtf.starts_with("#") {
+            let linevec: Vec<String> = linegtf
+                .split("\t")
+                .map(|x| x.to_string())
+                .collect::<Vec<String>>();
+            if linevec[2] == "transcript" {
+                let linecollect: String = linevec[8].split(";").collect::<Vec<_>>()[3]
+                    .replace(" ", "-")
+                    .split("-")
+                    .collect::<Vec<_>>()[2]
+                    .to_string()
+                    .replace("\"", "");
+                let genecollect: String = linevec[8].split(";").collect::<Vec<_>>()[0]
+                    .replace(" ", "-")
+                    .split("-")
+                    .collect::<Vec<_>>()[1]
+                    .to_string()
+                    .replace("\"", "");
+                gtf_vector.push(GenCodeTranscript {
+                    chrom: linevec[0].clone(),
+                    typeannotate: linevec[2].clone(),
+                    start: linevec[3].parse::<usize>().unwrap(),
+                    stop: linevec[4].parse::<usize>().unwrap(),
+                    geneid: genecollect,
+                    genename: linecollect,
+                })
+            }
+        }
     }
     Ok(gtf_vector)
 }
diff --git a/src/variantlinker.rs b/src/variantlinker.rs
diff --git a/src/varrefannot.rs b/src/varrefannot.rs
