hi,

i have a 4Gb BAM file with RNA-seq reads aligned with STAR to the hg38 version of the human genome and where, according to STAR, an important fraction of them (~25%) aligned to multiple loci. I"m interested in finding the genes that overlap these multimapping reads to have an idea of the origin of these reads. Could anyone suggest me a Rsamtools/GenomicAlignments/GenomicFeatures route to extract these multimapped alignments and genes overlapping them?

thanks!!

robert.

It's easy. For example,

library(GenomicAlignments)
ambiguousReads <- readGAlignmentPairs("/users/robert/project/example.bam", param = ScanBamParam(flag = scanBamFlag(isSecondaryAlignment = TRUE))) # Only import multi-mapping reads.

library(rtracklayer)
genesExonsTranscripts <- import("/users/robert/databases/hg38/GENCODEversion26.gtf") # GRanges with metadata.
genes <- subset(genesExonsTranscripts, type == "gene")

library(GenomicRanges)
genesAlignments <- countOverlaps(genes, ambiguousReads)
mcols(GENCODEgenes)[genesAlignments > 0, "gene_name"] # Gene symbols of genes with non-zero counts.

You may like to assign such alignments to a particular gene using RSEM.