#!/bin/sh
#PBS -N noUMI_varscan_pair
#PBS -j oe
#PBS -l ncpus=12
#PBS -l nodes=1
#PBS -l mem=20G

#inputpath=$1
#sample=$2
source  /home/liuxiangqiong/miniconda3/bin/activate base
hg19=/cgdata/Database/GATK/b37/human_g1k_v37_decoy.fasta
target=/cgdata/liuxiangqiong/work62pancancer/pipeline/v0/refdata/NanOnco_Plus_Panel_v2.0_Covered_b37_cg.parY2X.sort
#min_af=0.001
Vardict_Dir=/cgdata/liuxiangqiong/software/VarDict-1.5.7
outputdir=${inputpath}/1SV_vardict_pair
tumorbam=${bam_dir}/${tumor}_clean.bam
normalbam=${bam_dir}/${normal}_clean.bam
#pair
for i in `seq 1 8`
do
${Vardict_Dir}/vardict -G ${hg19} -f ${min_af} -N ${sample} -b "${tumorbam}|${normalbam}" -c 1 -S 2 -E 3 -g 4 ${target}-$i.bed |${Vardict_Dir}/testsomatic.R|${Vardict_Dir}/var2vcf_paired.pl -N "${tumor}|${normal}" -f ${min_af} > ${outputdir}/${sample}_pair.vd_$i.vcf&
done
wait
#bgzip
bcftools view  ${outputdir}/${sample}_pair.vd_1.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_1.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_2.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_2.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_3.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_3.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_4.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_4.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_5.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_5.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_6.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_6.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_7.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_7.vcf.gz
bcftools view  ${outputdir}/${sample}_pair.vd_8.vcf -Oz -o  ${outputdir}/${sample}_pair.vd_8.vcf.gz
#index
bcftools index ${outputdir}/${sample}_pair.vd_1.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_2.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_3.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_4.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_5.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_6.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_7.vcf.gz
bcftools index ${outputdir}/${sample}_pair.vd_8.vcf.gz
#
bcftools concat -a ${outputdir}/${sample}_pair.vd_1.vcf.gz  ${outputdir}/${sample}_pair.vd_2.vcf.gz  ${outputdir}/${sample}_pair.vd_3.vcf.gz ${outputdir}/${sample}_pair.vd_4.vcf.gz  ${outputdir}/${sample}_pair.vd_5.vcf.gz  ${outputdir}/${sample}_pair.vd_6.vcf.gz  ${outputdir}/${sample}_pair.vd_7.vcf.gz  ${outputdir}/${sample}_pair.vd_8.vcf.gz | bcftools sort - -O z -o ${outputdir}/${sample}_pair.vd.vcf
#标准化对齐
vt decompose ${outputdir}/${sample}_pair.vd.vcf |vt normalize - -r ${hg19} -o ${outputdir}/${sample}_pair.vd.normalized.vcf

#less ${outputdir}/${sample}_pair.vd.normalized.vcf | grep -vE "^#"  > ${outputdir}/${sample}_pair.vd.normalized.nohead.vcf

#3.3 将邻近的位点进行合并#3.3 将邻近的位点进行合并
#python3 /cgdata/liuxiangqiong/work62pancancer/Client/v0/script/complex_variant_new.py -i  ${outputdir}/${sample}_pair.vd.normalized.nohead.vcf -b ${bam_dir}/${tumor}_clean.bam -f ${hg19} -mode 2 -n 10 -p 0.5 -o ${outputdir}/${sample}_pair.vd.normalized.nohead.combind.vcf
#grep "^#" ${outputdir}/${sample}_pair.vd.normalized.vcf >${outputdir}/${sample}_pair.vd.normalized.head.combind.vcf
#less ${outputdir}/${sample}_pair.vd.normalized.nohead.combind.vcf >>${outputdir}/${sample}_pair.vd.normalized.head.combind.vcf
conda deactivate
#5.3对经过标准化的进行vep注释
source  /home/liuxiangqiong/miniconda3/bin/activate  ensemble-vep
PERL5LIB=/home/liuxiangqiong/miniconda3/envs/ensemble-vep/share/ensembl-vep-105.0-1/Bio/EnsEMBL/Variation/Utils/DB/HTS/Faidx
export PERL5LIB=$PERL5LIB
inputdir=${outputdir}/${sample}_pair.vd.normalized.vcf
outdir1=${outputdir}/${sample}_pair.vd.normalized.vepanno.vcf
vep --cache --fork 16 \
-i ${inputdir} \
-o ${outdir1} \
--assembly GRCh37 \
--cache_version 104 \
--dir_cache /home/liuxiangqiong/miniconda3/share/ensembl-vep-88.9-0/hg19 \
--refseq \
--canonical \
--offline \
--hgvs \
--fasta /cgdata/Database/GATK/b37/human_g1k_v37_decoy.fasta \
--hgvsg \
--max_af \
--check_existing \
--numbers \
--symbol \
--canonical \
--offline \
-vcf \
--use_given_ref \
--af_gnomad \
--check_existing


#5.4对vep的再进行annovar注释
outdir2=${outputdir}/${sample}_pair.vd.normalized.vepanno.annovar
perl /cgdata/CGTools/soft/src/annovar_20200608/table_annovar.pl ${outdir1} /cgdata/soft/src/annovar.archive/DB/ -buildver hg19 -out ${outdir2} -remove -protocol refGene,popfreq_max_20150413,gnomad_genome,avsnp150,snp138NonFlagged,dbnsfp33a,cosmic94,clinvar_20210501,intervar_20170202,icgc28 -operation g,f,f,f,f,f,f,f,f,f -nastring . -vcfinput --thread 40

#remove
rm ${outputdir}/${sample}_pair.vd_?.vcf*