[BioC] predictCoding() return empty Granges
Hervé Pagès
hpages at fhcrc.org
Fri Oct 5 08:01:26 CEST 2012
Hi Rebecca,
On 10/04/2012 02:44 PM, sun wrote:
> After I successfully adjusted the chromosomes of VCF and TranscriptDb
> objects to match the names of the BSgenome object, I ran predictCoding(),
> it returned nothing except a warning message, see below,
>
>> coding <- predictCoding(vcf, txdb, seqSource=Athaliana)
> *Warning message:
> In .setActiveSubjectSeq(query, subject) :
> circular sequence(s) in query 'chrM' ignored*
>> coding
> GRanges with 0 ranges and 1 metadata column:
> seqnames ranges strand | GENEID
> <Rle> <IRanges> <Rle> | <character>
> ---
> seqlengths:
>
> Here are the values of VCF, TranscriptDb and BSgenome Object,
>
>> vcf
> class: VCF
> dim: 551201 2
> genome: TAIR10
> exptData(1): header
> fixed(4): REF ALT QUAL FILTER
> info(18): DP DP4 ... PR VDB
> geno(6): GT GQ ... SP PL
> rownames(551201): Chr1:711 Chr1:956 ... ChrM:366919 ChrM:366920 *Qustion:
> After used renameSeqlevels() to change chr name of vcf (Chr to chr), the
> rownames here are still old chr name. not sure if this will be the problem
> for predictCoding(). *
> rowData values names(1): paramRangeID
> colnames(2): sample1_aln.sorted.bam sample2_aln.sorted.bam
> colData names(1): Samples
>
>> txdb
> TranscriptDb object:
> | Db type: TranscriptDb
> | Supporting package: GenomicFeatures
> | Data source: TAIR 10
> | Genus and Species: Arabodpsis
> | miRBase build ID: NA
> | transcript_nrow: 35386
> | exon_nrow: 207465
> | cds_nrow: 197160
> | Db created by: GenomicFeatures package from Bioconductor
> | Creation time: 2012-09-28 16:43:28 -0700 (Fri, 28 Sep 2012)
> | GenomicFeatures version at creation time: 1.9.37
> | RSQLite version at creation time: 0.11.1
> | DBSCHEMAVERSION: 1.0
>
>> Athaliana
> Arabidopsis genome
> |
> | organism: Arabidopsis thaliana (Arabidopsis)
> | provider: TAIR
> | provider version: 04232008
> | release date: NA
> | release name: dumped from ADB: Mar/14/08
> |
> | sequences (see '?seqnames'):
> | chr1 chr2 chr3 chr4 chr5 chrC chrM
> |
> | (use the '$' or '[[' operator to access a given sequence)
We provide 2 BSgenome packages for Arabidopsis: a very old one
BSgenome.Athaliana.TAIR.04232008
with chromosome names and lengths:
> seqlengths(Athaliana)
chr1 chr2 chr3 chr4 chr5 chrC chrM
30432563 19705359 23470805 18585042 26992728 154478 366924
and a somewhat less old one:
BSgenome.Athaliana.TAIR.TAIR9
with chromosome names and lengths:
> seqlengths(Athaliana)
Chr1 Chr2 Chr3 Chr4 Chr5 ChrM ChrC
30427671 19698289 23459830 18585056 26975502 366924 154478
They correspond to 2 different assemblies of course.
2 things:
(1) The code I showed you in the previous thread for renaming
the sequence names of your TranscriptDb object was assuming
that you were using the less old BSgenome package:
seqlevels(txdb) <- sub("^c", "C", seqlevels(txdb))
However it seems that you are using the very old one (where
chr names are "chr*"), and that your txdb chr names are
"Chr*", so you would actually need to do something like:
seqlevels(txdb) <- sub("^C", "c", seqlevels(txdb))
(2) However, and this is much more important than (1): it seems that
your VCF and TranscriptDb objects are based on the TAIR10 assembly.
So it makes no sense to use anything else but a BSgenome
package that contains the exact same assembly. Otherwise it's
very likely that most of the predictions returned by
predictCoding() will be wrong!
Unfortunately at the moment we don't provide a BSgenome package for
TAIR10 (nobody requested one so far). I'll try to make one in the
next couple of weeks. I'll also remove
BSgenome.Athaliana.TAIR.04232008 from the devel branch (BioC 2.12).
Cheers,
H.
>
>> sessionInfo()
> R version 2.15.1 (2012-06-22)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
> [1] C
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] BSgenome.Athaliana.TAIR.04232008_1.3.18
> BSgenome_1.26.0
> VariantAnnotation_1.4.0
> [4] Rsamtools_1.10.1
> Biostrings_2.26.0
> GenomicFeatures_1.10.0
> [7] AnnotationDbi_1.20.0
> Biobase_2.18.0
> GenomicRanges_1.10.0
> [10] IRanges_1.16.2
> BiocGenerics_0.4.0
> vimcom_0.9-2
> [13] setwidth_1.0-0
> colorout_0.9-9
>
> loaded via a namespace (and not attached):
> [1] DBI_0.2-5 RCurl_1.95-0.1.2 RSQLite_0.11.2
> XML_3.95-0.1 biomaRt_2.14.0 bitops_1.0-4.1
> [7] parallel_2.15.1 rtracklayer_1.18.0 stats4_2.15.1
> tools_2.15.1 zlibbioc_1.4.0
>
> any ideas? thank you.
>
> Rebecca
>
>
> On Thu, Oct 4, 2012 at 1:18 PM, Hervé Pagès <hpages at fhcrc.org> wrote:
>
>> Hi Rebecca,
>>
>>
>> On 10/04/2012 12:10 PM, sun wrote:
>>
>>> Hi All,
>>>
>>> I am going to use "coding <- predictCoding(vcf, txdb,
>>> seqSource=Athaliana)"
>>> to detect coding SNPs. The problem is that the chromosome names are not
>>> consistent among VCF, txdb and BSgenome. In vcf, the chromosome name is
>>> "Chr*", in txdb, the chr name is "Chr", but in BSgenome, the chr name is
>>> "chr*" .
>>>
>>> I know I can use renameSeqlevels() to adjust the seqlevels (chromosome
>>> names) of the VCF object to match that of the txdb annotation. But how can
>>> I adjust the chr name of BSgenome or TranscriptDB?
>>>
>>
>> In BioC 2.11 (released yesterday), you can rename the chromosomes of a
>> TranscriptDb object, so you could rename the chromosomes of your
>> VCF and TranscriptDb objects to match the names of the BSgenome object.
>>
>> E.g. for the TranscriptDb object:
>>
>> seqlevels(txdb) <- sub("^c", "C", seqlevels(txdb))
>>
>> Note that renaming the chromosomes of a TranscriptDb object is a new
>> feature and is not fully implemented yet. For example, if you use
>> select() on the object you'll still get the original names (those
>> stored in the db), and if you try to specify a chromosome name thru
>> the 'vals' arg of the transcripts(), exons() and cds() extractors,
>> you still need to use the original names. This will be addressed soon.
>>
>> Our plan is to also support renaming of the chromosomes of BSgenome
>> and SNPlocs objects very soon.
>>
>> Also, an additional level of convenience will be provided via the
>> seqnameStyle() getter and setter, so you'll be able to quickly rename
>> with something like:
>>
>> seqnameStyle(x) <- "UCSC"
>>
>> or
>>
>> seqnameStyle(vcf) <- seqnameStyle(txdb) <- seqnameStyle(genome)
>>
>> This will work on almost any 'x' object that contains chromosome
>> names (GRanges, GRangesList, GappedAlignments, TranscriptDb, VCF,
>> BSgenome, SNPlocs, etc...)
>>
>> Cheers,
>> H.
>>
>>
>>
>>
>>> Thanks,
>>>
>>> Rebecca
>>>
>>> [[alternative HTML version deleted]]
>>>
>>> ______________________________**_________________
>>> Bioconductor mailing list
>>> Bioconductor at r-project.org
>>> https://stat.ethz.ch/mailman/**listinfo/bioconductor<https://stat.ethz.ch/mailman/listinfo/bioconductor>
>>> Search the archives: http://news.gmane.org/gmane.**
>>> science.biology.informatics.**conductor<http://news.gmane.org/gmane.science.biology.informatics.conductor>
>>>
>>>
>> --
>> Hervé Pagès
>>
>> Program in Computational Biology
>> Division of Public Health Sciences
>> Fred Hutchinson Cancer Research Center
>> 1100 Fairview Ave. N, M1-B514
>> P.O. Box 19024
>> Seattle, WA 98109-1024
>>
>> E-mail: hpages at fhcrc.org
>> Phone: (206) 667-5791
>> Fax: (206) 667-1319
>>
>
> [[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
>
--
Hervé Pagès
Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B514
P.O. Box 19024
Seattle, WA 98109-1024
E-mail: hpages at fhcrc.org
Phone: (206) 667-5791
Fax: (206) 667-1319
More information about the Bioconductor
mailing list