This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Where to get homopolymer annotation for human genome ?

I got lots of SNV and Indel somatic mutations from whole genome sequencing data from tumor. And I want to annotate these somatic mutations with whether located in the homopolymer region of human genome. I have already search the ucsc repeatmasker track, and it does not contain homopolymer annotation. I want to know where to download this annotation or how can I calculate the homopolymer of human genome? Thanks !

snp homopolymer annotation
• 4,307 views • link
updated by bioinfo89 • bioinfo89 • written by lzy • lzy

Hi @ATpoint,

I would really appreciate if you could share the Rscript again. I couldn't find your Rscript.

Thanks!

• 0 views • link
written by bioinfo89 • bioinfo89

Sure, I just updated the link in the answer.

• 0 views • link
written by ATpoint ATpoint

2 answers

if your input is a VCF file, you can use GATK VariantAnnotator with -A HomopolymerRun

https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_annotator_HomopolymerRun.php

• 1 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thank you very much, I will try this method.

• 1 views • link
written by lzy • lzy

@Pierre Lindenbaum, -A HomopolymerRun isn't mentioned in the documentation on the GATK site or when you run "gatk VariantAnnotator --help" on the command line. Do you know what happened to this option or of another method to find all homopolymer regions in the human genome?

• 0 views • link
written by Charles-Alexandre Roy • Charles-Alexandre Roy

I wrote an R script for this. It finds homopolymer stretches with a minimum-length defined by the user from any BSgenome. Just copy the whole thing into R. An example on how to use it is at the bottom of the script. Just add the BSgenome you want and get a GRanges with all the polyX coordinates.

• 0 views • link
written by ATpoint ATpoint

Thank you very much, I have tried your R script for homopolymer finding, and got result as follow:

> PolyAll6_hg19 
GRanges object with 7760901 ranges and 0 metadata columns:
                  seqnames         ranges strand
                     <Rle>      <IRanges>  <Rle>
        [1]           chr1 [11541, 11546]      *
        [2]           chr1 [16506, 16511]      *
        [3]           chr1 [19736, 19741]      *
        [4]           chr1 [23703, 23709]      *
        [5]           chr1 [25273, 25278]      *
        ...            ...            ...    ...
  [7760897] chrUn_gl000249 [17529, 17534]      *
  [7760898] chrUn_gl000249 [22129, 22134]      *
  [7760899] chrUn_gl000249 [32852, 32857]      *
  [7760900] chrUn_gl000249 [35780, 35785]      *
  [7760901] chrUn_gl000249 [35937, 35942]      *
  -------
  seqinfo: 92 sequences from an unspecified genome; no seqlengths

But I don't know how to generate the bed file format for this result. Can you show me the example of bed file generating.

• 1 views • link
written by lzy • lzy

Sure, it is as easy as:

write.table(
    data.frame(seqnames(PolyAll6_hg19), start(PolyAll6_hg19)-1, end(PolyAll6_hg19),
    sep="\t", quote = F, col.names = F, row.names = F, file = "/path/to/output_file.bed"
)
• 1 views • link
written by ATpoint ATpoint

Hi @ATpoint,

Thank you very much for the Rscript, I used it for extracting the homopolymers from reference genome. I wanted to know if there is a way to use the same script to limit the homopolymers finding in the coding regions only?

• 0 views • link
written by bioinfo89 • bioinfo89

I would simply take the full output, then get the coding regions in GRanges format and use something like subsetByOverlaps() to get the homopolymers overlapping with it. Biomart can probably do the query for coding exons conveniently. Search its manual :)

• 0 views • link
written by ATpoint ATpoint

Log in to answer this question.

More posts like this