r1_trim=$1
r2_trim=$2
sample=$3
outdir=$4

curPath=$(dirname $(readlink -f "$0"))
subPath=$(dirname $(readlink -f $curPath))

no_IPC_target=${subPath}/DB/no_IPC_target.bed
IPC=${subPath}/DB/IPC.bed

ref=/cgdata/IVDRD/18gene_DB/qaseq_fusion/hg38.fa

mkdir -p ${outdir}/2.mapping;

bwa mem -t 16 ${ref} ${r1_trim} ${r2_trim} -R "@RG\tID:1\tLB:lib1\tPL:ILLUMINA\tSM:${sample}\tPU:unit1" > ${outdir}/2.mapping/${sample}.sam;

samtools view ${outdir}/2.mapping/${sample}.sam --threads 8 -q 30 -bS > ${outdir}/2.mapping/${sample}.pe.bam;

samtools view ${outdir}/2.mapping/${sample}.pe.bam > ${outdir}/2.mapping/${sample}.q30.sam;

sambamba sort ${outdir}/2.mapping/${sample}.pe.bam --tmpdir=${outdir}/tmp -o ${outdir}/2.mapping/${sample}.sort.bam -t 8;

samtools view -L $no_IPC_target ${outdir}/2.mapping/${sample}.sort.bam > ${outdir}/2.mapping/${sample}.sort.target.sam;

samtools view  -L $IPC ${outdir}/2.mapping/${sample}.sort.bam > ${outdir}/2.mapping/${sample}.sort.IPC.sam;

awk '{print $1}' ${outdir}/2.mapping/${sample}.sort.target.sam | sort | uniq > ${outdir}/2.mapping/${sample}.target.reads.txt;

awk '{print $1}' ${outdir}/2.mapping/${sample}.sort.IPC.sam | sort | uniq > ${outdir}/2.mapping/${sample}.IPC.reads.txt;

awk 'NR==FNR{a[$1]="OK";next}{print $0"\t"a[$1]}' ${outdir}/2.mapping/${sample}.IPC.reads.txt  ${outdir}/2.mapping/${sample}.q30.sam | awk '{if($NF!="OK"){print $0}}' > ${outdir}/2.mapping/${sample}.q30.tmp.sam;

awk 'NR==FNR{a[$1]="OK";next}{print $0"\t"a[$1]}' ${outdir}/2.mapping/${sample}.target.reads.txt ${outdir}/2.mapping/${sample}.q30.tmp.sam | awk '{if($NF!="OK"){print $0}}' > ${outdir}/2.mapping/${sample}.sort.nonspecific.sam;

awk '{print $1}' ${outdir}/2.mapping/${sample}.sort.nonspecific.sam | sort | uniq > ${outdir}/2.mapping/${sample}.nonspecific.reads.txt