#!/usr/bin/bash

################################################################
##       step1: recording all version informations            ##
################################################################

################# pipeline_version ########################
# pipeline: ${pipeline_version}

#################### soft_version ########################
# fastp: ${fastp_version}
# bwa: ${bwa_version}
# samtools: ${samtools_version}
# sentieon: ${sentieon_version}
# bcftools: ${bcftools_version}
# annovar: ${annovar_version}
# vep: ${vep_version}
# msisensor2: ${msisensor2_version}
# genefuse: ${genefuse_version}
# cnvkit: ${cnvkit_version}
# bamdst: ${bamdst_version}
# seqtk: ${seqtk_version}
# htslib: ${htslib_version}
# python2: ${python2_version}
# python3: ${python3_version}
# delly: ${delly_version}

#################### bed_version ########################
# bed: ${bed_version}

################# database_version ########################
# reference: ${reference_version}
# refgene: ${refgene_version}
# cytoBand: ${cytoBand_version}
# 1000g: ${g1000_version}
# esp: ${esp_version}
# exac: ${exac_version}
# cosmic: ${cosmic_version}
# clinvar: ${clinvar_version}
# intervar: ${intervar_version}
# dbnsfp: ${dbnsfp_version}
# avsnp: ${avsnp_version}
# private_database: ${private_database}
# pon_database: ${pon_database}



################################################################
##               step2: pipeline for analysis                 ##
################################################################

# Update with the location of the Sentieon software package and license file
export SENTIEON_LICENSE=${SENTIEON_LICENSE}


threads=${THREAD}

cd $outDir
echo "$sampleID -- beginning at: $(date)"

$python3 ${checkUMI} -f $fq1 -s ${sampleType}
################# fastp ########################
# fastp is used for remove lower quality bases
echo "fastp"
echo "$(date)"
$fastp -w ${threads} -i $fq1 -o ${sampleID}.R1_001.trim.fq.gz -I $fq2 -O ${sampleID}.R2_001.trim.fq.gz -l 35 -c --detect_adapter_for_pe -j ${sampleID}.jion -h ${sampleID}.html


################# UMI + alignment + dedup ########################
# extract umi info and alignment to reference
if [[ ${bedType} =~ ^(${umi_panel})$ ]];
then
# remove UMI from following panels
echo $'\nsentieon umi extract'
echo "$(date)"
$sentieon umi extract 21S+T,12M11S+T ${sampleID}.R1_001.trim.fq.gz ${sampleID}.R2_001.trim.fq.gz | \
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000  ${hg19fa} - || echo -n 'error') | \
$sentieon util sort -t ${threads} --sam2bam -i - -o ${sampleID}.sorted.bam

echo $'\nsentieon umi consensus'
echo "$(date)"
$sentieon umi consensus -i ${sampleID}.sorted.bam --input_format BAM -o ${sampleID}.umi_consensus.fq.gz

echo $'\nbwa' 
echo "$(date)"
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} ${sampleID}.umi_consensus.fq.gz || echo -n 'error' ) | \
$sentieon util sort --umi_post_process --sam2bam -i - -o ${sampleID}.sorted.dedup.bam

elif [[ ${bedType} =~ (${probe_panel}) && ${sampleType} =~ (ctdna) ]];
then
# remove UMI from following panels
echo $'\nsentieon umi'
echo "$(date)"
$sentieon umi extract -d 6M1S+T,6M1S+T ${sampleID}.R1_001.trim.fq.gz ${sampleID}.R2_001.trim.fq.gz | \
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000  ${hg19fa} - || echo -n 'error') | \
$sentieon util sort -t ${threads} --sam2bam -i - -o ${sampleID}.sorted.bam

echo $'\nsentieon umi consensus'
echo "$(date)"
$sentieon umi consensus -i ${sampleID}.sorted.bam --input_format BAM -o ${sampleID}.umi_consensus.fq.gz

echo $'\nbwa' 
echo "$(date)"
( $sentieon bwa mem -t ${threads} -p -C -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} ${sampleID}.umi_consensus.fq.gz || echo -n 'error' ) | \
$sentieon util sort --umi_post_process --sam2bam -i - -o ${sampleID}.sorted.dedup.bam

else
# following panels does not contain UMI
echo $'\nbwa' 
echo "$(date)"
( $sentieon bwa mem -t ${threads} -M -R "@RG\\tID:${sampleID}\\tSM:${sampleID}\\tPL:ILLUMINA\\tLB:lib" -K 10000000 ${hg19fa} ${sampleID}.R1_001.trim.fq.gz ${sampleID}.R2_001.trim.fq.gz || echo -n 'error' ) | \
$sentieon util sort -o ${sampleID}.sorted.bam -t ${threads} --sam2bam -i -

if [[ ${bedType} =~ ^(${probe_panel})$ ]]; then
echo $'\nRemove Duplicate'
echo "$(date)"
$sentieon driver -t ${threads} -i ${sampleID}.sorted.bam --algo LocusCollector --fun score_info ${sampleID}.score.txt
$sentieon driver -t ${threads} -i ${sampleID}.sorted.bam --algo Dedup --score_info ${sampleID}.score.txt --metrics ${sampleID}.dedup_metrics.txt ${sampleID}.sorted.dedup.bam
rm ${sampleID}.score.txt* ${sampleID}.dedup_metrics.txt
fi
fi


################## Variant Calling ######################
# TNscope is used for somatic SNVs & small INDELs calling
if [[ ${bedType} =~ ^(${pcr_panel})$ ]]; then

echo $'\nHaplotyper'
echo "$(date)"
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.bam \
--interval ${bedPath}/${bedType}.bed \
--algo Haplotyper --emit_conf=10 --call_conf=10 \
${sampleID}.germline.raw.vcf

echo $'\nTNscope'
echo "$(date)"

if [[ ${bedType} =~ (AT|CA) ]];
then
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.bam \
--interval ${bedPath}/${bedType}.bed --interval_padding 10 \
--algo TNscope \
--tumor_sample ${sampleID} \
--dbsnp ${dbsnp138vcf} \
--disable_detector sv \
--max_fisher_pv_active 0.05 \
--min_tumor_allele_frac 0.009 \
--sv_mask_ext 10 \
--prune_factor 0 \
--assemble_mode 2 \
--resample_depth 100000 \
--pon ${pon} \
${sampleID}.output_tnscope.pre_filter.vcf
else
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.bam \
--interval ${bedPath}/${bedType}.bed --interval_padding 10 \
--algo TNscope \
--tumor_sample ${sampleID} \
--dbsnp ${dbsnp138vcf} \
--disable_detector sv \
--max_fisher_pv_active 0.05 \
--min_tumor_allele_frac 0.009 \
--sv_mask_ext 10 \
--prune_factor 0 \
--assemble_mode 2 \
--resample_depth 100000 \
--trim_primer ${bedPath}/${bedType}.amplicon.bed \
--pon ${pon} \
${sampleID}.output_tnscope.pre_filter.vcf
fi

echo $'\nbcftools filter'
echo "$(date)"
$bcftools annotate -x "FILTER/triallelic_site" ${sampleID}.output_tnscope.pre_filter.vcf | \
$bcftools filter -m + -s "insignificant" -e "(PV>0.25 && PV2>0.25)" | \
$bcftools filter -m + -s "insignificant" -e "(INFO/STR == 1 && PV>0.05)" | \
$bcftools filter -m + -s "low_qual" -e "QUAL < 20" | \
$bcftools filter -m + -s "short_tandem_repeat" -e "RPA[0]>=10" | \
$bcftools filter -m + -s "noisy_region" -e "ECNT>5" | \
$sentieon util vcfconvert - ${sampleID}.somatic.raw.vcf
rm ${sampleID}.output_tnscope.pre_filter.vcf*

elif [[ ${bedType} =~ ^(${probe_panel})$ && ${sampleType} != "ctdna" ]]; then

echo $'\nHaplotyper'
echo "$(date)"
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam \
--interval ${bedPath}/${bedType}.bed \
--algo Haplotyper --emit_conf=10 --call_conf=10 \
${sampleID}.germline.raw.vcf

echo $'\nTNscope'
echo "$(date)"
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam \
--interval ${bedPath}/${bedType}.bed --interval_padding 10 \
--algo TNscope \
--tumor_sample ${sampleID} \
--dbsnp ${dbsnp138vcf} \
--disable_detector sv \
--max_fisher_pv_active 0.05 \
--min_tumor_allele_frac 0.005 \
--filter_t_alt_frac 0.005 \
--resample_depth 100000 \
--assemble_mode 4 \
--pon ${pon} \
${sampleID}.output_tnscope.pre_filter.vcf

echo $'\nbcftools filter'
echo "$(date)"
echo $'\nbcftools filter'
echo "$(date)"
$bcftools annotate -x "FILTER/triallelic_site" ${sampleID}.output_tnscope.pre_filter.vcf | \
$bcftools filter -m + -s "insignificant" -e "(PV>0.25 && PV2>0.25)" | \
$bcftools filter -m + -s "insignificant" -e "(INFO/STR == 1 && PV>0.05)" | \
$bcftools filter -m + -s "orientation_bias" -e "FMT/FOXOG[0] == 1" | \
$bcftools filter -m + -s "strand_bias" -e "SOR > 3" | \
$bcftools filter -m + -s "low_qual" -e "QUAL < 20" | \
$bcftools filter -m + -s "short_tandem_repeat" -e "RPA[0]>=10" | \
$bcftools filter -m + -s "noisy_region" -e "ECNT>5" | \
$bcftools filter -m + -s "read_pos_bias" -e "FMT/ReadPosRankSumPS[0] < -8" | \
$sentieon util vcfconvert - ${sampleID}.somatic.raw.vcf
rm ${sampleID}.output_tnscope.pre_filter.vcf*

else
echo $'\nHaplotyper'
echo "$(date)"
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam \
--interval ${bedPath}/${bedType}.bed \
--algo Haplotyper --emit_conf=10 --call_conf=10 \
${sampleID}.germline.raw.vcf

echo $'\nTNscope'
echo "$(date)"
$sentieon driver -r ${hg19fa} -t ${threads} -i ${sampleID}.sorted.dedup.bam --interval ${bedPath}/${bedType}.bed \
--algo TNscope \
--tumor_sample ${sampleID} \
--dbsnp ${dbsnp138vcf} \
--pcr_indel_model NONE \
--disable_detector sv \
--min_tumor_allele_frac 0.0003 \
--filter_t_alt_frac 0.0003 \
--max_error_per_read 3 \
--min_init_tumor_lod 3.0 \
--min_tumor_lod 3.0 \
--min_base_qual 40 \
--min_base_qual_asm 40 \
--resample_depth 100000 \
--assemble_mode 4 \
--pon ${pon} \
${sampleID}.output_tnscope.pre_filter.vcf

echo $'\nbcftools filter'
echo "$(date)"
$bcftools annotate -x "FILTER/triallelic_site" ${sampleID}.output_tnscope.pre_filter.vcf | \
$bcftools filter -m + -s "low_qual" -e "QUAL < 10" | \
$bcftools filter -m + -s "short_tandem_repeat" -e "RPA[0]>=10" | \
$bcftools filter -m + -s "read_pos_bias" -e "FMT/ReadPosRankSumPS[0] < -5" | \
$bcftools filter -m + -s "base_qual_bias" -e "FMT/BaseQRankSumPS[0] < -5" | \
$sentieon util vcfconvert - ${sampleID}.somatic.raw.vcf
rm ${sampleID}.output_tnscope.pre_filter.vcf*
fi


echo $'\nfilter repeat regions'
echo "$(date)"
grep "#" ${sampleID}.somatic.raw.vcf > ${sampleID}.somatic.filter.vcf
grep -v "#" ${sampleID}.somatic.raw.vcf | awk '{len=split($5,a,",");if(len<3 || (len>=3 && $7 !~ "panel_of_normals")) print $0}' >> ${sampleID}.somatic.filter.vcf


################# 1p19q ########################

if [[ ${bedType} =~ (AT) ]];
then
$samtools mpileup -l ${bedPath}/AT.force.bed -f ${hg19fa} -d 30000 -Q 0 -O -s ${sampleID}.sorted.bam > ${sampleID}.mpileup.snp.txt
$python3 ${dig_pileup} ${sampleID}.mpileup.snp.txt -o ${sampleID}.mpileup.snp.vcf
$python3 ${p1q19_paint} ${sampleID}.mpileup.snp.vcf
rm ${sampleID}.mpileup.snp.txt
fi

################# bcftools ########################
# bcftools is used for retrieving mutations of interest, or mutations with higher confidence

echo $'\nbcftools'
echo "$(date)"
#germline
$bcftools norm -m -both ${sampleID}.germline.raw.vcf -o ${sampleID}.germline.raw.split.vcf
$bcftools norm -f ${hg19fa} ${sampleID}.germline.raw.split.vcf -o ${sampleID}.germline.split.norm.vcf
rm ${sampleID}.germline.raw.vcf* ${sampleID}.germline.raw.split.vcf

# #somatic
$bcftools norm -m -both ${sampleID}.somatic.filter.vcf -o ${sampleID}.somatic.raw.split.vcf
$bcftools norm -f ${hg19fa} ${sampleID}.somatic.raw.split.vcf -o ${sampleID}.somatic.split.norm.vcf
rm ${sampleID}.somatic.filter.vcf ${sampleID}.somatic.raw.split.vcf 


###################### ANNOVAR ########################
# used for annotation

###################### germline ########################
echo $'\ngermline'
echo "$(date)"
mkdir anno
${annovar} ${sampleID}.germline.split.norm.vcf ${humandb} \
--outfile anno/${sampleID}.germline.anno --buildver hg19 \
--protocol refGene,cytoBand,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150 \
--operation g,r,f,f,f,f,f,f,f,f  --vcfinput --remove
rm anno/${sampleID}.germline.anno.avinput

if ([[ "$bedType" =~ (${germline_panel1}) ]]) || ([[ "$bedType" =~ (${germline_panel2}) ]]); then
$python3 ${germline_filter} anno/${sampleID}.germline.anno.hg19_multianno.vcf -o anno/${sampleID}.germline.inital.filter.vcf > anno/${sampleID}.germline.filter_reason.txt

# vep注释（主要为HGVS）
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
--use_given_ref --refseq --assembly GRCh37 --format vcf \
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
--hgvs --hgvsg --transcript_version \
-i anno/${sampleID}.germline.inital.filter.vcf -o anno/${sampleID}.germline.vep_anno.vcf
rm anno/${sampleID}.germline.inital.filter.vcf

# HGVS命名注释
$python3 ${anno_hgvs} anno/${sampleID}.germline.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.germline.hgvs_anno.vcf
rm anno/${sampleID}.germline.vep_anno.vcf*

# VEP修复后假阳性删除
$python3 ${snv_second_filter} anno/${sampleID}.germline.hgvs_anno.vcf -k germline -o anno/${sampleID}.germline.final.filter.vcf >>anno/${sampleID}.germline.filter_reason.txt
fi


###################### somatic ########################
echo $'\nsomatic'
echo "$(date)"

# 过滤 nopass 突变
grep "#" ${sampleID}.somatic.split.norm.vcf > anno/${sampleID}.somatic.split.filter.vcf
grep -v "#" ${sampleID}.somatic.split.norm.vcf | awk '{len=split($7,a,";");split($NF,b,":");if(len<2 || (len>=2 && b[3] >=0.01)) print $0}' >> anno/${sampleID}.somatic.split.filter.vcf

# annovar注释
${annovar} anno/${sampleID}.somatic.split.filter.vcf ${humandb} \
--outfile anno/${sampleID}.somatic.annovar --buildver hg19  \
--protocol refGene,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150,dbscsnv11  \
--operation g,f,f,f,f,f,f,f,f,f  --vcfinput --remove

# vep注释(主要为HGVS)
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
--use_given_ref --refseq --assembly GRCh37 --format vcf \
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
--hgvs --hgvsg --transcript_version \
-i anno/${sampleID}.somatic.annovar.hg19_multianno.vcf -o anno/${sampleID}.somatic.vep_anno.vcf

# HGVS命名注释
$python3 ${anno_hgvs} anno/${sampleID}.somatic.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.somatic.hgvs_anno.vcf

# 初步过滤(假阳性等)
$python3 ${snv_init_filter} anno/${sampleID}.somatic.hgvs_anno.vcf \
-o anno/${sampleID}.somatic.inital.filter.vcf -t ${sampleType} -m ${mode} > anno/${sampleID}.somatic.filter_reason.txt

# 合并同一单倍型上的复杂突变
$python3 ${merge_mnv} anno/${sampleID}.somatic.inital.filter.vcf $hg19fa -o anno/${sampleID}.somatic.merge_mnv.vcf

# 将合并后的mnp位置左移校准
$bcftools norm -m -both -f ${hg19fa} anno/${sampleID}.somatic.merge_mnv.vcf -o anno/${sampleID}.somatic.merged_norm.vcf

# 拆分MNV突变
$python3 ${find_mnp} anno/${sampleID}.somatic.merged_norm.vcf

if [ -f "anno/${sampleID}.somatic.mnp.vcf" ];
then
# annovar注释MNV
${annovar} anno/${sampleID}.somatic.mnp.vcf ${humandb} \
--outfile anno/${sampleID}.somatic.mnp.annovar --buildver hg19  \
--protocol refGene,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150,dbscsnv11  \
--operation g,f,f,f,f,f,f,f,f,f  --vcfinput --remove

# vep注释MNV(主要为HGVS)
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
--use_given_ref --refseq --assembly GRCh37 --format vcf \
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
--hgvs --hgvsg --transcript_version \
-i anno/${sampleID}.somatic.mnp.annovar.hg19_multianno.vcf -o anno/${sampleID}.somatic.mnp.vep_anno.vcf

# HGVS命名注释MNV
$python3 ${anno_hgvs} anno/${sampleID}.somatic.mnp.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.somatic.mnp.hgvs_anno.vcf

# 合并突变
$bgzip -f anno/${sampleID}.somatic.common.vcf
$tabix -p vcf anno/${sampleID}.somatic.common.vcf.gz

$bgzip -f anno/${sampleID}.somatic.mnp.hgvs_anno.vcf
$tabix -p vcf anno/${sampleID}.somatic.mnp.hgvs_anno.vcf.gz

$bcftools concat anno/${sampleID}.somatic.common.vcf.gz anno/${sampleID}.somatic.mnp.hgvs_anno.vcf.gz -a -d all -o anno/${sampleID}.somatic.merge.vcf

else
mv anno/${sampleID}.somatic.merged_norm.vcf anno/${sampleID}.somatic.merge.vcf
fi

# VEP修复后假阳性及良性变异删除
$python3 ${snv_second_filter} anno/${sampleID}.somatic.merge.vcf -k somatic -o anno/${sampleID}.somatic.final.filter.vcf >>anno/${sampleID}.somatic.filter_reason.txt

# 统计变异删除的原因
$python3 ${get_filter_reason} . -k somatic


###################### tmb ########################
if [[ ${bedType} =~ ^(${tmb_panel})$ ]]; then
echo $'\nTMB'
echo "$(date)"

# 过滤低频假阳性
grep "#" ${sampleID}.somatic.split.norm.vcf > anno/${sampleID}.tmb.filter_1.vcf
if [[ ${sampleType} =~ (ctdna) ]]; then
grep -v "#" ${sampleID}.somatic.split.norm.vcf | awk '{split($NF,a,":");if(a[3]>0.05 && a[3] <= 0.5) print $0}' >> anno/${sampleID}.tmb.filter_1.vcf
else
grep -v "#" ${sampleID}.somatic.split.norm.vcf | awk '{split($NF,a,":");if(a[3]>0.05) print $0}' >> anno/${sampleID}.tmb.filter_1.vcf
fi

# 过滤filter信息大于2的突变
grep "#" anno/${sampleID}.tmb.filter_1.vcf > anno/${sampleID}.tmb.filter_2.vcf
grep -v "#" anno/${sampleID}.tmb.filter_1.vcf | awk '{len=split($7,a,";");if(len<2) print $0}' >> anno/${sampleID}.tmb.filter_2.vcf

# 合并同一单倍型上的复杂突变
$python3 ${merge_mnv} anno/${sampleID}.tmb.filter_2.vcf $hg19fa -o anno/${sampleID}.tmb.merged.vcf

# 将合并后的mnp位置左移校准
$bcftools norm -m -both -f ${hg19fa} anno/${sampleID}.tmb.merged.vcf -o anno/${sampleID}.tmb.merged_norm.vcf
rm anno/${sampleID}.tmb.merged.vcf

# annovar 注释相关数据库
${annovar} anno/${sampleID}.tmb.merged_norm.vcf ${humandb} \
--outfile anno/${sampleID}.tmb.annovar --buildver hg19 \
--protocol refGene,cytoBand,1000g2015aug_all,esp6500siv2_all,exac03nonpsych,cosmic70,clinvar_20220624,intervar_20180118,dbnsfp35a,avsnp150 \
--operation g,r,f,f,f,f,f,f,f,f  --vcfinput --remove
rm anno/${sampleID}.tmb.merged_norm.vcf

# vep注释（主要为HGVS）
${vep} --dir ${vep_dir} --cache --offline --cache_version 98 \
--use_given_ref --refseq --assembly GRCh37 --format vcf \
--fa ${hg19fa} --force_overwrite --vcf --variant_class --gene_phenotype --vcf_info_field ANN \
--hgvs --hgvsg --transcript_version \
-i anno/${sampleID}.tmb.annovar.hg19_multianno.vcf -o anno/${sampleID}.tmb.vep_anno.vcf
rm anno/${sampleID}.tmb.annovar*

# HGVS命名注释
$python3 ${anno_hgvs} anno/${sampleID}.tmb.vep_anno.vcf ${clinic_transcripts} ${refFlat} -o anno/${sampleID}.tmb.hgvs_anno.vcf
rm anno/${sampleID}.tmb.vep_anno.vcf*

# 过滤其他假阳性，germline等
$python3 ${tmb_filter} anno/${sampleID}.tmb.hgvs_anno.vcf -t ${sampleType} -o anno/${sampleID}.tmb.final.filter.vcf -m ${mode} > anno/${sampleID}.tmb.filter_reason.txt
rm anno/${sampleID}.tmb.hgvs_anno.vcf

# 统计变异删除的原因
$python3 ${get_filter_reason} . -k tmb
rm anno/${sampleID}.tmb.filter_1.vcf anno/${sampleID}.tmb.filter_2.vcf
fi
rm ${sampleID}.somatic.split.norm.vcf


###################### concat vcf ########################
if ([[ "$bedType" =~ (${germline_panel1}) ]]) || ([[ "$bedType" =~ (QR) &&  ${sampleType} != "blood" && ${sampleType} != "dna" ]]);then
echo $'\nmerge_vcf'
echo "$(date)"

cp anno/${sampleID}.germline.final.filter.vcf anno/${sampleID}.germline.final.filter.concat.vcf
$bgzip -f anno/${sampleID}.germline.final.filter.concat.vcf
$tabix -p vcf anno/${sampleID}.germline.final.filter.concat.vcf.gz

cp anno/${sampleID}.somatic.final.filter.vcf anno/${sampleID}.somatic.final.filter.concat.vcf
$bgzip -f anno/${sampleID}.somatic.final.filter.concat.vcf
$tabix -p vcf anno/${sampleID}.somatic.final.filter.concat.vcf.gz

if [[ "$bedType" =~ (${tmb_panel}) ]];then
cp anno/${sampleID}.tmb.final.filter.vcf anno/${sampleID}.tmb.final.filter.concat.vcf
$bgzip -f anno/${sampleID}.tmb.final.filter.concat.vcf
$tabix -p vcf anno/${sampleID}.tmb.final.filter.concat.vcf.gz
$bcftools concat anno/${sampleID}.somatic.final.filter.concat.vcf.gz anno/${sampleID}.germline.final.filter.concat.vcf.gz anno/${sampleID}.tmb.final.filter.concat.vcf.gz -a -d all -o anno/${sampleID}.somatic.concat.vcf
rm anno/${sampleID}.tmb.final.filter.concat.vcf.gz*
else
$bcftools concat anno/${sampleID}.somatic.final.filter.concat.vcf.gz anno/${sampleID}.germline.final.filter.concat.vcf.gz -a -d all -o anno/${sampleID}.somatic.concat.vcf
fi
rm anno/${sampleID}.somatic.final.filter.concat.vcf.gz* anno/${sampleID}.germline.final.filter.concat.vcf.gz*
fi

###################### anno drugs ########################
echo $'\nanno_drugs'
echo "$(date)"

# 注释药物相关数据库
if ([[ "$bedType" =~ (${germline_panel1}) ]]) || ([[ "$bedType" =~ (QR) &&  ${sampleType} != "blood" && ${sampleType} != "dna" ]]);then
$python3 ${snv_out} anno/${sampleID}.somatic.concat.vcf -o anno/${sampleID}.somatic.xlsx -t ${sampleType} >> anno/${sampleID}.somatic.filter_reason.txt

elif [[ "$bedType" =~ (${germline_panel2}) && "$sampleType" =~ (blood|dna) ]];then
$python3 ${snv_out} anno/${sampleID}.germline.final.filter.vcf -o anno/${sampleID}.somatic.xlsx -t ${sampleType} >> anno/${sampleID}.somatic.filter_reason.txt

else
$python3 ${snv_out} anno/${sampleID}.somatic.final.filter.vcf -o anno/${sampleID}.somatic.xlsx -t ${sampleType} >> anno/${sampleID}.somatic.filter_reason.txt
fi

###################### add_pipe_tag ######################
# 添加流程pipe tag
$python3 ${add_pipe_tag} anno/${sampleID}.somatic.xlsx

#################### sanger_determine ####################
# 确定germline突变是否需要验证
$python3 ${sanger_determine} anno/${sampleID}.somatic.xlsx

#################### cis_trans ####################
# 确定T790M和G797S共突变的顺式或反式
$python3 ${cis_trans} anno/${sampleID}.somatic.xlsx

###################### core_filter ########################
if [[ "$bedType" =~ (${core_panel}) ]];then
$python3 ${core_filter} anno/${sampleID}.somatic.xlsx
fi

###################### qc_bamdst ########################
echo $'\nqc_bamdst'
echo "$(date)"

if [[ ${bedType} =~ ^(${pcr_panel})$ || ${bedType} =~ ^(${umi_panel})$ ]]; then
mkdir -p qc/bamdst
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 500 -o qc/bamdst ${sampleID}.sorted.bam
ls qc/bamdst/ | xargs -I{} mv qc/bamdst/{} qc/bamdst/${sampleID}.{}

elif [[ ${bedType} =~ ^(${probe_panel})$ && ${sampleType} != "ctdna" ]]; then
mkdir -p qc/bamdst
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 1000 -o qc/bamdst ${sampleID}.sorted.dedup.bam
ls qc/bamdst/ | xargs -I{} mv qc/bamdst/{} qc/bamdst/${sampleID}.{}

else
mkdir -p qc/bamdst/raw
if [[ "$sampleType" =~ (ctdna) ]]; then
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 1000 -o qc/bamdst/raw ${sampleID}.sorted.bam
else
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 500 -o qc/bamdst/raw ${sampleID}.sorted.bam
fi
ls qc/bamdst/raw/ | xargs -I{} mv qc/bamdst/raw/{} qc/bamdst/raw/${sampleID}.raw.{}

mkdir -p qc/bamdst/dedup
if [[ "$sampleType" =~ (ctdna) ]]; then
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 1000 -o qc/bamdst/dedup ${sampleID}.sorted.dedup.bam
else
$bamdst --uncover 30 -p ${bedPath}/${bedType}.bed --cutoffdepth 500 -o qc/bamdst/dedup ${sampleID}.sorted.dedup.bam
fi
ls qc/bamdst/dedup | xargs -I{} mv qc/bamdst/dedup/{} qc/bamdst/dedup/${sampleID}.dedup.{}
fi


###################### chemical ########################

echo $'\nchemical'
echo "$(date)"

if [[ ${bedType} =~ ^(DA)$ ]]; then
$samtools bedcov ${bedPath}/chemical_DA.bed ${sampleID}.sorted.bam >${sampleID}.chemical.depth.txt
else
$samtools bedcov ${bedPath}/chemical_DF.bed ${sampleID}.sorted.bam >${sampleID}.chemical.depth.txt
fi

if [[ "$bedType" =~ (${chemical_panel}) ]]; then
$python3 ${germline_chemical} anno/${sampleID}.germline.anno.hg19_multianno.txt >anno/${sampleID}.chemical.txt
$python3 ${check_sample_data} anno/${sampleID}.chemical.txt
rm ${sampleID}.chemical.depth.txt anno/${sampleID}.chemical.txt
fi


###################### Copy Number Aberration ######################
### hybrid capture: DG|DS|DH|DT|DL|DW|DR|DJ|DF|DQ|DE
### targeted amplicon seq: DP|HH|HT|HW|HR|HC|HF|ZJ|ZG|YW|QJ|YR|YZ|YW|YF
### 需要做cnv的panel：DF, DG, DJ, DQ, DR, DS, DW, DE

##################### CNVkit ########################
if [[  ${bedType} =~ ^(${cnv_panel1})$  ]];
then
echo "CNA starts at -- `date` "

if [[  ${bedType} =~ ^(${cnv_panel2})$  ]];
then
$python3 $cnvkit coverage ${sampleID}.sorted.dedup.bam ${cnv_target}/${bedType}.target.bed -p ${threads} -o ${sampleID}.targetcoverage.cnn
$python3 $cnvkit fix ${sampleID}.targetcoverage.cnn ${cnv_ref}/empty.cnn ${cnv_ref}/${bedType}.flat.cnn -o ${sampleID}.cnr
$python3 $cnvkit segment -p ${threads} ${sampleID}.cnr --drop-low-coverage -o ${sampleID}.cns

elif [[  ${bedType} =~ ^(${cnv_panel3})$  ]];
then
$python3 $cnvkit coverage ${sampleID}.sorted.bam ${cnv_target}/${bedType}.target.bed -p ${threads} -o ${sampleID}.targetcoverage.cnn
$python3 $cnvkit fix ${sampleID}.targetcoverage.cnn ${cnv_ref}/empty.cnn ${cnv_ref}/${bedType}.flat.cnn --no-edge -o ${sampleID}.cnr
$python3 $cnvkit segment -p ${threads} ${sampleID}.cnr --drop-low-coverage -o ${sampleID}.cns

else
$python3 $cnvkit coverage ${sampleID}.sorted.dedup.bam ${cnv_target}/${bedType}.target.bed -p ${threads} -o ${sampleID}.targetcoverage.cnn 
$python3 $cnvkit coverage ${sampleID}.sorted.dedup.bam ${cnv_anti}/${bedType}.antitarget.bed -p ${threads} -o ${sampleID}.antitargetcoverage.cnn 
$python3 $cnvkit fix ${sampleID}.targetcoverage.cnn ${sampleID}.antitargetcoverage.cnn ${cnv_ref}/${bedType}.flat.cnn -o ${sampleID}.cnr
$python3 $cnvkit segment -p ${threads} ${sampleID}.cnr --drop-low-coverage -o ${sampleID}.cns
fi

$python3 $cnvkit call ${sampleID}.cns --center --drop-low-coverage -o ${sampleID}.call.cns
$python3 $cnvkit export vcf ${sampleID}.call.cns -i "${sampleID}" -o ${sampleID}.cnv.vcf
$python3 $annotate_vcf -i ${sampleID}.cnv.vcf -r $refFlat -o ${sampleID}.annotated.cnv.csv

$python3 ${cna_out} ${sampleID}.annotated.cnv.csv

##################### THeTA ########################
# THeTA -- calculate sample purity (tumour/normal ratio)
$python3 $cnvkit export theta ${sampleID}.cns -r ${cnv_ref}/${bedType}.flat.cnn
${theta} ${sampleID}.interval_count --NUM_PROCESSES ${threads}


echo "CNA analysis completed at: `date` "
rm ${sampleID}.*coverage.cnn ${sampleID}.cnr ${sampleID}.cnv.vcf ${sampleID}.cns ${sampleID}.call.cns ${sampleID}.interval_count ${sampleID}.n2.graph.pdf ${sampleID}.n2.withBounds
fi


################### Microsatelite instability ######################
# ###################### MSIsensor2 ######################
# # Detect msi on below panel using msisensor

if [[ ${bedType} =~ ^(${msi_panel})$ ]];
then
echo $'\nmsisensor2'
echo "$(date)"
$msisensor2 msi -M ${msi_models} -t ${sampleID}.sorted.dedup.bam -o ${sampleID}.prefix
fi


########################## Fusion of genes ########################
##################### genefuse ########################
# Detect potential gene fusion on below panels.
if [[ ${bedType} =~ ^(${fusion_panel})$ ]]; then

echo $'\ngenefuse'
echo "$(date)"
$python3 ${make_fusion_gene} ${sampleID}
if [[ "$sampleType" =~ (ctdna) ]]; then
$genefuse -t ${threads} -r ${hg19fa} -f ${sampleID}.fusion_genes.csv -1 ${sampleID}.umi_consensus.fq.gz -h ${sampleID}.fusion.html -j ${sampleID}.fusion.json  1> ${sampleID}.fusion.txt
rm ${sampleID}.umi_consensus.fq.gz
else
$genefuse -t ${threads} -r ${hg19fa} -f ${sampleID}.fusion_genes.csv -1 ${sampleID}.R1_001.trim.fq.gz -2 ${sampleID}.R2_001.trim.fq.gz -h ${sampleID}.fusion.html -j ${sampleID}.fusion.json  1> ${sampleID}.fusion.txt
fi

echo $'\ndelly'
echo "$(date)"
$delly call -x ${delly_low_regions} -g ${hg19fa} -o ${sampleID}.bcf ${sampleID}.sorted.dedup.bam
$bcftools view ${sampleID}.bcf > ${sampleID}.sv.raw.vcf
rm ${sampleID}.bcf*
$python3 $delly_sv_filter ${sampleID}.sv.raw.vcf
rm ${sampleID}.delly*

$python3 $fusion_filter ${sampleID}.fusion.json -o ${sampleID}.fusion_filtered.json
$python3 ${fusion_out} ${sampleID}.fusion_filtered.json
$python3 ${check_fusion_result} ${sampleID}.fusion.xlsx
rm ${sampleID}.fusion.txt
fi


########################## Quality control #########################
##################### Fastp ########################
# This step is used for quality control by fastqc.

mkdir -p qc/fastp
echo $'\nQC_fastp'
echo "$(date)" 
mv ${sampleID}.jion ${sampleID}.html qc/fastp/

##################### samtools ########################
# This step is used for quality control by samtools.

mkdir -p qc/samtools
echo $'\nQC_samtools'
echo "$(date)"
$samtools stats -@ ${threads} ${sampleID}.sorted.bam >qc/samtools/${sampleID}.sorted.bam.stats.txt


##################### hla_type_neo_predict ########################

# This step is used for hla_type_neo_predict by hlahd and MuPeXI.
if [[ "$bedType" =~ ^(${hla_panel})$ ]]; then
echo $'\nHLA'
echo "$(date)"
cp anno/${sampleID}.tmb.final.filter.vcf anno/${sampleID}.tmb_copy.vcf
grep '#' anno/${sampleID}.tmb_copy.vcf > anno/${sampleID}.header
grep -v '#' anno/${sampleID}.tmb_copy.vcf > anno/${sampleID}.body

gunzip -k ${sampleID}.R1_001.trim.fq.gz
gunzip -k ${sampleID}.R2_001.trim.fq.gz
# get 300w reads, about 1GB to to hla typing
$seqtk sample -s100 ${sampleID}.R1_001.trim.fq 3000000 > ${sampleID}.R1_001.trim.sub.fq
$seqtk sample -s100 ${sampleID}.R2_001.trim.fq 3000000 > ${sampleID}.R2_001.trim.sub.fq

$python2 ${hla_optitype} -sampleID $sampleID -outDir $outDir  -results_TMB_txt anno/${sampleID}.body -fq1 ${sampleID}.R1_001.trim.sub.fq -fq2 ${sampleID}.R2_001.trim.sub.fq

# 拼接vcf
cat anno/${sampleID}.header anno/${sampleID}.select_variant_from_TMB.txt > anno/${sampleID}.mupexi.vcf

$python2 ${hla_optitype} -sampleID $sampleID -outDir $outDir -vcf anno/${sampleID}.mupexi.vcf -determine_tnb anno/${sampleID}.body
rm ${sampleID}.R*_001.trim.fq ${sampleID}.R*_001.trim.sub.fq \
anno/${sampleID}.body anno/${sampleID}.tmb_copy.vcf anno/${sampleID}.header anno/${sampleID}.select_variant_from_TMB.txt \
neo/${sampleID}_coverage_plot.pdf neo/${sampleID}.fasta neo/${sampleID}.g \
neo/${sampleID}_netMHCpan* neo/${sampleID}_pepmatch_* neo/${sampleID}_peptide_netMHCinput.txt neo/${sampleID}_reference_peptide_* \
neo/${sampleID}_vep_compatible_vcf.vcf neo/${sampleID}_vep.vep
fi


########################## merge_data #########################
$python3 ${single_sample_qc} . ${sampleID} ${sampleType}
$python3 ${parser_single_sample_data} . ${sampleID}
$python3 ${collect_single_sample_data} . ${sampleID}


rm ${sampleID}.R*_001.trim.fq.gz
if [[ ${bedType} =~ ^(${probe_panel})$ ]]; then
rm ${sampleID}.sorted.bam*
fi

echo "$sampleID -- ending at: $(date)"
echo "pipeline finished!" >finish.txt