Hi,

I have rna seq data that I have mapped with STAR which allows soft clipping of reads. I have imported the bam file with the following code:

param <- ScanBamParam(flag=scanBamFlag(isSecondaryAlignment = F), tag=c('NH'), what='seq')
reads <- readGAlignments('mapped.bam', use.names=T, param = param)

I know that the cigar string shows how many nucleotides have been soft clipped from either the 5' and 3' end of a read with the #S#M or #M#S. I'd like to extract just the nucleotides that have been clipped from the 5' and the 3' end for each read. I saw that there are some functions to operate on cigar strings and genomicAlignments, but didn't see any that looked like they would easily do this. Is there an easy way to extract this information? Thanks

Hi Jake,

Yep the cigar utilities in GenomicAlignments are pretty low-level and not particularly easy to use but they do the job:

library(GenomicAlignments)
cigar <- c("3S15M3I24S5H", "25M300N15M", "30M10S")
cigarRangesAlongQuerySpace(cigar, ops="S")
# IRangesList of length 3
# [[1]]
# IRanges of length 2
#     start end width
# [1]     1   3     3
# [2]    22  45    24
# 
# [[2]]
# IRanges of length 0
# 
# [[3]]
# IRanges of length 1
#     start end width
# [1]    31  40    10

Then use this IRangesList object to extract the soft-clipped sequences from the query sequences with extractAt() from the Biostrings package.

Cheers,

H.