#' Annotating pileup
#'
#' This function is used to get pileup data based on new transcript for the
#' purpose of the downstream outlier detection analysis.
#' @param pileup raw coverage data
#' @param exon exon
#' @param inputType type of intronic region contained in pileup, with choices
#' "whole_intron", "part_intron", or "only_exon"; the first is the default.
#' @param outputType type of intronic region that will be included in output,
#' with choices "whole_intron", "part_intron", or "only_exon"; the second is
#' the default.
#' @keywords
#' @export
#' @examples
annotate_pileup = function(pileup,exon,
inputType="whole_intron",
outputType="part_intron") {
# inputType = "whole_intron", "part_intron", "only_exon"
# output.ytpe = "whole_intron", "part_intron", "only_exon"
if (inputType=="whole_intron") {
if (outputType=="whole_intron") {
ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=NULL) ;
output = pileup;
} else if (outputType=="part_intron") {
intron.len = ceiling(len.intron.hy(exon=exon)*0.5);
ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
output = pileup[ep.new$coverage.col,] ; # Area to be included.
} else if (outputType=="only_exon") {
intron.len = 0;
ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
output = pileup[ep.new$coverage.col,] ; # Area to be included.
} else {
stop(outputType," is not an option for outputType.")
}
} else if (inputType=="part_intron") {
if (outputType=="whole_intron") {
stop(outputType," is not an option when inputType=part_intron.")
} else if (outputType=="part_intron") {
intron.len = ceiling(len.intron.hy(exon=exon)*0.5);
ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
output = pileup ; # Area to be included.
##### If the dimension of the given pileup data is not equal to the expected one
##### from the intron.len calculated, should display an error. Fix this!
} else if (outputType=="only_exon") {
intron.len.temp = ceiling(len.intron.hy(exon=exon)*0.5);
ep.new.temp = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len.temp) ;
exonic.region = c()
for (i in 1:nrow(ep.new.temp$ep)) {
exonic.region = c(exonic.region,(ep.new.temp$epl[i]:ep.new.temp$epr[i]));
}
output = pileup[exonic.region,];
rm(intron.len.temp,ep.new.temp,exonic.region);
intron.len = 0;
ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
} else {
stop(outputType," is not an option for outputType.")
}
} else if (inputType=="only_exon") {
if (outputType=="whole_intron") {
stop(outputType," is not an option when inputType==only_exon.")
} else if (outputType=="part_intron") {
stop(outputType," is not an option when inputType==only_exon.")
} else if (outputType=="only_exon") {
intron.len = 0;
ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
output = pileup ;
} else {
stop(outputType," is not an option for outputType.")
}
} else {
stop(inputType," is not an option for inputType.")
}
# Save exon position as a matrix for the new pileup
epm0 = ep.new$ep; epm = ep.new$ep;
# Change the directions
seqdir = strsplit(exon,":")[[1]][3] ;
if (seqdir=="-") {
output = apply(output,2,rev); # log counts
seqlen = epm0[nrow(epm0),4]
epm[,1] = seqlen - rev(epm0[,4]) + 1;
epm[,2] = seqlen - rev(epm0[,3]) + 1;
epm[,3] = seqlen - rev(epm0[,2]) + 1;
epm[,4] = seqlen - rev(epm0[,1]) + 1;
}
# pileup annotation info (dai)
dai = as.list(NULL)
dai$epm = epm
dai$intron.len = intron.len
return(list(new.pileup=output,dai=dai))
}
##########################################################
# annotate.pileup.v1 = function(pileup,exon,pileup.type="whole_intron",
# outputType="part_intron",intron.len=NULL) {
# # pileup.type = "whole_intron", "part_intron", "only_exon"
# # output.ytpe = "whole_intron", "part_intron", "only_exon"
#
# if (pileup.type=="whole_intron") {
#
# if (outputType=="whole_intron") {
# ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=NULL) ;
# output = pileup;
# } else if (outputType=="part_intron") {
# if (is.null(intron.len)) {
# intron.len = ceiling(len.intron.hy(exon=exon)*0.5);
# }
# ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
# output = pileup[ep.new$coverage.col,] ; # Area to be included.
# } else if (outputType=="only_exon") {
# intron.len = 0;
# ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
# output = pileup[ep.new$coverage.col,] ; # Area to be included.
# } else {
# stop(outputType,"is not an option for outputType.")
# }
#
# } else if (pileup.type=="part_intron") {
#
# if (outputType=="whole_intron") {
# stop(outputType,"is not an option when pileup.type=part_intron.")
# } else if (outputType=="part_intron") {
# if (is.null(intron.len)) {
# intron.len = ceiling(len.intron.hy(exon=exon)*0.5);
# }
# ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
# output = pileup ; # Area to be included.
# ##### If the dimension of the given pileup data is not equal to the expected one
# ##### from the intron.len calculated, should display an error. Fix this!
# } else if (outputType=="only_exon") {
# intron.len.temp = ceiling(len.intron.hy(exon=exon)*0.5);
# ep.new.temp = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len.temp) ;
#
# exonic.region = c()
# for (i in 1:nrow(ep.new.temp$ep)) {
# exonic.region = c(exonic.region,(ep.new.temp$epl[i]:ep.new.temp$epr[i]));
# }
# output = pileup[exonic.region,];
# rm(intron.len.temp,ep.new.temp,exonic.region);
#
# len.intron = 0;
# ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=len.intron) ;
# } else {
# stop(outputType,"is not an option for outputType.")
# }
#
# } else if (pileup.type=="only_exon") {
#
# if (outputType=="whole_intron") {
# stop(outputType,"is not an option when pileup.type==only_exon.")
# } else if (outputType=="part_intron") {
# stop(outputType,"is not an option when pileup.type==only_exon.")
# } else if (outputType=="only_exon") {
# intron.len = 0;
# ep.new = find.exon.hy(exon,is.intron=TRUE,num.intron=intron.len) ;
# output = pileup ;
# } else {
# stop(outputType,"is not an option for outputType.")
# }
#
# } else {
# stop(pileup.type,"is not an option for outputType.")
# }
#
# # Save exon position as a matrix for the new pileup
# epm0 = ep.new$ep; epm = ep.new$ep;
#
# # Change the directions
# seqdir = strsplit(exon,":")[[1]][3] ;
# if (seqdir=="-") {
# output = apply(output,2,rev); # log counts
# seqlen = epm0[nrow(epm0),4]
# epm[,1] = seqlen - rev(epm0[,4]) + 1;
# epm[,2] = seqlen - rev(epm0[,3]) + 1;
# epm[,3] = seqlen - rev(epm0[,2]) + 1;
# epm[,4] = seqlen - rev(epm0[,1]) + 1;
# }
#
# # pileup annotation info (dai)
# dai = as.list(NULL)
# dai$epm = epm
# dai$intron.len = intron.len
# return(list(out.pileup=output,dai=dai))
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.