stereoPipeline_pt_v7.1.sh

#!/bin/bash
set -e

if [[ $# -lt 12 ]];then
    echo "usage: sh $0 -splitCount -maskFile -rnaFq1 -rnaFq2 -adtFq1 -adtFq2 -proteinList -refIndex -genomeFile -speciesName -tissueType -annotationFile -outDir -imageRecordFile -imageCompressedFile -doCellBin -rRNAremove -pidStart -threads -sif
    -splitCount : count of splited stereochip mask file, usually 16 for Q4 fq data and 1 for Q40 fq data
    -maskFile : stereochip mask file
    -rnaFq1 : RNA fastq file path of read1, if there are more than one fastq file, please separate them with comma, e.g:lane1_read_1.fq.gz,lane2_read_1.fq.gz
    -rnaFq2 : RNA fastq file path of read2, if there are more than one fastq file, please separate them with comma, not requested for Q4 fastq data, e.g:lane1_read_2.fq.gz,lane2_read_2.fq.gz
    -adtFq1 : ADT fastq file path of read1, if there are more than one fastq file, please separate them with comma, e.g:lane1_read_1.fq.gz,lane2_read_1.fq.gz
    -adtFq2 : ADT fastq file path of read2, if there are more than one fastq file, please separate them with comma, not requested for Q4 fastq data, e.g:lane1_read_2.fq.gz,lane2_read_2.fq.gz
    -proteinList : protein list file which contain protein id, sequences and names.
    -refIndex : reference genome indexed folder, please build index before SAW analysis run
    -annotationFile :  annotations file in gff or gtf format, the file must contain gene and exon annotations
    -speciesName : specie of the sample
    -tissueType : tissue type of the sample
    -outDir : output directory path
    -imageRecordFile : image file(*.ipr) generated by ImageStudio software, not requested
    -imageCompressedFile : image file(*.tar.gz) generated by ImageStudio software, not requested
    -doCellBin : [Y/N]
    -rRNAremove : [Y/N]
    -pidStart : PID start position. 21 for sequencing transcriptome and ADT libraries together. 0 for sequencing separately.
    -threads : the number of threads to be used in running the pipeline
    -sif : the file format of the visual software
    "
    exit
fi

while [[ -n "$1" ]]
do
    case "$1" in
        -splitCount) splitCnt="$2"
            shift ;;
        -maskFile) maskFile="$2"
            shift ;;
        -rnaFq1) rnaRead1="$2"
            shift ;;
        -rnaFq2) rnaRead2="$2"
            shift ;;
        -adtFq1) adtRead1="$2"
            shift ;;
        -adtFq2) adtRead2="$2"
            shift ;;
        -proteinList) proteinList="$2"
            shift ;;
        -refIndex) GDir="$2"
            shift ;;
        -speciesName) refName="$2"
            shift ;;
        -tissueType) tissueType="$2"
            shift ;;
        -annotationFile) annoFile="$2"
            shift ;;
        -outDir) outDir="$2"
            shift ;;
        -imageRecordFile) iprFile="$2"
            shift ;;
        -imageCompressedFile) imageTarFile="$2"
            shift ;;
        -doCellBin) doCell="$2"
            shift ;;
        -rRNAremove) rRNAremove="$2"
            shift ;;
        -pidStart) pidStart="$2"
            shift ;;
        -threads) threads="$2"
            shift ;;
        -sif) sif="$2"
            shift ;;
    esac
        shift
done


#software check
if [ `command -v singularity` ]
then
    singularityPath=`command -v singularity`
    echo `date` " singularity check: pass, and singularity path is ${singularityPath}"
else
    echo `date` " singularity check: singularity does not exits, please verify that you have installed singularity and exported it to your system PATH variable"
    exit
fi

if [[ -n $sif ]]
then
    echo `date` " singularity image file check: file exist and SIF path is ${sif}"
else
    echo `date` " singularity image file check: file does not exist, please double check your SIF file is in the current directory or the path given by the option -s is valid."
fi


if [[ ! -d $outDir ]];then
    mkdir -p $outDir
fi

# Get basic information
maskname=$(basename $maskFile)
SN=${maskname%%.*}

maskDIR=$(dirname $maskFile)
annoDIR=$(dirname $annoFile)
refDIR=$(dirname $GDir)
proteinListDIR=$(dirname $proteinList)
if [[ $iprFile ]] && [[ $imageTarFile ]];then
    iprDIR=$(dirname $iprFile)
    imgTarDIR=$(dirname $imageTarFile)
fi

# Prepare output directories
arr_result=( ${outDir}/00T.mapping ${outDir}/00P.mapping ${outDir}/01T.merge ${outDir}/02T.count ${outDir}/03.calibration ${outDir}/05T.tissuecut ${outDir}/05P.tissuecut ${outDir}/06T.spatialcluster ${outDir}/06P.spatialcluster ${outDir}/07T.saturation ${outDir}/07P.saturation ${outDir}/08.multiomics ${outDir}/09.report ${outDir}/visualization)
for each in "${arr_result[@]}";
do
    if [[ ! -d $each ]];then
        mkdir -p $each
    fi
done

if [[ $iprFile ]] && [[ $imageTarFile ]];then
    mkdir -p ${outDir}/04.register
fi

if [[ $doCell == "Y" ]]; then
    mkdir -p ${outDir}/051T.cellcut
    mkdir -p ${outDir}/061T.cellcluster
    mkdir -p ${outDir}/051P.cellcut
    mkdir -p ${outDir}/061P.cellcluster
fi

#Run splitMask or CIDCount for preparation
#  ulimit -c 100000000000
echo `date` "=> splitMask, compute CID count and predict the memory of mapping start......"
rnaRead1List=(`echo $rnaRead1 | tr ',' ' '`)
fqbases=()
starBams=()
bcStat=()
bcLogFinalOut=()
bcReadsCounts=()
GSize=(`du -sh --block-size=G ${GDir}/Genome | cut -f 1`)

if [[ ! -n "$rnaRead2" ]]; then
    fqType="Q4"
    arr_result=( ${outDir}/00T.mapping/splitBin ${outDir}/00T.mapping/mergeList )
    for each in "${arr_result[@]}";do
        if [[ ! -d $each ]];then mkdir -p $each; fi
    done
    export SINGULARITY_BIND=$outDir,$maskDIR,$annoDIR,$refDIR
    /usr/bin/time -v singularity exec ${sif} splitMask \
        ${maskFile} ${outDir}/00T.mapping/splitBin $threads $splitCnt 2_25
    for ((i=1;i<=$splitCnt;i++)); do
        if [[ $(echo ${#i}) == '1' ]];then a=0$i; else a=$i;fi
        echo $rnaRead1 | sed 's/,/\n/g' | grep _$i.fq.gz > ${outDir}/00T.mapping/mergeList/$a.${SN}.Q4.fq.list
    done
else
    fqType="Q40"
    rnaRead2List=(`echo $rnaRead2 | tr ',' ' '`)
    fqNumber=`echo ${#rnaRead1List[@]}`
    export SINGULARITY_BIND=$outDir,$maskDIR,$annoDIR,$refDIR
    /usr/bin/time -v singularity exec ${sif} CIDCount \
        -i ${maskFile} \
        -s ${refName} \
        -g ${GSize} > ${outDir}/00T.mapping/CIDCount
fi
echo "Your rna sequencing reads are in ${fqType} format."

## Run SAW mapping to perform CID mapping and STAR alignment.
echo `date` "=> CID mapping, adapter filtering and RNA alignment start......"
if [[ $fqType == 'Q40' ]]; then
    for ((i=0;i<=`expr $(echo $fqNumber) - 1`;i++)); do
        fqname=$(basename ${rnaRead1List[i]})
        fqdir=$(dirname ${rnaRead1List[i]})
        fqbase=${fqname%%.*}
        fqbases[i]=$fqbase
        bcPara=${outDir}/00T.mapping/${fqbase}.bcPara
        barcodeReadsCount=${outDir}/00T.mapping/${fqbase}.barcodeReadsCount.txt
        echo  " ~~~ mapping - $fqname ~~~"
        echo "in=${maskFile}" > $bcPara
        echo "in1=${rnaRead1List[i]}" >> $bcPara
        echo "in2=${rnaRead2List[i]}" >> $bcPara
        echo "barcodeReadsCount=${barcodeReadsCount}" >> $bcPara
        echo "barcodeStart=0" >> $bcPara
        echo "barcodeLen=25" >> $bcPara
        echo "umiStart=25" >> $bcPara
        echo "umiLen=10" >> $bcPara
        echo "mismatch=1" >> $bcPara
        echo "bcNum=`head -1 ${outDir}/00T.mapping/CIDCount`" >> $bcPara
        echo "polyAnum=15" >> $bcPara
        echo "mismatchInPolyA=2" >> $bcPara
        if [[ $rRNAremove == "Y" ]]; then
            echo "rRNAremove" >> $bcPara
        fi
        rnaRead1DIR=$(dirname ${rnaRead1List[i]})
        rnaRead2DIR=$(dirname ${rnaRead2List[i]})
        export SINGULARITY_BIND=$rnaRead1DIR,$rnaRead2DIR,$outDir,$maskDIR,$annoDIR,$refDIR
        /usr/bin/time -v singularity exec ${sif} mapping \
            --outSAMattributes spatial \
            --outSAMtype BAM SortedByCoordinate \
            --genomeDir ${GDir} \
            --runThreadN ${threads} \
            --outFileNamePrefix ${outDir}/00T.mapping/${fqbase}. \
            --sysShell /bin/bash \
            --stParaFile ${bcPara} \
            --readNameSeparator \" \" \
            --limitBAMsortRAM 63168332971 \
            --limitOutSJcollapsed 10000000 \
            --limitIObufferSize=280000000 \
            --outBAMsortingBinsN 50 \
            --outSAMmultNmax 1 \
            > ${outDir}/00T.mapping/${fqbase}.run.log

        starBam=${outDir}/00T.mapping/${fqbase}.Aligned.sortedByCoord.out.bam
        starBams[i]=$starBam
        bcStat[i]=${outDir}/00T.mapping/${fqbase}.CIDMap.stat
        bcFinalOut[i]=${outDir}/00T.mapping/${fqbase}.Log.final.out
        bcReadsCounts[i]=$barcodeReadsCount
    done
elif [[ $fqType == 'Q4' ]]; then
    for ((i=1;i<=$splitCnt;i++)); do
        if [[ $(echo ${#i}) == '1' ]];then a=0$i; else a=$i;fi
        /usr/bin/time -v singularity exec ${sif} CIDCount \
            -i $(ls ${outDir}/00T.mapping/splitBin/${a}.${SN}.barcodeToPos.bin) \
            -s ${refName} \
            -g ${GSize} > ${outDir}/00T.mapping/CIDCount
        fqbase=$a.${SN}.Q4
        bcPara=${outDir}/00T.mapping/${fqbase}.bcPara
        barcodeReadsCount=${outDir}/00T.mapping/${fqbase}.barcodeReadsCount.txt
        echo  " ~~~ mapping Q4 FASTQ Data - $a ~~~"
        echo "in=$(ls ${outDir}/00T.mapping/splitBin/${a}.${SN}.barcodeToPos.bin)" > $bcPara
        rnaRead1List=${outDir}/00T.mapping/mergeList/$a.${SN}.Q4.fq.list
        echo "in1=${rnaRead1List}" >> $bcPara
        rnaRead1DIR=$(dirname $(cat $rnaRead1List)|tr '\n' ',')
        echo "barcodeReadsCount=${barcodeReadsCount}" >> $bcPara
        echo "barcodeStart=0" >> $bcPara
        echo "barcodeLen=24" >> $bcPara
        echo "umiStart=25" >> $bcPara
        echo "umiLen=10" >> $bcPara
        echo "mismatch=1" >> $bcPara
        echo "bcNum=`head -1 ${outDir}/00T.mapping/CIDCount`" >> $bcPara
        echo "polyAnum=15" >> $bcPara
        echo "mismatchInPolyA=2" >> $bcPara
        if [[ $rRNAremove == "Y" ]]; then
            echo "rRNAremove" >> $bcPara
        fi
        export SINGULARITY_BIND=$rnaRead1DIR,$outDir,$maskDIR,$annoDIR,$refDIR
        /usr/bin/time -v singularity exec ${sif} mapping \
            --outSAMattributes spatial \
            --outSAMtype BAM SortedByCoordinate \
            --genomeDir ${GDir} \
            --runThreadN ${threads} \
            --outFileNamePrefix ${outDir}/00T.mapping/${fqbase}. \
            --sysShell /bin/bash \
            --stParaFile ${bcPara} \
            --readNameSeparator \" \" \
            --limitBAMsortRAM 63168332971 \
            --limitOutSJcollapsed 10000000 \
            --limitIObufferSize=280000000 \
            --outBAMsortingBinsN 50 \
            --outSAMmultNmax 1 \
            > ${outDir}/00T.mapping/${fqbase}.run.log

        starBams[i]=${outDir}/00T.mapping/${fqbase}.Aligned.sortedByCoord.out.bam
        bcStat[i]=${outDir}/00T.mapping/${fqbase}.CIDMap.stat
        bcFinalOut[i]=${outDir}/00T.mapping/${fqbase}.Log.final.out
        bcReadsCounts[i]=$barcodeReadsCount
    done
fi

if [[ $(echo ${#bcReadsCounts[*]}) == '1' ]]; then
    bcReadsCountsStr=$bcReadsCounts
    starBamsStr=${starBams[0]}
    bcFinalOutStr=${bcFinalOut[0]}
    bcStatStr=${bcStat[0]}
else
    bcReadsCountsStr=$( IFS=','; echo "${bcReadsCounts[*]}" )
    starBamsStr=$( IFS=','; echo "${starBams[*]}" )
    bcFinalOutStr=$( IFS=','; echo "${bcFinalOut[*]}" )
    bcStatStr=$( IFS=','; echo "${bcStat[*]}" )
fi


## RUN SAW merge to merge barcodeReadsCount file
echo `date` "=> merge barcode reads count tables start......"
barcodeReadsCounts=${outDir}/01T.merge/${SN}.merge.barcodeReadsCount.txt
export SINGULARITY_BIND=$outDir,$maskDIR
if [[ $fqType == 'Q4' ]] && [[ $(echo ${#bcReadsCounts[*]}) > '1' ]]; then
    echo 'Q4'
    /usr/bin/time -v singularity exec ${sif} merge \
        ${maskFile} \
        $bcReadsCountsStr \
        $barcodeReadsCounts
elif [[ $fqType == 'Q40' ]]; then
    echo 'Q40'
    if [[ $(echo ${#bcReadsCounts[*]}) == '1' ]]
    then
        cp $bcReadsCountsStr $barcodeReadsCounts
    else
        /usr/bin/time -v singularity exec ${sif} merge \
            ${maskFile} \
            $bcReadsCountsStr \
            $barcodeReadsCounts
    fi
fi

## Run SAW count to do annotation, deduplication, and generate gene expression matrix
echo `date` "=> annotation, deduplication, and generate gene expression matrix start......"
export SINGULARITY_BIND=$outDir,$annoDIR,$refDIR
export HDF5_USE_FILE_LOCKING=FALSE
/usr/bin/time -v singularity exec ${sif} count \
    -i ${starBamsStr} \
    -o ${outDir}/02T.count/${SN}.Aligned.sortedByCoord.out.merge.q10.dedup.target.bam \
    -a ${annoFile} \
    -s ${outDir}/02T.count/${SN}.Aligned.sortedByCoord.out.merge.q10.dedup.target.bam.summary.stat \
    -e ${outDir}/02T.count/${SN}.raw.gef \
    --sat_file ${outDir}/02T.count/${SN}_raw_barcode_gene_exp.txt \
    --umi_on \
    --umi_len 10 \
    --save_lq \
    --save_dup \
    --sn ${SN} \
    -c ${threads}

## Run SAW mapping-SP to perform CID mapping, MID filtration and PID alignment
echo `date` "=> CID mapping and PID alignment start......"
echo `date` "=> Please make sure your protein list only records actual input proteins"
echo $adtRead1|sed 's/,/\n/g'>${outDir}/00P.mapping/adtFq1.list
adtRead1List=${outDir}/00P.mapping/adtFq1.list
adtRead1DIR=$(dirname $(cat $adtRead1List)|tr '\n' ',')

if [[ ! -n "$adtRead2"  ]]; then
    adtFqType="Q4"
    echo "Your sequencing reads are in ${adtFqType} format."
    find ${outDir}/00T.mapping/splitBin/*${SN}.barcodeToPos.bin>${outDir}/00P.mapping/splitMask.txt
    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$adtRead1DIR,$outDir,$maskDIR,${proteinListDIR}
    /usr/bin/time -v singularity exec ${sif} mapping-SP \
        --thread ${threads} \
        --cidStart 0 --cidLen 24 --cidMismatch 1 \
        --midStart 25 --midLen 10 \
        --pidStart ${pidStart} --pidLen 15 --pidMismatch 1 \
        --mask ${outDir}/00P.mapping/splitMask.txt \
        --fastq ${outDir}/00P.mapping/adtFq1.list \
        --reference ${proteinList} \
        --output ${outDir}/00P.mapping \
        --sn ${SN}
else
    adtFqType="Q40"
    echo "Your sequencing reads are in ${adtFqType} format."
    echo "~~~ mapping Q40 ${adtRead1} ~~~"
    echo $adtRead2|sed 's/,/\n/g'>${outDir}/00P.mapping/adtFq2.list
    adtRead2List=${outDir}/00P.mapping/adtFq2.list
    adtRead2DIR=$(dirname $(cat $adtRead2List)|tr '\n' ',')
    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$adtRead1DIR,$adtRead2DIR,$outDir,$maskDIR,${proteinListDIR}
    /usr/bin/time -v singularity exec ${sif} mapping-SP \
        --thread ${threads} \
        --cidStart 0 --cidLen 25 --cidMismatch 1 \
        --midStart 25 --midLen 10 \
        --pidStart ${pidStart} --pidLen 15 --pidMismatch 1 \
        --mask ${maskFile} \
        --fastq ${outDir}/00P.mapping/adtFq1.list \
        --fastq2 ${outDir}/00P.mapping/adtFq2.list \
        --reference ${proteinList} \
        --output ${outDir}/00P.mapping \
        --sn ${SN}
fi
gzip -f ${outDir}/00P.mapping/${SN}.protein.raw.gem

## Run SAW calibration to do multi-GEF alignment
echo `date` "=> calibration start......."
export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} calibration \
    -i ${outDir}/00P.mapping/${SN}.protein.raw.gef,${outDir}/02T.count/${SN}.raw.gef \
    -o ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef,${outDir}/03.calibration/${SN}.calibrated.raw.gef \
    -O Proteomics,Transcriptomics

## Convert calibrated GEF to GEM
export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} cellCut view \
    -s ${SN} \
    -i ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
    -o ${outDir}/03.calibration/${SN}.calibrated.gem
gzip -f ${outDir}/03.calibration/${SN}.calibrated.gem

export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} cellCut view \
    -s ${SN} \
    -i ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef \
    -o ${outDir}/03.calibration/${SN}.protein.calibrated.gem
gzip -f ${outDir}/03.calibration/${SN}.protein.calibrated.gem

## Run SAW register to stitch microscope tile images to a panoramic image, perform tissue and cell (optional, depends on -doCellBin) segmentation, and register the panoramic image and the segmentated images with the gene expression matrix.
if [[ -f $imageTarFile ]] && [[ -f $iprFile ]]  && [[ $doCell == "Y" ]]; then
    # Run SAW register (stitch, tissue segmentation, cell segmentation) + SAW imageTools
    echo `date` "=> image processing and registration start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    imgTarDIR=$(dirname $imageTarFile)
    iprDIR=$(dirname $iprFile)
    export SINGULARITY_BIND=$outDir,$imgTarDIR,$iprDIR
    out_iprFile=$(find ${outDir}/04.register -maxdepth 1 -name \*.ipr | head -1)
    /usr/bin/time -v singularity exec ${sif} register \
        -i ${imageTarFile} \
        -c ${iprFile} \
        -v ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
        -o ${outDir}/04.register \
        -w True --core ${threads}

    out_iprFile=$(find ${outDir}/04.register -maxdepth 1 -name \*.ipr | head -1)
    /usr/bin/time -v singularity exec ${sif} imageTools ipr2img \
        -i ${imageTarFile} \
        -c ${out_iprFile} \
        -d tissue cell \
        -r True \
        -o ${outDir}/04.register

    registerTif=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif)
    regTifStr=$(echo $registerTif | tr ' ' ',')
    regGroup=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif -exec sh -c 'for f do basename -- "$f" _fov_stitched_transformed.tif;done' sh {} +)
    regGroupStr=$(echo $regGroup | sed 's/ \|$/\/Image,/g' | sed 's/.$//')
    echo $regTifStr
    echo $regGroupStr
    /usr/bin/time -v singularity exec ${sif} imageTools img2rpi \
        -i ${regTifStr} \
        -g ${regGroupStr} \
        -b 1 10 50 100 \
        -o ${outDir}/04.register/fov_stitched_transformed.rpi
elif [[ -f $imageTarFile ]] && [[ -f $iprFile ]]  && [[ $doCell == "N" ]]; then
    # Run SAW rapidRegister (stitch, tissue segmentation) + SAW imageTools
    echo `date` "=> image processing and registration start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    imgTarDIR=$(dirname $imageTarFile)
    iprDIR=$(dirname $iprFile)

    export SINGULARITY_BIND=$outDir,$imgTarDIR,$iprDIR
    /usr/bin/time -v singularity exec ${sif} register \
        -i ${imageTarFile} \
        -c ${iprFile} \
        -v ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
        -o ${outDir}/04.register \
        -w False --core ${threads}
    out_iprFile=$(find ${outDir}/04.register -maxdepth 1 -name \*.ipr | head -1)
    /usr/bin/time -v singularity exec ${sif} imageTools ipr2img \
        -i ${imageTarFile} \
        -c ${out_iprFile} \
        -d tissue \
        -r True \
        -o ${outDir}/04.register
    registerTif=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif)
    regTifStr=$(echo $registerTif | tr ' ' ',')
    regGroup=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif -exec sh -c 'for f do basename -- "$f" _fov_stitched_transformed.tif;done' sh {} +)
    regGroupStr=$(echo $regGroup | sed 's/ \|$/\/Image,/g' | sed 's/.$//')
    echo $regTifStr
    echo $regGroupStr
    /usr/bin/time -v singularity exec ${sif} imageTools img2rpi \
        -i ${regTifStr} \
        -g ${regGroupStr} \
        -b 1 10 50 100 \
        -o ${outDir}/04.register/fov_stitched_transformed.rpi
fi


# Run SAW tissueCut
if [[ -f $imageTarFile ]] && [[ -f $iprFile ]]; then
    # Run tissueCut to get the spatial gene expression profile of the tissue-covered region
    nucleusLayer=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif -exec sh -c 'for f do basename -- "$f" _fov_stitched_transformed.tif;done' sh {} + | grep -v IF | awk -F_ '{print$1}')
    tissueMaskFile=$(find ${outDir}/04.register -maxdepth 1 -name ${nucleusLayer}_${SN}_tissue_cut.tif)

    echo `date` "=> tissueCut (Transcriptomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$outDir
    /usr/bin/time -v singularity exec ${sif} tissueCut \
        -i ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
        --dnbfile ${barcodeReadsCounts} \
        -s ${tissueMaskFile} \
        --sn ${SN} \
        -O Transcriptomics \
        -d -t 8 \
        -o ${outDir}/05T.tissuecut

    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$outDir
    /usr/bin/time -v singularity exec ${sif} cellCut bgef \
        -i ${outDir}/05T.tissuecut/${SN}.tissue.gef \
        -o ${outDir}/05T.tissuecut/${SN}.${nucleusLayer}.gef \
        -O Transcriptomics \
        -b 1,5,10,20,50,100,150,200

    echo `date` "=> tissueCut (Proteomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$outDir
    validCidReads=${outDir}/00P.mapping/${SN}_valid_cid_reads.tsv
    /usr/bin/time -v singularity exec ${sif} tissueCut \
        -i ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef \
        --dnbfile ${validCidReads} \
        -s ${tissueMaskFile} \
        --sn ${SN} \
        -O Proteomics \
        -d -t 8 \
        -o ${outDir}/05P.tissuecut

    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$outDir
    /usr/bin/time -v singularity exec ${sif} cellCut bgef \
        -i ${outDir}/05P.tissuecut/${SN}.protein.tissue.gef \
        -o ${outDir}/05P.tissuecut/${SN}.protein.${nucleusLayer}.gef \
        -O Proteomics \
        -b 1,5,10,20,50,100,150,200
else
    # Run SAW tissueCut based on the gene expression matrix directly
    echo `date` "=> tissueCut (Transcriptomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    export SINGULARITY_BIND=$outDir,$annoDIR,$refDIR
    echo `date` "=> there is no image, tissueCut (Transcriptomics) based on the gene expression matrix start......."
    /usr/bin/time -v singularity exec ${sif} tissueCut \
        -i ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
        --dnbfile ${barcodeReadsCounts} \
        --sn ${SN} \
        -O Transcriptomics \
        -d -t 8 \
        -o ${outDir}/05T.tissuecut

    validCidReads=${outDir}/00P.mapping/${SN}_valid_cid_reads.tsv
    echo `date` "=> tissueCut (Proteomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    /usr/bin/time -v singularity exec ${sif} tissueCut \
        -i ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef \
        --dnbfile ${validCidReads} \
        --sn ${SN} \
        -O Proteomics \
        -d -t 8 \
        -o ${outDir}/05P.tissuecut
fi

## Complete raw GEF to visual GEF
export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} cellCut bgef \
    -i ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
    -o ${outDir}/03.calibration/${SN}.gef \
    -O Transcriptomics \
    -b 1,5,10,20,50,100,150,200

export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} cellCut bgef \
    -i ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef \
    -o ${outDir}/03.calibration/${SN}.protein.gef \
    -O Proteomics \
    -b 1,5,10,20,50,100,150,200

## Convert GEF to GEM
export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} cellCut view \
    -s ${SN} \
    -i ${outDir}/05T.tissuecut/${SN}.tissue.gef \
    -o ${outDir}/05T.tissuecut/${SN}.tissue.gem
gzip -f ${outDir}/05T.tissuecut/${SN}.tissue.gem

export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} cellCut view \
    -s ${SN} \
    -i ${outDir}/05P.tissuecut/${SN}.protein.tissue.gef \
    -o ${outDir}/05P.tissuecut/${SN}.protein.tissue.gem
gzip -f ${outDir}/05P.tissuecut/${SN}.protein.tissue.gem


# Run spatialCluster for Transcriptomics and Proteomics
binSize=200
resolution=1.0
export SINGULARITY_BIND=$outDir
export HDF5_USE_FILE_LOCKING=FALSE
mkdir -p ${outDir}/tmp
export NUMBA_CACHE_DIR=${outDir}/tmp
export MPLCONFIGDIR=${outDir}/tmp

echo `date` "=> spatialCluster (Transcriptomics) start......."
/usr/bin/time -v singularity exec ${sif} spatialCluster \
    -i ${outDir}/05T.tissuecut/${SN}.tissue.gef \
    -o ${outDir}/06T.spatialcluster/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad \
    -s ${binSize} \
    -r ${resolution}

echo `date` "=> spatialCluster (Proteomics) start......."
export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
/usr/bin/time -v singularity exec ${sif} spatialCluster-SP \
    -i ${outDir}/05P.tissuecut/${SN}.protein.tissue.gef \
    -o ${outDir}/06P.spatialcluster/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad \
    -s ${binSize}

# Run cellCut, cellCluster and cellChunk for Transcriptomics
export SINGULARITY_BIND=$outDir
if [[ -f $imageTarFile ]] && [[ -f $iprFile ]]  && [[ $doCell == "Y" ]]; then
    echo `date` "=> cellCut (Transcriptomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE

    nucleusLayer=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif -exec sh -c 'for f do basename -- "$f" _fov_stitched_transformed.tif;done' sh {} + | grep -v IF | awk -F_ '{print$1}')
    nucleusMask=$(find ${outDir}/04.register -maxdepth 1 -name ${nucleusLayer}_${SN}_mask.tif)

    /usr/bin/time -v singularity exec ${sif} cellCut cgef \
        -i ${outDir}/03.calibration/${SN}.calibrated.raw.gef \
        -m ${nucleusMask} \
        -o ${outDir}/051T.cellcut/${SN}.cellbin.gef

    echo `date` "=> cellCorrect (Transcriptomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    /usr/bin/time -v singularity exec ${sif} cellCorrect \
            -i ${outDir}/03.calibration/${SN}.calibrated.raw.gef  \
            -m ${nucleusMask} \
            -o ${outDir}/051T.cellcut
    gzip -f ${outDir}/051T.cellcut/${SN}.adjusted.gem

    ## Write the cellCorrect mask image into SN.rpi
    cellCorrectMask=$(find ${outDir}/051T.cellcut -maxdepth 1 -name ${nucleusLayer}\*_mask_edm_dis\*.tif)
    /usr/bin/time -v singularity exec ${sif} imageTools img2rpi \
        -i ${cellCorrectMask} \
        -g ${nucleusLayer}/CellMask_adjusted \
        -b 2 10 50 100 150 \
        -o ${outDir}/04.register/${SN}.rpi

    ## convert cellbin GEF to GEM
    # /usr/bin/time -v singularity exec ${sif} cellCut view \
    #     -i ${outDir}/051T.cellcut/${SN}.cellbin.gef \
    #     -d ${outDir}/03.calibration/${SN}.gef \
    #     -o ${outDir}/051T.cellcut/${SN}.cellbin.gem \
    #     -s ${SN}

    echo `date` "=> cellCluster (Transcriptomics) start......."
    mkdir -p ${outDir}/tmp
    export NUMBA_CACHE_DIR=${outDir}/tmp
    export MPLCONFIGDIR=${outDir}/tmp
    /usr/bin/time -v singularity exec ${sif} cellCluster \
        -i ${outDir}/051T.cellcut/${SN}.adjusted.cellbin.gef \
        -o ${outDir}/061T.cellcluster/${SN}.adjusted.cell.cluster.h5ad

    /usr/bin/time -v singularity exec ${sif} cellCluster \
        -i ${outDir}/051T.cellcut/${SN}.cellbin.gef \
        -o ${outDir}/061T.cellcluster/${SN}.cell.cluster.h5ad

    echo `date` "=> cellChunk (Transcriptomics) start......."
    # Write rendering data into cellbin.gef
    /usr/bin/time -v singularity exec ${sif} cellChunk \
        -i ${outDir}/051T.cellcut/${SN}.adjusted.cellbin.gef \
        -o ${outDir}/051T.cellcut/

    /usr/bin/time -v singularity exec ${sif} cellChunk \
        -i ${outDir}/051T.cellcut/${SN}.cellbin.gef \
        -o ${outDir}/051T.cellcut/

fi

# Run cellCut and cellCluster for Proteomics
export SINGULARITY_BIND=$outDir
if [[ -f $imageTarFile ]] && [[ -f $iprFile ]]  && [[ $doCell == "Y" ]]; then
    echo `date` "=> cellCut (Proteomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    nucleusLayer=$(find ${outDir}/04.register -maxdepth 1 -name \*fov_stitched_transformed.tif -exec sh -c 'for f do basename -- "$f" _fov_stitched_transformed.tif;done' sh {} + | grep -v IF | awk -F_ '{print$1}')
    nucleusMask=$(find ${outDir}/04.register -maxdepth 1 -name ${nucleusLayer}_${SN}_mask.tif)
    /usr/bin/time -v singularity exec ${sif} cellCut cgef \
        -i ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef \
        -m ${nucleusMask} \
        -o ${outDir}/051P.cellcut/${SN}.protein.cellbin.gef

    echo `date` "=> cellCorrect (Proteomics) start......."
    export HDF5_USE_FILE_LOCKING=FALSE
    /usr/bin/time -v singularity exec ${sif} cellCorrect \
            -i ${outDir}/03.calibration/${SN}.protein.calibrated.raw.gef \
            -m ${nucleusMask} \
            -o ${outDir}/051P.cellcut
    mv ${outDir}/051P.cellcut/${SN}.adjusted.cellbin.gef ${outDir}/051P.cellcut/${SN}.protein.adjusted.cellbin.gef
    mv ${outDir}/051P.cellcut/${SN}.adjusted.gem ${outDir}/051P.cellcut/${SN}.protein.adjusted.gem
    gzip -f ${outDir}/051P.cellcut/${SN}.protein.adjusted.gem


    ## convert cellbin GEF to GEM
    # /usr/bin/time -v singularity exec ${sif} cellCut view \
    #     -i ${outDir}/051T.cellcut/${SN}.cellbin.gef \
    #     -d ${outDir}/02T.count/${SN}.gef \
    #     -o ${outDir}/051T.cellcut/${SN}.cellbin.gem \
    #     -s ${SN}

    echo `date` "=> cellCluster (Proteomics) start......."
    mkdir -p ${outDir}/tmp
    export NUMBA_CACHE_DIR=${outDir}/tmp
    export MPLCONFIGDIR=${outDir}/tmp
    export HDF5_USE_FILE_LOCKING=FALSE
    /usr/bin/time -v singularity exec ${sif} cellCluster-SP \
        -i ${outDir}/051P.cellcut/${SN}.protein.adjusted.cellbin.gef \
        -o ${outDir}/061P.cellcluster/${SN}.protein.adjusted.cell.cluster.h5ad

    /usr/bin/time -v singularity exec ${sif} cellCluster-SP \
        -i ${outDir}/051P.cellcut/${SN}.protein.cellbin.gef \
        -o ${outDir}/061P.cellcluster/${SN}.protein.cell.cluster.h5ad

    echo `date` "=> cellChunk (Proteomics) start......."
    # Write rendering data into cellbin.gef
    /usr/bin/time -v singularity exec ${sif} cellChunk \
        -i ${outDir}/051P.cellcut/${SN}.protein.adjusted.cellbin.gef \
        -o ${outDir}/051P.cellcut/

    /usr/bin/time -v singularity exec ${sif} cellChunk \
        -i ${outDir}/051P.cellcut/${SN}.protein.cellbin.gef \
        -o ${outDir}/051P.cellcut/
fi

# Run saturation
export SINGULARITY_BIND=$outDir
echo `date` "=> saturation (Transcriptomics) start ......"
export HDF5_USE_FILE_LOCKING=FALSE
bcStatStr=$(find ${outDir}/00T.mapping -name \*stat | tr '\n' ',' | sed 's/.$//')
/usr/bin/time -v singularity exec ${sif} saturation \
    -i ${outDir}/02T.count/${SN}_raw_barcode_gene_exp.txt \
    --tissue ${outDir}/05T.tissuecut/${SN}.tissue.gef \
    --bcstat ${bcStatStr} \
    --summary ${outDir}/02T.count/${SN}.Aligned.sortedByCoord.out.merge.q10.dedup.target.bam.summary.stat \
    -o ${outDir}/07T.saturation
## plot_200x200_saturation.png also named as ${SN}.saturation.bin200.png in some case.

export SINGULARITY_BIND=$outDir
echo `date` "=> saturation (Proteomics) start ......"
export HDF5_USE_FILE_LOCKING=FALSE
/usr/bin/time -v singularity exec ${sif} saturation \
    -i ${outDir}/00P.mapping/${SN}_cid_pid_mid_reads.tsv \
    --tissue ${outDir}/05P.tissuecut/${SN}.protein.tissue.gef \
    --bcstat ${outDir}/00P.mapping/${SN}_map.stat \
    --protein \
    -o ${outDir}/07P.saturation

## Run multiomicsAnalysis
echo `date` "=> multiomicsAnalysis start ......"
echo `date` "=> Please make sure your protein list only records actual input proteins"
export HDF5_USE_FILE_LOCKING=FALSE
export NUMBA_CACHE_DIR=${outDir}/tmp
export MPLCONFIGDIR=${outDir}/tmp
export NUMBA_DEFAULT_NUM_THREADS=10
export SINGULARITY_BIND=$outDir,${proteinListDIR}
/usr/bin/time -v singularity exec ${sif} multiomicsAnalysis \
    -r ${outDir}/05T.tissuecut/${SN}.tissue.gem.gz \
    -p ${outDir}/05P.tissuecut/${SN}.protein.tissue.gem.gz \
    -b 50 \
    -pl ${proteinList} \
    -o ${outDir}/08.multiomics


# Run report to generate HTML report
echo `date` "=> report generation start......"
export HDF5_USE_FILE_LOCKING=FALSE
export SINGULARITY_BIND=$outDir
out_iprFile=$(find ${outDir}/04.register -maxdepth 1 -name \*.ipr | head -1)

bcStatStr=$(find ${outDir}/00T.mapping -name \*.stat | tr '\n' ',' | sed 's/.$//')
bcFinalOutStr=$(find ${outDir}/00T.mapping -name \*.final.out | tr '\n' ',' | sed 's/.$//')
pipever=$(basename ${sif} .sif)

if [[ -n ${out_iprFile} ]] && [[ -e ${out_iprFile} ]] && [[ $doCell == 'Y' ]]; then
    /usr/bin/time -v singularity exec ${sif} report-PT \
        --RNAMapStat  ${bcStatStr}  \
        --ProteinMapStat ${outDir}/00P.mapping/${SN}_map.stat \
        -a ${bcFinalOutStr} \
        -g ${outDir}/02T.count/${SN}.Aligned.sortedByCoord.out.merge.q10.dedup.target.bam.summary.stat \
        --RNATissueCutStat ${outDir}/05T.tissuecut/tissuecut.stat \
        --ProteinTissueCutStat ${outDir}/05P.tissuecut/tissuecut.stat \
        --RNAVisGef ${outDir}/03.calibration/${SN}.gef \
        --ProteinVisGef ${outDir}/03.calibration/${SN}.protein.gef \
        --RNASquareClusterFile ${outDir}/06T.spatialcluster/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad \
        --ProteinSquareClusterFile ${outDir}/06P.spatialcluster/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad \
        --rpi ${outDir}/04.register/${SN}.rpi \
        --saturation ${outDir}/07T.saturation/plot_200x200_saturation.png \
        --sn ${SN} \
        --RNACellGef ${outDir}/051T.cellcut/${SN}.adjusted.cellbin.gef \
        --ProteinCellGef ${outDir}/051P.cellcut/${SN}.protein.adjusted.cellbin.gef \
        --RNACellCluster ${outDir}/061T.cellcluster/${SN}.adjusted.cell.cluster.h5ad  \
        --ProteinCellCluster ${outDir}/061P.cellcluster/${SN}.protein.adjusted.cell.cluster.h5ad \
        --iprFile ${out_iprFile} \
        --species ${refName} \
        --tissue ${tissueType} \
        --rna_tissue_fig ${outDir}/05T.tissuecut/tissue_fig \
        --protein_tissue_fig ${outDir}/05P.tissuecut/tissue_fig \
        --reference ${refName} \
        --pipelineVersion ${pipever} \
        --adt_fastq_name ${adtRead1} \
        --multimomics_spatial ${outDir}/08.multiomics/${SN}_50_spatial_leiden_totalVI_02.png \
        --multimomics_umap ${outDir}/08.multiomics/${SN}_50_UMAP_leiden_totalVI_01.png \
        --multimomics_heatmap ${outDir}/08.multiomics/${SN}_50_matrixplot_Protein_totalVI_04.png \
        --multimomics_bubble  ${outDir}/08.multiomics/${SN}_50_dotplot_RNA_totalVI_03.png \
        -o ${outDir}/09.report

elif [[ -n ${out_iprFile} ]] && [[ -e ${out_iprFile} ]] && [[ $doCell == 'N' ]]; then
    /usr/bin/time -v singularity exec ${sif} report-PT \
        --RNAMapStat  ${bcStatStr}  \
        --ProteinMapStat ${outDir}/00P.mapping/${SN}_map.stat \
        -a ${bcFinalOutStr} \
        -g ${outDir}/02T.count/${SN}.Aligned.sortedByCoord.out.merge.q10.dedup.target.bam.summary.stat \
        --RNATissueCutStat ${outDir}/05T.tissuecut/tissuecut.stat \
        --ProteinTissueCutStat ${outDir}/05P.tissuecut/tissuecut.stat \
        --RNAVisGef ${outDir}/03.calibration/${SN}.gef \
        --ProteinVisGef ${outDir}/03.calibration/${SN}.protein.gef \
        --RNASquareClusterFile ${outDir}/06T.spatialcluster/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad \
        --ProteinSquareClusterFile ${outDir}/06P.spatialcluster/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad \
        --rpi ${outDir}/04.register/${SN}.rpi \
        --saturation ${outDir}/07T.saturation/plot_200x200_saturation.png \
        --sn ${SN} \
        --iprFile ${out_iprFile} \
        --species ${refName} \
        --tissue ${tissueType} \
        --rna_tissue_fig ${outDir}/05T.tissuecut/tissue_fig \
        --protein_tissue_fig ${outDir}/05P.tissuecut/tissue_fig \
        --reference ${refName} \
        --pipelineVersion ${pipever} \
        --adt_fastq_name ${adtRead1} \
        --multimomics_spatial ${outDir}/08.multiomics/${SN}_50_spatial_leiden_totalVI_02.png \
        --multimomics_umap ${outDir}/08.multiomics/${SN}_50_UMAP_leiden_totalVI_01.png \
        --multimomics_heatmap ${outDir}/08.multiomics/${SN}_50_matrixplot_Protein_totalVI_04.png \
        --multimomics_bubble  ${outDir}/08.multiomics/${SN}_50_dotplot_RNA_totalVI_03.png \
        -o ${outDir}/09.report
else
    /usr/bin/time -v singularity exec ${sif} report-PT \
        --RNAMapStat  ${bcStatStr}  \
        --ProteinMapStat ${outDir}/00P.mapping/${SN}_map.stat \
        -a ${bcFinalOutStr} \
        -g ${outDir}/02T.count/${SN}.Aligned.sortedByCoord.out.merge.q10.dedup.target.bam.summary.stat \
        --RNATissueCutStat ${outDir}/05T.tissuecut/tissuecut.stat \
        --ProteinTissueCutStat ${outDir}/05P.tissuecut/tissuecut.stat \
        --RNAVisGef ${outDir}/03.calibration/${SN}.gef \
        --ProteinVisGef ${outDir}/03.calibration/${SN}.protein.gef \
        --RNASquareClusterFile ${outDir}/06T.spatialcluster/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad \
        --ProteinSquareClusterFile ${outDir}/06P.spatialcluster/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad \
        --saturation ${outDir}/07T.saturation/plot_200x200_saturation.png \
        --sn ${SN} \
        --species ${refName} \
        --tissue ${tissueType} \
        --rna_tissue_fig ${outDir}/05T.tissuecut/tissue_fig \
        --protein_tissue_fig ${outDir}/05P.tissuecut/tissue_fig \
        --reference ${refName} \
        --pipelineVersion ${pipever} \
        --adt_fastq_name ${adtRead1} \
        --multimomics_spatial ${outDir}/08.multiomics/${SN}_50_spatial_leiden_totalVI_02.png \
        --multimomics_umap ${outDir}/08.multiomics/${SN}_50_UMAP_leiden_totalVI_01.png \
        --multimomics_heatmap ${outDir}/08.multiomics/${SN}_50_matrixplot_Protein_totalVI_04.png \
        --multimomics_bubble  ${outDir}/08.multiomics/${SN}_50_dotplot_RNA_totalVI_03.png \
        -o ${outDir}/09.report
fi

cd ${outDir}
singularity exec ${sif} python3 << CODE
import os
import zipfile
reportdir = "./09.report/AnalysisReport"
file_news = reportdir +'.zip'
z = zipfile.ZipFile(file_news,'w',zipfile.ZIP_DEFLATED)
for dirpath, dirnames, filenames in os.walk(reportdir):
	fpath = dirpath.replace(reportdir,'')
	fpath = fpath and fpath + os.sep or ''
	for filename in filenames:
	  z.write(os.path.join(dirpath, filename),fpath+filename)
	print('zip success')
z.close()
CODE
mv ${outDir}/09.report/AnalysisReport.zip ${outDir}/09.report/${SN}.AnalysisReport.zip

# Organize files for visualization required by StereoMap
## 03.calibration
if [[ -f ${outDir}/03.calibration/${SN}.gef ]] || [[ -f ${outDir}/03.calibration/${SN}.protein.gef ]]
then
    ln -s ${outDir}/03.calibration/${SN}.gef ${outDir}/visualization/${SN}.gef
    ln -s ${outDir}/03.calibration/${SN}.protein.gef ${outDir}/visualization/${SN}.protein.gef
fi

## 04.register
if [[ -f ${outDir}/04.register/fov_stitched_transformed.rpi ]] || [[ -f ${outDir}/04.register/${SN}.rpi ]] || [[ -f out_iprFile ]]
then
    ln -s ${outDir}/04.register/fov_stitched_transformed.rpi ${outDir}/visualization/fov_stitched_transformed.rpi
    ln -s ${outDir}/04.register/${SN}.rpi ${outDir}/visualization/${SN}.rpi
    ln -s $out_iprFile ${outDir}/visualization/$(basename "$out_iprFile")
fi

## 051T.cellcut & 051P.cellcut
if [[ -f ${outDir}/051T.cellcut/${SN}.adjusted.cellbin.gef ]] || [[ -f ${outDir}/051P.cellcut/${SN}.protein.adjusted.cellbin.gef ]]
then
    ln -s ${outDir}/051T.cellcut/${SN}.adjusted.cellbin.gef ${outDir}/visualization/${SN}.adjusted.cellbin.gef
    ln -s ${outDir}/051P.cellcut/${SN}.protein.adjusted.cellbin.gef ${outDir}/visualization/${SN}.protein.adjusted.cellbin.gef
fi

## 06T.spatialcluster && 06P.spatialcluster
if [[ -f ${outDir}/06T.spatialcluster/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad ]] || [[ -f ${outDir}/06P.spatialcluster/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad ]]
then
    ln -s ${outDir}/06T.spatialcluster/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad ${outDir}/visualization/${SN}.bin${binSize}_${resolution}.spatial.cluster.h5ad
    ln -s ${outDir}/06P.spatialcluster/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad ${outDir}/visualization/${SN}.protein.bin${binSize}_0.1.spatial.cluster.h5ad
fi
echo `date` " All done! "