BED files

Question

[adarsh_pp]

Hi,

I am trying to compare, bed files of multiple human exome (hg38) capture kits using Bedtools Intersect to get differences in terms of total number of genes, omim genes, cds, exons and total number of bases.

For this I need to get an updated reference bed file of each parameter for comparison.

Please let me know how can I get these files.

Thank you

[barslmn]

I am guessing you want the annotation files. https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/

[Ram]

You deleted your previous post where Pierre had added a comment asking you to validate all your previous posts. Did you go back and validate them?

[Ram]

Also, how is this question different from the one you asked 8 weeks ago here: Calculate the percentage of genomic region covered from the BED file.

[adarsh_pp]

Sorry for the inconvenience happened here. Both of the questions are different and require different approaches. That's why I posted it as a second question.

[Ram]

Thank you for the clarification. Please do go back and validate your older posts.

Answer 1

[Alex Reynolds]

Generically via bedmap:

bedmap --echo --echo-map --count exome.bed annotations.bed > answer.bed

The annotations.bed file can be GENCODE genes, OMIM genes, etc. which have CDSs, exons, etc.

For example, GENCODE gene counts can be obtained via convert2bed:

wget -qO- https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_44/gencode.v44.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

Then:

bedmap --echo --echo-map-id-uniq --count exome.bed genes.bed > answer.bed

The last column of answer.bed will have the number of genes that overlap an entry in exome.bed.

You can repeat this for other sets and modify awk statements etc. to get specific subsets of annotations.

OMIM annotations may need an intermediate map step to go from MIM numbers to ENSG to bedmap back to BED entries.