Convert Paired End Reads BAM file to GRanges

[ayushraman]

Hi,

I would like to know ways to convert the Paired end bam file into GRanges (R Data Structure). I have readBamGappedAlignments function but I don't know if it works for paired end reads or not ? Any help/advice/hints would be appreciated..

thanks

[jurassicMarc]

Hi ayushraman

I've used the readGappedAlignments function in the Rsamtools package to plot gene coverage across genes. I had paired end Illumina data, and I used the function to read the data into a GRanges object as follows:

Firstly, I created a GRanges object of genome position (with coordinated extracted from the CCDS database)

library(ggplot2)
library(ggbio) 
library(IRanges)
library(GenomicRanges)
library(Rsamtools)

gene_coords <- GRanges (seqnames = chr, ranges = IRanges(start = cds_start, end = cds_end))

I then used the readGappedAlignments function to read in paired-end data from the BAM file

bam_fields <- c("pos", "qwidth") #I believe these are the minimal parameters required
				
param <- ScanBamParam(which = gene_coords, what = bam_fields)
				
bam_file <- "path/to/paired-end-data-file.bam"
				
bam <- readGappedAlignments(bam_file, param = param)

I converted the input data to a GRanges object

gr <- granges(bam)

Once I had that I could plot the coverage (Providing the gene had any coverage) using the ggbio's autoplot function:

#Plot a gene if the gene region has any coverage
gene_plot <- if(length(width(gr)) >= 1){
        #exons contains the start and stop regions for each exon in the gene
        exons <- data.frame(xmin = exon_starts, xmax = exon_ends, ymin = 0, ymax = Inf)
					
	cov_plot <- autoplot (
		gr, 
		stat = "coverage",
		geom = "area",
		colour = "black",
		fill = "green",
		xlab = "Position",
		ylab = "Coverage",
		main = "Coverage across gene"
	
        #Marks the exon regions out so that you can see how much of the sequencing is on target				
	) + geom_rect(
		data = exons, 
		aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax), 
		fill = "yellow", 
		colour = "black", alpha = 0.5
					
	)
}

#plot the graph				
gene_plot

This method for reading in data is relatively quick, not that memory intensive and quite simple.

Hope this helps!