Snakefile-features

###########################################################################################################################
# EXTRACTING FEATURES
###########################################################################################################################

##########
## SET UP 
##########

configfile: "config.yaml"

all_rna = [
    "functional-protein-exon2", 
    "functional-protein-exon3", 
    "functional-lncrna-exon1", 
    "functional-lncrna-exon2", 
    "functional-short-ncrna", 
    "protein-exon2-negative-control", 
    "protein-exon3-negative-control", 
    "lncrna-exon1-negative-control", 
    "lncrna-exon2-negative-control", 
    "short-ncrna-negative-control"
]

xxx = [
    "protein-exon2", 
    "protein-exon3", 
    "lncrna-exon1", 
    "lncrna-exon2", 
    "short-ncrna"
]

rna_types = [
    "protein", 
    "lncrna", 
    "short-ncrna"
]

###################
## RUN ALL FEATURES 
###################
rule all_features:
    input:
        expand("data/{sample_final}-features.csv", sample_final=rna_types)


#############
## Intrinsic features 
#############
# GC percentage: divide number of guanine and cytosine nucleotides by total sequence length
#
# Sequence Low complexity Density: Find regions of low sequence complexity with dustmasker (XXX reference and version!). 
# Density is calculated as the ratio between nucleotides in the low complexity regions and total nucleotides on sequence
                                  
rule run_intrinsic_seq_features:
    input:
        "data/datasets/{sample}-dataset.csv",
		"data/datasets/{sample}-seq.fa"
    
    output:
        "data/intrinsic/{sample}-intrinsic.csv"
    
    params:
        sample=all_rna
    
    shell:
        "scripts/B0-intrinsic-seq-features.sh {input} > {output}"


###############
## Conservation 
###############
# The maximum phyloP score for each sequence is extracted from its corresponding bigWig files using ```bigWigSummary``` (Haeussler et al., 2019). 
# # Parameters used for running bigWigSummary
# -type=max reports the maximum score
# $chr $start $end are the chromosome coordiantes for the sequence extracted from the input dataset (initial_data)
# 1 causes the score across the whole sequence to be reported

rule run_seq_conservation_features:
    input:
        initial_data="data/datasets/{sample}-dataset.csv",
        bigWigSummary_exe=config["bigWigSummary_exe"],
        _100w_phyloP_bw=config["_100w_phyloP_bw"],
		zoonomia_phyloP_bw=config["zoonomia_phyloP_bw"] 

    output:
        "data/conservation/{sample}-conservation.csv"
    
    params:
        sample=all_rna
    
    shell:
        "scripts/B1-seq-conservation-features.sh {input} > {output}"


################
## Transcriptome
################
# Read counts were calculated using ```samtools view -c``` and maximum read depth was calculated using ```samtools depth -r``` (Li et al., 2009).

rule run_transcriptome:
    input:
        "data/datasets/{sample}-dataset.csv", 
        encode_folder=config["encode_folder"]

    output:
        "data/transcriptome/{sample}-transcriptome.csv"
    
    params:
        sample=all_rna
    
    shell:
        "scripts/B2-transcriptome-expression-features.sh {input} > {output}"


##################        
## Genomic repeats 
##################
# Parameters used to run mmseqs
# Distances to the nearest repetitive elements were calculated from a file of non-redundant hits of repetitive DNA elements in the human genome obtained from Dfam v3.1 (Hubley et al., 2016). Once this was converted to bed file format, the distance to the nearest hit was calculated using ```bedtools closest```  (Quinlan, 2014). 
# # Parameters used for running bedtools:
# bedtools closest reports the either overlapping or the nearest sequence in input 
# -a $input_bed is the sequences of interest, converted to a bed file
# -b dfam-hg38-nrph.bed is the Dfam non-redundent repetitive DNA elements
# -io ignores overlaps, as some functional ncRNAs are considered repetitive DNA elements
# -D ref reports distance with respect to the reference genome
# -iu ignores upstream to report closest downstream element
# -id ignores downstream to report closest upstream element


rule run_genomic_repeat:
    input:
        initial_data="data/datasets/{sample}-dataset.csv",
        initial_fasta="data/datasets/{sample}-seq.fa", 
        human_genome_mmseqs=config["human_genome_mmseqs"],
        dfam_hits=config["dfam_hits"]   

    output:
        "data/repeats/{sample}-repeats.csv"

    params:
        sample=all_rna

    shell:
        "scripts/B3-genomic-repeat-associated-features.sh {input} > {output}"


###########################        
## Protein and RNA specific 
###########################
rule run_protein_and_rna_specific_features:
    input:
        initial_data="data/datasets/{sample}-dataset.csv",
        initial_fasta="data/datasets/{sample}-seq.fa",
        interaction_database=config["interaction_database"],
        bigBedToBed_exe=config["bigBedToBed_exe"]
        
        
    output:
        "data/specific/{sample}-specific.csv"

    params:
        sample=all_rna

    shell:
        "scripts/B4-protein-and-rna-specific-features.sh {input} > {output}"


#######################        
## Population variation
#######################
# SNP frequencies, which were used to caluculate the average minor allele frequency, were calculated using ```vcftools``` (Danecek et al., 2011). All chromosome SNP data from the Genome Aggregation Database (gnomAD) v3 was downloaded as a local copy, with the 
# SNPs for each region extracted using ```tabix -f -h``` (Karczewski et al. 2020; Li et al., 2009).

rule run_populatio_variation:
    input:
        initial_data="data/datasets/{sample}-dataset.csv",
        gnomad_directory=config["gnomad_directory"]

    output:
        "data/population/{sample}-population.csv"

    params:
        sample=all_rna

    shell:
        "scripts/B5-population-variation-features.sh {input} > {output}"


########################        
## Histone modifications 
########################
rule run_histone_features:
    input:
        expand("data/datasets/{sample}-dataset.csv", sample=all_rna, histone_name=all_histones)#,
        #histone_marks_directory=config["histone_marks_directory"]

    output:
        "data/histone_feature/{histone_name}/{sample}-histone-feature.csv"
        
    params:
        sample=all_rna,
        histone_name=all_histones

    shell:
        "scripts/histone_processing.sh {input} {histone_name} {output}"


##########################       
## Chromatin accessibility 
##########################
rule run_chromatin_accessibility_features:
    input:
        expand("data/datasets/{sample}-dataset.csv", sample=all_rna)#,
        #histone_marks_directory=config["chromatin_accessibility_directory"]

    output:
        "data/chrm_acc_feature/{sample}-chrm_acc-feature.csv"
        
    params:
        sample=all_rna

    shell:
        "scripts/chrm_acc_processing.sh {input} {output}"


################################        
## Methylome
################################
rule run_methylome_features:
    input:
        expand("data/datasets/{sample}-dataset.csv", sample=all_rna)#,
        #histone_marks_directory=config["methylome_directory"]

    output:
        "data/methylome_feature/{sample}-methylome-feature.csv"
        
    params:
        sample=all_rna

    shell:
        "scripts/methylome_processing.sh {input} {output}"


################################        
## Dinucleotide Frequencies
################################
rule run_dinucleotide_features:
    input:
        initial_data="data/datasets/{sample}-dataset.csv"

    output:
        "data/dinucleotide_feature/{sample}-dinucleotide-feature.csv"
        
    params:
        sample=all_rna

    shell:
        "scripts/B0.1_Dinucleotide_frequencies.sh {input} {output}"


#### FORMAT FINAL FEATURES DATASET #### 
rule concatenate_features:
    input:
        "data/datasets/{sample}-dataset.csv",
        "data/intrinsic/{sample}-intrinsic.csv",
        "data/conservation/{sample}-conservation.csv",
        #"data/transcriptome/{sample}-transcriptome.csv",
        "data/repeats/{sample}-repeats.csv",
        "data/specific/{sample}-specific.csv", 
        "data/population/{sample}-population.csv",
        #expand("data/histone_feature/{histone_name}_{sample}-histone-feature.csv", sample=all_rna, histone_name=all_histones),
        #"data/chrm_acc_feature/chrm_acc_{sample}-chrm_acc-feature.csv",
        #"data/methylome_feature/{sample}-methylome-feature.csv",
        "data/dinucleotide_feature/{sample}-dinucleotide-feature.csv"
    output:
        "data/{sample}-dataset-features.csv"
    
    params:
        sample=all_rna, histone_name=all_histones

    shell:
        "paste -d ',' {input} > {output}"


# Removes column "Distance" and header from negative control files
# Joins functional file with its corresponding negative control
# Removes columns "ID","Chromosome","End", "Sequence", and "GeneId"/"DistanceGene"       
rule join_positive_negative:
    input:
        positive="data/functional-{sample_dataset}-dataset-features.csv",
        negative="data/{sample_dataset}-negative-control-dataset-features.csv"
        
    output:
        "data/{sample_dataset}-features.csv"

    params:
        sample_dataset=all_rna
    
    shell:
        """
        awk 'NR>1' {input.negative} > {output}.tmp &&
        cat {input.positive} {output}.tmp > {output} &&
        cut -d',' -f2,8- {output} > {output}.final &&                             
        mv {output}.final {output} &&
        rm {output}.tmp
        """

# Merges exons datasets 
rule join_protein_exons:
	input:
		exon2="data/protein-exon2-features.csv",
		exon3="data/protein-exon3-features.csv"
    
	output: 
		"data/protein-features.csv"
    
	shell:
		"""
    	cat {input.exon2} {input.exon3} > {output}
		"""

rule join_lncrna_exons:
	input:
		exon1="data/lncrna-exon1-features.csv",
		exon2="data/lncrna-exon2-features.csv"
    
	output: 
		"data/lncrna-features.csv"
    
	shell:
		"""
		cat {input.exon1} {input.exon2} > {output}
		"""