main(){
    PYTHON_PATH=python3.4
    RSCRIPT_PATH=Rscript
    READEMPTION_PATH=../bin/reademption
    SEGEMEHL_PATH=../bin/segemehl.x
    LACK_PATH=../bin/lack.x
    FASTQ_DUMP_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
    EDGER_SALMONELLA_FOLDER=edgeR_Salmonella
    EDGER_HUMAN_FOLDER=edgeR_Human
    PLOTTING_FOLDER=correlation_plots

    create_folders
    set_up_analysis_folder
    link_reference_sequences
    link_annotations
    download_sra
    convert_to_fastq
    convert_to_quality_trimmed_fasta
    run_read_alignment
    run_gene_quantification
    change_reademption_folders_permission
    run_salmonella_edger
    run_human_edger
}

create_folders(){
    mkdir -p bin $SRA_FOLDER $FASTQ_FOLDER $EDGER_SALMONELLA_FOLDER \
	$EDGER_HUMAN_FOLDER
}

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

download_sra(){
    BASE_URL="ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX/"
    for LIB in \
	GSM1967950,SRX147/SRX1470939/SRR2981489/SRR2981489.sra,HeLa-S3:WT:00_h:R1 \
	GSM1967932,SRX147/SRX1470921/SRR2981471/SRR2981471.sra,HeLa-S3:DeltaryhB:00_h:R1 \
	GSM1967923,SRX147/SRX1470912/SRR2981462/SRR2981462.sra,HeLa-S3:Deltaomr:00_h:R1 \
	GSM1967941,SRX147/SRX1470930/SRR2981480/SRR2981480.sra,HeLa-S3:Deltayrl:00_h:R1 \
	GSM1967959,SRX147/SRX1470948/SRR2981498/SRR2981498.sra,HeLa-S3:mock:04_h:R1 \
	GSM1967953,SRX147/SRX1470942/SRR2981492/SRR2981492.sra,HeLa-S3:WT:04_h:R1 \
	GSM1967935,SRX147/SRX1470924/SRR2981474/SRR2981474.sra,HeLa-S3:DeltaryhB:04_h:R1 \
	GSM1967926,SRX147/SRX1470915/SRR2981465/SRR2981465.sra,HeLa-S3:Deltaomr:04_h:R1 \
	GSM1967944,SRX147/SRX1470933/SRR2981483/SRR2981483.sra,HeLa-S3:Deltayrl:04_h:R1 \
	GSM1967962,SRX147/SRX1470951/SRR2981501/SRR2981501.sra,HeLa-S3:mock:16_h:R1 \
	GSM1967956,SRX147/SRX1470945/SRR2981495/SRR2981495.sra,HeLa-S3:WT:16_h:R1 \
	GSM1967938,SRX147/SRX1470927/SRR2981477/SRR2981477.sra,HeLa-S3:DeltaryhB:16_h:R1 \
	GSM1967929,SRX147/SRX1470918/SRR2981468/SRR2981468.sra,HeLa-S3:Deltaomr:16_h:R1 \
	GSM1967947,SRX147/SRX1470936/SRR2981486/SRR2981486.sra,HeLa-S3:Deltayrl:16_h:R1 \
	GSM1967951,SRX147/SRX1470940/SRR2981490/SRR2981490.sra,HeLa-S3:WT:00_h:R2 \
	GSM1967933,SRX147/SRX1470922/SRR2981472/SRR2981472.sra,HeLa-S3:DeltaryhB:00_h:R2 \
	GSM1967924,SRX147/SRX1470913/SRR2981463/SRR2981463.sra,HeLa-S3:Deltaomr:00_h:R2 \
	GSM1967942,SRX147/SRX1470931/SRR2981481/SRR2981481.sra,HeLa-S3:Deltayrl:00_h:R2 \
	GSM1967960,SRX147/SRX1470949/SRR2981499/SRR2981499.sra,HeLa-S3:mock:04_h:R2 \
	GSM1967954,SRX147/SRX1470943/SRR2981493/SRR2981493.sra,HeLa-S3:WT:04_h:R2 \
	GSM1967936,SRX147/SRX1470925/SRR2981475/SRR2981475.sra,HeLa-S3:DeltaryhB:04_h:R2 \
	GSM1967927,SRX147/SRX1470916/SRR2981466/SRR2981466.sra,HeLa-S3:Deltaomr:04_h:R2 \
	GSM1967945,SRX147/SRX1470934/SRR2981484/SRR2981484.sra,HeLa-S3:Deltayrl:04_h:R2 \
	GSM1967963,SRX147/SRX1470952/SRR2981502/SRR2981502.sra,HeLa-S3:mock:16_h:R2 \
	GSM1967957,SRX147/SRX1470946/SRR2981496/SRR2981496.sra,HeLa-S3:WT:16_h:R2 \
	GSM1967939,SRX147/SRX1470928/SRR2981478/SRR2981478.sra,HeLa-S3:DeltaryhB:16_h:R2 \
	GSM1967930,SRX147/SRX1470919/SRR2981469/SRR2981469.sra,HeLa-S3:Deltaomr:16_h:R2 \
	GSM1967948,SRX147/SRX1470937/SRR2981487/SRR2981487.sra,HeLa-S3:Deltayrl:16_h:R2 \
	GSM1967952,SRX147/SRX1470941/SRR2981491/SRR2981491.sra,HeLa-S3:WT:00_h:R3 \
	GSM1967934,SRX147/SRX1470923/SRR2981473/SRR2981473.sra,HeLa-S3:DeltaryhB:00_h:R3 \
	GSM1967925,SRX147/SRX1470914/SRR2981464/SRR2981464.sra,HeLa-S3:Deltaomr:00_h:R3 \
	GSM1967943,SRX147/SRX1470932/SRR2981482/SRR2981482.sra,HeLa-S3:Deltayrl:00_h:R3 \
	GSM1967961,SRX147/SRX1470950/SRR2981500/SRR2981500.sra,HeLa-S3:mock:04_h:R3 \
	GSM1967955,SRX147/SRX1470944/SRR2981494/SRR2981494.sra,HeLa-S3:WT:04_h:R3 \
	GSM1967937,SRX147/SRX1470926/SRR2981476/SRR2981476.sra,HeLa-S3:DeltaryhB:04_h:R3 \
	GSM1967928,SRX147/SRX1470917/SRR2981467/SRR2981467.sra,HeLa-S3:Deltaomr:04_h:R3 \
	GSM1967946,SRX147/SRX1470935/SRR2981485/SRR2981485.sra,HeLa-S3:Deltayrl:04_h:R3 \
	GSM1967964,SRX147/SRX1470953/SRR2981503/SRR2981503.sra,HeLa-S3:mock:16_h:R3 \
	GSM1967958,SRX147/SRX1470947/SRR2981497/SRR2981497.sra,HeLa-S3:WT:16_h:R3 \
	GSM1967940,SRX147/SRX1470929/SRR2981479/SRR2981479.sra,HeLa-S3:DeltaryhB:16_h:R3 \
	GSM1967931,SRX147/SRX1470920/SRR2981470/SRR2981470.sra,HeLa-S3:Deltaomr:16_h:R3 \
	GSM1967949,SRX147/SRX1470938/SRR2981488/SRR2981488.sra,HeLa-S3:Deltayrl:16_h:R3
    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_BIN \
	    --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/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
}

run_salmonella_edger(){
    # Extract gene names and locus tag colums
    grep -Pv "\trRNA\t|\ttRNA\t" ../reference_sequences_and_annotations/output/Salmonella.csv \
    	| cut -f 10,11 \
    	> tmp_names_and_locus_tags

    # Extract Salmonella countings for non-mock libraries
    grep -P "^Orienta|\tNC_01" $READEMPTION_FOLDER/output/gene_quanti/gene_quanti_combined/gene_wise_quantifications_combined.csv \
    	| grep -v "^anti-sense" \
    	| grep -Pv "\trRNA\t|\ttRNA\t" \
    	| cut -f 1-46 \
    	> tmp

    paste tmp tmp_names_and_locus_tags \
    	> $EDGER_SALMONELLA_FOLDER/Salmonella_gene_wise_read_countings.csv

    rm tmp tmp_names_and_locus_tags

    $RSCRIPT_PATH bin/edgeR_Salmonella.R \
    	$EDGER_SALMONELLA_FOLDER/Salmonella_gene_wise_read_countings.csv \
    	$EDGER_SALMONELLA_FOLDER/
}

run_human_edger(){
    grep -Pv "\tNC_01" $READEMPTION_FOLDER/output/gene_quanti/gene_quanti_combined/gene_wise_quantifications_combined.csv \
        | grep -v "^anti-sense" \
        | grep -Ev "=rRNA|=Mt_rRNA|=tRNA|=Mt_tRNA" \
    	> $EDGER_HUMAN_FOLDER/Human_gene_wise_read_countings.csv
    
    $RSCRIPT_PATH bin/edgeR_Human.R \
        $EDGER_HUMAN_FOLDER/Human_gene_wise_read_countings.csv \
        $EDGER_HUMAN_FOLDER/
}

main
