#!/usr/bin/sh

genome=/home/chenfan/mnt/sdc/RNAseq/hs37d5.fasta

refgene=/home/chenfan/mnt/sdc/RNAseq/hg19_Refgene.gtf

picard=/home/chenfan/mnt/sdc/RNAseq/picard.jar

gatk=/home/chenfan/mnt/sdc/RNAseq/GenomeAnalysisTK.jar

knownindel=/home/chenfan/VQSR/data/Mills_and_1000G_gold_standard.indels.b37.vcf

knownsnp=/home/chenfan/VQSR/data/1000G_phase1.snps.high_confidence.b37.vcf

result_path=/home/chenfan/mnt/sdc/RNAseq/result

#------------------------------STAR 2-pass alignment steps-----------------------------------------

echo "------------------------------STAR 2-pass alignment steps-----------------------------------------"

in_dir=${result_path}/1_bamtofastq_output

out_dir=${result_path}/2_star_output

cd ${out_dir}/1index

STAR --runMode genomeGenerate --genomeDir ${out_dir}/1index --genomeFastaFiles $genome --sjdbGTFfile $refgene --runThreadN 30

cd ${out_dir}/1pass

STAR --genomeDir ${out_dir}/1index --readFilesIn ${in_dir}/1.fastq  ${in_dir}/2.fastq --runThreadN 30

cd ${out_dir}/2index

STAR --runMode genomeGenerate --genomeDir ${out_dir}/2index --genomeFastaFiles $genome --sjdbFileChrStartEnd ${out_dir}/1pass/SJ.out.tab --sjdbOverhang 154 --runThreadN 30

cd ${out_dir}/2pass

STAR --genomeDir ${out_dir}/2index --readFilesIn ${in_dir}/1.fastq ${in_dir}/2.fastq  --runThreadN 30

#------------------------------Picard add read groups,sort,mark duplicates,and create index-----------------------------------------

echo "------------------------------Picard add read groups,sort,mark duplicates,and create index-----------------------------------------"

in_dir=${out_dir}

out_dir=${result_path}/3_picard_output

cd ${out_dir}

java -jar $picard AddOrReplaceReadGroups I=${in_dir}/2pass/Aligned.out.sam O=${out_dir}/sort_addRG_align.bam SO=coordinate RGID=hap1 RGLB=rnaseq RGPL=illumina RGPU=runname RGSM=C20

java -jar $picard MarkDuplicates I=${out_dir}/sort_addRG_align.bam O=${out_dir}/markdup_sort_addRG_align.bam  CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT M=${out_dir}/output.metrics

#-------------------------------SplitNCigarReads---------------------------------------

echo "------------------------------SplitNCigarReads-----------------------------------------"

in_dir=${out_dir}

out_dir=${result_path}/4_splitncigarReads_output

cd ${out_dir}

java -jar $gatk -T SplitNCigarReads -R $genome -I ${in_dir}/markdup_sort_addRG_align.bam -o ${out_dir}/split_markdup_sort_addRG_align.bam -rf ReassignOneMappingQuality -RMQF 255 -RMQT 60 -U ALLOW_N_CIGAR_READS

#-------------------------------Indel Realignment ---------------------------------------

echo "------------------------------Indel Realignment -----------------------------------------"

in_dir=${out_dir}

out_dir=${result_path}/5_indelrealign_output

cd ${out_dir}

java -jar $gatk \

-T RealignerTargetCreator \

-nt 32 \

-R $genome \

-I ${in_dir}/split_markdup_sort_addRG_align.bam \

-o ${out_dir}/realign_interval.list \

-known $knownindel

java -jar $gatk \

-T IndelRealigner \

-R $genome \

-I ${in_dir}/split_markdup_sort_addRG_align.bam \

-known $knownindel \

-o ${out_dir}/realign_split_markdup_sort_addRG_align.bam \

-targetIntervals ${out_dir}/realign_interval.list

#-------------------------------Base Recalibration---------------------------------------

echo "------------------------------Base Recalibration-----------------------------------------"

in_dir=${out_dir}

out_dir=${result_path}/6_bqsr_output

cd ${out_dir}

java -jar $gatk \

-T BaseRecalibrator \

-nct 32 \

-R $genome \

-I ${in_dir}/realign_split_markdup_sort_addRG_align.bam \

-knownSites $knownsnp \

-knownSites $knownindel \

-o ${out_dir}/recal_data.table

java -jar $gatk \

-T PrintReads \

-nct 32 \

-R $genome \

-I ${in_dir}/realign_split_markdup_sort_addRG_align.bam \

-BQSR ${out_dir}/recal_data.table \

-o ${out_dir}/bqsr_realign_split_markdup_sort_addRG_align.bam

#-------------------------------Variant calling---------------------------------------

echo "------------------------------Variant calling-----------------------------------------"

in_dir=${out_dir}

out_dir=${result_path}/7_callvariants_output

cd ${out_dir}

java -jar $gatk \

-T HaplotypeCaller \

-nct 32 \

-R $genome \

-ploidy 1 \

-I ${in_dir}/bqsr_realign_split_markdup_sort_addRG_align.bam \

-dontUseSoftClippedBases \

-stand_call_conf 20.0 \

-o ${out_dir}/Hap1_rnaseq.vcf

#-------------------------------Variant filtering---------------------------------------

echo "------------------------------Variant filtering-----------------------------------------"

in_dir=${out_dir}

out_dir=${result_path}/8_filtervariants_output

cd ${out_dir}

java -jar $gatk \

-T VariantFiltration \

-R $genome \

-V ${in_dir}/Hap1_rnaseq.vcf \

-window 35 \

-cluster 3 \

--filterName FS -filter "FS>30.0" \

--filterName QD -filter "QD<2.0" \

-o ${out_dir}/filtered_Hap1_rnaseq.vcf

java -jar $gatk \

-T SelectVariants \

-ef \

-R $genome \

-V ${out_dir}/filtered_Hap1_rnaseq.vcf \

-o ${out_dir}/ef_Hap1_C20_rnaseq.vcf