#! /bin/bash
# mapping_bwa-mem_CF
# last edit: 2021-06
# authors: Christelle Fraïsse

#-# usage: 
#	while read line; do
#	qsub mapping_bwa-mem_CF.sh ${line}
#	done </share/nas-isem_e/genomique_integrative/Cione_WGS/NamesSampleTag.txt
#-# prog required: bwa v0.7.5a ; bamutil v1.0.14 ; picard v1.119 ; samtools v1.9 ; GATK v3.4-0-g7e26428
#-# files required: NamesSampleTag.txt [produced with merge_fasta.sh]
#-# description: Mapping fastq files to the assembly.


## define global variables
toBWA=/home/fchristelle/bin/bwa-0.7.5a
toBAM=/home/fchristelle/bin/bamUtil
toJAVA=/home/fchristelle/bin/java
toSAM=/home/fchristelle/bin/samtools-1.9
toGATK=/share/nas-isem_e/genomique_integrative/Cione_WGS
refgen=/share/nas-isem_e/genomique_integrative/Cione_WGS/RefCiona/GCA_009617815.1_Ciona_intestinalis_HT_Hoya_T-line_assembly_2019_genomic.fna

## your picard sample
f=$1
pfx=`echo $f | cut -d'_' -f 1`
lib=`echo $f | cut -d'_' -f 3`
infolder=/share/nas-isem_e/genomique_integrative/Ciona_raw_data/WGS_Trios/Sample_${pfx}
 
# create new folder
outfolder=/share/nas-isem_e/genomique_integrative/Cione_WGS/assembly_2019/ciona_bwa-mapping_${pfx}
mkdir -p ${outfolder}
 
# common
fq1=${infolder}/${f}_R1.fastq.gz
fq2=${infolder}/${f}_R2.fastq.gz

# marking secondary hits with -M to comply with Picard
# using 'nthr' processors in parallel
# mem is more demanding and needs more than 3Gb per thread
nthr=10
 
rgstring="@RG\tID:ciona-2019_WGS_$pfx\tCN:IBL\tDT:2020-01-06\tLB:lib-$lib\tPL:ILLUMINA-2019\tSM:$f"
 
echo "# mapping paired reads with **bwa mem**"


## default numeric settings are left unchanged as:
#  -k INT     minimum seed length [19]
#  -w INT     band width for banded alignment [100]
#  -d INT     off-diagonal X-dropoff [100]
#  -r FLOAT   look for internal seeds inside a seed longer than {-k} * FLOAT [1.5]
#  -c INT     skip seeds with more than INT occurrences [10000]
#  -A INT     score for a sequence match [1]
#  -B INT     penalty for a mismatch [4]
#  -O INT     gap open penalty [6]
#  -E INT     gap extension penalty; a gap of size k cost {-O} + {-E} * k [1]
#  -L INT     penalty for clipping [5]
#  -U INT     penalty for an unpaired read pair [17]
#  -T INT     minimum score to output [30]
# '-H':	Use hard clipping ’H’ in the SAM output. This option may dramatically reduce the redundancy of output when mapping long contig or BAC sequences. 
# '-M': mark shorter split hits as secondary (for Picard/GATK compatibility)
 
bwape=${pfx}_bwa_mem
 
# store the full command line to include it in the next part
cmd="${toBWA}/bwa mem -R \"${rgstring}\" \
	-M \
	-t ${nthr} \
	${refgen} \
	${fq1} ${fq2} > ${outfolder}/${bwape}.sam"
 
# display the full command with expanded variables
echo "# the command will be:"
echo ${cmd}
 
# execute the command
eval ${cmd} 2>${outfolder}/${pfx}_bwa_mem.err


########################### post process results ##############

#### add @PG group to results with **bamUtil 'bam polishBam' **
# $pgstring is called quoted to be replaced by its building variables
 
bwaver=$(expr "$(bwa 2>&1)" : '.*Version:\ \(.*\)')
pgstring=@PG$'\t'ID:BWA$'\t'PN:bwa$'\t'VN:${bwaver}$
 
echo "# now adding a PG field to the data"
 
# clean up initial file after success to save space
${toBAM}/bam polishBam --PG "${pgstring}" \
	-i ${outfolder}/${bwape}.sam \
	-o ${outfolder}/${bwape}_pg.sam \
	-l ${outfolder}/polishBam_${bwape}.log && rm -f ${outfolder}/${bwape}.sam \
 
##### convert to sorted bam, index & cleanup
java -Xmx4g -jar ${toJAVA}/SortSam.jar \
	I=${outfolder}/${bwape}_pg.sam \
	O=${outfolder}/${bwape}.bam \
	SO=coordinate \
	CREATE_INDEX=true \
	VALIDATION_STRINGENCY=LENIENT \
	
##### mark duplicates with Picard & index
java -Xmx4g -jar ${toJAVA}/MarkDuplicates.jar \
	I=${outfolder}/${bwape}.bam \
	O=${outfolder}/${bwape}_mdup.bam \
	ASSUME_SORTED=TRUE \
	REMOVE_DUPLICATES=FALSE \
	CREATE_INDEX=TRUE \
	METRICS_FILE=${outfolder}/${bwape}_duplicate_metrics.txt \
	VALIDATION_STRINGENCY=LENIENT \
	
##### get basic stats from the resulting bam file with Samtools
${toSAM}/samtools flagstat ${outfolder}/${bwape}_mdup.bam \
	>${outfolder}/${bwape}_mdup_flagstat.txt
	
##### get basic stats from the resulting bam file with GTAK
## Default numeric settings are left unchanged as:
# --maxBaseQuality 	127
# --maxMappingQuality	2147483647
# --minBaseQuality 	-1
# --minMappingQuality 	-1
# --partitionType 	NA [default coverage partition by sample, but it can be any combination of sample, readgroup and/or library]
# --countType 		NA [How should overlapping reads from the same fragment be handled?]
java -Xmx4g -jar ${toGATK}/GenomeAnalysisTK.jar \
	-R ${refgen} \
	-T DepthOfCoverage \
	-o ${outfolder}/${bwape}_mdup_depthofcoverage \
	-I ${outfolder}/${bwape}_mdup.bam \
	--includeDeletions \
	--omitDepthOutputAtEachBase \
	--omitIntervalStatistics \
	--omitLocusTable \

##### Cleanup
rm -f ${outfolder}/${bwape}_pg.sam
rm -f ${outfolder}/${bwape}.bam
rm -f ${outfolder}/${bwape}.bai
    

echo "# finished mapping ${pfx} reads with BWA mem"
