Merge pull request #8 from Juke34/agat

Agat module for normalization of GFF3 files

Author: Juke34

modules/agat.nf

@@ -1,18 +1,17 @@
-process convert_sp_gxf2gxf {
+process normalize_gxf {
     label 'agat'
-    tag "$sample_id"
-    publishDir "${params.outdir}/bamutil_clipoverlap", mode: 'copy'
+    publishDir "${params.outdir}/agat_gff3", mode: 'copy'
 
     input:
         path(gxf)
 
     output:
-        path ("*.gff"), emit: gff
+        path ("*.gff3"), emit: gff
 
     script:
-        gff_file = genome_fasta.baseName.replaceAll(/\..+(\.gz)?$/, '')
+        base_name = gxf.baseName.replaceAll(/\..+(\.gz)?$/, '')
         """
-        agat_convert_sp_gxf2gxf.pl --gxf ${gxf} -o ${gff_file}.gff3
+        agat config --expose --tabix
+        agat_convert_sp_gxf2gxf.pl --gxf ${gxf} -o ${base_name}_normalized.gff3
         """
-
 }

nextflow.config

@@ -54,7 +54,8 @@ profiles {
         params.aline_profiles = "${baseDir}/config/ressources/base_aline.config" 
         params.aligner        = "STAR" 
         params.reads          = "${baseDir}/data/chr21/chr21_small_R1.fastq.gz "
-        params.genome         = "${baseDir}/data/chr21/chr21_small.fasta.gz" 
+        params.genome         = "${baseDir}/data/chr21/chr21_small.fasta.gz"
+        params.annotation     = "${baseDir}/data/chr21/chr21_small_filtered.gff3"
         params.library_type   = "ISR" 
         params.read_type      = "short_single"
     }
@@ -63,6 +64,7 @@ profiles {
         params.aligner        = "STAR" 
         params.reads          = "${baseDir}/data/chr21/"
         params.genome         = "${baseDir}/data/chr21/chr21_small.fasta.gz" 
+        params.annotation     = "${baseDir}/data/chr21/chr21_small_filtered.gff3"
         params.library_type   = "ISR" 
         params.read_type      = "short_paired"
     }

rain.nf

@@ -12,6 +12,7 @@ import java.nio.file.*
 // Input/output params
 params.reads = "/path/to/reads_{1,2}.fastq.gz/or/folder"
 params.genome = "/path/to/genome.fa"
+params.annotation = "/path/to/annotations.gff3"
 params.outdir = "result"
 params.reads_extension = ".fastq.gz" // Extension used to detect reads in folder
 params.paired_reads_pattern = "_{1,2}"
@@ -80,6 +81,9 @@ def helpMSG() {
     --reads                     path to the illumina read file (fastq or fastq.gz) (default: $params.reads)
     --genome                    path to the genome (default: $params.genome)
 
+        Annotation input:
+    --annotation                path to a GFF3 file with annotations of genomic features
+
         Output:
     --output                    path to the output directory (default: $params.outdir)
 
@@ -133,6 +137,7 @@ include {samtools_index; samtools_fasta_index} from './modules/samtools.nf'
 include {reditools2} from "./modules/reditools2.nf"
 include {jacusa2} from "./modules/jacusa2.nf"
 include {sapin} from "./modules/sapin.nf"
+include {normalize_gxf} from "./modules/agat.nf"
 
 //*************************************************
 // STEP 3 - Deal with parameters
@@ -248,4 +253,5 @@ workflow rain {
         samtools_fasta_index(genome)
         jacusa2(samtools_index.out.tuple_sample_bam_bamindex, samtools_fasta_index.out.tuple_fasta_fastaindex)
         sapin(bamutil_clipoverlap.out.tuple_sample_clipoverbam, genome)
+        normalize_gxf(params.annotation)
 }

src/stats/utils.py

@@ -0,0 +1,23 @@
+import Bio
+from typing import Generator
+
+def group_by_overlap(record: Bio.SeqRecord) -> Generator[list[Bio.SeqFeature], None, None]:
+    """Create an iterator that yields groups of features with overlapping genomic positions within a record"""
+
+    feature: Bio.SeqFeature = record.features[0]
+    start: Bio.SeqFeature.ExactPosition = feature.location.start
+    end: Bio.SeqFeature.ExactPosition = feature.location.end
+
+    group: list[Bio.SeqFeature] = []
+
+    for feature in record.features:
+        if start <= feature.location.start <= end:
+            end = feature.location.end
+            group.append(feature)
+        else:
+            yield group
+            start = feature.location.start
+            end = feature.location.end
+            group = [feature]
+
+    yield group