[Bioc-sig-seq] complement and reverseComplement from a character vector

[Hervé Pagès]

Hi Nicolas,

Added to Biostrings 2.13.21.

Note that in order to keep it simple I've modified it so the input must
be the sizes of the contigs, not the contigs themselves. That way it
remains independent on how the contigs are stored and there is no need
to make it a generic:

   N50 <- function(csizes)
   {
     sorted_csizes <- sort(csizes)
     tmp <- cumsum(sorted_csizes)
     N50 <- sorted_csizes[which(tmp >= max(tmp)/2)[1]]
     return(N50)
   }

Can you please send me the documentation (off-list) with the exact definition
of the N50 value (or a reference to it) and an example?

Thanks!
H.


Nicolas Delhomme wrote:
> Dear Hervé,
> 
> Sorry for answering so late. Thanks for integrating these methods. I 
> agree that making them endomorphic is probably the best and most 
> intuitive solution.
> 
> I've been playing with Biostrings lately on "de-novo" assembly data and 
> one important value that one wants to calculate is the so-called "N50" 
> value. It is the size of the contig, which addition to the cumulative 
> contig size (ordered by decreasing sizes) makes it bigger than half the 
> total size of all contigs. Basically it gives an idea of the contig size 
> distribution within the data (although, one still needs to know the 
> contig size range). I could not find any function doing this, so I have 
> written a function to calculate it and I think it could be an 
> interesting addition to the Biostrings package. Here is the code:
> 
> setGeneric("N50",function(x) standardGeneric("N50"))
> setMethod("N50","XStringSet",function(x){
>   # the order
>   x.order<-order(width(x),decreasing=TRUE)
>   width(x[x.order])[which(cumsum(width(x[x.order]))>=sum(width(x)/2))[[1]]]
>   })
> 
> An example (using the non endomorphic functions...):
> 
> my.contig<-DNAStringSet(
>     sapply(
>         sample(c(100:10000),10),
>         function(size){
>             paste(sample(c("A","C","G","T"),size,replace=TRUE),collapse="")
>     })
> )
> N50(my.contig)
> 
> 
> If you find it valuable, I could contribute the documentation.
> 
> Best wishes,
> 
> ---------------------------------------------------------------
> Nicolas Delhomme
> 
> High Throughput Functional Genomics Center
> 
> European Molecular Biology Laboratory
> 
> Tel: +49 6221 387 8426
> Email: nicolas.delhomme at embl.de
> Meyerhofstrasse 1 - Postfach 10.2209
> 69102 Heidelberg, Germany
> ---------------------------------------------------------------
> 
> P.S.
> 
>  > sessionInfo()
> R version 2.9.0 (2009-04-17)
> i386-apple-darwin8.11.1
> 
> locale:
> en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8
> 
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
> 
> other attached packages:
> [1] Biostrings_2.12.0 IRanges_1.2.0     Biobase_2.4.0
> 
> loaded via a namespace (and not attached):
> [1] tools_2.9.0
> 
> 
> 
> On 2 Jun 2009, at 20:36, Hervé Pagès wrote:
> 
>> Hi Nicolas,
>>
>> Nicolas Delhomme wrote:
>>> Hi all,
>>> As described in the doc for reverse, complement and 
>>> reverseComplement, x can be a character vector.
>>> When I tried, it turned out that these functions are not implemented 
>>> for complement and reverseComplement.
>>> There are of some use for me, so I just wrote them up:
>>> setMethod("complement", "character",
>>>    function(x, ...){
>>>            if(length(x)==1) x<-DNAString(x)
>>>        else x<-DNAStringSet(x)
>>>        complement(x)
>>>    }
>>> )
>>> setMethod("reverseComplement", "character",
>>>    function(x, ...){
>>>            if(length(x)==1) x<-DNAString(x)
>>>            else x<- DNAStringSet(x)
>>>            reverseComplement(x)
>>>    }
>>> )
>>> I just post them in case there would be of use for someone else. I 
>>> recognize that it does not save much compared to first converting the 
>>> character vector into a DNAString or DNAStringSet. At least, for me, 
>>> it allows to skip some "input" evaluation test checking whether I 
>>> start with a character vector or a DNAString.
>>
>> Thanks for the feedback!
>>
>> The reason these method were not defined is that when the input is 
>> character
>> it's not clear whether it should be treated as DNA or RNA input. However
>> choosing to treat it as DNA is probably what the user will want 99% of 
>> the
>> time so they could indeed be implemented as in your code above. So if the
>> input contains the "U" letter, they will simply fail (without trying to
>> be smart).
>>
>> Note that there is no method for BString/BStringSet objects for 
>> exactly the
>> same reason.
>>
>> A question subsists though: should these methods return a 
>> DNAString/DNAStringSet
>> object or should the result be turned back into an ordinary character 
>> vector
>> before it's returned to the user? The latter would make these methods
>> "endomorphisms" (i.e. the output has the same type as the input) which is
>> more consistent with what the other methods do but I'm not against making
>> an exception when the input is character (not a big deal as long as 
>> this is
>> clearly documented). Then if the user really want this result to be 
>> character
>> then s/he can always apply as.character() to it.
>>
>> If there are no objections, I will add these methods to the Biostrings 
>> package.
>>
>> Cheers,
>> H.
>>
>>> Best,
>>> ---------------------------------------------------------------
>>> Nicolas Delhomme
>>> High Throughput Functional Genomics Center
>>> European Molecular Biology Laboratory
>>> Tel: +49 6221 387 8426
>>> Email: nicolas.delhomme at embl.de
>>> Meyerhofstrasse 1 - Postfach 10.2209
>>> 69102 Heidelberg, Germany
>>> _______________________________________________
>>> Bioc-sig-sequencing mailing list
>>> Bioc-sig-sequencing at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/bioc-sig-sequencing
>>
>> -- 
>> Hervé Pagès
>>
>> Program in Computational Biology
>> Division of Public Health Sciences
>> Fred Hutchinson Cancer Research Center
>> 1100 Fairview Ave. N, M2-B876
>> P.O. Box 19024
>> Seattle, WA 98109-1024
>>
>> E-mail: hpages at fhcrc.org
>> Phone:  (206) 667-5791
>> Fax:    (206) 667-1319
> 

-- 
Hervé Pagès

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M2-B876
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: hpages at fhcrc.org
Phone:  (206) 667-5791
Fax:    (206) 667-1319