I am aligning RNA-seq data to a custom reference using subjunc():

align.res <- subjunc(index="Repvec", file1=fwdname, readfile2=revname, output_file=bamname, nthreads = 6, sortReadsByCoordinates=TRUE)

With the last parameter, the BAM file is sorted and the .bai file is created.

Everthing is fine, when I load the BAM file into an genome viewer; however, when I want to do a quick check of the alignment with samtools idxstat, the 'mapped read-segments' and 'unmapped read-segments' columns are empty:

samtools idxstats sample1.bam
# mCherry_1495    1495    0       0
# GFP__1468       1468    0       0
# BFP2_1494       1494    0       0
# AMPR_2739       2739    0       0
# *               0       0       0

After recreating the index from the BAM file with samtools index, the columns are populated:

samtools index sample1.bam
samtools idxstats sample1.bam
# mCherry_1495    1495    40412   0
# GFP__1468       1468    525     0
# BFP2_1494       1494    23843   0
# AMPR_2739       2739    3       0
# *               0       0       40681972

It would be helpful, if subjunc() (and align()) would add also the information required by samtools idxstats to the .bai file.

The environment is as follows:

R version 4.3.0 (2023-04-21 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 11 x64 (build 22621)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] Rsubread_2.14.2

loaded via a namespace (and not attached):
[1] compiler_4.3.0 Matrix_1.5-4.1 tools_4.3.0    grid_4.3.0     lattice_0.21-8

subjunc stores alignment statistics in align.res, although admitedly that gives overall statistics for the whole file rather than per sequence.