Question: [R] Coverage Per Gene
gravatar for Jetse
6.6 years ago by
Jetse0 wrote:

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])

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?

R bam coverage ucsc • 2.4k views
ADD COMMENTlink modified 4.7 years ago by Biostar ♦♦ 20 • written 6.6 years ago by Jetse0
gravatar for Sean Davis
6.6 years ago by
Sean Davis25k
National Institutes of Health, Bethesda, MD
Sean Davis25k wrote:

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

ADD COMMENTlink written 6.6 years ago by Sean Davis25k
gravatar for Jetse
6.6 years ago by
Jetse0 wrote:

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

topDataframe <-[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 =,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.

ADD COMMENTlink modified 6.6 years ago • written 6.6 years ago by Jetse0
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1845 users visited in the last hour