#!/usr/bin/env bash
set -euo pipefail
mkdir -p /tmp/rustar-mre-32 && cd /tmp/rustar-mre-32

BASE=https://raw.githubusercontent.com/nf-core/test-datasets/626c8fab639062eade4b10747e919341cbf9b41a
curl -fsLO $BASE/reference/genome.fasta
curl -fsL  $BASE/reference/genes_with_empty_tid.gtf.gz | gunzip -c > genes.gtf
curl -fsLO $BASE/testdata/GSE110004/SRR6357072_1.fastq.gz
curl -fsLO $BASE/testdata/GSE110004/SRR6357072_2.fastq.gz

RUSTAR=ghcr.io/scverse/rustar-aligner:dev
STAR=community.wave.seqera.io/library/htslib_samtools_star_gawk:ae438e9a604351a4

mkdir -p idx-rustar idx-star
docker run --rm -v $PWD:/w -w /w $RUSTAR rustar-aligner --runMode genomeGenerate \
    --genomeDir idx-rustar --genomeFastaFiles genome.fasta --sjdbGTFfile genes.gtf \
    --sjdbOverhang 100 --genomeSAindexNbases 7
docker run --rm -v $PWD:/w -w /w $STAR STAR --runMode genomeGenerate \
    --genomeDir idx-star --genomeFastaFiles genome.fasta --sjdbGTFfile genes.gtf \
    --sjdbOverhang 100 --genomeSAindexNbases 7

COMMON=(--readFilesIn SRR6357072_1.fastq.gz SRR6357072_2.fastq.gz --readFilesCommand zcat
        --runThreadN 4 --sjdbGTFfile genes.gtf --twopassMode Basic --runRNGseed 0
        --outSAMtype BAM Unsorted --outSAMattributes NH HI AS NM MD
        --quantMode TranscriptomeSAM --quantTranscriptomeSAMoutput BanSingleEnd
        --outSAMattrRGline ID:WT_REP2 SM:WT_REP2)

docker run --rm -v $PWD:/w -w /w $RUSTAR rustar-aligner \
    --genomeDir idx-rustar "${COMMON[@]}" --outFileNamePrefix RUS.
docker run --rm -v $PWD:/w -w /w $STAR STAR \
    --genomeDir idx-star "${COMMON[@]}" --outFileNamePrefix STAR.

for label in rustar STAR; do
    if [ $label = rustar ]; then G=RUS./Aligned.out.bam; T=RUS./Aligned.toTranscriptome.out.bam
    else                         G=STAR.Aligned.out.bam; T=STAR.Aligned.toTranscriptome.out.bam; fi
    g_hdr=$(docker run --rm -v $PWD:/w $STAR samtools view -H /w/$G | grep -c '^@RG')
    g_rec=$(docker run --rm -v $PWD:/w $STAR samtools view /w/$G | grep -c 'RG:Z:')
    t_hdr=$(docker run --rm -v $PWD:/w $STAR samtools view -H /w/$T | grep -c '^@RG')
    t_rec=$(docker run --rm -v $PWD:/w $STAR samtools view /w/$T | grep -c 'RG:Z:')
    t_total=$(docker run --rm -v $PWD:/w $STAR samtools view -c /w/$T)
    printf "%-6s genome: @RG=%d, RG:Z:=%d | transcriptome: @RG=%d, RG:Z:=%d / %d\n" \
        $label $g_hdr $g_rec $t_hdr $t_rec $t_total
done