prod_reproduction/202307-FeBY2/bak/run_FeBY2_cnv_call_L2.sh.bak

set -e

function usage() {
    echo "Usage: $0 [ -s <sample> ] [ -o <out_dir> ] [ -l <sam_list> ] [ -x <script_dir> ] [ -c <FeBY2_CNV_outdir> ]"
    echo "-s Sample Name"
    echo "-o Output dir"
    echo "-l sample list"
    echo "-x script_dir"
    echo "-c FeBY2 CNV outdir"
}

if [ $# -eq 0  ]
then
    usage
    exit
fi
while getopts ":s:o:l:x:c:h" opt; do
    case $opt in
        s)
            sample=$OPTARG
            ;;
        o)
            out_dir=$OPTARG
            ;;
        l)
            sam_list=$OPTARG
            ;;
        c)
            cnv_dir=$OPTARG
            ;;
        x)
            feby2_script=$OPTARG
            ;;
        h)
            usage
            exit
            ;;
        \?)
            echo "Invalid option: -$OPTARG" > /dev/stderr
            usage
            exit
            ;;
    esac
done

if [ ! -f $cnv_dir/ref.list ]
then
    echo "Referenc Sample List:$cnv_dir/ref.list does NOT exist!"
    echo "Please Check."
    cp $sam_list ref.list
fi
if [ -d $cnv_dir ]
then
    echo $cnv_dir.
else
    echo $cnv_dir " NOT exist!"
    exit
fi
cd $cnv_dir
echo "out_dir:" $out_dir
. $feby2_script/conf/feby2.cfg 

sex=`echo ${sample}_?_QC_plot.pdf |sed "s/${sample}_//;s/_QC_plot.pdf//" `
echo sex is $sex.
test -d ${sample}_cnv || mkdir ${sample}_cnv
test -e ${sample}_cnv/$sample.normals.txt && rm ${sample}_cnv/$sample.normals.txt

for sam in ` cat ref.list |grep -v $sample `
do
echo $sam.pcov | cat  >> ${sample}_cnv/$sample.normals.txt
done

$gatkb CombineReadCounts -inputList ${sample}_cnv/$sample.normals.txt -O ${sample}_cnv/$sample.combined-normals.tsv
$gatkb CreatePanelOfNormals -I ${sample}_cnv/$sample.combined-normals.tsv -O ${sample}_cnv/$sample.normals.pon -noQC --disableSpark --minimumTargetFactorPercentileThreshold 0.2

ref_cov=${sample}_cnv/$sample.normals.pon
echo $ref_cov
$gatkb NormalizeSomaticReadCounts -I ${sample}.pcov -PON ${ref_cov} -PTN ${sample}_cnv/${sample}.ptn.tsv -TN ${sample}_cnv/${sample}.tn.tsv
$gatkb PerformSegmentation -TN ${sample}_cnv/${sample}.tn.tsv -O ${sample}_cnv/${sample}.seg -LOG
$gatkb CallSegments -TN ${sample}_cnv/${sample}.tn.tsv -S ${sample}_cnv/${sample}.seg -O ${sample}_cnv/${sample}.called
$gatkb PlotSegmentedCopyRatio -TN ${sample}_cnv/${sample}.tn.tsv -PTN ${sample}_cnv/${sample}.ptn.tsv -S ${sample}_cnv/${sample}.seg -O ${sample}_cnv -pre ${sample} -SD ${b37_dict} -LOG

perl $exe_dir/seg_merge.pl ${sample}_cnv/${sample}.seg
awk 'NR>1 && ($6 > 1.38  || $6 < 0.62) ' ${sample}_cnv/${sample}.seg.xls |awk -v sex=$sex '{print $1"\t"sex"\t"$2"\t"$3"\t"$4"\t"$5"\t"$6}' > ${sample}-cnv.input.tsv
echo -e "Sample\tGender\tChr\tStart\tEnd\tLength\tCopy-ratio\tBand_region" |cat - ${sample}-cnv.input.tsv > ${sample}-${sex}-cnv.input.xls && rm ${sample}-cnv.input.tsv
$anno_cnv -name ${sample} -outdir . -clean ${sample}-${sex}-cnv.input.xls
perl $exe_dir/trim_cnv_result.pl $sample && mv $sample.cnv.anno.xls $sample-?-cnv.input.xls ${sample}_cnv/

## plot
cp ${sample}_cnv/$sample.ptn.tsv .
cp ${sample}_cnv/$sample.tn.tsv .
cp ${sample}_cnv/$sample.seg .
cp $out_dir/${sample}/conest/${sample}.ps.tab .
perl $exe_dir/pre_plot.pl ${sample}
/cgdata/guocq/soft/R-3.6.0/bin/Rscript $exe_dir/plot_panel_CR.r ${sample}
convert -density 300 $sample.CR_VAF.pdf $sample.CR_VAF.png && rm -r $sample.ptn.tsv $sample.tn.tsv $sample.seg ${sample}.ps.tab $sample.CR_VAF.pdf ${sample}_plot.tsv