On Tue, Mar 1, 2011 at 4:20 AM, Pages, Herve <hpages@fhcrc.org> wrote:

> Hi Ravi,
>
> Sean's suggestion to map the probes directly to the transcriptome can
> easily be achieved with the Biostrings/BSgenome/GenomicFeatures
> infrastructure.
>
> The code below uses vwhichPDict (a variant of matchPDict) and the danRer6
> genome from UCSC (release name: Sanger Institute Zv8). Note that UCSC just
> released danRer7 (Sanger Institute Zv9) but we don't have a BSgenome
> package
> for it yet.
>
> Also I don't think we support the Nimblegen Zebrafish 12 x135K Expression
> chip so I'm using the zebrafish chip from Affymetrix instead.
>
> Depending on your machine and the quality of your internet connection,
> this code should run in 8-15 minutes. Less than 800 Mb of RAM is used.
>
>  ## Load the probes:
>  library(zebrafishprobe)
>  library(Biostrings)
>  probes <- DNAStringSet(zebrafishprobe)
>
>  ## Compute the transcriptome (32992 transcripts):
>  library(GenomicFeatures)
>  txdb <- makeTranscriptDbFromUCSC(genome="danRer6", tablename="ensGene")
>  library(BSgenome.Drerio.UCSC.danRer6)
>  txseqs <- extractTranscriptsFromGenome(Drerio, txdb)
>  ## Note that saving 'txseqs' with write.XStringSet() generates
>  ## a FASTA file that is 51M.
>
>  ## Map 'probes' to 'txseqs' (allowing 2 mismatches):
>  pdict2 <- PDict(probes, max.mismatch=2)
>  m2 <- vwhichPDict(pdict2, txseqs, max.mismatch=2)
>  makeMap <- function(m, probeids, txnames)
>  {
>    probeid <- probeids[unlist(m)]
>    txname <- rep.int(txnames, elementLengths(m))
>    ans <- unique(data.frame(probeid, txname))
>    rownames(ans) <- NULL
>    ans
>  }
>  map <- makeMap(m2, zebrafishprobe$Probe.Set.Name, names(txseqs))
>
> Note that the longest step is not the call to vwhichPDict() but the
> extraction of the transcriptome.
>
> However, the final result doesn't look that good. A quick look at the
> mapping shows that only 62.4% of the probe ids are mapped and that some
> probe ids are mapped to more than 1 transcript:
>
>
Nice example, Herve!

Mapping to multiple transcripts is to be expected, so one needs to come up
with ways of dealing with the situation.  Missing probes is also typical.
 These can be dealt with by mapping to the unmatched probes to other
transcript databases.  For example, one could map to Ensembl first, RefSeq
second, all zebrafish mRNAs next, all zebrafish ESTs next.


>  > head(map)
>             probeid             txname
>  1 Dr.19494.1.S1_at ENSDART00000048610
>  2 Dr.19494.1.S1_at ENSDART00000103479
>  3 Dr.19494.1.S1_at ENSDART00000103478
>  4 Dr.10343.1.S1_at ENSDART00000006449
>  5 Dr.12057.1.A1_at ENSDART00000054814
>  6 Dr.22271.1.A1_at ENSDART00000054814
>  > length(unique(zebrafishprobe$Probe.Set.Name))
>  [1] 15617
>  > length(unique(map$probeid))
>  [1] 9746
>  > any(duplicated(map$probeid))
>  [1] TRUE
>
> But maybe you will have more luck with the Nimblegen Zebrafish 12 x135K
> Expression chip.
>
>
Keep in mind that with 60mer probes, there is potentially a need to allow
gaps in the alignment and that consecutive matches of even 40 or 50 base
pairs are likely to generate a signal on a microarray.

Sean


> One more thing. While testing the above code, I discovered 2 problems
> that were preventing it to work: one with the
> extractTranscriptsFromGenome()
> function and one with the vwhichPDict() function. Those problems are now
> fixed in the release and devel versions of GenomicFeatures (v 1.2.4 &
> 1.3.13)
> and Biostrings (v 2.18.3 & 2.19.12). These new versions will become
> available
> via biocLite() in the next 24-36 hours.
>
> Hope this helps.
> H.
>
>
> ----- Original Message -----
> From: "Ravi Karra" <ravi.karra@gmail.com>
> To: "Sean Davis" <sdavis2@mail.nih.gov>
> Cc: bioconductor@r-project.org
> Sent: Saturday, February 26, 2011 2:56:02 PM
> Subject: Re: [BioC] Mapping microarray probes to the genome using
> findOverlaps
>
> Hi Sean,
> I could do that, but am not sure how to.  The annotated zebrafish genome is
> the Tubingen strain and the some of the probes on the array are from the EK
> and AB strains.   This means that I need to allow for SNP's in the
> alignment.  I originally tried to align the probes to the Ensembl
> Transcripts using matchPDict but when I allowed for 2 mismatches
> (max.mismatch = 2) across the probe sequences, my computer never stopped
> running the program!  I found TopHat to be much faster (8 min) and TopHat
> allows for a few nucleotide wobble by default!.
> Do you have a suggestion(s) on another way to align the array to the
> Ensembl Transcripts?
> Thanks,
> Ravi
>
>
>
>
> On Feb 26, 2011, at 5:45 PM, Sean Davis wrote:
>
> >
> >
> > On Sat, Feb 26, 2011 at 5:29 PM, Ravi Karra <ravi.karra@gmail.com>
> wrote:
> > Hello,
> >
> > I am still trying to map probes on the Nimblegen Zebrafish 12 x135K
> Expression to the Zv9 version of the zebrafish genome available from
> Ensembl!  I am very reluctantly pursuing an alignment approach to annotation
> as the original annotation provided with the array is quite outdated.
> >
> > I performed a gapped alignment using the individual probe sequences
> (60-mers) from the array using TopHat and loaded the results into
> Bioconductor as a GappedAlignments object.  I have made a TranscriptDb
> object using the Zv9 genome from Ensembl.  Next, I plan to use findOverlaps
> for the annotation. What is the best way to get the overlaps (by exon, cds,
> or by transcript)?  I am a little concerned that using transcriptsByOverlaps
> might be a bit too broad and result in mapping reads to multiple genes (for
> example transcripts in the genome that have overlapping genomic
> coordinates). By contrast, mapping with exonsByOverlaps and cdsByOverlaps
> might be too restrictive and miss information in the UTR's.  My gut feeling
> is that annotating by cds is the "in-between" approach.  What is the
> recommended approach for RNA-seq?  As you can tell, I am quite new to this!
> >
> >
> > Hi, Ravi.  Sorry to answer a question with more questions, but why not
> just map the probes against Ensembl Transcripts or refseq?  What is the
> advantage of mapping to the genome and then going back to the transcripts?
> >
> > Sean
>
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@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]]

