main(){
    PYTHON_PATH=python3.4
    PERL=perl
    RSCRIPT_PATH=Rscript
    READEMPTION_PATH=../bin/reademption
    SEGEMEHL_PATH=../bin/segemehl.x
    LACK_PATH=../bin/lack.x
    FASTQ_TRIMMER_PATH=../bin/fastq_quality_trimmer
    FASTQ_TO_FASTA_PATH=../bin/fastq_to_fasta
    FASTQ_DUMP_PATH=../bin/fastq-dump

    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
    PLOT_FOLDER=plots
    DIFFEXP_PLOTS=${PLOT_FOLDER}/differential_expression
    INF_VS_MOCK_PLOTS=${PLOT_FOLDER}/infected_vs_mock

    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
    plot_salmonella_srna_kinetics
    plot_salmonella_regulon_correlation 
    run_human_edger
    generate_infected_vs_mock_plots
    generate_diff_exp_plots
    run_interspecies_correlation
}

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

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/SRX/"
    for LIB in \
	GSM1924772,SRX141/SRX1410954/SRR2936886/SRR2936886.sra,HeLa-S3:WT:00_h:R1 \
	GSM1924754,SRX141/SRX1410936/SRR2936868/SRR2936868.sra,HeLa-S3:DeltapinT:00_h:R1 \
	GSM1924790,SRX141/SRX1410972/SRR2936904/SRR2936904.sra,HeLa-S3:mock:02_h:R1 \
	GSM1924775,SRX141/SRX1410957/SRR2936889/SRR2936889.sra,HeLa-S3:WT:02_h:R1 \
	GSM1924757,SRX141/SRX1410939/SRR2936871/SRR2936871.sra,HeLa-S3:DeltapinT:02_h:R1 \
	GSM1924793,SRX141/SRX1410976/SRR2936907/SRR2936907.sra,HeLa-S3:mock:04_h:R1 \
	GSM1924778,SRX141/SRX1410960/SRR2936892/SRR2936892.sra,HeLa-S3:WT:04_h:R1 \
	GSM1924760,SRX141/SRX1410942/SRR2936874/SRR2936874.sra,HeLa-S3:DeltapinT:04_h:R1 \
	GSM1924796,SRX141/SRX1410979/SRR2936910/SRR2936910.sra,HeLa-S3:mock:08_h:R1 \
	GSM1924781,SRX141/SRX1410963/SRR2936895/SRR2936895.sra,HeLa-S3:WT:08_h:R1 \
	GSM1924763,SRX141/SRX1410945/SRR2936877/SRR2936877.sra,HeLa-S3:DeltapinT:08_h:R1 \
	GSM1924799,SRX141/SRX1410982/SRR2936913/SRR2936913.sra,HeLa-S3:mock:16_h:R1 \
	GSM1924784,SRX141/SRX1410966/SRR2936898/SRR2936898.sra,HeLa-S3:WT:16_h:R1 \
	GSM1924766,SRX141/SRX1410948/SRR2936880/SRR2936880.sra,HeLa-S3:DeltapinT:16_h:R1 \
	GSM1924802,SRX141/SRX1410985/SRR2936916/SRR2936916.sra,HeLa-S3:mock:24_h:R1 \
	GSM1924787,SRX141/SRX1410969/SRR2936901/SRR2936901.sra,HeLa-S3:WT:24_h:R1 \
	GSM1924769,SRX141/SRX1410951/SRR2936883/SRR2936883.sra,HeLa-S3:DeltapinT:24_h:R1 \
	GSM1924773,SRX141/SRX1410955/SRR2936887/SRR2936887.sra,HeLa-S3:WT:00_h:R2 \
	GSM1924755,SRX141/SRX1410937/SRR2936869/SRR2936869.sra,HeLa-S3:DeltapinT:00_h:R2 \
	GSM1924791,SRX141/SRX1410973/SRR2936905/SRR2936905.sra,HeLa-S3:mock:02_h:R2 \
	GSM1924776,SRX141/SRX1410958/SRR2936890/SRR2936890.sra,HeLa-S3:WT:02_h:R2 \
	GSM1924758,SRX141/SRX1410940/SRR2936872/SRR2936872.sra,HeLa-S3:DeltapinT:02_h:R2 \
	GSM1924794,SRX141/SRX1410977/SRR2936908/SRR2936908.sra,HeLa-S3:mock:04_h:R2 \
	GSM1924779,SRX141/SRX1410961/SRR2936893/SRR2936893.sra,HeLa-S3:WT:04_h:R2 \
	GSM1924761,SRX141/SRX1410943/SRR2936875/SRR2936875.sra,HeLa-S3:DeltapinT:04_h:R2 \
	GSM1924797,SRX141/SRX1410980/SRR2936911/SRR2936911.sra,HeLa-S3:mock:08_h:R2 \
	GSM1924782,SRX141/SRX1410964/SRR2936896/SRR2936896.sra,HeLa-S3:WT:08_h:R2 \
	GSM1924764,SRX141/SRX1410946/SRR2936878/SRR2936878.sra,HeLa-S3:DeltapinT:08_h:R2 \
	GSM1924800,SRX141/SRX1410983/SRR2936914/SRR2936914.sra,HeLa-S3:mock:16_h:R2 \
	GSM1924785,SRX141/SRX1410967/SRR2936899/SRR2936899.sra,HeLa-S3:WT:16_h:R2 \
	GSM1924767,SRX141/SRX1410949/SRR2936881/SRR2936881.sra,HeLa-S3:DeltapinT:16_h:R2 \
	GSM1924803,SRX141/SRX1410986/SRR2936917/SRR2936917.sra,HeLa-S3:mock:24_h:R2 \
	GSM1924788,SRX141/SRX1410970/SRR2936902/SRR2936902.sra,HeLa-S3:WT:24_h:R2 \
	GSM1924770,SRX141/SRX1410952/SRR2936884/SRR2936884.sra,HeLa-S3:DeltapinT:24_h:R2 \
	GSM1924774,SRX141/SRX1410956/SRR2936888/SRR2936888.sra,HeLa-S3:WT:00_h:R3 \
	GSM1924756,SRX141/SRX1410938/SRR2936870/SRR2936870.sra,HeLa-S3:DeltapinT:00_h:R3 \
	GSM1924792,SRX141/SRX1410974/SRR2936906/SRR2936906.sra,HeLa-S3:mock:02_h:R3 \
	GSM1924777,SRX141/SRX1410959/SRR2936891/SRR2936891.sra,HeLa-S3:WT:02_h:R3 \
	GSM1924759,SRX141/SRX1410941/SRR2936873/SRR2936873.sra,HeLa-S3:DeltapinT:02_h:R3 \
	GSM1924795,SRX141/SRX1410978/SRR2936909/SRR2936909.sra,HeLa-S3:mock:04_h:R3 \
	GSM1924780,SRX141/SRX1410962/SRR2936894/SRR2936894.sra,HeLa-S3:WT:04_h:R3 \
	GSM1924762,SRX141/SRX1410944/SRR2936876/SRR2936876.sra,HeLa-S3:DeltapinT:04_h:R3 \
	GSM1924798,SRX141/SRX1410981/SRR2936912/SRR2936912.sra,HeLa-S3:mock:08_h:R3 \
	GSM1924783,SRX141/SRX1410965/SRR2936897/SRR2936897.sra,HeLa-S3:WT:08_h:R3 \
	GSM1924765,SRX141/SRX1410947/SRR2936879/SRR2936879.sra,HeLa-S3:DeltapinT:08_h:R3 \
	GSM1924801,SRX141/SRX1410984/SRR2936915/SRR2936915.sra,HeLa-S3:mock:16_h:R3 \
	GSM1924786,SRX141/SRX1410968/SRR2936900/SRR2936900.sra,HeLa-S3:WT:16_h:R3 \
	GSM1924768,SRX141/SRX1410950/SRR2936882/SRR2936882.sra,HeLa-S3:DeltapinT:16_h:R3 \
	GSM1924804,SRX141/SRX1410987/SRR2936918/SRR2936918.sra,HeLa-S3:mock:24_h:R3 \
	GSM1924789,SRX141/SRX1410971/SRR2936903/SRR2936903.sra,HeLa-S3:WT:24_h:R3 \
	GSM1924771,SRX141/SRX1410953/SRR2936885/SRR2936885.sra,HeLa-S3:DeltapinT: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/Human.fa \
    # 	| cut -f 1 -d" " | sed "s/>//" | tr "\n" "," | sed "s/,$//")

    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 \
    	--segemehl_bin=${SEGEMEHL_PATH} \
        --lack_bin=${LACK_PATH} \
    	--split \
    	-r \
        -p ${PROCESSORS_MAPPING} \
        -a $ACCURACY \
        -l $MIN_MAPPING_LEN \
    	--poly_a_clipping \
    	--progress \
    	--adapter AGATCGGAAGAGCACACGTCTGAAC \
    	$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_WT.R \
    	$EDGER_SALMONELLA_FOLDER/Salmonella_gene_wise_read_countings.csv \
    	$EDGER_SALMONELLA_FOLDER/

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

plot_salmonella_srna_kinetics(){
    $PYTHON_PATH \
	../bin/plot_srna_kinetics.py \
	--output_pdf $PLOT_FOLDER/sRNA_kinetics_HeLa_time_course.pdf \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_02_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_04_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_08_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_16_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_24_h_vs_WT_00_h.csv

    $PYTHON_PATH \
	../bin/plot_srna_kinetics.py \
	--output_pdf $PLOT_FOLDER/sRNA_kinetics_HeLa_time_course_pinT_only.pdf \
	--mark_only_pinT \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_02_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_04_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_08_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_16_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_24_h_vs_WT_00_h.csv
}

plot_salmonella_regulon_correlation(){
    $PYTHON_PATH bin/calculate_correlations_based_on_edger_results.py \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_02_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_04_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_08_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_16_h_vs_WT_00_h.csv \
	$EDGER_SALMONELLA_FOLDER/edgeR-comparison-WT_24_h_vs_WT_00_h.csv \
    	> $PLOT_FOLDER/pearson_corr_coeff_of_all_gene_vs_pinT.csv

    $PYTHON_PATH bin/generate_corr_values_boxplots.py \
        --output_pdf $PLOT_FOLDER/Regulon_correlation_with_pinT_boxplots.pdf \
    	$PLOT_FOLDER/pearson_corr_coeff_of_all_gene_vs_pinT.csv \
    	input/Gene_name_translation_table.csv \
        input/regulons/$REGULON_FOLDER 
    
    python3 bin/generate_phop_heat_map.py \
        --output_pdf $PLOT_FOLDER/Heatmap_edgeR_foldchanges_PhoP_regulon.pdf \
    	input/regulons/PhoP_regulon.csv \
    	$PLOT_FOLDER/pearson_corr_coeff_of_all_gene_vs_pinT.csv
}

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-10,14-28,32-61 \
    	> $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/
}

generate_infected_vs_mock_plots(){
    $RSCRIPT_PATH \
	bin/plots_infected_vs_mock.R \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_mock-02_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_mock-04_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_mock-08_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_mock-16_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_mock-24_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-WT_vs_mock-02_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-WT_vs_mock-04_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-WT_vs_mock-08_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-WT_vs_mock-16_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-WT_vs_mock-24_h.csv \
	$INF_VS_MOCK_PLOTS
}

generate_diff_exp_plots(){
    $RSCRIPT_PATH bin/generate_heatmaps.R \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_WT-02_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_WT-04_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_WT-08_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_WT-16_h.csv \
	$EDGER_HUMAN_FOLDER/RUV-edgeR-DeltapinT_vs_WT-24_h.csv \
	$DIFFEXP_PLOTS
}

run_interspecies_correlation(){
    
    $PERL_PATH ./bin/RUV2table.pl \
	$EDGER_HUMAN_FOLDER/RUV-Normalized-WT_Mutant_02_h-24_h.csv \
	> data/gene_wise_quantifications_hsa_2h-24h.tsv

    $PERL_PATH ./bin/RUV2table.pl \
	$EDGER_SALMONELLA_FOLDER/RUV-Normalized-WT_and_DeltapinT_00_h-16_h.csv \
	> data/gene_wise_quantifications_bac_0h-16h.tsv

    ## extract mRNA gene expression values for bac and hsa
    $PERL_PATH ./bin/splitter.pl mRNA           ./data/gene_wise_quantifications_bac_0h-16h.tsv > gene_wise_quantifications_bac_mRNA.tsv
    $PERL_PATH ./bin/splitter.pl protein_coding ./data/gene_wise_quantifications_hsa_2h-24h.tsv > gene_wise_quantifications_hsa_mRNA.tsv

    ## calculate z-scores over all samples for each gene (remove all genes with SD < 2)
    $PERL_PATH ./bin/normalizer.pl gene_wise_quantifications_hsa_mRNA.tsv    > gene_wise_quantifications_hsa_mRNA_norm.tsv
    $PERL_PATH ./bin/normalizer.pl gene_wise_quantifications_bac_mRNA.tsv    > gene_wise_quantifications_bac_mRNA_norm.tsv

    ## correlate all bac genes with all hsa genes
    rsquared=0.4225 # equally r = 0.65
    pvalue=0.0001
    rm -f correlation__bac_vs_hsa.csv correlation__bac_vs_hsa_sig.csv

    for roi in $(cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print $F[0] unless ($.==1)')
    do
	cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print if ($.==1 || $F[0]=~m/'$roi'/)' > regulator.csv
	cat gene_wise_quantifications_hsa_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print unless ($F[0]=~m/'$roi'/)'      > regulated.csv

	echo 'raw_regulator <- read.table(file="regulator.csv", header=F, skip = 1, sep="\t")' >  r.script
	echo 'raw_targets   <- read.table(file="regulated.csv", header=F, skip = 1, sep="\t")' >> r.script
	echo 'target_names    <- raw_targets[,1]' >> r.script
	echo 'regulator_name  <- raw_regulator[,1]' >> r.script
	echo 'regulator <- data.matrix(raw_regulator[,3:length(raw_regulator[1,])])' >> r.script
	echo 'targets   <- data.matrix(raw_targets[,3:length(raw_targets[1,])])' >> r.script
	echo 'colnames(regulator) <- NULL' >> r.script
	echo 'colnames(targets)   <- NULL' >> r.script
	echo '' >> r.script
	echo 'get_corr <- function(x, r=regulator[1,]){' >> r.script
	echo '  cor.res <- cor.test(x, r, alternative = "two.sided", method ="pearson", conf.level = 0.95)' >> r.script
	echo '  c(cor.res$estimate, cor.res$p.value)' >> r.script
	echo '}' >> r.script
	echo '' >> r.script
	echo 'res_corr_m   <- t(apply(targets, MARGIN=1, FUN=get_corr ))' >> r.script
	echo 'res_corr     <- data.frame(gene = target_names, regulator = regulator_name, R = res_corr_m[,1], Rsquared = (res_corr_m[,1])^2, pval = res_corr_m[,2], padj = p.adjust(res_corr_m[,2], method = "fdr"))' >> r.script   
	echo "res_corr_sig <- subset(res_corr, (Rsquared > $rsquared & padj < $pvalue))" >> r.script
	echo 'if (file.exists("correlation__bac_vs_hsa.csv")){' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_hsa.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_hsa.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	echo 'if (file.exists("correlation__bac_vs_hsa_sig.csv")){' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_hsa_sig.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_hsa_sig.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	$RSCRIPT_PATH BATCH --no-save --no-restore r.script
	rm -f r.script regulator.csv regulated.csv
    done

    ## check for GO enrichment in each correlation clusters considering only positive correlation
    wget http://purl.obolibrary.org/obo/go/go-basic.obo
    obo=$(readlink -f go-basic.obo)
    wget http://geneontology.org/gene-associations/gene_association.goa_ref_human.gz
    anno=$(readlink -f gene_association.goa_ref_human.gz)

    for roi in $(cat correlation__bac_vs_hsa_sig.csv | $PERL_PATH -F"\t" -lane '$h{$F[1]}++;END{print join ("\n", grep {!m/regulator/} sort keys %h)}')
    do
	
	cat correlation__bac_vs_hsa.csv | $PERL_PATH -F"\t" -lane 'if($F[1] eq "'$roi'"){if($F[0]=~m/(ENSG\d+)/){print "$1"}elsif($.>1){print STDERR "Warning: $F[0]";}}' > population_ensembl.txt
	cat population_ensembl.txt | $PERL_PATH ./bin/ensemble2uniprot.pl > population_uniprot.txt
	path="Ontologizer/Cluster_${roi}"
	mkdir -p $path
	cat correlation__bac_vs_hsa_sig.csv | $PERL_PATH -F"\t" -lane 'if($F[1] eq "'$roi'" && $F[2]>0){if($F[0]=~m/(ENSG\d+)/){print "$1"}elsif($.>1){print STDERR "Warning: $F[0]";}}' | $PERL_PATH ./bin/ensemble2uniprot.pl > ${roi}.txt
	java -jar  ../bin/Ontologizer.jar --mtc Westfall-Young-Single-Step --calculation Term-For-Term --association $anno --go $obo --ignore --annotation --outdir $path --population population_uniprot.txt --resamplingsteps 5000 --studyset ${roi}.txt -d 0.1 --dot 2> /dev/null
	cat $path/table-${roi}-Term-For-Term-Westfall-Young-Single-Step.txt | $PERL_PATH -F"\t" -lane 'if($.==1){print "%$F[0]\t$F[6]"}else{print "$F[0]\t$F[6]" if($F[6]<0.05)}' > $path/table-${roi}-significant4revigo.csv
	cat $path/table-${roi}-Term-For-Term-Westfall-Young-Single-Step.txt | $PERL_PATH -F"\t" -lane 'print join(",",$F[0],$F[8],$F[6]) if ($F[6]<0.05)' > $path/simREVIGO__ALL_cluster_${roi}.csv
	
	mv ${roi}.txt $path
	rm -f population_ensembl.txt population_uniprot.txt
	
    done

    ## make correlation for all bac genes agains all bac genes
    rsquared=0.4225 # equally r = 0.65
    pvalue=0.0001
    rm -f correlation__bac_vs_bac.csv correlation__bac_vs_bac_sig.csv

    for roi in $(cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print $F[0] unless ($.==1)')
    do
	cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print if ($.==1 || $F[0]=~m/'$roi'/)' > regulator.csv
	cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print unless ($F[0]=~m/'$roi'/)'      > regulated.csv

	echo 'raw_regulator <- read.table(file="regulator.csv", header=F, skip = 1, sep="\t")' >  r.script
	echo 'raw_targets   <- read.table(file="regulated.csv", header=F, skip = 1, sep="\t")' >> r.script
	echo 'target_names    <- raw_targets[,1]' >> r.script
	echo 'regulator_name  <- raw_regulator[,1]' >> r.script
	echo 'regulator <- data.matrix(raw_regulator[,3:length(raw_regulator[1,])])' >> r.script
	echo 'targets   <- data.matrix(raw_targets[,3:length(raw_targets[1,])])' >> r.script
	echo 'colnames(regulator) <- NULL' >> r.script
	echo 'colnames(targets)   <- NULL' >> r.script
	echo '' >> r.script
	echo 'get_corr <- function(x, r=regulator[1,]){' >> r.script
	echo '  cor.res <- cor.test(x, r, alternative = "two.sided", method ="pearson", conf.level = 0.95)' >> r.script
	echo '  c(cor.res$estimate, cor.res$p.value)' >> r.script
	echo '}' >> r.script
	echo '' >> r.script
	echo 'res_corr_m   <- t(apply(targets, MARGIN=1, FUN=get_corr ))' >> r.script
	echo 'res_corr     <- data.frame(gene = target_names, regulator = regulator_name, R = res_corr_m[,1], Rsquared = (res_corr_m[,1])^2, pval = res_corr_m[,2], padj = p.adjust(res_corr_m[,2], method = "fdr"))' >> r.script   
	echo "res_corr_sig <- subset(res_corr, (Rsquared > $rsquared & padj < $pvalue))" >> r.script
	echo 'if (file.exists("correlation__bac_vs_bac.csv")){' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_bac.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_bac.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	echo 'if (file.exists("correlation__bac_vs_bac_sig.csv")){' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_bac_sig.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_bac_sig.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	$RSCRIPT_PATH r.script
	rm -f r.script regulator.csv regulated.csv
    done

    ## Plot heatmap of correlation data
    # hsa: all genes with R^2 > 0.4225 && p-value adjusted < 0.0001 && positive correlated
    # bac: all secreted genes with R^2 > 0.4225 && p-value adjusted < 0.0001 && positive correlated
    # linker gene: sseC (== YP_005181244)

    cat correlation__bac_vs_bac_sig.csv | $PERL_PATH -F"\t" -lane 'print if $F[1] eq "YP_005181244"' | $PERL_PATH -F"\t" -lane 'print if ($F[-1] < 0.0001 && $F[3] > '$rsquared' && $F[2] > 0)'  | awk '{print $1}' | grep -f data/secreted.grep > bac_genes_to_plot.list
    cat correlation__bac_vs_hsa_sig.csv | $PERL_PATH -F"\t" -lane 'print if $F[1] eq "YP_005181244"' | $PERL_PATH -F"\t" -lane 'print if ($F[-1] < 0.0001 && $F[3] > '$rsquared' && $F[2] > 0)'  | awk '{print $1}' > hsa_genes_to_plot.list

    cat gene_wise_quantifications_bac_mRNA_norm.tsv | head -n 1 | $PERL_PATH -F"\t" -lane '$F[1]="group"; print join("\t", @F)'                     >  data_to_heatmap.csv
    cat gene_wise_quantifications_bac_mRNA_norm.tsv | grep -f bac_genes_to_plot.list | $PERL_PATH -F"\t" -lane '$F[1]="bac"; print join("\t", @F)'  >> data_to_heatmap.csv
    cat gene_wise_quantifications_hsa_mRNA_norm.tsv | grep -f hsa_genes_to_plot.list | $PERL_PATH -F"\t" -lane '$F[1]="hsa"; print join("\t", @F)'  >> data_to_heatmap.csv
    cat gene_wise_quantifications_bac_mRNA_norm.tsv | grep  YP_005181244             | $PERL_PATH -F"\t" -lane '$F[1]="reg"; print join("\t", @F)'  >> data_to_heatmap.csv

    $PERL_PATH bin/translate_plotdata.pl data/bac_genbank2genename.dict data/hsa_ensemble2genename.dict data_to_heatmap.csv > data_to_heatmap_translated.csv

    $RSCRIPT_PATH bin/sseC_correlationheatmap.R
    rm -f sseC_correlationheatmap.Rout
    ## extract mRNA gene expression values for bac and hsa
    $PERL_PATH ./bin/splitter.pl mRNA           ./data/gene_wise_quantifications_bac_0h-16h.tsv > gene_wise_quantifications_bac_mRNA.tsv
    $PERL_PATH ./bin/splitter.pl protein_coding ./data/gene_wise_quantifications_hsa_2h-24h.tsv > gene_wise_quantifications_hsa_mRNA.tsv

    ## calculate z-scores over all samples for each gene (remove all genes with SD < 2)
    $PERL_PATH ./bin/normalizer.pl gene_wise_quantifications_hsa_mRNA.tsv    > gene_wise_quantifications_hsa_mRNA_norm.tsv
    $PERL_PATH ./bin/normalizer.pl gene_wise_quantifications_bac_mRNA.tsv    > gene_wise_quantifications_bac_mRNA_norm.tsv

    ## correlate all bac genes with all hsa genes
    rsquared=0.4225 # equally r = 0.65
    pvalue=0.0001
    rm -f correlation__bac_vs_hsa.csv correlation__bac_vs_hsa_sig.csv

    for roi in $(cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print $F[0] unless ($.==1)')
    do
	cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print if ($.==1 || $F[0]=~m/'$roi'/)' > regulator.csv
	cat gene_wise_quantifications_hsa_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print unless ($F[0]=~m/'$roi'/)'      > regulated.csv

	echo 'raw_regulator <- read.table(file="regulator.csv", header=F, skip = 1, sep="\t")' >  r.script
	echo 'raw_targets   <- read.table(file="regulated.csv", header=F, skip = 1, sep="\t")' >> r.script
	echo 'target_names    <- raw_targets[,1]' >> r.script
	echo 'regulator_name  <- raw_regulator[,1]' >> r.script
	echo 'regulator <- data.matrix(raw_regulator[,3:length(raw_regulator[1,])])' >> r.script
	echo 'targets   <- data.matrix(raw_targets[,3:length(raw_targets[1,])])' >> r.script
	echo 'colnames(regulator) <- NULL' >> r.script
	echo 'colnames(targets)   <- NULL' >> r.script
	echo '' >> r.script
	echo 'get_corr <- function(x, r=regulator[1,]){' >> r.script
	echo '  cor.res <- cor.test(x, r, alternative = "two.sided", method ="pearson", conf.level = 0.95)' >> r.script
	echo '  c(cor.res$estimate, cor.res$p.value)' >> r.script
	echo '}' >> r.script
	echo '' >> r.script
	echo 'res_corr_m   <- t(apply(targets, MARGIN=1, FUN=get_corr ))' >> r.script
	echo 'res_corr     <- data.frame(gene = target_names, regulator = regulator_name, R = res_corr_m[,1], Rsquared = (res_corr_m[,1])^2, pval = res_corr_m[,2], padj = p.adjust(res_corr_m[,2], method = "fdr"))' >> r.script   
	echo "res_corr_sig <- subset(res_corr, (Rsquared > $rsquared & padj < $pvalue))" >> r.script
	echo 'if (file.exists("correlation__bac_vs_hsa.csv")){' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_hsa.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_hsa.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	echo 'if (file.exists("correlation__bac_vs_hsa_sig.csv")){' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_hsa_sig.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_hsa_sig.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	$RSCRIPT_PATH r.script
	rm -f r.script* regulator.csv regulated.csv
    done

    ## check for GO enrichment in each correlation clusters considering only positive correlation
    wget http://purl.obolibrary.org/obo/go/go-basic.obo
    obo=$(readlink -f go-basic.obo)
    wget http://geneontology.org/gene-associations/gene_association.goa_ref_human.gz
    anno=$(readlink -f gene_association.goa_ref_human.gz)

    for roi in $(cat correlation__bac_vs_hsa_sig.csv | $PERL_PATH -F"\t" -lane '$h{$F[1]}++;END{print join ("\n", grep {!m/regulator/} sort keys %h)}')
    do
	
	cat correlation__bac_vs_hsa.csv | $PERL_PATH -F"\t" -lane 'if($F[1] eq "'$roi'"){if($F[0]=~m/(ENSG\d+)/){print "$1"}elsif($.>1){print STDERR "Warning: $F[0]";}}' > population_ensembl.txt
	cat population_ensembl.txt | $PERL_PATH ./bin/ensemble2uniprot.pl > population_uniprot.txt
	path="Ontologizer/Cluster_${roi}"
	mkdir -p $path
	cat correlation__bac_vs_hsa_sig.csv | $PERL_PATH -F"\t" -lane 'if($F[1] eq "'$roi'" && $F[2]>0){if($F[0]=~m/(ENSG\d+)/){print "$1"}elsif($.>1){print STDERR "Warning: $F[0]";}}' | $PERL_PATH ./bin/ensemble2uniprot.pl > ${roi}.txt
	java -jar  ../bin/Ontologizer.jar --mtc Westfall-Young-Single-Step --calculation Term-For-Term --association $anno --go $obo --ignore --annotation --outdir $path --population population_uniprot.txt --resamplingsteps 5000 --studyset ${roi}.txt -d 0.1 --dot 2> /dev/null
	cat $path/table-${roi}-Term-For-Term-Westfall-Young-Single-Step.txt | $PERL_PATH -F"\t" -lane 'if($.==1){print "%$F[0]\t$F[6]"}else{print "$F[0]\t$F[6]" if($F[6]<0.05)}' > $path/table-${roi}-significant4revigo.csv
	cat $path/table-${roi}-Term-For-Term-Westfall-Young-Single-Step.txt | $PERL_PATH -F"\t" -lane 'print join(",",$F[0],$F[8],$F[6]) if ($F[6]<0.05)' > $path/simREVIGO__ALL_cluster_${roi}.csv
	
	mv ${roi}.txt $path
	rm -f population_ensembl.txt population_uniprot.txt
	
    done

    ## make correlation for all bac genes agains all bac genes
    rsquared=0.4225 # equally r = 0.65
    pvalue=0.0001
    rm -f correlation__bac_vs_bac.csv correlation__bac_vs_bac_sig.csv

    for roi in $(cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print $F[0] unless ($.==1)')
    do
	cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print if ($.==1 || $F[0]=~m/'$roi'/)' > regulator.csv
	cat gene_wise_quantifications_bac_mRNA_norm.tsv | $PERL_PATH -F"\t" -lane 'print unless ($F[0]=~m/'$roi'/)'      > regulated.csv

	echo 'raw_regulator <- read.table(file="regulator.csv", header=F, skip = 1, sep="\t")' >  r.script
	echo 'raw_targets   <- read.table(file="regulated.csv", header=F, skip = 1, sep="\t")' >> r.script
	echo 'target_names    <- raw_targets[,1]' >> r.script
	echo 'regulator_name  <- raw_regulator[,1]' >> r.script
	echo 'regulator <- data.matrix(raw_regulator[,3:length(raw_regulator[1,])])' >> r.script
	echo 'targets   <- data.matrix(raw_targets[,3:length(raw_targets[1,])])' >> r.script
	echo 'colnames(regulator) <- NULL' >> r.script
	echo 'colnames(targets)   <- NULL' >> r.script
	echo '' >> r.script
	echo 'get_corr <- function(x, r=regulator[1,]){' >> r.script
	echo '  cor.res <- cor.test(x, r, alternative = "two.sided", method ="pearson", conf.level = 0.95)' >> r.script
	echo '  c(cor.res$estimate, cor.res$p.value)' >> r.script
	echo '}' >> r.script
	echo '' >> r.script
	echo 'res_corr_m   <- t(apply(targets, MARGIN=1, FUN=get_corr ))' >> r.script
	echo 'res_corr     <- data.frame(gene = target_names, regulator = regulator_name, R = res_corr_m[,1], Rsquared = (res_corr_m[,1])^2, pval = res_corr_m[,2], padj = p.adjust(res_corr_m[,2], method = "fdr"))' >> r.script   
	echo "res_corr_sig <- subset(res_corr, (Rsquared > $rsquared & padj < $pvalue))" >> r.script
	echo 'if (file.exists("correlation__bac_vs_bac.csv")){' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_bac.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr, file = "correlation__bac_vs_bac.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	echo 'if (file.exists("correlation__bac_vs_bac_sig.csv")){' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_bac_sig.csv", col.names = FALSE, append = TRUE,  quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} else{' >> r.script
	echo 'write.table(res_corr_sig, file = "correlation__bac_vs_bac_sig.csv", col.names = TRUE,  append = FALSE, quote = FALSE, sep = "\t", row.names = FALSE)' >> r.script
	echo '} ' >> r.script
	$RSCRIPT_PATH r.script
	rm -f r.script* regulator.csv regulated.csv
    done

    ## Plot heatmap of correlation data
    # hsa: all genes with R^2 > 0.4225 && p-value adjusted < 0.0001 && positive correlated
    # bac: all secreted genes with R^2 > 0.4225 && p-value adjusted < 0.0001 && positive correlated
    # linker gene: sseC (== YP_005181244)

    cat correlation__bac_vs_bac_sig.csv | $PERL_PATH -F"\t" -lane 'print if $F[1] eq "YP_005181244"' | $PERL_PATH -F"\t" -lane 'print if ($F[-1] < 0.0001 && $F[3] > '$rsquared' && $F[2] > 0)'  | awk '{print $1}' | grep -f data/secreted.grep > bac_genes_to_plot.list
    cat correlation__bac_vs_hsa_sig.csv | $PERL_PATH -F"\t" -lane 'print if $F[1] eq "YP_005181244"' | $PERL_PATH -F"\t" -lane 'print if ($F[-1] < 0.0001 && $F[3] > '$rsquared' && $F[2] > 0)'  | awk '{print $1}' > hsa_genes_to_plot.list

    cat gene_wise_quantifications_bac_mRNA_norm.tsv | head -n 1 | $PERL_PATH -F"\t" -lane '$F[1]="group"; print join("\t", @F)'                     >  data_to_heatmap.csv
    cat gene_wise_quantifications_bac_mRNA_norm.tsv | grep -f bac_genes_to_plot.list | $PERL_PATH -F"\t" -lane '$F[1]="bac"; print join("\t", @F)'  >> data_to_heatmap.csv
    cat gene_wise_quantifications_hsa_mRNA_norm.tsv | grep -f hsa_genes_to_plot.list | $PERL_PATH -F"\t" -lane '$F[1]="hsa"; print join("\t", @F)'  >> data_to_heatmap.csv
    cat gene_wise_quantifications_bac_mRNA_norm.tsv | grep  YP_005181244             | $PERL_PATH -F"\t" -lane '$F[1]="reg"; print join("\t", @F)'  >> data_to_heatmap.csv

    $PERL_PATH bin/translate_plotdata.pl data/bac_genbank2genename.dict data/hsa_ensemble2genename.dict data_to_heatmap.csv > data_to_heatmap_translated.csv

    $RSCRIPT_PATH bin/sseC_correlationheatmap.R
    rm -f sseC_correlationheatmap.Rout
}

main
