[BioC] ChIPpeakAnno

[Bahar]

Paul Shannon <pshannon at ...> writes:

> 
> I offer below a short example of using rGADEM with MotifDb which you may
find useful.  
> 
> The example motivates and then demonstrates "seeded" search, in which a
candidate TF is provided to
> rGADEM.  (The rGADEM vignette illustrates the simpler task of de novo
motif search.
> 
>  - Paul
> 
> On Feb 27, 2013, at 9:03 AM, Zhu, Lihua (Julie) wrote:
> 
> > Jose,
> > 
> > Thanks for your kind word!
> > 
> > You might want to check out rGADEM package for de nova motif discovery.
> > 
> > Best regards,
> > 
> > Julie
> > 
> 
> library(rGADEM)
> library(MotifDb)
> library(BSgenome.Hsapiens.UCSC.hg19)
> path <- system.file("extdata/Test_100.fasta", package="rGADEM")
> sequences <- readDNAStringSet(path, "fasta")
> length(sequences)   # 49
> 
>    # you can use the UCSC genome browser to BLAT the first 4 sequences,
one at a time,
>    # against human hg19, and see what TFs bind there
>    #   http://genome.ucsc.edu/cgi-bin/hgBlat?command=start
> 
> as.character(substr(sequences[1:4], 1, 80))
> 
>    # Once you blat, click through to view the blat match in the genome
browser,
>    # and enable the "Transcription Factor ChIP-seq from ENCODE" track.
>    # The FOXA1 Tf has a high score in this region of chr1
>    # for most of these sequences
> 
>    # obtain the reported motif for this TF from the Bioc MotifDb package
> 
> query(MotifDb, 'foxa1')
> 
>    # MotifDb object of length 1
>    # | Created from downloaded public sources: 2012-Nov-01
>    # | 1 position frequency matrices from 1 source:
>    # |        JASPAR_CORE:    1
>    # | 1 organism/s
>    # |           Hsapiens:    1
>    # Hsapiens-JASPAR_CORE-FOXA1-MA0148.1 
> 
>    # now use rGADEM to do a seeded search:
> 
> pwm.foxa1 <- as.list(MotifDb["Hsapiens-JASPAR_CORE-FOXA1-MA0148.1"])
> gadem.foxa1 <-GADEM(sequences,verbose=1,genome=Hsapiens,Spwm=pwm.foxa1,
seed=TRUE)
> consensus(gadem.foxa1)
> 
>    # [1] "nTGTTTACwyw"      "GmwGrrrsswGGvAGn"
> 
>    # how do the 49 sequences match to these two motifs?
> mean(sapply(gadem.foxa1 <at> motifList[[1]] <at> alignList, function(x) x
<at> pval))
> mean(sapply(gadem.foxa1 <at> motifList[[2]] <at> alignList, function(x) x
<at> pval))
> 
> On Feb 27, 2013, at 9:03 AM, Zhu, Lihua (Julie) wrote:
> 
> > Jose,
> > 
> > Thanks for your kind word!
> > 
> > You might want to check out rGADEM package for de nova motif discovery.
> > 
> > Best regards,
> > 
> > Julie
> > 
> 
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at ...
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
http://news.gmane.org/gmane.science.biology.informatics.conductor
> 
> 

Hi Paul,
Thank you very much for your example rGADEM and motIV analysis.
I am working with a set of sequences in fasta format and looking for
potential motif denovo. Once identified, I wanted to verify if those match
with available databases (JASPAR etc.) However I fail to get already
reported motifs using the same sequences. I copied the script below for your
input and suggestions please.
library(rGADEM)
pwd<-"" #INPUT FILES- BedFiles, FASTA, etc.
path<- system.file("extdata/mysequences.fasta",package="rGADEM")
FastaFile<-paste(pwd,path,sep="")
Sequences <- readDNAStringSet(FastaFile, "fasta")

gadem<-GADEM(Sequences,
             seed=1,
             verbose=1,
             numWordGroup=3,
             numTop3mer=20,
             numTop4mer=40,
             numTop5mer=60,
             numGeneration=5,
             pValue=0.0002,
             eValue=0.0,
             maskR=1,
             maxSpaceWidth=10,
             useChIPscore=0,
             numEM=40,
             fEM=0.5,
             fullScan=1,
             slideWinPWM=6,
             stopCriterion=1,
             numBackgSets=10,
             weightType=0,
             bFileName="NULL",
             nmotifs=25,
             extTrim=1,
             minSpaceWidth=0,
             minSites =-1)

motifs <- getPWM(gadem)
print(motifs)
Thank you very much for your comments.
Bahar