main(){
    PYTHON_PATH=python3.4
    READEMPTION_PATH=../bin/reademption
    SEGEMEHL_PATH=../bin/segemehl.x
    LACK_PATH=../bin/lack.x
    FASTQ_DUMP_PATH=../bin/fastq-dump
    FASTQ_TRIMMER_PATH=../bin/fastq_quality_trimmer
    FASTQ_TO_FASTA_PATH=../bin/fastq_to_fasta
    SRA_FOLDER=Libs_as_SRA
    FASTQ_FOLDER=Libs_as_FASTQ
    READEMPTION_FOLDER=READemption_analysis
    ACCURACY=90
    MIN_MAPPING_LEN=20
    PROCESSORS_MAPPING=40
    PROCESSORS_GENE_QUANTI=5
    PLOTTING_FOLDER=correlation_plots

    create_folders
    set_up_analysis_folder
    link_reference_sequences
    link_annotations
    download_sra_files
    convert_to_fastq
    convert_to_quality_trimmed_fasta
    run_read_alignment
    run_coverage_file_generation
    run_gene_quantification
    change_reademption_folders_permission
    generate_correlation_plots
}

create_folders(){
    mkdir -p bin $SRA_FOLDER $FASTQ_FOLDER $PLOTTING_FOLDER
}

set_up_analysis_folder(){
    $PYTHON_PATH $READEMPTION_PATH create $READEMPTION_FOLDER
}

download_sra_files(){
    BASE_URL="ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX/"

    for LIB in \
	GSM1655401,SRX986/SRX986623/SRR1964702/SRR1964702.sra,HeLa-S3:WT:16_h:total_RNA \
	GSM1655402,SRX986/SRX986624/SRR1964703/SRR1964703.sra,HeLa-S3:WT:16_h:rRNA_free
    do
	GEO_ACCESSION=$(echo $LIB | cut -f 1 -d",")
	FILE_PATH=$(echo $LIB | cut -f 2 -d",")
	LIB_NAME=$(echo $LIB | cut -f 3 -d",")

	wget \
	    -O $SRA_FOLDER/${GEO_ACCESSION}_${LIB_NAME}.sra \
	    -c \
	    ${BASE_URL}/${FILE_PATH}
    done
}

convert_to_fastq(){
    for SRA_FILE in $(ls $SRA_FOLDER)
    do
	COMPRESSED_FASTQ_FILE=$(echo $SRA_FILE \
	    | sed -e "s/GSM......._//" -e "s/.sra/.fq.bz2/")
	$FASTQ_DUMP_PATH \
	    --stdout \
	    $SRA_FOLDER/$SRA_FILE \
	    | bzip2 \
	    --stdout - \
	    > ${FASTQ_FOLDER}/$COMPRESSED_FASTQ_FILE &
    done
    wait
}

convert_to_quality_trimmed_fasta(){
    for FASTQ_FILE in $(ls $FASTQ_FOLDER)
    do
	FASTA_FILE=$(echo $FASTQ_FILE | sed "s/.fq.bz2/.fa.bz2/")
	bunzip2 --stdout $FASTQ_FOLDER/$FASTQ_FILE \
	    | $FASTQ_TRIMMER_PATH -t 20 -Q33 -l 1 \
	    | $FASTQ_TO_FASTA_PATH -Q33 \
	    | bzip2 --stdout - > $READEMPTION_FOLDER/input/reads/$FASTA_FILE &
    done
    wait
}

link_reference_sequences(){
    ln -s ../../../../reference_sequences_and_annotations/output/hs37d5.fa \
	$READEMPTION_FOLDER/input/reference_sequences/
    ln -s ../../../../reference_sequences_and_annotations/output/NC_016810.fa \
	$READEMPTION_FOLDER/input/reference_sequences/
    ln -s ../../../../reference_sequences_and_annotations/output/NC_017718.fa \
	$READEMPTION_FOLDER/input/reference_sequences/
    ln -s ../../../../reference_sequences_and_annotations/output/NC_017719.fa \
	$READEMPTION_FOLDER/input/reference_sequences/
    ln -s ../../../../reference_sequences_and_annotations/output/NC_017720.fa \
	$READEMPTION_FOLDER/input/reference_sequences/
}

link_annotations(){
    ln -s ../../../../reference_sequences_and_annotations/output/Human.gff \
	$READEMPTION_FOLDER/input/annotations/
    ln -s ../../../../reference_sequences_and_annotations/output/Salmonella.gff \
	$READEMPTION_FOLDER/input/annotations/
}

run_read_alignment(){
    SEQ_ID_STRING=$(grep ">" $READEMPTION_FOLDER/input/reference_sequences/hs37d5.fa \
	| cut -f 1 -d" " | sed "s/>//" | tr "\n" "," | sed "s/,$//")
    CROSSMAPPING_PARAMETERS="-x Salmonella:NC_016810.1,NC_017718.1,NC_017719.1,NC_017720.1;Human:$SEQ_ID_STRING"

    $PYTHON_PATH $READEMPTION_PATH \
        align \
	--split \
	-r \
        -p ${PROCESSORS_MAPPING} \
        -a $ACCURACY \
        -l $MIN_MAPPING_LEN \
	--poly_a_clipping \
	--progress \
	$CROSSMAPPING_PARAMETERS \
        $READEMPTION_FOLDER
}

run_gene_quantification(){
    $PYTHON_PATH $READEMPTION_PATH \
        gene_quanti \
        -p ${PROCESSORS_GENE_QUANTI} \
    	--unique_only \
    	--no_count_splitting_by_gene_no \
    	--no_count_split_by_alignment_no \
        $READEMPTION_FOLDER
}

run_coverage_file_generation(){
    $PYTHON_PATH $READEMPTION_PATH \
    coverage \
        -p ${PROCESSORS_MAPPING} \
        $READEMPTION_FOLDER
}

change_reademption_folders_permission(){
    chmod -R a-w $READEMPTION_FOLDER
}

generate_correlation_plots(){
    grep -Ev "NC_01|^anti-sense" \
    	READemption_analysis/output/gene_quanti/gene_quanti_combined/gene_wise_quantifications_combined_tnoar.csv \
    	> $PLOTTING_FOLDER/gene_wise_quantifications_combined_tnoar_human.csv 

    grep -E "NC_01|^Orien" READemption_analysis/output/gene_quanti/gene_quanti_combined/gene_wise_quantifications_combined_tnoar.csv \
    	| grep -v "^anti-sense" \
    	> $PLOTTING_FOLDER/gene_wise_quantifications_combined_tnoar_Salmonella.csv

    $PYTHON_PATH bin/plot_correllation.py \
    	$PLOTTING_FOLDER/gene_wise_quantifications_combined_tnoar_Salmonella.csv \
    	"Salmonella" \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Salmonella_scatterplot.pdf \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Salmonella_histogram.csv \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Salmonella_rRNA_list.csv \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Salmonella_top_changed_non_rRNAs.csv

    $PYTHON_PATH bin/plot_correllation.py \
    	$PLOTTING_FOLDER/gene_wise_quantifications_combined_tnoar_human.csv \
    	"Human" \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Human_scatterplot.pdf \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Human_histogram.csv \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Human_rRNA_list.csv \
    	$PLOTTING_FOLDER/RiboZero_untreated_vs_treated_Human_top_changed_non_rRNAs.csv
}

main
