#!/bin/bash

#purpose:recalibration of reads quality
##Directories
workshop="/home/zhanghl/data_3/workshop/practice1_ovary"

##resources
human_b37_bundle="/home/zhanghl/supporting_files/Homo_sapiens_glk_v37/human_b37_bundle"

human_fasta=${human_b37_bundle}/human_g1k_v37_decoy.fasta
hapmap=${human_b37_bundle}/hapmap_3.3.b37.vcf
omni_1000G=${human_b37_bundle}/1000G_omni2.5.b37.vcf
phase1_1000G=${human_b37_bundle}/1000G_phase1.snps.high_confidence.b37.vcf
snp=${human_b37_bundle}/dbsnp_138.b37.vcf
mills=${human_b37_bundle}/Mills_and_1000G_gold_standard.indels.b37.vcf
gold_indels=${human_b37_bundle}/1000G_phase1.indels.b37.vcf

##Softwares
softwares="/home/zhanghl/supporting_softwares"

export  samtools="${softwares}/samtools-0.1.19/samtools"
export  STAR="${softwares}/STAR-master/bin/Linux_x86_64_static/STAR"
export  tophat2="${softwares}/tophat-2.0.12/tophat2"
export  Trimmomatic="${softwares}/Trimmomatic-0.33/trimmomatic-0.33.jar"
export  VarScan="${softwares}/Varscan.v2.4.1/VarScan.v2.4.1.jar"
export  GenomeAnalysisTK="${softwares}/GenomeAnalysisTK/GenomeAnalysisTK.jar"
export  HTseq="${softwares}/HTSeq-0.6.1p1/HTSeq/scripts/count.py"
export  bowtie2="${softwares}/bowtie2-2.2.3/bowtie2"
export  qc="${softwares}/qc_trim/QC.sh"
export  fastq_dump="${software}/sratoolkit.2.3.5-2-centos_linux64/bin/fastq-dump"
export  featurecount="${softwares}/featureCounts"
export  bwa="${softwares}/bwa-0.7.10/bwa"

##JAVA Setting
export JAVA_HOME=${softwares}/jdk1.8.0_20
export PATH=$JAVA_HOME/bin:$PATH
export CLASSPATH=.:$JAVA_HOME/lib/dt.jar:$JAVA_HOME/lib/tools.jar
java_single_task_Mem=5g
java_combine_task_Mem=47g
export picard=${softwares}/picard/picard.jar

###################################################################    
#$GenomeAnalysisTK recalibration
recalibration_function(){
    local sample_name=$1
    local input=$2
    local output=$3
    
#Analyze patterns of covariation in the sequence dataset
    java -jar $GenomeAnalysisTK \
    -T BaseRecalibrator \
    -R $human_fasta \
    -I $input/${sample_name}.sorted.dedup.bam \
    -knownSites $snp \
    -knownSites $gold_indels \
    -o ${output}/${sample_name}.recal_data.table
    
#Do a second pass to analyze covariation remaining after recalibration
    java -jar $GenomeAnalysisTK \
    -T BaseRecalibrator \
    -R $human_fasta \
    -I ${input}/${sample_name}.sorted.dedup.bam \
    -knownSites $snp \
    -knownSites $gold_indels \
    -BQSR ${output}/${sample_name}.recal_data.table \
    -o ${output}/${sample_name}.post_recal_data.table

#Generate before/after plots
    java -jar $GenomeAnalysisTK \
    -T AnalyzeCovariates \
    -R $human_fasta \
    -before ${output}/${sample_name}.recal_data.table \
    -after ${output}/${sample_name}.post_recal_data.table \
    -plots ${output}/${sample_name}.recalibration_plots.pdf 
    
#Apply the recalibration to your sequence data
    java -jar $GenomeAnalysisTK \
    -T PrintReads \
    -R $human_fasta \
    -I ${input}/${sample_name}.sorted.dedup.bam \
    -BQSR ${output}/${sample_name}.recal_data.table \
    -o    ${output}/${sample_name}.sorted.dedup.recal.bam 
}


for sample in `cat /your_path/sample_list`;do  
 recalibration_function \
 $sample \
/your_path1 \
/your_path2
done