#' saturation analyses for CDR3 sequences
#'
#' By 1 2 3 4 5 6 7 8 9 10 20 30 40 50 60 70 80 90 percent
#' Or user defined the percent and times, eg: 10% for 10 times.
#' To counts the reads/distinct cdr3aa number/QC/diversity
#'
#' @param IRsample cdr3aa column is needed, other(V,D,J,cdr3nt,type) are optional.
#' @param number range from 1~99
#' @param times how many times should be sampled.
#' @return saturation_matrix 18 rows and 8 columns
#' @import ggplot2
#' @export
#' @keywords plot_cdr3_stat
#' @examples
#' #' plot_cdr3_stat(stat_matrix,'cdr3aa-distribution-stat.pdf')
stat_cdr3aa_saturation <- function(IRsample,by,times){
require(ggplot2);
all_reads=nrow(IRsample)
all_cdr3=length(unique(IRsample$cdr3aa))
base_reads=floor(all_reads/100)
if(! hasArg(by)){
len=lapply((c(1:9,1:9*10)), function(i){
sub_IRsample=IRsample[sample(all_reads, base_reads*i),]
tmp=sub_IRsample[,4]
tmp=tmp[tmp!='N/A']
return(c(
length(unique(sub_IRsample$cdr3aa)),
QC_IR(sub_IRsample),
stat_cdr3_diversity(tmp,'d50'),
stat_cdr3_diversity(tmp,'shannon'),
stat_cdr3_diversity(tmp,'simpson'),
stat_cdr3_diversity(tmp,'gini')
))
})
saturation=do.call(rbind,len)
saturation=cbind(paste0(c(1:9,1:9*10),"%"),saturation)
#saturation=rbind(c('100%',all_cdr3,QC_IR(IRsample),,),saturation)
colnames(saturation)=c('percent','unique cdr3aa number','all reads','not qualified','not cdr3 aa','d50','shannon','simpson','gini')
return(saturation)
}else{
len=lapply(rep(1,times), function(i){
sub_IRsample=IRsample[sample(all_reads, base_reads*by ),]
tmp=sub_IRsample[,4]
tmp=tmp[tmp!='N/A']
return(c(
length(unique(sub_IRsample$cdr3aa)),
QC_IR(sub_IRsample),
stat_cdr3_diversity(tmp,'d50'),
stat_cdr3_diversity(tmp,'shannon'),
stat_cdr3_diversity(tmp,'simpson'),
stat_cdr3_diversity(tmp,'gini')
))
})
saturation=do.call(rbind,len)
saturation=cbind(paste0( rep(by,times) ,"%"),saturation)
#saturation=rbind(c('100%',all_cdr3,QC_IR(IRsample),,),saturation)
colnames(saturation)=c('percent','unique cdr3aa number','all reads','not qualified','not cdr3 aa','d50','shannon','simpson','gini')
return(saturation)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.