[BioC] alignShortReads

[Martin Morgan]

On 08/21/2013 02:54 PM, Sonja Althammer wrote:
> Hi!
> I am having trouble aligning a DNAStringSet to a BSgenome. The example in
> the help works for me....
> Could anyone explain to me the error I am getting, please and how I can
> make it work? Actually I would like to align these reads to a specific
> gene... Is there maybe a better way?
> Thanks!
> Sonja
>
>
>> alignShortReads(reads, Hsapiens, seqNames="chr1")
> Error in .Call2("ACtree2_build", tb, pp_exclude, base_codes, nodebuf_ptr,

alignShortReads is from R453Plus1Toolbox.

DNAStringSet can contain ambiguity letters such as 'M', to indicate that the 
nucleotide could be either A or C

 > IUPAC_CODE_MAP
      A      C      G      T      M      R      W      S      Y      K      V
    "A"    "C"    "G"    "T"   "AC"   "AG"   "AT"   "CG"   "CT"   "GT"  "ACG"
      H      D      B      N
  "ACT"  "AGT"  "CGT" "ACGT"

but the aligner in alignShortReads expects only

 > DNA_BASES
[1] "A" "C" "G" "T"

If your sequences do contain just bases, it would make sense to run 
Biostrings::matchPDict directly on the sequence of the gene of interest. You 
could obtained the sequence of your gene(s) of interest by, e.g., using the 
function ?select and the Homo.sapiens package to get coordinates of the gene, 
?getSeq to extract the sequence from the Hsapiens BSgenome object you already 
have, and matchPDict to do the alignment.

Probably in a way that requires more care, I did the following to arrive at the 
sequence of the transcripts of BRCA1

     library(BSgenome.Hsapiens.UCSC.hg19)
     library(Homo.sapiens)

Discover what I can do:

     ## 'keytypes' I can query with...
     keytypes(Homo.sapiens)
     ## discover the correct way to specify a gene symbol I'm interested in
     grep("BRCA", keys(Homo.sapiens, "SYMBOL"))
     ## discover what info I can extract
     cols(Homo.sapiens)

Now create ranges for transcripts

     ## select genomic coordinates for BRCA1
     xx = select(Homo.sapiens, "BRCA1",
            c("TXID", "TXCHROM", "TXSTRAND", "TXSTART", "TXEND"),
            "SYMBOL")
     ## turn coordinates into a GRanges
     yy = with(xx, GRanges(TXCHROM, IRanges(TXSTART, TXEND), TXSTRAND))

or alternatively map to the central 'Entrez' identifier and...

     ## figure out Entrez ID corresponding to symbol
     eid = select(Homo.sapiens, "BRCA1", "ENTREZID", "SYMBOL")$ENTREZID
     ## extract transcript coordinates from TxDb
     yy = transcriptsBy(TxDb.Hsapiens.UCSC.hg19.knownGene, "gene")[eid]

What's a gene? Maybe it's just the range of transcripts:

     gn = range(yy)

Finally retrieve the relevant sequence(s)

     seq = getSeq(Hsapiens, yy)  ## or getSeq(Hsapiens, gn)

Judging by the package you're using, you might have longer 454 reads, where the 
ungapped high fidelity matches of matchPDict are not appropriate. If you're 
interested in more general alignments, or if the sequences you're trying to 
align do have ambiguity letters, there is Biostrings::pairwiseAlignment.  It 
might still be that this general purpose aligner is not appropriate for the task 
you're trying to accomplish, so some detail (and perhaps input from others on 
the list!) might be needed.

Hope that provides some guidance.

Martin

> :
>    non base DNA letter found in Trusted Band for pattern 248
>> class(reads)
> [1] "DNAStringSet"
> attr(,"package")
> [1] "Biostrings"
>
>> class(Hsapiens)
> [1] "BSgenome"
> attr(,"package")
> [1] "BSgenome"
>
> 	[[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
>


-- 
Computational Biology / Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N.
PO Box 19024 Seattle, WA 98109

Location: Arnold Building M1 B861
Phone: (206) 667-2793