#!/bin/bash

#purpose：mapping reads to reference genome to get the bam file
#set env
workshop="/home/zhanghl/data_3/workshop/practice1_ovary"
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

#preprocessing steps
preprocessing_function(){
    local sample_name=$1
    local input=$2
    local output=$3
#mapping
    $bwa mem \
    -t 10 \
    -M  \
    -R "@RG\tID:HKV2KALXX\tSM:${sample_name}\tPL:illumina\tLB:${sample_name}" \
    $human_fasta \
    ${input}/${sample_name}_1.clean.fq.gz ${input}/${sample_name}_2.clean.fq.gz > ${output}/${sample_name}.sam
    echo "$sample_name mapping step is done" >> ${output}/run.log    

#转换成BAM文件，sort, mark duplicates
#sortbam
    java -jar $picard SortSam \
     INPUT=${output}/${sample_name}.sam \
     OUTPUT=${output}/${sample_name}.sorted.bam \
     SORT_ORDER=coordinate
     echo "$sample_name sortsam step is done" >> ${output}/run.log
#MarkDuplicates
    java -jar $picard MarkDuplicates \
      INPUT=${output}/${sample_name}.sorted.bam \
      OUTPUT=${output}/${sample_name}.sorted.dedup.bam \
      METRICS_FILE=${output}/${sample_name}.metrics.txt
      echo "$sample_name MarkDuplicates step is done" >> ${output}/run.log
#Index bam
    java -jar $picard BuildBamIndex \
      INPUT=${output}/${sample_name}.sorted.dedup.bam
    echo "$sample_name index bam step is done" >> ${output}/run.log     
}

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