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?