[BioC] Error in locateVariants

[Valerie Obenchain]

Thanks for sending the files. Please keep the conversation 'on-list' by 
cc'ing bioconductor at r-project.org.

There was a bug in locateVariants() related to 'features' not having 
txid or cdsid. That is now fixed in devel (1.7.38) and release (1.6.7). 
Thanks for reporting this.

A couple things to note:

(1) The second argument to readVcf() is 'genome'. Here you are setting 
the genome to the string 'Pinfestans.fa' string which is probably not 
what you intended to do.

     vcf1<-readVcf("CGTACG_filteredSNPS.vcf","Pinfestans.fa")

(2) When 'features' is a GRangesList it must be an appropriate 
extraction to match the 'region' argument (explained on the man page). 
You are interested in CodingVariants so the GRangesList must be 
equivalent to a cdsBy(txdb, "transcript") extraction. The 
'myGRangesList' object is of length 1 with 801 different seqlevels in 
the single list element.

myGRangesList<-GRangesList(myGRanges)
 > length(myGRangesList)
[1] 1
 > length(seqlevels(myGRangesList))
[1] 801

The idea is that cds regions from the same transcript (i.e., same 
seqlevel) are grouped in list elements. Assuming that the seqlevels in 
this file are transcript names, you could subset by 'type' and split on 
seqname to get cds by transcripts.

 > cds <- myGRanges[myGRanges$type == "CDS"]
 > cdsby <- split(cds, seqnames(cds))
 > length(cdsby)
[1] 801

An alternative approach to is to make a TranscriptDb from your gff file 
with makeTranscriptDbFromGFF() in GenomicFeatures. Either the txdb or 
cdsby extraction can be used as 'features'.

     txdb <- makeTranscriptDbFromGFF("annotation.gff3")
     cdsby <- cdsBy(txdb, "transcript")

Oops. Looks like this gff file has something strange going on. I'll send 
the report to Marc.

 >  txdb <- 
makeTranscriptDbFromGFF("PI_T30-4_FINAL_CALLGENES_3.annotation.gff3")
Error in .parse_attrCol(attrCol, file, colnames) :
   Some attributes do not conform to 'tag=value' format

At any rate, makeTranscriptDbFromGFF() can be a useful function to 
convert gff to txdb for well behaved files.

Valerie


On 07/24/2013 01:34 PM, Valerie Obenchain wrote:
> Hi Jolly,
>
> I need a reproducible example in order to help. Are 'rd' and
> 'myGRangesList' too large to send? Also provide the output of
> sessionInfo().
>
> Valerie
>
>
> On 07/24/2013 11:38 AM, Jolly Shrivastava wrote:
>> Hello Valerie,
>>
>> I am running VariantAnnotation package to locate the position of SNPs
>> using
>> locateVariants.
>>
>> The commands that I gave to extract vcf file and the gff file are as
>> follows:
>>
>> ###Reading the vcf file
>> vcf1<-readVcf("CGTACG_filteredSNPS.vcf","Pinfestans.fa")
>> rd<-rowData(vcf1)
>>
>> ####imported and converted the gff3 file to GRangesList object
>> subject<-import.gff("annotation.gff3")
>> myGRanges<-as(subject,"GRanges")
>> myGRangesList<-GRangesList(myGRanges)
>>
>>
>> ### made sure that both of them have same sequence levels
>> rd<-keepSeqlevels(rd,intersect(seqlevels(rd),seqlevels(myGRangesList))
>> myGRangesList<-keepSeqlevels(myGRangesList,intersect(seqlevels(rd),seqlevels(myGRangesList))
>>
>>
>> ###Ran the locateVariants command
>> variants_found<-locateVariants(rd,myGRangesList,CodingVariants())
>>
>> ###Got the following error
>> Error in DataFrame(...) : different row counts implied by arguments
>>
>> Can you please point to what I am doing wrong?
>>
>> Regards
>> Jolly
>>
>>
>
> _______________________________________________
> 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