Most reads counted as NOT_NH_1 without apparent reason?

[pedrolmoura]

Hi! I've used rMATS-turbo for a while, always on my own STAR aligned files, and it has worked very well throughout my use cases.

However, I had to access some pre-aligned literature .bam files, and they are unreadable to rMATS for some reason I have not yet been able to identify. These are paired-end .bam files which are fully readable by featureCounts and IGV. Here is an example of the rMATS bam outcome log:

Example.bam
USED: 0
NOT_PAIRED: 11232424
NOT_NH_1: 174044596
NOT_EXPECTED_CIGAR: 0
NOT_EXPECTED_READ_LENGTH: 0
NOT_EXPECTED_STRAND: 0
EXON_NOT_MATCHED_TO_ANNOTATION: 0
JUNCTION_NOT_MATCHED_TO_ANNOTATION: 0
CLIPPED: 0
TOTAL_FOR_BAM: 185277020

For this particular .bam file, featureCounts returns 70 million successfully assigned reads, and according to samtools only ~1 million are flagged as multimapped. I looked into other issues, but all I found related to NOT_NH_1 being a potential cause for rMATS crashes did not help, since this particular case doesn't seem to be a case of annotation or read length mismatches. My current options are --allow-clipping, --novelSS and --variable-read-lengths.

Thank you in advance and hope this is a straightforward issue.

[EricKutschera]

Since samtools indicates that not many reads are multimapped but rMATS says that most reads are NOT_NH_1, my guess is that the NH tag is just not set for most reads

You could try checking the bam headers to see how the reads were aligned. That might explain why the alignments don't have the NH tag set as expected:

samtools view -h /path/to/file.bam | grep '@PG'

To let rMATS use the alignments you could filter out secondary alignments and then set the NH tag to 1. This samtools command filters secondary alignments and removes any NH tag. Then the sed command appends NH:i:1 to any line that doesn't start with @ (comment lines)

samtools view -h -F SECONDARY --remove-tag NH /path/to/file.bam | sed '/^[^@]/s/\(.*\)/\1\tNH:i:1/' > /path/to/updated.sam

Finally convert the updated sam files to bam files:

samtools view -o /path/to/updated.bam /path/to/updated.sam

Here's a related post:
#298 (comment)