Rsubread: .bai file lacks data for `samtools idxstats`
Entering edit mode
Last seen 4 months ago

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
Rsubread • 628 views
Entering edit mode
Last seen 11 hours ago
WEHI, Melbourne, Australia

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

Entering edit mode

Thanks, Gordon, for pointing out the align.res file. However, I would be really interested in the per-sequence statistics.


Login before adding your answer.

Traffic: 436 users visited in the last hour
Help About
Access RSS

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6