[R] Coverage Per Gene

Question

[Jetse]

Hello everyone,

To calculate the coverage of a bam file, I used A: Tools To Calculate Average Coverage For A Bam File? answer. I have an GRangesList object with all genes and their ranges:

txdb <- makeTranscriptDbFromUCSC( genome='hg19', tablename='ensGene' )
tx_by_gene <- transcriptsBy( txdb, 'gene')

Another object I have is an GRanges object with my bam file:

alignment <- readBamGappedAlignments( fileName )

I also have a vector with highly differential expressed genes.

topGenes <- c("ENSG00000258724","ENSG00000259141",...)

Now I am doing this:

subset <- subsetByOverlaps(alignment,tx_by_gene[topGenes])
mean(coverage(subset))

But then I get the coverage per chromosome with only using the topGenes, and I want the coverage per gene. How do I do this?

Answer 1

[Sean Davis]

You might want to look at summarizeOverlaps() instead of subsetByOverlaps.

Answer 2

[Jetse]

I found another solution, with summarizeOverlaps() I get per chromosome too. Now I write to a bed file:

topDataframe <- as.data.frame(tx_by_gene[topGenes])


topBed <- data.frame(chrom = gsub( "chr([0-9(MT|X|Y)])" , "\\1" ,topOnlyDataframe$seqnames),
                        chromStart = topOnlyDataframe$start,
                        chromEnd = topOnlyDataframe$end,
                        name = topOnlyDataframe$element,
                        score = rep.int(0,nrow(topOnlyDataframe)),
                        strand = topOnlyDataframe$strand)

Then write this to a bed file:

write.table(topBed, file="/data/jetse/wntSignalling/topGenesOnly.BED",sep="\t", row.names = FALSE, col.names = FALSE, quote = FALSE)

Then use bedtools to get the coverage.