Reannotation.Snakemake

# =====================================
# Reannotation.Snakemake
# =====================================
#
# Purpose:
#   This Snakemake workflow performs genome reannotation using RNA-seq data,
#   combining evidence from multiple sources including Trinity assemblies,
#   PASA updates, and protein predictions.
#
# Workflow Steps:
#   1. RNA-seq Mapping (STAR)
#   2. Transcriptome Assembly
#      - De novo assembly (Trinity)
#      - Genome-guided assembly (Trinity)
#      - Reference-guided assembly (Cufflinks)
#   3. PASA Analysis
#      - Transcript alignment and assembly
#      - Annotation comparison and update
#   4. Protein Prediction
#      - TransDecoder
#      - CPAT
#      - BLAST/PFAM validation
#   5. lncRNA Identification
#
# Input Requirements:
#   - Genome sequence (FASTA)
#   - Initial annotation (GFF3)
#   - RNA-seq reads (FASTQ)
#   - Protein databases (UniProt, PFAM)
#
# Output Files:
#   - Updated annotation (GFF3)
#   - Novel transcripts
#   - lncRNA annotations
#   - Protein predictions
#
# Dependencies:
#   - STAR (RNA-seq aligner)
#   - Trinity (transcriptome assembler)
#   - Cufflinks
#   - PASA
#   - TransDecoder
#   - CPAT
#   - BLAST
#   - PFAM
#   - MySQL (for PASA)
#
# Usage:
#   snakemake -s Reannotation.Snakemake \
#     --cluster 'sbatch {params.cluster}' \
#     --stats snakemake.stats \
#     -j <jobs>
#
# Author: htafer
# Last Updated: 2025-07-28
# =====================================

import math
import os
import os.path
from snakemake.utils import validate
import fnmatch

# Helper Functions
#################################

def getFileNumber(file_name, entryNumber):
    """Count FASTA entries and calculate optimal split number."""
    with open(file_name) as f:
        total = sum(1 for line in f if line.startswith(">"))
    return math.floor(2 * total / entryNumber + 1)

def getBamFile(dir, species):
    """Get list of BAM files for a species."""
    pattern = "*.*.bam"
    species_dir = os.path.join(dir, species + ".reads")
    return [f for f in os.listdir(species_dir) 
            if fnmatch.fnmatch(f, pattern)]

def uniquify(seq, idfun=None):
    """Remove duplicates from a sequence while maintaining order."""
    seen = set()
    if idfun is None:
        idfun = lambda x: x
    
    return [x for x in seq if not (idfun(x) in seen or seen.add(idfun(x)))]

# Configuration
#################################

# Load config file
configfile: "config.yaml"

# Directory Structure
WORKDIR = config.get("workdir", "/media/work/genomes/exoDer/annotation/reAnnotation")
COMPUTEDIR = config.get("computedir", "/global/lv70539/htafer")
READS = config.get("reads_dir", "/global/lv70539/htafer/reads")
LOGDIR = "logs"

# Computing Resources
THREADS = config.get("threads", 16)
MAX_MEMORY = config.get("max_memory", "64G")

# Project Information
ID = config.get("project_id", "exoDer")
GENOME = ID + ".fasta"
ANNOT = ID + ".Protein.gff3"

# Input Files
BAMFILES = os.path.join(READS, "{samples}.fastq")
BAMS, = glob_wildcards(BAMFILES)
BAMS = sorted(BAMS)  # Ensure consistent order

# Database Paths
BLAST = config.get("blast_db", "/home/htafer/share/database/UniProt90pSaccharomyceta.fasta")
PFAM = config.get("pfam_db", "/home/htafer/bin/PfamScan/db/")

# Analysis Parameters
INTRON = config.get("max_intron_length", 2000)
MIN_TRANSCRIPT_LENGTH = config.get("min_transcript_length", 200)
MIN_PROTEIN_LENGTH = config.get("min_protein_length", 100)

# Cluster Configuration
CLUSTER_PARTITION = config.get("cluster_partition", "mem_0064")
CLUSTER_QOS = config.get("cluster_qos", "normal_0064")

# Tool Parameters
TRINITY_MEMORY = config.get("trinity_memory", "200G")
STAR_MEMORY = config.get("star_memory", "64G")
CUFFLINKS_MEMORY = config.get("cufflinks_memory", "32G")

    ##################################                                 
    #FILES                           #                                 
    ##################################                                 

#GENOME           
ID="exoDer";
REF=ID+".fasta"

#BAMFILES         
BAMFILES="/global/lv70539/htafer/reads/{samples}.fastq"
BAMS ,= glob_wildcards(BAMFILES);

ANNOT=ID+".Protein.gff3"
BLAST="/home/htafer/share/database/UniProt90pSaccharomyceta.fasta"
PFAM="/home/htafer/bin/PfamScan/db/"

##################################                                     
#BIOLOGY                         #                                     
##################################                                     

#BIOLOGY          
INTRON=2000



#Which rules are run locally # depends on which computer the pipeline is running                         
#localrules: all, clean, STAR, STARIdx, bamToFastq, mergeBam, cufflinks, composeMerge, mergeAssemblies, normalization, trinityAlignment, trinityDeNovo, PASA, PASAtrainingset, PASAhints, mergeFastq, Annotcompare
localrules: all, clean, composeMerge, mergeAssemblies, PASAtrainingset, PASAhints, mergeFastq, STARIdx

# Workflow Targets
#################################

# Rule: all
# Purpose: Define all output files required from the workflow
rule all:
    input: 
        annotation = ID + ".AnnotationUpdate",             # Updated genome annotation
        pasa = ID + ".PASA",                             # PASA analysis results
        introns = ID + ".introns.gff",                   # Intron annotations
        transdecoder = ID + ".PASAtrainingset",          # TransDecoder predictions
        new_proteins = ID + ".newProtein.gff3",          # Novel protein annotations
        other_transcripts = ID + ".noKnownFunc.gff3",    # Uncharacterized transcripts
        validation = "logs/validate.done"                # Input validation flag

# Rule: validate_inputs
# Purpose: Check that all required input files and tools exist
rule validate_inputs:
    output:
        touch("logs/validate.done")
    run:
        # Check required directories
        for dir in [WORKDIR, COMPUTEDIR, READS]:
            if not os.path.exists(dir):
                raise ValueError(f"Required directory not found: {dir}")
        
        # Check required input files
        required_files = [
            os.path.join(WORKDIR, GENOME),
            os.path.join(WORKDIR, ANNOT),
            BLAST,
            os.path.join(PFAM, "Pfam-A.hmm")
        ]
        
        missing_files = [f for f in required_files if not os.path.exists(f)]
        if missing_files:
            raise ValueError("Missing required files:\n" + "\n".join(missing_files))
        
        # Check RNA-seq input files
        missing_reads = []
        for sample in BAMS:
            fastq = os.path.join(READS, f"{sample}.fastq")
            if not os.path.exists(fastq):
                missing_reads.append(fastq)
            elif os.path.getsize(fastq) == 0:
                missing_reads.append(f"{fastq} (empty)")
        
        if missing_reads:
            raise ValueError("Missing or empty read files:\n" + "\n".join(missing_reads))
        
        # Check required tools are in PATH
        required_tools = ["STAR", "trinity", "cufflinks", "blastx", "pfam_scan.pl"]
        missing_tools = []
        for tool in required_tools:
            if not shutil.which(tool):
                missing_tools.append(tool)
        
        if missing_tools:
            raise ValueError("Required tools not found in PATH:\n" + "\n".join(missing_tools))

##################################################################################                              
#                                                                                #
#                  Get The Rest                                                  #
#                                                                                #
##################################################################################                              

rule otherTranscript:
     input: lnc=ID+".lncRNA.gff3", coding=ID+".allProtein.gff3", pasaAssemblies=ID+".PASA/"+ID+".pasa_assemblies.gff3"
     output: ID+".noKnownFunc.gff3"
     shell:"""    
     mergealign.pl < {input.pasaAssemblies} | cut -f 1-9 >  merged.gff 
     bedtools intersect -s -v -a merged.gff -b {input.lnc} | bedtools intersect -s -v -a - -b {input.coding} > {output}                                    
     """
#################################################################################                              
#                                                                               #
#                  Get lncRNAs                                                  #
#                                                                               #
#################################################################################                              
########################            
####    lncRNA      ####            
########################            

rule lncRNA:
     input: ids="./CPAT/cpat.results", blast="BLAST/blast.coding.candidates", pfam="PFAM/"+ID+".pfam.coding.candidates", pasaAssemblies=ID+".PASA/"+ID+".pasa_assemblies.gff3"
     output: lnc=ID+".lncRNA.gff3"
     threads: THREADS
     params:
     shell:"""                      
     cat {WORKDIR}/{input.blast} {WORKDIR}/{input.pfam} | sort -u > coding.candidates                               
     cut -f 1 {WORKDIR}/{input.ids}  > id.cpat.results                      
     cat id.cpat.results | perl -lane 'my $count=`grep -cw $F[0] coding.candidates`; chomp($count); if($count==0){{print $F[0]}};' | sort -u  | fgrep -w -f - {WORKDIR}/{input.pasaAssemblies} > {WORKDIR}/{output.lnc}                     
     """
########################            
####      BLAST     ####            
########################            

rule blast:
     input: cpat="./CPAT/"+ID+".cpat.fasta"
     output: blast="BLAST/blast.coding.candidates"
     threads: THREADS
     params:
     shell:"""                      
     cd BLAST                       
     formatdb -i {BLAST} -p T -n uniprot90S                                 
     blastx -query {WORKDIR}/{input.cpat} -db uniprot90S -num_threads {threads} -outfmt 7 -out blast.out            
     cat blast.out | perl -lane 'if($F[7]>$F[6] && $F[10]<0.001){{print;}}' | cut -f 1 | sort -u | sed -r 's/asmbl/>asmbl/' > {WORKDIR}/{output.blast}      
     """

########################            
####     PFAM       ####            
########################            

rule pfam:
     input: "./CPAT/"+ID+".cpat.fasta"
     output: "PFAM/"+ID+".pfam.coding.candidates"
     threads: THREADS
     shell: """                     
     mkdir -p PFAM                  
     cd PFAM                        
     pfam_scan.pl -fasta {WORKDIR}/{input} -dir {PFAM} -cpu {threads} -outfile pfam.out                             
     cat pfam.out | grep -v '#' | perl -lane 'if($F[2]>$F[1] && $F[12]<0.001){{print;}}' | cut -f 1 | sort -u | sed -r 's/asmbl/>asmbl/' > {WORKDIR}/{output}                                       
                                    
     """



################                    
#### CPAT   ####                    
################                    

rule cpat:
     input: gene="noOverlap.fasta", fa="assemblies.fa", nonc="nonconding.fasta", coding="codingseq.fasta"
     output: fasta="./CPAT/"+ID+".cpat.fasta", results="CPAT/cpat.results"
     threads: THREADS
     shell:"""                      
     mkdir -p CPAT                  
     cd CPAT                        
     make_hexamer_tab.py -c {WORKDIR}/{input.coding} -n {WORKDIR}/{input.nonc} > hexamer.tab                        
     make_logitModel.py -c {WORKDIR}/{input.coding} -n {WORKDIR}/{input.nonc} -x hexamer.tab -o {ID}                
     cpat.py -d {ID}.logit.RData -x hexamer.tab -g {WORKDIR}/{input.gene} -o {ID}.cpat                              
     cat {ID}.cpat |  perl -lane 'if($F[5]<0.01){{print;}}' | sort -k 2,2gr | sed -r 's/ASMBL/asmbl/' > {WORKDIR}/{output.results}                          
     parallel -j {threads} 'grep -P "{{}}$" {WORKDIR}/{input.fa} -A 1 ' ::: `cat {WORKDIR}/{output.results}  | cut -f 1 | grep asmbl` > {WORKDIR}/{output.fasta}                                    
     """


rule prepareFiles:
     input: fa="assemblies.fa", intron=ID+".introns.gff", updateAnn=ID+".allProtein.gff3"
     output: coding="codingseq.fasta", nonc="nonconding.fasta"
     params:
     threads: THREADS
     shell:"""                      
     #We have to prepare a training data set in order to train CPAT         
     #Start with the protein dataset
     getAnnoFasta.pl {input.updateAnn} --seqfile {REF}                      
     cat {ID}.allProtein3.cdsexons | sed -r 's/.cds[0-9]//g' | perl -lane 'BEGIN{{my $previousId="";}} if($F[0]=~/>/){{if($F[0] eq $previousId){{next;}}else{{$previousId=$F[0]; print "\n",$F[0]}}}}else{{printf $F[0]}}' > {output.coding}
     #Now generate the non coding dataset with the intron                   
     cat {input.intron} | sed -r 's/;src=E//g' > temp.introns               
     parallel -j {threads} 'grep -P "{{}}$" {input.fa} -A 1 ' ::: `cat temp.introns | cut -f 9 | sed -r 's/.+Target=([^$]+)/\\1/' | sort -u` > {output.nonc}
     """



################                    
#   INTERSECT  #                    
################                    
rule intersect:
     input: transAnn=ID+".PASAtrainingset/"+ID+".assemblies.fasta.transdecoder.genome.gff3", pasaAssemblies=ID+".PASA/"+ID+".pasa_assemblies.gff3", pasaSequence=ID+".PASA/"+ID+".assemblies.fasta", updateAnn=ID+".allProtein.gff3"
     output: fasta="noOverlap.fasta", fa="assemblies.fa"
     params:
     threads: THREADS
     shell:"""                      
     #put the gene together         
     mergealign.pl < {input.pasaAssemblies} | cut -f 1-9 >  merged.gff      
     #check which gene overlap neither with the updated protein annotation nor with the pasa protein annotation     
     intersectBed -v -s -a merged.gff -b {input.updateAnn} > no.gff3        
     intersectBed -v -s -a no.gff3 -b {input.transAnn} > noOverlap.gff      
     perl ~/bin/rmenterdb.pl < {input.pasaSequence} > {output.fa}           
     #Fetch the sequences           
     parallel -j {threads} 'grep -P "{{}}$" {output.fa} -A 1 ' ::: `cat noOverlap.gff | cut -f 9 | sed -r 's/.+Target=([^$]+)/\\1/' | sort -u` > {output.fasta}                                     
     """



#################################################################################                                   
#                                                                               #                                   
#                  Get New Proteins                                             #                                   
#                                                                               #                                   
#################################################################################                  

rule newProtein:
       input: transAnn=ID+".PASAtrainingset/"+ID+".assemblies.fasta.transdecoder.genome.gff3", updateAnn=ID+".AnnotationUpdate/"+ID+".update.gff3"
       output:newProtein=ID+".newProtein.gff3", allProtein=ID+".allProtein.gff3"
       shell:"""       
       grep -P "\tCDS\t" {input.transAnn} | mergealign.pl > {ID}.transdecoder.CDS.gff3                 
       grep -P "\tCDS\t" {input.updateAnn} | mergealign.pl > {ID}.update.CDS.gff3                      
       bedtools intersect -v -a {ID}.transdecoder.CDS.gff3 -b {ID}.update.CDS.gff3 | cut -f 9 | sed -r 's/.+Parent=//' | fgrep -f - {input.transAnn} | grep -P "\tCDS\t" > {ID}.newProtein.CDS.gff3                            
       cut -f 9 {ID}.newProtein.CDS.gff3 | sed -r 's/.+Parent=//' | fgrep -f - {input.transAnn} | perl -lane 'if($F[2]=~/mRNA/){{$anno=$F[8]; $anno=~s/ID=.+Parent=/ID=/; $anno=~/=(asmbl[^g]+g[^;]+)/; my  $id=$1; $id=~s/g\./m./; $anno=~s/ORF/$id/; print "$F[0]\t$F[1]\tgene\t$F[3]\t$F[4]\t$F[5]\t$F[6]\t$F[7]\t$anno"; $anno=$F[8]; $anno=~s/ORF/$id/; print join("\t",@F[0..7]),"\t$ann\
o";}}else{{print;}}' > {output.newProtein}                     
       cat {input.updateAnn} {output.newProtein}  > {output.allProtein}                                
       """


#################################################################################                                   
#                                                                               #                                   
#                  Annotation Update                                            #                                   
#                                                                               #                                   
#################################################################################                                   

rule Annotcompare:
      input: genome=ID+".fasta", annot=ANNOT, dn="./reads/"+ID+".trinityDN", gg="./reads/"+ID+".trinityGG",dirIn=ID+".PASA"
      output: dir=ID+".AnnotationUpdate", updateGFF3=ID+".AnnotationUpdate/"+ID+".update.gff3"
      params: cluster="--partition=mem_0064 --qos=normal_0064"
      threads: THREADS
      shell: """                    
      #mysql start                  
      cd $HOME/bin/mysql-5.6.22/    
      $HOME/bin/mysql-5.6.22/scripts/mysql_install_db --default-files=./my-new.cnf --datadir=$HOME/bin/mysql-5.6.22/data                                    
      cp $HOME/my-new.cnf $HOME/bin/mysql-5.6.22/                           
      $HOME/bin/mysql-5.6.22/bin/mysqld_safe --defaults-file=./my-new.cnf --skip-grant-tables &                     
      sleep 10                      
      prefix=`date --rfc-3339=ns  | md5sum | head -c 16`                    
      mkdir -p {COMPUTEDIR}/PASA${{prefix}}                                 
      cd {COMPUTEDIR}/PASA${{prefix}}                                       
      #check compatibility          
      #$PASAHOME/misc_utilities/pasa_gff3_validator.pl {WORKDIR}/{input.annot}                                      
      cat {WORKDIR}/{input.dn}/Trinity.fasta {WORKDIR}/{input.gg}/Trinity-GG.fasta > {COMPUTEDIR}/PASA${{prefix}}/transcripts.fasta                         
      #Prepare annotCompare.config  
      DBNAME={ID}                   
      cat $PASAHOME/pasa_conf/pasa.annotationCompare.Template.txt | sed -r "s/DBNAME/${{DBNAME}}/" > ./annotCompare.config                                  
      $PASAHOME/scripts/Launch_PASA_pipeline.pl -c ./annotCompare.config -A -L --annots_gff3 {WORKDIR}/{input.annot} -g {WORKDIR}/{input.genome} -t ./transcripts.fasta                             
      mv {COMPUTEDIR}/PASA${{prefix}} {WORKDIR}/{output.dir}                
      cd {WORKDIR}/{output.dir}     
      mv *gene_structures_post_PASA_updates.*.gff3 {output.updateGFF3}      
      #mysql shutdown               
      cd $HOME/bin/mysql-5.6.22/    
      $HOME/bin/mysql-5.6.22/bin/mysqladmin --defaults-file=./my-new.cnf shutdown                                   
      """

##################################  
##                               #  
##      PASA  FIRST PASS         #  
##                               #  
##################################  

rule PASA:
     input: genome=ID+".fasta", dn="./reads/"+ID+".trinityDN", gg="./reads/"+ID+".trinityGG", cf=ID+".cufflinks.gtf"
     output: dir=ID+".PASA"
     params: cluster="--partition=mem_0064 --qos=normal_0064"
     threads: THREADS
     shell:"""                      
     #mysql start                   
     cd $HOME/bin/mysql-5.6.22/     
     $HOME/bin/mysql-5.6.22/scripts/mysql_install_db --default-files=./my-new.cnf --datadir=$HOME/bin/mysql-5.6.22/data                                     
     cp $HOME/my-new.cnf $HOME/bin/mysql-5.6.22/                            
     $HOME/bin/mysql-5.6.22/bin/mysqld_safe --defaults-file=./my-new.cnf --skip-grant-tables &                      
     sleep 10                       
     prefix=`date --rfc-3339=ns  | md5sum | head -c 16`                     
     mkdir -p {COMPUTEDIR}/PASA${{prefix}}                                  
     cd {COMPUTEDIR}/PASA${{prefix}}
     #Prepare files                 
     cat {WORKDIR}/{input.dn}/Trinity.fasta {WORKDIR}/{input.gg}/Trinity-GG.fasta > {COMPUTEDIR}/PASA${{prefix}}/transcripts.fasta                          
     # NOT genome-guided, only de novo                                      
     $PASAHOME/misc_utilities/accession_extractor.pl < {WORKDIR}/{input.dn}/Trinity.fasta > {COMPUTEDIR}/PASA${{prefix}}/tdn.accs                           
     #Prepare alignAssembly.config  
     DBNAME={ID}                    
     cat $PASAHOME/pasa_conf/pasa.alignAssembly.Template.txt | sed -r "s/DBNAME/${{DBNAME}}/" > ./alignAssembly.config                                      
     $PASAHOME/scripts/Launch_PASA_pipeline.pl -c ./alignAssembly.config -C -r -R -g {WORKDIR}/{input.genome} --ALIGNERS blat,gmap -t ./transcripts.fasta --transcribed_is_aligned_orient --TDN ./tdn.accs --cufflinks_gtf {WORKDIR}/{input.cf} -I {INTRON} --stringent_alignment_overlap 30.0 --CPU {threads}              
     $PASAHOME/scripts/build_comprehensive_transcriptome.dbi -c ./alignAssembly.config -t ./transcripts.fasta --min_per_ID 95 --min_per_aligned 30          
     mv {COMPUTEDIR}/PASA${{prefix}} {WORKDIR}/{output.dir}                 
     #shutdown mysql                
     cd $HOME/bin/mysql-5.6.22/     
     $HOME/bin/mysql-5.6.22/bin/mysqladmin --defaults-file=./my-new.cnf shutdown                                    
     """


rule PASAtrainingset:
     input: ID+".PASA"
     output: dir=ID+".PASAtrainingset"
     params: cluster="--partition=mem_0064 --qos=normal_0064"
     threads: THREADS
     shell:"""                      
     cd {WORKDIR}/{input}           
     $PASAHOME/scripts/pasa_asmbls_to_training_set.dbi --pasa_transcripts_fasta {WORKDIR}/{input}/{ID}.assemblies.fasta --pasa_transcripts_gff3 {WORKDIR}/{input}/{ID}.pasa_assemblies.gff3         
     mkdir {WORKDIR}/{output}       
     mv {ID}.assemblies.fasta.transdecoder.* {WORKDIR}/{output}             
     """

rule PASAhints:
     input: ID+".PASA"
     output: ID+".introns.gff"
     params: cluster="-cwd -V"
     threads: THREADS
     shell:"""                      
     perl {WORKDIR}/introncalc2.0.pl < {WORKDIR}/{input}/{ID}.pasa_assemblies.gff3 > {WORKDIR}/{output}             
     """

##################################                                     
##                               #                                     
##      TRINITY                  #                                     
##                               #                                     
##################################                                     


rule mergeFastq:
     input: expand("./reads/{samples}.fastq", samples=BAMS)
     output: "./reads/"+ID+".merged.fq"
     params: cluster="--partition=mem_0064 --qos=normal_0064"
     threads: THREADS
     shell:"""    
     cat {WORKDIR}/{input} > {WORKDIR}/{output}                        
     """

rule normalization:
     input: "./reads/"+ID+".merged.fq"
     output: "./reads/"+ID+".trinityIn"
     params: cluster="--partition=mem_0256 --qos=normal_0256"
     threads: THREADS
     shell:"""    
     prefix=`date --rfc-3339=ns  | md5sum | head -c 16`                
     mkdir -p {COMPUTEDIR}/trinityIn${{prefix}}                        
     $HOME/bin/trinityrnaseq-2.0.6/util/insilico_read_normalization.pl --seqType fq --single {input} --SS_lib_type F --JM 200G --max_cov 30 --CPU {threads} --output {COMPUTEDIR}/trinityIn${{prefix}}                                     
     mv {COMPUTEDIR}/trinityIn${{prefix}} {WORKDIR}/{output}           
     """

rule trinityDeNovo:
     input: "./reads/"+ID+".trinityIn"
     output: "./reads/"+ID+".trinityDN"
     params: cluster="--partition=mem_0256 --qos=normal_0256"
     threads: THREADS
     shell:"""    
     prefix=`date --rfc-3339=ns  | md5sum | head -c 16`                
     mkdir -p {COMPUTEDIR}/trinity${{prefix}}                          
     fastQName=`ls {WORKDIR}/{input}/*.ok  | sed -r 's/.ok//'`         
     Trinity --seqType fq --single ${{fastQName}}  --SS_lib_type F --CPU {threads} --output {COMPUTEDIR}/trinity${{prefix}}  --max_memory 240G             
     mv {COMPUTEDIR}/trinity${{prefix}} {WORKDIR}/{output}             
     """


rule trinityAlignment:
     input: "./reads/"+ID+".merged.bam"
     output: dire="./reads/"+ID+".trinityGG"
     params: cluster="--partition=mem_0128 --qos=normal_0128"
     threads: THREADS
     shell:"""    
     #prefix=`date --rfc-3339=ns  | md5sum | head -c 16`               
     prefix="cc3498d43dd4cabb" 
     mkdir -p {COMPUTEDIR}/trinity${{prefix}}                          
     Trinity --genome_guided_bam {input} --genome_guided_max_intron {INTRON} --SS_lib_type F --max_memory 120G --CPU {threads} --output {COMPUTEDIR}/trinity${{prefix}}                            
     mv {COMPUTEDIR}/trinity${{prefix}} {WORKDIR}/{output}             
     """

##################################                                     
##                               #                                     
##      CUFFLINKS                #                                     
##                               #                                     
##################################                                     


rule mergeAssemblies:
     input: "./reads/assemblies.txt"
     output: dir="./reads", file=ID+".cufflinks.gtf"
     params: cluster="--partition=mem_0064 --qos=normal_0064"
     threads: THREADS
     shell:"""    
     cuffmerge -o {WORKDIR}/{output.dir} -s {REF} {input} -p {threads} 
     mv {WORKDIR}/{output.dir}/merged.gtf {output.file}                
     """


rule composeMerge:
     input: expand("./reads/{samples}.cufflinks/transcripts.gtf", samples=BAMS)
     output: txt="./reads/assemblies.txt"
     shell:"""    
     ls {WORKDIR}/{input} > {output.txt}                               
     """

rule cufflinks:
     input: "./reads/{samples}.mapped.sam"
     output: dir="./reads/{samples}.cufflinks", file="./reads/{samples}.cufflinks/transcripts.gtf"
     params: cluster="--partition=mem_0064 --qos=normal_0064"
     threads: THREADS
     shell:"""    
     prefix=`date --rfc-3339=ns  | md5sum | head -c 16`                
     mkdir -p {COMPUTEDIR}/${{prefix}}                                 
     cufflinks -o {COMPUTEDIR}/${{prefix}} -p {threads} -u -I {INTRON} --max-bundle-length 10000 --library-type ff-firststrand --min-intron-length 30 {WORKDIR}/{input}                            
     mv {COMPUTEDIR}/${{prefix}}/* {WORKDIR}/{output.dir}              
     """

##################################                                     
##                               #                                     
##     Spliced RNAseq mapping    #                                     
##                               #                                     
##################################                                     


rule mergeBam:
     input: expand("./reads/{samples}.mapped.bam", samples=BAMS)
     output: "./reads/"+ID+".merged.bam"
     params: cluster="--partition=mem_0064 --qos=normal_0064"
     threads: THREADS
     shell:"""    
     prefix=`date --rfc-3339=ns  | md5sum | head -c 16`                
     mkdir -p {COMPUTEDIR}/${{prefix}}                                 
     samtools merge {COMPUTEDIR}/${{prefix}}/merged.bam {input}        
     samtools sort {COMPUTEDIR}/${{prefix}}/merged.bam {COMPUTEDIR}/${{prefix}}/merged.sorted                  
     #mergeBam INPUT={WORKDIR}/{input} SORT_ORDER=coordinate OUTPUT={WORKDIR}/{output}                         
     mv {COMPUTEDIR}/${{prefix}}/merged.sorted.bam {WORKDIR}/{output}  
     rm -rf {COMPUTEDIR}/${{prefix}}                                   
     """


# RNA-seq Mapping Rules
#################################

# Rule: STAR
# Purpose: Map RNA-seq reads to the genome using STAR
rule STAR:
    input: 
        reads = "./reads/{samples}.fastq",
        idx = ID + ".idx"
    output: 
        bam = "./reads/{samples}.mapped.bam",
        sam = "./reads/{samples}.mapped.sam"
    params:
        # Cluster parameters
        cluster = "--partition={partition} --qos={qos}".format(
            partition=CLUSTER_PARTITION,
            qos=CLUSTER_QOS
        ),
        # STAR parameters
        chim_segment_min = config.get("star", {}).get("chim_segment_min", 12),
        overhang_min = config.get("star", {}).get("sjdb_overhang", 100),
        max_intron = INTRON
    threads: THREADS
    resources:
        mem_mb = lambda wildcards, attempt: attempt * int(STAR_MEMORY.rstrip("G")) * 1024
    log:
        "logs/mapping/{samples}_star.log"
    shell: """
        (# Create temp directory
        prefix=`date --rfc-3339=ns | md5sum | head -c 16`
        tempdir="{COMPUTEDIR}/${{prefix}}"
        mkdir -p "$tempdir"
        cd "$tempdir"
        
        # Run STAR alignment
        STAR --genomeDir {WORKDIR}/{input.idx} \
             --readFilesIn {WORKDIR}/{input.reads} \
             --runThreadN {threads} \
             --twopassMode Basic \
             --outReadsUnmapped None \
             --chimSegmentMin {params.chim_segment_min} \
             --chimJunctionOverhangMin {params.chim_segment_min} \
             --alignSJDBoverhangMin 10 \
             --alignIntronMax {params.max_intron} \
             --chimSegmentReadGapMax 3 \
             --alignSJstitchMismatchNmax 5 -1 5 5 \
             --outSAMtype BAM SortedByCoordinate \
             2> {log}

        # Process alignment files
        samtools sort Aligned.sortedByCoord.out.bam -o Aligned.sorted.bam
        samtools view -h Aligned.sorted.bam | \
            awk 'BEGIN {{OFS="\t"}} {{
                if ($0 ~ /^@/) {{
                    print
                    next
                }}
                split($6,C,/[0-9]*/)
                split($6,L,/[SMDIN]/)
                if (C[2]=="S") {{
                    $10=substr($10,L[1]+1)
                    $11=substr($11,L[1]+1)
                }}
                if (C[length(C)]=="S") {{
                    L1=length($10)-L[length(L)-1]
                    $10=substr($10,1,L1)
                    $11=substr($11,1,L1)
                }}
                gsub(/[0-9]*S/,"",$6)
                print
            }}' > Aligned.sorted.sam
        
        # Move results to final location
        mv Aligned.sorted.sam {WORKDIR}/{output.sam}
        mv Aligned.sorted.bam {WORKDIR}/{output.bam}
        
        # Cleanup
        rm -rf "$tempdir") 2>> {log}
        """

# Rule: STARIdx
# Purpose: Generate STAR genome index
rule STARIdx:
    input:
        genome = ID + ".fasta",
        gff = ID + ".Protein.gff"
    output:
        dir = directory(ID + ".idx")
    params:
        overhang = config.get("star", {}).get("sjdb_overhang", 100)
    threads: THREADS
    resources:
        mem_mb = lambda wildcards, attempt: attempt * int(STAR_MEMORY.rstrip("G")) * 1024
    log:
        "logs/mapping/star_index.log"
    shell: """
        (# Create temp directory
        prefix=`date --rfc-3339=ns | md5sum | head -c 16`
        tempdir="{COMPUTEDIR}/${{prefix}}"
        mkdir -p "$tempdir"
        
        # Generate STAR index
        STAR --runMode genomeGenerate \
             --genomeDir "$tempdir" \
             --genomeFastaFiles {WORKDIR}/{input.genome} \
             --sjdbGTFfile {input.gff} \
             --sjdbOverhang {params.overhang} \
             --runThreadN {threads} \
             2> {log}
        
        # Move results to final location
        mv "$tempdir" {WORKDIR}/{output.dir}) 2>> {log}
        """