Re: [Bioc-sig-seq] making target from fasta file

Wed, 04 Jun 2008 14:02:01 -0700 

Hi Herve
The read.DNAStringSet works now, I think the problem was my downloaded fasta file. 
Here is the result and the session Info:


Fly_RM
[1] "/Users/jdhahbi/JosephDhahbi/SOLEXA/Fly/R/export.files/s-2/*.RData/RepeatMasker.fa" 
fly_rm=read.DNAStringSet(Fly_RM, format="fasta")

Read 1126586 items

fly_rm

  A DNAStringSet instance of length 154993
width seq names [1] 433 AATTCGCGTCCGCTTACCCATGTGCCTGTGGATGCCGAACAGGAGGCGCCGTTGACGGCGAATG...GCTACCCGTCTCTAAGCTTGCAGTTTTGGATTTAAGTGAATCGGTTATTCACGGGGTCGGGGA dm3_rmsk_NINJA_I ... [2] 177 TGTCGCGGATCGAACGGTGCAATCGATAGGCGTAATCAGTATTTCCAGATAGTGATAAGATTTG...AAGCTTAGCGTGCATTGTCGATCGAGAGTTTGGAGGGCAAACTGCGGTAAGATAAGATTAAAT dm3_rmsk_NINJA_LT... [3] 1086 ATACGATGGCTGTACCTCATGGTAGTCAGACCATACTAACCTATGAGGCAATAGCCTGGGGTGA...ACTGCATGGAGGAGGAGAGCTCTGCACATATTATATGTGACTGCATGGCGCTTTCCATCAGGA dm3_rmsk_Baggins1... [4] 62 TCTGCATAGTCTCCGTGTCCATCAACGCCAACGAACGGTTAAGGTGCAGCATAAACACTTCT dm3_rmsk_DMLTR5 r... [5] 346 AAATAATCCTGGAGGAGGCAAGCCAGCCCTCACAGCGACATTTTATTGTTTTTGGAAACAAAAA...GCACAACCGCGCCCATAGGGCAGCCCAGGAAAATGTTAAGCAAGTCCTTCAGGACTCTATACT dm3_rmsk_Gypsy_I ... [6] 243 CCTGGCTCAGTTATTCCGAGAGCACTTCCCACGGAGAGCAGAGAAGGAGATGAGCTGACGGCAG...CGAGAGTGTGAGAGGTGCGGAGTCCTTTGCCGAACGCCAATCCGAGTATAGAGTCCTGAAGGA dm3_rmsk_DMRT1C r... [7] 149 CGGATGCTATAGGCCGCACCCAGTTTCTCAAAGCTGGCGACTTCAACGCATGGTCAAAAAGCTG...AGAGGCAAGATGGTGCTGGAGGCATTCGCGACGTTGGACCTGGCTCAGTTATTCCGAGAGCAC dm3_rmsk_DMRT1C r... [8] 230 AACCAGAATCACTGCACGGCTGTCCAAGACCTGCCCAAGACGGTGCGCGAACACAGAGTGGAGC...ACCGATGTTTTTTCGGACATTGGGTTTGTTAGGGCATGGGTGGGCAGATGGTGGGTGTACAGC dm3_rmsk_DMRT1C r... [9] 210 CAGCAGCAGTGGTGGCTATCAAGAGCATCCGTCAAACGCCGTATGGAGGGAAGTCTGCTATAAT...TGGAGCCACGCCAAAGATGCTTCAAATGTCTGGAGGAAGGCCACATAGCGGGCCACTGTAGAA dm3_rmsk_DMRT1C r... 
     ...   ... ...
[154985] 164 GAGATGAGATGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGA...AGAAGAGAAGAGAAGAGAAGAGAAGAGAAGATAAGAGAAGAGAATAGAATAGAAGAGAATAGA dm3_rmsk_(GAGAA)n... [154986] 169 AGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAA...GAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGACTAGACGAGACGATAAGAGAATAGATGAGAA dm3_rmsk_(GAGAA)n... [154987] 157 CAACACAACACAACACAACACAACACAACACAACACAACACAACACAACACAACACAACACAAC...ACAACACAACACAACACAACACAACACAACACAACACAACACAACACAACATTTCACAACACA dm3_rmsk_(CACAA)n... [154988] 153 AGAGAGCAGAGAGAAGAGAGAAGAGAGAAGAGAGAAGAGAGAAGAGAGAAGAGAGAAGAGAGAA...ACGAGAGACGAGAGACGAGAGTAGAGAGACGAGAGAAGAGAGACGAGAGAAGAGAGAAGAGAG dm3_rmsk_(GA)n ra... [154989] 119 ATCGACTCCCTGGAAATCCCAAATGGGGATGCAGAAAGTATGGCGGCTATCCTCATGAACATGCTGGGCAGACTTTGCGACCCATCCATGCCAAGAAAAAAGAAGCCAAAATGGAAGCC dm3_rmsk_DMRT1B r... [154990] 52 ATGGAAGCCGGTCCAAAAAGATACCGTGTTTCGGCAGGTGTTCCCCAAGGAT dm3_rmsk_DMRT1B r... [154991] 155 TCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTC...TCTTCTCTTCTCTTCTCTTCTCATCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCTTCTCT dm3_rmsk_(TTCTC)n... [154992] 150 AGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAG...AGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGAAGAGACGAGA dm3_rmsk_(GAGAA)n... [154993] 147 ACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAA...GACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGACAAGAC dm3_rmsk_(CAAGA)n...
A quick question: when I use "" with the file name as you did I get the folloiwng error: 
fly_rm=read.DNAStringSet("Fly_RM", format="fasta")

Error in file(file, "r") : cannot open the connection
In addition: Warning message:
In file(file, "r") : cannot open file 'Fly_RM': No such file or directory Error in XStringSet("DNAString", x, start = start, end = end, width = width, : error in evaluating the argument 'x' in selecting a method for function 'XStringSet' 



sessionInfo()

R version 2.7.0 (2008-04-22)
i386-apple-darwin8.10.1

locale:
en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] tools stats graphics grDevices utils datasets methods base 

other attached packages:
[1] BiostringsCinterfaceDemo_0.1.2 BSgenome.Dmelanogaster.UCSC.dm3_1.3.7 BSgenome_1.8.3 Biostrings_2.8.9 
[5] Biobase_2.0.1





Thanks









On Tue, 03 Jun 2008 19:36:47 -0700
 Herve Pages <[EMAIL PROTECTED]> wrote:

Joseph Dhahbi, P.h.D. wrote:


I downloaded RepeatMasker from the Table Browser:
http://genome.ucsc.edu/cgi-bin/hgTables?command=start


Thanks! So I downloaded it and loaded in R with:
> allrepeats <- read.DNAStringSet("dm3rm", format="fasta") 
  Read 1126586 items
  > allrepeats
    A DNAStringSet instance of length 154993
width seq names [1] 433 AATTCGCGTCCGCTTACCCAT...GTTATTCACGGGGTCGGGGA dm3_rmsk_NINJA_I ... [2] 177 TGTCGCGGATCGAACGGTGCA...GCGGTAAGATAAGATTAAAT dm3_rmsk_NINJA_LT... [3] 1086 ATACGATGGCTGTACCTCATG...CATGGCGCTTTCCATCAGGA dm3_rmsk_Baggins1... [4] 62 TCTGCATAGTCTCCGTGTCCA...GGTGCAGCATAAACACTTCT dm3_rmsk_DMLTR5 r... [5] 346 AAATAATCCTGGAGGAGGCAA...GTCCTTCAGGACTCTATACT dm3_rmsk_Gypsy_I ... [6] 243 CCTGGCTCAGTTATTCCGAGA...GAGTATAGAGTCCTGAAGGA dm3_rmsk_DMRT1C r... [7] 149 CGGATGCTATAGGCCGCACCC...CTCAGTTATTCCGAGAGCAC dm3_rmsk_DMRT1C r... [8] 230 AACCAGAATCACTGCACGGCT...GCAGATGGTGGGTGTACAGC dm3_rmsk_DMRT1C r... [9] 210 CAGCAGCAGTGGTGGCTATCA...CATAGCGGGCCACTGTAGAA dm3_rmsk_DMRT1C r... 
       ...   ... ...
[154985] 164 GAGATGAGATGAGAAGAGAAG...ATAGAATAGAAGAGAATAGA dm3_rmsk_(GAGAA)n... [154986] 169 AGAAGAGAAGAGAAGAGAAGA...GATAAGAGAATAGATGAGAA dm3_rmsk_(GAGAA)n... [154987] 157 CAACACAACACAACACAACAC...ACACAACATTTCACAACACA dm3_rmsk_(CACAA)n... [154988] 153 AGAGAGCAGAGAGAAGAGAGA...CGAGAGAAGAGAGAAGAGAG dm3_rmsk_(GA)n ra... [154989] 119 ATCGACTCCCTGGAAATCCCA...AAGAAGCCAAAATGGAAGCC dm3_rmsk_DMRT1B r... [154990] 52 ATGGAAGCCGGTCCAAAAAGA...GGCAGGTGTTCCCCAAGGAT dm3_rmsk_DMRT1B r... [154991] 155 TCTCTTCTCTTCTCTTCTCTT...TCTCTTCTCTTCTCTTCTCT dm3_rmsk_(TTCTC)n... [154992] 150 AGAGAAGAGAAGAGAAGAGAA...AGAGAAGAGAAGAGACGAGA dm3_rmsk_(GAGAA)n... [154993] 147 ACAAGACAAGACAAGACAAGA...AAGACAAGACAAGACAAGAC dm3_rmsk_(CAAGA)n... 

It works fine for me. Here is my sessionInfo():

  > sessionInfo()
  R version 2.7.0 (2008-04-22)
  x86_64-unknown-linux-gnu

  locale:

LC_CTYPE=en_US.UTF-8;LC_NUMERIC=C;LC_TIME=en_US.UTF-8;LC_COLLATE=en_US.UTF-8;LC_MONETARY=C;LC_MESSAGES=en_US.UTF-8;LC_PAPER=en_US.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.UTF-8;LC_IDENTIFICATION=C

  attached base packages:
[1] stats graphics grDevices utils datasets methods base 

  other attached packages:
  [1] Biostrings_2.8.12

Could you send yours please?


I will try your suggestion
Note that if you only want to count the number of matches of your PDict object in the repeat regions of the *entire* genome (without actually caring about where the matches are actually occuring), then you could do something like this (use read.XStringViews instead of read.DNAStringSet):
> allrepeats <- read.XStringViews("dm3rm", format="fasta", subjectClass="DNAString", collapse="-") 

and then call countPDict(pdict, subject(allrepeats)).
The purpose of using 'collapse="-"' when reading the file is to separate the RepeatMasker chunks with this letter when they are all put together in the big 'subject(allrepeats)' sequence. This is in order to avoid the matches 
that would span across several chunks.

Cheers,
H.



Thank you for your help

On Tue, 03 Jun 2008 17:14:10 -0700
 Herve Pages <[EMAIL PROTECTED]> wrote:

Hi Joseph,

Joseph Dhahbi, P.h.D. wrote:


Hi
I downloaded the drosophila RepeatMasker from UCSC GB as a text file which is in fasta format and looks like this:
dm3_rmsk_NINJA_I range=chr4:2-434 5'pad=0 3'pad=0 strand=+ repeatMasking=none 
AATTCGCGTCCGCTTA......
dm3_rmsk_NINJA_LTR range=chr4:435-611 5'pad=0 3'pad=0 strand=+ repeatMasking=none 
TGTCGCGGATC....
dm3_rmsk_Baggins1 range=chr4:638-1723 5'pad=0 3'pad=0 strand=- repeatMasking=none 
ATACGATGG......
I made the input dictionary and I would like to make the RepeatMasker sequences as the target. When I used ‘read.DNAStringSet’ it recognized only the first sequence of the fasts file. Ho do I merge all of the sequences in and make them as a target.
If your file is really FASTA then read.DNAStringSet() should extract all the records and return a DNAStringSet object where each element corresponds to a record in the original file. So it seems like you've hit a bug in the read.DNAStringSet() function. Can you please provide the URL to the file 
you downloaded so we can try to reproduce?
Anyway, what you are trying to achieve can be done in an easier (and more efficient) way. You don't need to download the RepeatMasker sequences for this; just use the BSgenome.Dmelanogaster.UCSC.dm3 package. The RepeatMasker information is already included in it as part of the built-in masks provided 
for each chromosome:

  > library(BSgenome.Dmelanogaster.UCSC.dm3)
  > Dmelanogaster
  Fly genome
  |
  | organism: Drosophila melanogaster
  | provider: UCSC
  | provider version: dm3
  | release date: Apr. 2006
  | release name: BDGP Release 5
  |
  | single sequences (see '?seqnames'):
| chr2L chr2R chr3L chr3R chr4 chrX chrU | chrM chr2LHet chr2RHet chr3LHet chr3RHet chrXHet chrYHet 
  |   chrUextra
  |
  | multiple sequences (see '?mseqnames'):
  |   upstream1000  upstream2000  upstream5000
  |
| (use the '$' or '[[' operator to access a given sequence) 
  > chr2L <- Dmelanogaster$chr2L
  > chr2L
23011544-letter "MaskedDNAString" instance (# for masking) seq: CGACAATGCACGACAGAGGAAGCAGAACAGATATTT...GCATATTTGCAAATTTTGATGAACCCCCCTTTCAAA 
  masks:
maskedwidth maskedratio active names 1 200 8.691290e-06 FALSE assembly gaps 2 1966561 8.545976e-02 FALSE RepeatMasker 3 61603 2.677048e-03 FALSE Tandem Repeats Finder [period<=12] 
  all masks together:
    maskedwidth maskedratio
        1988181  0.08639929
  all active masks together:
    maskedwidth maskedratio
              0           0
Note that the built-in masks are always inactive by default. To activate 
a mask do:
> active(masks(chr2L))[2] <- TRUE # activate the RepeatMasker mask
Now only the parts of chr2L that are NOT repeat regions are visible. 
To invert this, use gaps():

  > chr2Lrepeats <- gaps(chr2L)
  > chr2Lrepeats
23011544-letter "MaskedDNAString" instance (# for masking) seq: #GACAATGCACGACAGAGGAAGCAGAACAGATATTT...GCATATTTGCAAATTTT################### 
  masks:
    maskedwidth maskedratio active
  1    21044983   0.9145402   TRUE
Then use matchPDict() (or countPDict()) in the usual way.
The GenomeSearching vignette in the BSgenome package has more information about masking (some sections are still incomplete but 
will be completed soon).

Hope this helps,
H.



Thank you for your help


Regards,
Joseph

Joseph M. Dhahbi, PhD
Childrens Hospital Oakland Research Institute
5700 Martin Luther King Jr. Way
Oakland, CA 94609
USA
Ph.(510)428-3885 EXT.5743
Cell.(702)335-0795
Fax (510)450-7910
[EMAIL PROTECTED]
The email message (and any attachments) is for the sole...{{dropped:3}} 

_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing







Regards,
Joseph

Joseph M. Dhahbi, PhD
Childrens Hospital Oakland Research Institute
5700 Martin Luther King Jr. Way
Oakland, CA 94609
USA
Ph.(510)428-3885 EXT.5743
Cell.(702)335-0795
Fax (510)450-7910
[EMAIL PROTECTED]
The email message (and any attachments) is for the sole use of the intended recipient(s) and may contain confidential information. Any unauthorized review, use, disclosure or distribution is prohibited. If you are not the intended recipient, please contact the sender by reply email and destroy all copies of the original message (and any attachments). Thank You. 








Regards,
Joseph

Joseph M. Dhahbi, PhD
Childrens Hospital Oakland Research Institute
5700 Martin Luther King Jr. Way
Oakland, CA 94609
USA
Ph.(510)428-3885 EXT.5743
Cell.(702)335-0795
Fax (510)450-7910
[EMAIL PROTECTED]
The email message (and any attachments) is for the sole...{{dropped:3}}

_______________________________________________
Bioc-sig-sequencing mailing list
[email protected]
https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing

Reply via email to