#!/bin/bash
BAM="$1"
sampleID="$2"
REFERENCE="/dssg06/InternalResearch06/zhangjp/MT/MT.fa"
BED="/dssg/home/zhangjp/soft/mutserve/install/MT.bed"
if [ "$#" -ne 2 ]; then
    echo "Usage: $0 <BAM file> <Sample ID>"
    exit 1
fi
sambamba sort -t 10 -n $BAM  -o "${sampleID}.sort.name.bam" --tmpdir  /dssg06/InternalResearch06/zhangjp/MT/tmp_test
samtools fastq -1  "${sampleID}_R1.MT.fastq" -2 "${sampleID}_R2.MT.fastq" -0 /dev/null -s /dev/null -n -F 0x900 "${sampleID}.sort.name.bam" -@ 10

#bedtools bamtofastq -i  "${sampleID}.sort.name.bam" -fq "${sampleID}_R1.MT.fastq" -fq2 "${sampleID}_R2.MT.fastq"
# 比对到线粒体参考基因组
# 将 SAM 文件转换为 BAM 文件
/dssg06/InternalResearch06/zhangjp/soft/bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -t 10 "$REFERENCE" "${sampleID}_R1.MT.fastq" "${sampleID}_R2.MT.fastq" | samtools view -bS - > "${sampleID}.bam"
# 对 BAM 文件进行排序
sambamba sort -t 10 -o "${sampleID}.sorted.bam" "${sampleID}.bam" --tmpdir  /dssg06/InternalResearch06/zhangjp/MT/tmp_test
# 为排序后的 BAM 文件创建索引
sambamba markdup  --overflow-list-size=600000 -t 10 "${sampleID}.sorted.bam" "${sampleID}.dedup.bam" --remove-duplicates
sambamba depth region -t 8 -L  "$BED"   "${sampleID}.dedup.bam" > "${sampleID}.region_depth.txt"
### calling
 /dssg/home/zhangjp/soft/openlogic-openjdk-11.0.26+4-linux-x64/bin/java -jar /dssg/home/zhangjp/soft/mutserve/install/mutserve.jar call "${sampleID}.dedup.bam"  --output "${sampleID}.vcf" --reference  "${REFERENCE}"  --threads 4 --mode mtdna --deletions --insertions  --alignQ 30  --baseQ 20
###annotate
/dssg/home/zhangjp/soft/openlogic-openjdk-11.0.26+4-linux-x64/bin/java -jar /dssg/home/zhangjp/soft/mutserve/install/mutserve.jar annotate --input "${sampleID}.txt"  --output  "${sampleID}.ann.txt" --annotation /dssg/home/zhangjp/soft/mutserve/install//rCRS_annotation_2020-08-20.txt
###Feature
python /dssg/home/zhangjp/soft/mutserve/install/MT.Feature.py "${sampleID}.ann.txt"   "${sampleID}.region_depth.txt"  $sampleID
###Feature
### 清理中间文件
rm "${sampleID}.bam"
rm "${sampleID}.sorted.bam"
rm "${sampleID}.sort.name.bam"
rm "${sampleID}.dedup.bam"
rm "${sampleID}_R1.MT.fastq"
rm "${sampleID}_R2.MT.fastq"