Hi Steve Thank you for your answer. Then I ran this script, but there is error in the end. And I don't know how to solve it. library(Biostrings) library(GenomicRanges) library(BSgenome.Amellifera.UCSC.apiMel2) bee.txdb<- makeTranscriptDbFromUCSC(genome="apiMel2", tablename="ensGene") tx <- transcripts(bee.txdb, columns=c("tx_id", "tx_name", "gene_id")) intragenic_seq<-getSeq(Amellifera,gaps(tx)) #####Error in .starfreeStrand(strand(names)) : cannot mix "*" with other strand values Yating > Hi, > > On Tue, Nov 29, 2011 at 6:25 AM, Yating Cheng <yating.cheng at="" charite.de=""> > wrote: >> Dear All, >> >> Does anyone know that how to extract intragenic sequences from the >> genome. >> >> Like in the genescan, it is mentioned that the predictions are based on >> transcriptional, translational and donor/acceptor splicing signals as >> well >> as the length and compositional distributions of exons, introns and >> intergenic regions. >> >> But I am not sure which function I should use. > > Use the GenomicFeatures package to build a gene annotation database > for your organism. By getting a bit more intimate with that package > (good place to start is by reading its vignette, along with that for > GenomicRanges) you will figure out how to identify where the > intergenic regions lie on your genome. > > Once you have the intergenic regions of your genome stored in a > GRanges object, its easy to use those ranges to get the sequence from > your genome using the appropriate BSgenome.* package for your > organism. The Biostrings::getSeq function could do that for you, among > other things. > > HTH, > -steve >> >> >> Thank you very much. >> >> >> >> >> Yating Cheng >> Molecular Medicine Master's Program >> Charit?-Universit?tsmedizin Berlin >> >> _______________________________________________ >> 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 > > > > -- > Steve Lianoglou > Graduate Student: Computational Systems Biology > ?| Memorial Sloan-Kettering Cancer Center > ?| Weill Medical College of Cornell University > Contact Info: http://cbio.mskcc.org/~lianos/contact >

Hi, On Sun, Dec 4, 2011 at 4:16 PM, Yating Cheng <yating.cheng at="" charite.de=""> wrote: > Hi Steve > > Thank you for your answer. Then I ran this script, but there is error in > the end. And I don't know how to solve it. > > library(Biostrings) > library(GenomicRanges) > > > library(BSgenome.Amellifera.UCSC.apiMel2) > bee.txdb<- makeTranscriptDbFromUCSC(genome="apiMel2", > ?tablename="ensGene") > ?tx <- transcripts(bee.txdb, columns=c("tx_id", "tx_name", "gene_id")) > > > intragenic_seq<-getSeq(Amellifera,gaps(tx)) > > #####Error in .starfreeStrand(strand(names)) : > ?cannot mix "*" with other strand values When you call `gaps()` you get ranges (usually as long as your chromosome) with strands assigned as '*' -- you have to filter those out. Running with your example, you'd do: R> intergenic <- gaps(tx) R> intergenic <- intergenic[strand(intergenic) != '*'] R> intergenic_seq <- getSeq(Amellifera, intergenic) On a related note -- I think I find that behavior from gaps (returning full chromosome-length ranges assigned to the '*' strand) more annoying than useful ... especially when none of my ranges defined in the thing I am calling gaps on has ranges assigned to `*` ... but a tale for a different thread, perhaps ... HTH, -steve -- Steve Lianoglou Graduate Student: Computational Systems Biology ?| Memorial Sloan-Kettering Cancer Center ?| Weill Medical College of Cornell University Contact Info: http://cbio.mskcc.org/~lianos/contact