How to calculate the coverage of list of genes in whole exome data

Question

[anii]

Hi all,

I have bam file of whole exome human data. Now I want to check coverage of a list of genes in my data. Basically I want get output like this:

Gene    Percentage of coding region covered
A1BG               100 %
A1BG-AS1           100%
A1CF                99.94%
NEK8                100%
FANCI              100%
A2ML1             100%

[GenoMax]

If this data is using a specific kit then the manufacturer may have the BED files available for the regions that are covered by the kit. You should use that file.

[anii]

Yes, they used the Agilent kit and i have downloaded the bed file from the agilent site. should i use this bed file or should i get the file from the manufacture?

[GenoMax]

Use the agilent file. Make sure it is for the genome build you are using otherwise things would be totally off.

[anii]

Okay , Thakyou GenoMax

Answer 1

[Shred]

Create a BED file of your genes of interests and then use mosdepth as described https://github.com/brentp/mosdepth#exome-example

[anii]

shred can you please elaborate the steps to create the BED file of genes of interests.

[Trivas]

use the UCSC table browser. You can upload a list of identifiers (refseq IDs, I believe) and it will spit out a bed file of the corresponding coordinates.

[Shred]

A BED file is a tab separated file which defines a region. It has 3 columns, chromosome, start, end.

You could generate using the method pointed by @Trivas or using BioMart https://www.ensembl.org/biomart/martview

https://en.wikipedia.org/wiki/BED_(file_format)

[anii]

I have get this output by using mosdepth. It give me coverage at 1X, 5X,10X,15X ,20X , 30X. In some regions genes are repeated as you can see in the example. Now how Can i convert the coverage in the form of percentage?

example:

#chrom  start   end region  id  1X  5X  10X 15X 20X 30X

chr1    69090   70008   OR4F5   918 918 918 918 918 918
chr1    134772  140566  LOC729737   1985    0   0   0   0   0
chr1    182387  184878  DDX11L17    2464    1214    860 550 412 273
chr1    187890  187958  MIR6859-1   68  68  68  68  68  68
chr1    1020119 1056114 AGRN    22418   15841   14836   14026   13286   12229
chr1    1020119 1056114 AGRN    22418   15841   14836   14026   13286   12229
chr1    1033992 1056114 AGRN    17233   14700   13867   13225   12695   11715

[Shred]

Never post screenshots here. Tabular data could be pasted using the code format. Which is the BED file for these regions?

[anii]

I have used the bed file(hg38) from ucsc

[anii]

I have get this output by using mosdepth. It give me coverage at 1X, 5X,10X,15X ,20X , 30X. In some regions genes are repeated as you can see in the example. Now how Can i convert the coverage in the form of percentage?

example:

chrom          start         end             region              1X 5X 10X 15X 20X 30X
chr1             69090    70008        OR4F5               918 918 918 918 918 918 
chr1          134772    140566       LOC729737        1985 0 0 0 0 0 
chr1       182387       184878        DDX11L17         2464 1214 860 550 412 273 
chr1       187890         187958        MIR6859-1       68 68 68 68 68 68 
chr1         1020119   1056114         AGRN              22418    15841 14836 14026 13286 12229 
chr1       1020119     1056114        AGRN              22418    15841 14836 14026 13286 12229 
chr1      1033992      1056114        AGRN                17233 14700 13867 13225 12695 11715

[Shred]

This may be useful https://github.com/brentp/mosdepth/issues/166#issuecomment-993250933

Answer 2

[amy__]

https://genviz.org/module-03-genvisr/0003/04/01/gencov_GenVisR/ Is one way