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

    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
}

download_sra(){
    BASE_URL="ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/"
    for LIB in \
	GSM1655351,SRX/SRX986/SRX986573/SRR1964650/SRR1964650.sra,HeLa-S3:WT_LB:00_h:R1 \
	GSM1655352,SRX/SRX986/SRX986574/SRR1964651/SRR1964651.sra,HeLa-S3:WT_input:00_h:R1 \
	GSM1655353,SRX/SRX986/SRX986575/SRR1964652/SRR1964652.sra,HeLa-S3:mock_GFP_minus:04_h:R1 \
	GSM1655354,SRX/SRX986/SRX986576/SRR1964653/SRR1964653.sra,HeLa-S3:WT_GFP_plus:04_h:R1 \
	GSM1655355,SRX/SRX986/SRX986577/SRR1964654/SRR1964654.sra,HeLa-S3:WT_GFP_minus:04_h:R1 \
	GSM1655356,SRX/SRX986/SRX986578/SRR1964655/SRR1964655.sra,HeLa-S3:mock_GFP_minus:24_h:R1 \
	GSM1655357,SRX/SRX986/SRX986579/SRR1964656/SRR1964656.sra,HeLa-S3:WT_GFP_plus:24_h:R1 \
	GSM1655358,SRX/SRX986/SRX986580/SRR1964657/SRR1964657.sra,HeLa-S3:WT_GFP_minus:24_h:R1 \
	GSM1655359,SRX/SRX986/SRX986581/SRR1964658/SRR1964658.sra,HeLa-S3:WT_LB:00_h:R2 \
	GSM1655360,SRX/SRX986/SRX986582/SRR1964659/SRR1964659.sra,HeLa-S3:WT_input:00_h:R2 \
	GSM1655361,SRX/SRX986/SRX986583/SRR1964660/SRR1964660.sra,HeLa-S3:mock_GFP_minus:04_h:R2 \
	GSM1655362,SRX/SRX986/SRX986584/SRR1964661/SRR1964661.sra,HeLa-S3:WT_GFP_plus:04_h:R2 \
	GSM1655363,SRX/SRX986/SRX986585/SRR1964662/SRR1964662.sra,HeLa-S3:WT_GFP_minus:04_h:R2 \
	GSM1655364,SRX/SRX986/SRX986586/SRR1964663/SRR1964663.sra,HeLa-S3:mock_GFP_minus:24_h:R2 \
	GSM1655365,SRX/SRX986/SRX986587/SRR1964664/SRR1964664.sra,HeLa-S3:WT_GFP_plus:24_h:R2 \
	GSM1655366,SRX/SRX986/SRX986588/SRR1964665/SRR1964665.sra,HeLa-S3:WT_GFP_minus:24_h:R2 \
	GSM1655367,SRX/SRX986/SRX986589/SRR1964666/SRR1964666.sra,HeLa-S3:WT_LB:00_h:R3 \
	GSM1655368,SRX/SRX986/SRX986590/SRR1964667/SRR1964667.sra,HeLa-S3:WT_input:00_h:R3 \
	GSM1655369,SRX/SRX986/SRX986591/SRR1964668/SRR1964668.sra,HeLa-S3:mock_GFP_minus:04_h:R3 \
	GSM1655370,SRX/SRX986/SRX986592/SRR1964669/SRR1964669.sra,HeLa-S3:WT_GFP_plus:04_h:R3 \
	GSM1655371,SRX/SRX986/SRX986593/SRR1964670/SRR1964670.sra,HeLa-S3:WT_GFP_minus:04_h:R3 \
	GSM1655372,SRX/SRX986/SRX986594/SRR1964671/SRR1964671.sra,HeLa-S3:mock_GFP_minus:24_h:R3 \
	GSM1655373,SRX/SRX986/SRX986595/SRR1964672/SRR1964672.sra,HeLa-S3:WT_GFP_plus:24_h:R3 \
	GSM1655374,SRX/SRX986/SRX986596/SRR1964673/SRR1964673.sra,HeLa-S3:WT_GFP_minus:24_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_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
}

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-10,17-28 \
    	> 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" \
	| cut -f 1-22,29-34 \
    	> $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
