[BioC] Repeat masker sequences as GRanges object
Michael Lawrence
lawrence.michael at gene.com
Fri Sep 12 22:47:47 CEST 2014
On Fri, Sep 12, 2014 at 11:08 AM, Hervé Pagès <hpages at fhcrc.org> wrote:
> Hi,
>
>
> On 09/12/2014 06:30 AM, James W. MacDonald wrote:
>
>> Hi Hermann,
>>
>> How about this:
>>
>> library(AnnotationHub)
>>> hub <- AnnotationHub()
>>> hub$goldenpath.hg19.database.rmsk_0.0.1.RData
>>>
>> GRanges with 5298130 ranges and 2 metadata columns:
>> seqnames ranges strand |
>> name
>> <Rle> <IRanges> <Rle> |
>> <character>
>> [1] chr1 [16777161, 16777470] + |
>> AluSp
>> [2] chr1 [25165801, 25166089] - |
>> AluY
>> [3] chr1 [33553607, 33554646] + |
>> L2b
>> [4] chr1 [50330064, 50332153] + |
>> L1PA10
>> [5] chr1 [58720068, 58720973] - |
>> L1PA2
>> ... ... ... ... ...
>> ...
>> [5298126] chr21_gl000210_random [25379, 25875] + |
>> MER74B
>> [5298127] chr21_gl000210_random [26438, 26596] - |
>> MIRc
>> [5298128] chr21_gl000210_random [26882, 27022] - |
>> MIRc
>> [5298129] chr21_gl000210_random [27297, 27447] + |
>> HAL1-2a_MD
>> [5298130] chr21_gl000210_random [27469, 27682] + |
>> HAL1-2a_MD
>> score
>> <numeric>
>> [1] 2147
>> [2] 2626
>> [3] 626
>> [4] 12545
>> [5] 8050
>> ... ...
>> [5298126] 1674
>> [5298127] 308
>> [5298128] 475
>> [5298129] 371
>> [5298130] 370
>> ---
>> seqlengths:
>> chr1 chr2 ... chr18_gl000207_random
>> 249250621 243199373 ... 4262
>>
>> This is a GRanges of all features from UCSC's Repeat Masker table.
>>
>
> Nice to see the "name" and "score" metadata cols but I wonder why the
> original UCSC names for these cols (which are "repName" and "swScore")
> were not preserved. Also, other UCSC cols in the rmsk table at UCSC
> might be of interest (e.g. "repClass" and "repFamily").
>
> FWIW, here is one way to get these cols:
>
> local_file <- tempfile()
>
> download.file("http://hgdownload.soe.ucsc.edu/
> goldenPath/hg19/database/rmsk.txt.gz", local_file)
>
> ## Get the col names from
> ## http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/rmsk.sql
> COLNAMES <- c("bin", "swScore", "milliDiv", "milliDel", "milliIns",
> "genoName", "genoStart", "genoEnd", "genoLeft",
> "strand", "repName", "repClass", "repFamily",
> "repStart", "repEnd", "repLeft", "id")
>
> library(GenomicRanges)
> df <- read.table(local_file, col.names=COLNAMES)
> rmsk <- makeGRangesFromDataFrame(df, keep.extra.columns=TRUE,
> seqnames.field="genoName",
> start.field="genoStart",
> end.field="genoEnd",
> strand.field="strand",
> starts.in.df.are.0based=TRUE)
>
> Then:
>
> > head(rmsk)
> GRanges with 6 ranges and 13 metadata columns:
> seqnames ranges strand | bin swScore milliDiv
> milliDel
> <Rle> <IRanges> <Rle> | <integer> <integer> <integer>
> <integer>
> [1] chr1 [10001, 10468] + | 585 1504 13
> 4
> [2] chr1 [10469, 11447] - | 585 3612 114
> 270
> [3] chr1 [11504, 11675] - | 585 437 235
> 186
> [4] chr1 [11678, 11780] - | 585 239 294
> 19
> [5] chr1 [15265, 15355] - | 585 318 230
> 38
> [6] chr1 [16713, 16749] + | 585 203 162
> 0
> milliIns genoLeft repName repClass repFamily repStart
> <integer> <integer> <factor> <factor> <factor>
> <integer>
> [1] 13 -249240153 (CCCTAA)n Simple_repeat Simple_repeat 1
> [2] 13 -249239174 TAR1 Satellite telo -399
> [3] 35 -249238946 L1MC LINE L1 -2236
> [4] 10 -249238841 MER5B DNA hAT-Charlie -74
> [5] 0 -249235266 MIR3 SINE MIR -119
> [6] 0 -249233872 (TGG)n Simple_repeat Simple_repeat 1
> repEnd repLeft id
> <integer> <integer> <integer>
> [1] 463 0 1
> [2] 1712 483 2
> [3] 5646 5449 3
> [4] 104 1 4
> [5] 143 49 5
> [6] 37 0 6
>
> I also tried with rtracklayer but got the following error:
>
> > library(rtracklayer)
> > session <- browserSession()
> > genome(session) <- "hg19"
> > query <- ucscTableQuery(session, "RepeatMasker", table="rmsk")
> > system.time(rmsk <- getTable(query))
> Error in scan(file, what, nmax, sep, dec, quote, skip, nlines,
> na.strings, :
> line 3813855 did not have 17 elements
> Timing stopped at: 551.027 7.24 1304.386
>
> Could be due to my flaky internet connection though...
>
>
This is just hitting up against the row count limit of the table browser.
rtracklayer uses the table browser as an abstraction, but it has its
limitations.
Anyway, it would be great if the AnnotationHub track could be improved.
Those extra fields are indeed useful.
> Cheers,
> H.
>
> > sessionInfo()
> R version 3.1.1 (2014-07-10)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
> [9] LC_ADDRESS=C LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] parallel stats graphics grDevices utils datasets methods
> [8] base
>
> other attached packages:
> [1] rtracklayer_1.25.16 GenomicRanges_1.17.40 GenomeInfoDb_1.1.19
> [4] IRanges_1.99.28 S4Vectors_0.2.3 BiocGenerics_0.11.5
>
> loaded via a namespace (and not attached):
> [1] BatchJobs_1.3 BBmisc_1.7 BiocParallel_0.99.19
> [4] Biostrings_2.33.14 bitops_1.0-6 brew_1.0-6
> [7] checkmate_1.4 codetools_0.2-9 DBI_0.3.0
> [10] digest_0.6.4 fail_1.2 foreach_1.4.2
> [13] GenomicAlignments_1.1.29 iterators_1.0.7 RCurl_1.95-4.3
> [16] Rsamtools_1.17.33 RSQLite_0.11.4 sendmailR_1.1-2
> [19] stats4_3.1.1 stringr_0.6.2 tools_3.1.1
> [22] XML_3.98-1.1 XVector_0.5.8 zlibbioc_1.11.1
>
>
>> Best,
>>
>> Jim
>>
>>
>>
>>
>> On Thu, Sep 11, 2014 at 3:16 AM, Hermann Norpois <hnorpois at gmail.com>
>> wrote:
>>
>> Hello,
>>>
>>> I would like to have repeat sequences as GRanges object
>>> I started with ...
>>>
>>> library (BSgenome.Hsapiens.UCSC.hg19)
>>> ch1 <- Hsapiens$chr1
>>> active (masks (ch1))
>>> AGAPS AMB RM TRF
>>> TRUE TRUE FALSE FALSE
>>> active (masks(ch1))["RM"] <- TRUE
>>> active (masks (ch1))
>>> AGAPS AMB RM TRF
>>> TRUE TRUE TRUE FALSE
>>>
>>> Can anyboldy give me a hint how to continue.
>>>
>>> Thanks
>>> hermann
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> Bioconductor mailing list
>>> Bioconductor at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>>> Search the archives:
>>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>>
>>>
>>
>>
>>
> --
> Hervé Pagès
>
> Program in Computational Biology
> Division of Public Health Sciences
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N, M1-B514
> P.O. Box 19024
> Seattle, WA 98109-1024
>
> E-mail: hpages at fhcrc.org
> Phone: (206) 667-5791
> Fax: (206) 667-1319
>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives: http://news.gmane.org/gmane.
> science.biology.informatics.conductor
>
[[alternative HTML version deleted]]
More information about the Bioconductor
mailing list