Commit 1bfd6fd2 by Eric Coissac

Add a Kit-Ligand dataset, and correct the genotyping function

parent 2fea466f
......@@ -178,8 +178,11 @@ genotyping = function(sequencing,sequencing.error,phetero) {
probas = matrix(0.0,nrow = nrow(sequencing),
ncol = nrow(pnucs))
ACGT.cols = match(c('A','C','G','T'),names(sequencing))
for (i in seq_len(nrow(sequencing))) {
obs = as.numeric(sequencing[i,4:7])
obs = as.numeric(sequencing[i,ACGT.cols])
for (j in seq_len(nrow(pnucs))) {
theo= as.numeric(pnucs[j,1:4])
probas[i,j] = dmultinomialACGT(obs[1],obs[2],obs[3],obs[4],
......@@ -206,11 +209,13 @@ genotyping = function(sequencing,sequencing.error,phetero) {
best.genotype = apply(probas, MARGIN = 1, which.max)
second.genotype = apply(probas, MARGIN = 1, function(x) which(x==second(x))[1])
g.good = rownames(pnucs)[best.genotype]
haplotype1 = substr(as.character(g.good),1,1)
haplotype2 = substr(as.character(g.good),2,2)
g.alt = rownames(pnucs)[second.genotype]
lod = apply(probas, MARGIN = 1, max)
......@@ -218,7 +223,7 @@ genotyping = function(sequencing,sequencing.error,phetero) {
true = apply(sequencing[,1:2],MARGIN = 1,function(x) paste0(sort(x),collapse = ""))
sequencing = cbind(sequencing,probas,g.good,g.alt,lod,score = exp(score),ok=true==g.good)
sequencing = cbind(sequencing,probas,g.good,haplotype1,haplotype2,g.alt,lod,score = exp(score),ok=true==g.good)
return(sequencing)
}
This source diff could not be displayed because it is too large. You can view the blob instead.
kitlg.gene = read.table("RawData/ROO-D6-5104.3_124450000-124600000.bases.txt",
header=FALSE)
names(kitlg.gene)=c("position","A","C","G","T")
usethis::use_data(kitlg.gene,overwrite = TRUE)
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment