[BioC] Finding coding SNPs with predictCoding
Valerie Obenchain
vobencha at fhcrc.org
Tue May 15 22:17:52 CEST 2012
Hi Thomas, Alex,
Marc has been working on a makeTranscriptDbFromGFF() function which is
now in GenomicFeatures v 1.9.11 in devel. This addresses the last of
Thomas' points below wrt the user providing annotations in GFF or GTF
format. The function creates a TxDb which can be given to
locateVariants() and predictCoding(). Let us know how it goes.
Valerie
On 04/05/2012 01:47 PM, Valerie Obenchain wrote:
> Hi Thomas, Alex,
>
> Changes below have been made in release 1.2.4 and devel 1.3.4.
>
> predictCoding :
>
> - Now returns a proteinID that is the triplet number in the protein
>
> - The refAA and varAA are the 1-letter amino acid code but not the
> 3-letter code (i.e., 'G' but not 'Gly'). We don't currently have a
> function that returns the 3-letter code. Is this of interest/use?
>
> - Over the next dev cycle I will implement methods for annotations
> (subject argument) to be gff/gtf. Currently genomes (seqSource
> argument) can be a BSgenome or fasta file.
>
>
> locateVariants:
>
> To add flexibility for finding variants in a particular region the
> function now dispatches on a 'region' argument (e.g.,
> CodingVariants(), IntronVariants(), etc.). I've included a
> SpliceSiteVariants() that counts ranges that overlap with any portion
> of the first 2 or last 2 nucleotides of an intron. Called as,
>
> locateVariants(query, subject, range=SpliceSiteVariants())
>
> Thanks again for the suggestions.
>
> Valerie
>
>
> On 03/05/2012 01:25 AM, Alex Gutteridge wrote:
>> Hi Valerie,
>>
>> I'd exactly echo what Thomas wrote (thanks Thomas!). My specific use
>> case is mapping coding SNPs to PDB protein structures, so something
>> that encodes 'His88 (CAT) -> Gln88 (CAA)' is all I need.
>>
>> Alex Gutteridge.
>>
>> On 03.03.2012 00:58, Thomas Girke wrote:
>>> Hi Valerie,
>>>
>>> Based on my experience the position in the complete protein (rather
>>> than
>>> CDS) sequence would be the most important piece of information to
>>> record
>>> here. For instance, if a SNP changes the 88th triplet CAT (His) in the
>>> ORF of a transcript to CAA (Gln) then you want to record it like this:
>>> His88 (CAT) -> Gln88 (CAA). This way the user can map this change to a
>>> protein structure or inject it into the corresponding protein sequence
>>> without any additional remapping or coding.
>>>
>>> Another feature to consider are SNPs affecting splice sites (commonly
>>> first and last two nucleotides of an intron).
>>>
>>> If possible it would be also very useful to support users who want to
>>> work with their custom genomes and annotations files provided as FASTA
>>> and GFF/GTF files, respectively.
>>>
>>> Best,
>>>
>>> Thomas
>>>
>>>
>>> On Fri, Mar 02, 2012 at 11:31:57PM +0000, Valerie Obenchain wrote:
>>>> Slight change to this -
>>>>
>>>> I'm now returning the following new columns,
>>>>
>>>> \item{\code{seqTxLoc}}{
>>>> Location in transcript-based coordinates of the first
>>>> nucleotide in
>>>> the codon sequence to be translated. This position
>>>> corresponds to the
>>>> first nucleotide in both the \code{refSeq} and \code{varSeq}
>>>> columns.
>>>> }
>>>> \item{\code{varTxLoc}}{
>>>> Location in transcript-based coordinates of the first
>>>> nucleotide in
>>>> the variant. This value will be the same as \code{seqTxLoc}
>>>> when the
>>>> variant starts exactly at the beginning of a codon.
>>>> }
>>>> \item{\code{varCdsLoc}}{
>>>> Location in cds-based coordinates of the first nucleotide in
>>>> the variant. This position is relative to the start of the
>>>> cds region
>>>> defined in the \code{subject} annotation.
>>>> }
>>>> \item{\code{subjStrand}}{
>>>> The strand of the \code{subject} the variant matched.
>>>> \code{predictCoding}
>>>> determines which variants fall in a coding region by finding
>>>> the
>>>> overlaps
>>>> between the \code{query} and \code{subject}. The
>>>> \code{query} may be
>>>> un-stranded \sQuote{*} but the \code{subject} annotation will
>>>> have a strand.
>>>> }
>>>>
>>>>
>>>> You are interested in 'protein coordinates'. Does 'varCdsLoc'
>>>> described
>>>> above meet the need or are you looking for the actual codon number in
>>>> the coding sequence? I am interested in hearing more about what you
>>>> are
>>>> doing with the protein coordinates, how you are using them. It would
>>>> help us better design future functions.
>>>>
>>>> Thanks,
>>>> Valerie
>>>>
>>>> On 03/02/2012 01:11 AM, Alex Gutteridge wrote:
>>>> > Thanks Valerie - much appreciated!
>>>> >
>>>> > On 01.03.2012 21:30, Valerie Obenchain wrote:
>>>> >> A 'txLoc' column has been added to the output of predictCoding.
>>>> >> Available in devel version 1.1.57.
>>>> >>
>>>> >> Valerie
>>>> >>
>>>> >>
>>>> >> On 02/28/2012 08:20 AM, Valerie Obenchain wrote:
>>>> >>> Good suggestion. Yes, predictCoding is does this internally. I'll
>>>> >>> post back here when this has been added.
>>>> >>>
>>>> >>> Valerie
>>>> >>>
>>>> >>>
>>>> >>>
>>>> >>> On 02/28/2012 01:49 AM, Alex Gutteridge wrote:
>>>> >>>> Hi Valerie,
>>>> >>>>
>>>> >>>> Thanks everything works great now. One small feature request -
>>>> >>>> would it be hard to output the protein sequence position of the
>>>> >>>> coding SNPs? At the moment once I've run predictCoding I'm
>>>> >>>> re-extracting the cds and working out the position of each coding
>>>> >>>> SNP so I can see where in the protein sequence it is, but it
>>>> seems
>>>> >>>> like this is probably just replicating what predictCoding must be
>>>> >>>> doing internally anyway?
>>>> >>>>
>>>> >>>> Alex Gutteridge
>>>> >>>
>>>> >>>
>>>> >>> On 02/24/2012 10:39 AM, Valerie Obenchain wrote:
>>>> >>>> Hi Alex,
>>>> >>>>
>>>> >>>> Thanks for the bug report. The cdsID was taken from an overlap
>>>> >>>> between the query and GRangesList of cds by transcripts. This
>>>> gave
>>>> >>>> the correct transcript number but (incorrectly) took the first
>>>> cds
>>>> >>>> number in the list by default. Now fixed in devel 1.1.55.
>>>> >>>>
>>>> >>>> I've also updated the man page.
>>>> >>>>
>>>> >>>> Valerie
>>>> >>>>
>>>> >>>>
>>>> >>>>
>>>> >>>> On 02/24/2012 02:08 AM, Alex Gutteridge wrote:
>>>> >>>>> On 22.02.2012 18:58, Herv? Pag?s wrote:
>>>> >>>>>> Hi Alex,
>>>> >>>>>>
>>>> >>>>>> On 02/22/2012 03:56 AM, Alex Gutteridge wrote:
>>>> >>>>>
>>>> >>>>> [...]
>>>> >>>>>
>>>> >>>>>>> But the predictCoding call gives this error:
>>>> >>>>>>>
>>>> >>>>>>> Error in .setSeqNames(x, value) :
>>>> >>>>>>> The replacement value for isActiveSeq must be a logical
>>>> vector,
>>>> >>>>>>> with
>>>> >>>>>>> names that match the seqlevels of the object
>>>> >>>>>>
>>>> >>>>>> The error message doesn't help much but I think the pb is
>>>> that you
>>>> >>>>>> didn't rename chMT properly. Try to do this:
>>>> >>>>>>
>>>> >>>>>> seqlevels(snps) <- gsub("chrMT", "chrM", seqlevels(snps))
>>>> >>>>>>
>>>> >>>>>> before you start the for(eg in entrez.ids){..} loop again.
>>>> >>>>>>
>>>> >>>>>> Cheers,
>>>> >>>>>> H.
>>>> >>>>>
>>>> >>>>> Thanks Herv? that nailed it. I'm having some difficulty
>>>> joining up
>>>> >>>>> the output of predictCoding() with the query SNPs though. If
>>>> >>>>> someone could point out where the disconnect in my thinking is I
>>>> >>>>> would appreciate it!
>>>> >>>>>
>>>> >>>>> Here's my (now edited down) script:
>>>> >>>>>
>>>> >>>>> library(BSgenome.Hsapiens.UCSC.hg19)
>>>> >>>>> library(VariantAnnotation)
>>>> >>>>> library(SNPlocs.Hsapiens.dbSNP.20110815)
>>>> >>>>> library(TxDb.Hsapiens.UCSC.hg19.knownGene)
>>>> >>>>>
>>>> >>>>> entrez.ids = c('6335')
>>>> >>>>> txdb19 = TxDb.Hsapiens.UCSC.hg19.knownGene
>>>> >>>>>
>>>> >>>>> snps = getSNPlocs(c("ch1","ch2"),as.GRanges=T)
>>>> >>>>> seqlevels(snps) <- gsub("ch", "chr", seqlevels(snps))
>>>> >>>>> seqlevels(snps) <- gsub("chrMT", "chrM", seqlevels(snps))
>>>> >>>>>
>>>> >>>>> gene.list = cdsBy(txdb19, by="gene")
>>>> >>>>> vsd.list = gene.list[entrez.ids]
>>>> >>>>> cds.list = cdsBy(txdb19,by="tx")
>>>> >>>>>
>>>> >>>>> eg = entrez.ids[1]
>>>> >>>>>
>>>> >>>>> snp.idx = unique(queryHits(findOverlaps(snps, vsd.list[[eg]])))
>>>> >>>>> eg.snps = snps[snp.idx]
>>>> >>>>> iupac = values(eg.snps)[,"alleles_as_ambig"]
>>>> >>>>> eg.snps.exp = rep(eg.snps, nchar(IUPAC_CODE_MAP[iupac]))
>>>> >>>>> variant.alleles =
>>>> >>>>>
>>>> DNAStringSet(strsplit(paste(IUPAC_CODE_MAP[iupac],collapse=""),"")[[1]])
>>>>
>>>> >>>>>
>>>> >>>>>
>>>> >>>>>
>>>> >>>>> aa =
>>>> >>>>>
>>>> predictCoding(eg.snps.exp,txdb19,seqSource=Hsapiens,varAllele=variant.alleles)
>>>> >>>>>
>>>> >>>>> #####
>>>> >>>>>
>>>> >>>>> Then if I query the predictCoding results in aa in an
>>>> interactive
>>>> >>>>> session I get the following (see inline comments for what I
>>>> think
>>>> >>>>> should be happening, but I must be misinterpreting what queryID
>>>> >>>>> means)
>>>> >>>>>
>>>> >>>>> The docs for predictCoding() contain a small typo
>>>> >>>>> (s/queryHits/queryID), but otherwise seem clear?
>>>> >>>>>
>>>> >>>>> Columns include ?queryID?, ?consequence?, ?refSeq?, ?varSeq?,
>>>> >>>>> ?refAA?, ?varAA?, ?txID?, ?geneID?, and ?cdsID?.
>>>> >>>>>
>>>> >>>>> ?queryHits? The ?queryHits? column provides a map back
>>>> to the
>>>> >>>>> variants in the original ?query?. If the ?query? was a
>>>> >>>>> ?VCF?
>>>> >>>>> object this index corresponds to the row in the
>>>> >>>>> ?GRanges? in
>>>> >>>>> the ?rowData? slot. If ?query? was an expanded
>>>> ?GRanges?,
>>>> >>>>> ?RangedData? or ?RangesList? the index corresponds to
>>>> >>>>> the row
>>>> >>>>> in the expanded object.
>>>> >>>>>
>>>> >>>>> #####
>>>> >>>>>
>>>> >>>>>> aa[1,]
>>>> >>>>> DataFrame with 1 row and 9 columns
>>>> >>>>> queryID consequence refSeq
>>>> varSeq refAA
>>>> >>>>> <integer> <factor> <DNAStringSet> <DNAStringSet> <AAStringSet>
>>>> >>>>> 1 1 nonsynonymous CTC
>>>> ATC L
>>>> >>>>> varAA txID geneID cdsID
>>>> >>>>> <AAStringSet> <character> <factor> <integer>
>>>> >>>>> 1 I 10921 6335 33668
>>>> >>>>>> #So the first SNP (queryID: 1) is nonsynonymous and maps to tx
>>>> >>>>>> '10921' and cds '33668'.
>>>> >>>>>> #If I look at the first query SNP I get this:
>>>> >>>>>> eg.snps.exp[aa[1,'queryID'],]
>>>> >>>>> GRanges with 1 range and 2 elementMetadata values:
>>>> >>>>> seqnames ranges strand | RefSNP_id
>>>> >>>>> alleles_as_ambig
>>>> >>>>> <Rle> <IRanges> <Rle> | <character> <character>
>>>> >>>>> [1] chr2 [167055370, 167055370] * | 111558968
>>>> >>>>> R
>>>> >>>>> ---
>>>> >>>>> seqlengths:
>>>> >>>>> chr1 chr2 chr3 chr4 chr5 chr6 ... chr20 chr21 chr22
>>>> chrX
>>>> >>>>> chrY chrM
>>>> >>>>> NA NA NA NA NA NA ... NA NA
>>>> NA NA
>>>> >>>>> NA NA
>>>> >>>>>> #So SNP 1 is at 167055370 on chr2
>>>> >>>>>> #But if I check tx '10921' I see that the cds overlapping
>>>> >>>>>> 167055370 is actually '33651'
>>>> >>>>>> #And cds '33668' is at the other end of the tx:
>>>> >>>>>> cds.list[[aa[1,'txID']]]
>>>> >>>>> GRanges with 26 ranges and 3 elementMetadata values:
>>>> >>>>> seqnames ranges strand | cds_id
>>>> >>>>> cds_name
>>>> >>>>> <Rle> <IRanges> <Rle> | <integer> <character>
>>>> >>>>> [1] chr2 [167168009, 167168266] - | 33668 <NA>
>>>> >>>>> [2] chr2 [167163466, 167163584] - | 33667 <NA>
>>>> >>>>> [3] chr2 [167163020, 167163109] - | 33666 <NA>
>>>> >>>>> [4] chr2 [167162302, 167162430] - | 33647 <NA>
>>>> >>>>> [5] chr2 [167160748, 167160839] - | 33646 <NA>
>>>> >>>>> [6] chr2 [167159600, 167159812] - | 33645 <NA>
>>>> >>>>> [7] chr2 [167151109, 167151172] - | 33644 <NA>
>>>> >>>>> [8] chr2 [167149741, 167149882] - | 33643 <NA>
>>>> >>>>> [9] chr2 [167144947, 167145153] - | 33642 <NA>
>>>> >>>>> ... ... ... ... ... ...
>>>> >>>>> ...
>>>> >>>>> [18] chr2 [167099012, 167099166] - | 33659 <NA>
>>>> >>>>> [19] chr2 [167094604, 167094777] - | 33658 <NA>
>>>> >>>>> [20] chr2 [167089850, 167089972] - | 33657 <NA>
>>>> >>>>> [21] chr2 [167085201, 167085482] - | 33656 <NA>
>>>> >>>>> [22] chr2 [167084180, 167084233] - | 33655 <NA>
>>>> >>>>> [23] chr2 [167083077, 167083214] - | 33654 <NA>
>>>> >>>>> [24] chr2 [167060870, 167060974] - | 33653 <NA>
>>>> >>>>> [25] chr2 [167060465, 167060735] - | 33652 <NA>
>>>> >>>>> [26] chr2 [167055182, 167056374] - | 33651 <NA>
>>>> >>>>> exon_rank
>>>> >>>>> <integer>
>>>> >>>>> [1] 2
>>>> >>>>> [2] 3
>>>> >>>>> [3] 4
>>>> >>>>> [4] 5
>>>> >>>>> [5] 6
>>>> >>>>> [6] 7
>>>> >>>>> [7] 8
>>>> >>>>> [8] 9
>>>> >>>>> [9] 10
>>>> >>>>> ... ...
>>>> >>>>> [18] 19
>>>> >>>>> [19] 20
>>>> >>>>> [20] 21
>>>> >>>>> [21] 22
>>>> >>>>> [22] 23
>>>> >>>>> [23] 24
>>>> >>>>> [24] 25
>>>> >>>>> [25] 26
>>>> >>>>> [26] 27
>>>> >>>>> ---
>>>> >>>>> seqlengths:
>>>> >>>>> chr1 chr2 ...
>>>> >>>>> chr18_gl000207_random
>>>> >>>>> 249250621 243199373 ...
>>>> >>>>> 4262
>>>> >>>>>
>>>> >>>>>
>>>> >>>>
>>>> >>>> _______________________________________________
>>>> >>>> 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
>>>> >>>
>>>> >>> _______________________________________________
>>>> >>> 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
>>>> >
>>>>
>>>> _______________________________________________
>>>> 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
>>
>
More information about the Bioconductor
mailing list