From 00a2a30c1ed9a6db205f0d6a59caab9616977a97 Mon Sep 17 00:00:00 2001 From: shackett Date: Sun, 16 Sep 2012 17:56:53 -0400 Subject: [PATCH 1/6] temporary mod to allow running without pipeline --- fit-hmm.R | 41 ++++++++++++++++++++++++++++++++++++----- 1 file changed, 36 insertions(+), 5 deletions(-) diff --git a/fit-hmm.R b/fit-hmm.R index f966dde..08a93f6 100755 --- a/fit-hmm.R +++ b/fit-hmm.R @@ -3,24 +3,55 @@ ## options(warn=2,error=recover); options(error=quote(q("yes"))) -args <- commandArgs() -dollar0 <- substring(args[grep("^--file=", args)], 8) +### manual param filling +setwd("/Users/Sean/Desktop/Andolfatto/9_12_toy") + +dollar0 <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/msg/blah" +toy_indivs <- c("indivA12_AATAAG", "indivE2_CAGCCG", "indivE4_TAGGAG") +direct <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/hmm_data" +outdir <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/hmm_fit" + +args <- + +-d $outdir -i $indiv -s $sex -o $Routdir -p $deltapar1 -q $deltapar2 -r $rfac -c $chroms -x $sexchroms -y $chroms2plot -z $priors -t $theta" + +opts <- list(s = "male", i = toy_indivs[1], d = direct, o = outdir, p = 0.1, q = 0.1, r = 0.000001, c = "2L,2R,3L,3R,4", x = "X", y = "all", z = "0,.5,.5", t = 1, H = "both", e = 0.8, j = 30, l = 20) + + + + + + + +#args <- commandArgs() +#dollar0 <- substring(args[grep("^--file=", args)], 8) source(sprintf("%s/ded.R", dirname(dollar0))) source(sprintf("%s/hmmlib.R", dirname(dollar0))) library("R.methodsS3", lib.loc = dirname(dollar0)) library("R.oo", lib.loc = dirname(dollar0)) -opts <- getopts() + +#opts <- getopts() indivs <- unlist(strsplit(opts$i,split=",")) sex <- opts$s dir <- opts$d outdir <- opts$o deltapar1 <- as.numeric(opts$p) deltapar2 <- as.numeric(opts$q) -rfac <- as.numeric(opts$r) + priors <- unlist(strsplit(opts$z,split=",")) theta <- as.numeric(opts$t) +HMMtype <- opts$H + +if(HMMtype %in% c("ML", "both")){ + rfac <- as.numeric(opts$r) + } +if(HMMtype %in% c("viterbi", "both")){ + HMM_seqpair <- as.numeric(opts$j) + HMM_diffthresh <- as.numeric(opts$l) + HMM_decay <- as.numeric(opts$e) + } stopifnot(!is.null(indivs), !is.null(dir), !is.null(outdir), length(indivs) == 1) @@ -42,7 +73,7 @@ plotPadding <- 10^(ceiling(log10(aveSpace))-2) alleles <- c("A","C","G","T") - +indiv <- indivs for(indiv in indivs) { cat(indiv, "\n") ## if(opts$c == "all") From fe151a8479eaffb723971f54c70a814244694001 Mon Sep 17 00:00:00 2001 From: shackett Date: Mon, 17 Sep 2012 12:12:08 -0400 Subject: [PATCH 2/6] added the ability to use either or both of the ML or viterbi HMMs added viterbi paramters passing added the ability to run to either or both HMMs changes plotting to make a haplotype plot for each HMM and to label breakpoints and mismatch frequency for each HMM --- .DS_Store | Bin 0 -> 15364 bytes fit-hmm.R | 115 ++++++++++++----- hmmlib.R | 274 ++++++++++++++++++++++++++++++++++++++++- instructions/.DS_Store | Bin 0 -> 6148 bytes 4 files changed, 353 insertions(+), 36 deletions(-) create mode 100644 .DS_Store create mode 100644 instructions/.DS_Store diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..4f34841f59480dbe274b7907564198f7ad96f074 GIT binary patch literal 15364 zcmeHNU2Ggz6~4!I?Adr?j~&~MH;ueWL}REOoHS|T6dM1gO_M0De`==<-OTPx;!$UJ zHaoMngGuWaB#KZ)i+G3-;zxi4Eg+Df1&OCVKty?>0a1i_03q?zBKlHZfbZOyb!TSR zt$_;KbVi!FckVss+;h(T&i$JaLeN;nI3Ps75JJJtsivYBFPT=ngf6b)(XCMVMkGtiM~z*-j2M-_Mn7%k_qaenSn7OH>5j8SD>o_*@=L70_ z%kUPnTJP%WfuqBR2Qq_)hOZ7}R_{D=WFRw){MD4Aq_oGB?mhcIa__8u_0h-H9*0i2 zh&AVHo1d&dKPk2T!P{xQPM>5F46)L)p#u+ILUTPYq2+}<_5rdi6Qm0ptsPThgQqbz?D1pIxUM6Xk7`loRXZZ{UXjN3TSOQ+Lribj6(@1 zZM2xHrDxkBEo=?wOVG>#2YVzRZG&J>ghz$!qVz{uRHSz;@S-HjC#_itS+gNuZuz}? zZh^ku(5@IJ+tGu%cmyaI1nL+OCj3{DK6k}p%zNCkVecvH+&-GjD8202FnyP#FUgjQ z;*vyd1--=45cHQsZyf#H+hIA)kmG`%@ft(A?BTzJdU#tvI}Y0_5+RKEdYo)sm@^e{ zR&iT5DkB`-0mmHlGT>K^t|_D(=S946ocr+|(^}#=+;#^fyvcU}SS};V-;BB1zCoNe zEQIXV0tqhK&O%lZ8fBrMj`DhO4(lWQ+avsAVmz=&haO{MUX0fBA@_JhiV7(!m_rk%c?W;HSGxXVu&&{*S+m&_d9S#pI%)(c*pP% zpV7bdAw8ebKf$XV49vl9rRI?SQ>YA?5}V$*2Eze>5_boJqsAQp<;~!Zz?7IfTE(om zALD3st$4W}(g~P~xJH#Lg6e!WjmW#WQQSrN6fp<)^P0k+^9s?CmLsa^$q*u2#S%se zo^7jbZ7@$oxLUx)V}a|hAfjzZTf%KmjSO@c9bBpygS zka!^CfvtGkj&%;N?(&#&cwLpnsKo1E9&h}ecD)e|uRbN9(0|uH!g)`HLc>B(Y#ztY*8m;4edutO4jV^O{X8?n(?xUr=cja zc=w*PI_t6udO&#FK2@E41pjPkZpvr9eLzjm*?P{jH6!J-4INO^^Gm+_v{OOh;Vw13 zkar3NDCY)c!^hI;hq>9)y6c@SnMJ5{9c|xze55non;zfDL`6BrT%(2oRrgi{aqmdK zpWS=>uYM-AcVsB2Rm3iAT_?n8aRK4ADxMP8#B<_@;(76ccu~A0UKX#2KZ#exYvOO> zb@2ufwUI_0)JZ#OANA7*=^!1UVLC>4(FC2OQ*@fn&;xXy9;Pf^Ad4=OOCD9|QTjAJ zNuQz5)7R+h^bPtZeUF}}pVH6h=X9NZMK9Cu>Cf~kz2uxkUnc`0lZFRK%bV2Bbux4>JkBD!e4PBbdGc zl<&NCm)bW6l$tB?=%_+UYvI^XB=zn+eSMlDk;lxp?ZayWV#BY!eV|X(8f64Lp;?*(LN8E`@?=nf9N@GJq<#W8U8Be83HmHO zMW3Um=^6SWeTlwIU!kwkx9M4Wj=l@rzCb^wpU{i+65XI*#}IoehS)#UA~v=KX_J#9 zANu*nVc}!%I)}WGNotFnRE^qYI&US!wa^qAJBW{N3glF|PVJw2WHHnJr~7~2fXb_% zzyD8e{QsX%l{^yka*8Qd3e^Tl*)oWX^>{1lH=eu_sb zKgA=hE@7BvB%odMQ 1){ + #Run Viterbi-HMM forward and backwards and output the state sequence + dual.hmm <- new.HMM_pathFR(phi, d, prob, zero, HMM_decay, r, Int) + datapos <- data$pos + + #Compare recombination events called by foward and reverse HMM runs and determine a consensus, likely set of recombinations + conc.hap <- new.HMM_path.consensus(dual.hmm, contig, K, data$pos) + + Pr.z.given.y_vit <- t(sapply(1:L, function(x){ifelse(1:K == conc.hap$Ancestry[x], 1, 0)})) + + if(length(unique(conc.hap$Ancestry)) != 1){ + #gradient based on ambiguous breakpoint boundaries + + for(recon in 1:length(conc.hap$recoZ[,1])){ + reco <- conc.hap$recoZ[recon,] + interpol_p <- 1 - (datapos[reco$mark_start:reco$mark_end] - datapos[reco$mark_start])/(datapos[reco$mark_end] - datapos[reco$mark_start]) + + Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$starthap] <- interpol_p + Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$endhap] <- 1-interpol_p + }} + } + colnames(Pr.z.given.y_vit) <- paste("viterbi Pr(", ancestries, "|y)") + vit_bp_calls[[indiv]][[contig]] <- conc.hap$recoZ + } + + if("ML" %in% HMMtype){data <- cbind(data, Pr.z.given.y_ML)} + if("viterbi" %in% HMMtype){data <- cbind(data, Pr.z.given.y_vit)} + + data <- cbind(data, Pr.z.given.y_ML, Pr.z.given.y_vit) attr(data, "badpos") <- badpos dataa[[contig]] <- data @@ -259,16 +302,20 @@ for(indiv in indivs) { plotfile <- file.path(outdir, indiv, paste(indiv, "hmmprob.pdf", sep="-")) if(file.exists(plotfile)) { cat("plot already exists\n") ; next } pdf(file=plotfile, width=7, height=1.5) + par(mar=c(2,2.5,0.5,0.5),bg="transparent",cex.main=.68,cex.lab=.8,font.lab=2,cex.axis=.38,mgp=c(1,.000001,0),xaxs="i") + for(hmmRunning in HMMtype){ + plot(0, 0, xlab="", ylab="", col="transparent", xlim=c(1,sum(as.numeric(contigLengths$length)) + plotPadding*(length(plot.contigs)+1)), ylim=c(-1.01,1.01), axes=F) axis(side=2,at=c(-1,0,1),labels=c("","",""),col="gray38"); mtext(c("par2","par1"),side=2,line=.68,at=c(-1,1),font=2,cex=.8,col=c("blue","red"),las=2); - box(col="gray68"); + mtext(paste("HMM:", hmmRunning, collapse = " "), side = 2, line = -0.8, at = -1.4, font=2, cex=0.6, las=2) + box(col="gray68"); current_start <- plotPadding; - for(contig in plot.contigs) { + for(contig in plot.contigs) { mtext(side=1,at=current_start,contig,font=2,cex=.8,line=1,xpd=T,adj=0) current_end <- current_start + contigLengths[contigLengths$chr == contig,"length"] - 1; @@ -289,31 +336,31 @@ for(indiv in indivs) { if (sum(names(dataa) %in% contig)!=0) { contig_data <- dataa[[contig]]; x <- contig_data$pos - y <- contig_data[,paste("Pr(", ancestries, "|y)")] + y <- contig_data[,paste(hmmRunning, "Pr(", ancestries, "|y)")] ### divvy up homozygous and heterozygous blocks - byBlocks <- breakpoint.width(x, y[,par1homo_col], y[,par2homo_col], indiv=indiv, contig=contig, conf1=.05 ,conf2=.95); + byBlocks <- breakpoint.width(x, y[,par1homo_col], y[,par2homo_col], indiv=indiv, contig=contig, conf1=.05 ,conf2=.95, hmmRunning); if (is.null(byBlocks[["bps"]])==F) { breakpoints <- rbind(breakpoints,byBlocks[["bps"]]); } ### plot like.par1 <- contig_data[contig_data$read_allele==contig_data$par1ref,]$pos; like.par2 <- contig_data[contig_data$read_allele==contig_data$par2ref,]$pos; plot.posterior(x+current_start, y, ancestries, like.par1+current_start, like.par2+current_start, bounds=c(1,contigLengths[contigLengths$chr==contig,]$length)+current_start-1, subtract=current_start, tickwidth=5*10^7) - + ### report mismatch fraction for homozygous regions if (nrow(byBlocks[["blocks"]])>0) { ## plot fraction of par1/(par1+par2) among informative markers (between -1 and 1) - matchMismatch <- rbind(matchMismatch, reportCounts(contig_data, as.vector(byBlocks[["blocks"]][,"V1"]), as.vector(byBlocks[["blocks"]][,"V4"]), as.numeric(as.vector(byBlocks[["blocks"]][,"V7"])), as.numeric(as.vector(byBlocks[["blocks"]][,"V8"])))) + matchMismatch <- rbind(matchMismatch, cbind(reportCounts(contig_data, as.vector(byBlocks[["blocks"]][,"V1"]), as.vector(byBlocks[["blocks"]][,"V4"]), as.numeric(as.vector(byBlocks[["blocks"]][,"V7"])), as.numeric(as.vector(byBlocks[["blocks"]][,"V8"]))), hmmtype = hmmRunning)) } } current_start <- current_start + contigLengths[contigLengths$chr == contig,"length"] + plotPadding; if (contig != plot.contigs[length(plot.contigs)]) { - abline(v=current_start-(plotPadding/2),col="gray68",lwd=1) + abline(v=current_start-(plotPadding/2),col="gray68",lwd=1) } } - + } etc <- "" main <- sprintf("%s (%s): delta=(%.0e, %.0e)", indiv, sex, deltapar1, deltapar2) dev.off() diff --git a/hmmlib.R b/hmmlib.R index cdaff8c..b24692e 100755 --- a/hmmlib.R +++ b/hmmlib.R @@ -238,11 +238,11 @@ findBreak <- function(indiv,contig,y,z1,z2,conf1,conf2,parentage) { list(blocks=as.data.frame(blocks),bps=as.data.frame(widths)); } -breakpoint.width <- function(x, y1, y2, indiv=NA, contig=NA, conf1=.05 ,conf2=.95) { +breakpoint.width <- function(x, y1, y2, indiv=NA, contig=NA, conf1=.05 ,conf2=.95, hmmtype) { width1 <- findBreak(indiv,contig,y1,y1>=conf1,y1<=conf2,conf1,conf2,'homozygous_par1') width2 <- findBreak(indiv,contig,y2,y2>=conf1,y2<=conf2,conf1,conf2,'homozygous_par2') - list(blocks=rbind(width1[["blocks"]],width2[["blocks"]]),bps=rbind(width1[["bps"]],width2[["bps"]])) + list(blocks=cbind(rbind(width1[["blocks"]],width2[["blocks"]]), hmmtype),bps=cbind(rbind(width1[["bps"]],width2[["bps"]]), hmmtype)) } @@ -814,3 +814,273 @@ cleanupReadPileup <- function(x, ref) { x <- mapply(gsub, pattern=list("[\\.,]"), replacement=ref, x=x) return(x); } + +#Write a transition matrix + +TransMat <- function(drI, a, nstates, zerocells){ + tMat <- matrix(data = 0, ncol = nstates*2, nrow = nstates*2) + diag(tMat[1:nstates,1:nstates]) <- 1-drI + diag(tMat[(nstates+1):(2*nstates), 1:nstates]) <- drI + tMat[1:nstates, (nstates+1):(2*nstates)] <- a; diag(tMat[1:nstates, (nstates+1):(2*nstates)]) <- 0 + diag(tMat[(nstates+1):(2*nstates), (nstates+1):(2*nstates)]) <- a + if(!is.null(zerocells)){ + tMat[zerocells[,1],zerocells[,2]] <- 0} + tMat + } + +HMMeval = function(phi, d, prob, zero, a, r, Int){ + +nstates <- length(phi) +npos <- length(d) +phis <- as.numeric(c(phi, rep(0, times = nstates))) +drI <- d*r*Int + +#delta - state-probability matrix +deltad <- matrix(data = NA, ncol = nstates*2, nrow = npos) + +#psi - state matrix +psid <- matrix(data = NA, ncol = nstates*2, nrow = npos) + +#initialization + +deltad[1,] <- log(phis*c(rep(prob[1,], times = 2)), base = 2) +psid[1,] <- rep(0, times = 6) + +#recursion: + +for (i in 2:npos){ + + Ad <- TransMat(drI[i], a, nstates, zero) + betad <- matrix(data = 0, ncol = nstates*2, nrow = nstates*2) + for (j in 1:(nstates*2)){ + betad[j,] <- deltad[i-1,j]+log(Ad[j,], base = 2)+log(c(rep(prob[i,], times = 2)), base = 2) + } + + for (j in 1:nstates){ + if(diag(betad)[j] == max(diag(betad)[1:nstates])){ + deltad[i,j] <- diag(betad)[j] + psid[i,j] <- j + }else{ + deltad[i,j] <- diag(betad[c((nstates+1):(nstates*2)),c(1:nstates)])[j] + psid[i,j] <- j+nstates + } + } + for (j in (nstates+1):(2*nstates)){ + deltad[i,j] <- max(betad[,j]) + psid[i,j] <- c(1:(nstates*2))[betad[,j] == max(betad[,j])][1] + } + } + +#termination: +finald = max(deltad[npos,][1:nstates]) +finpsid = c(1:nstates)[deltad[npos,][1:nstates] == finald] + +#path (state sequence) backtracking +path = rep(NA, times = npos) +path[npos] = finpsid +for (i in (npos-1):1){ + path[i] <- psid[i+1,path[i+1]] + } +path +} + +### remove recombination events that occur close together in f or r hmmpath - <20 markers +### pair forward and reverse HMM-run recombination events ### +#freco <- forwardtrans +#rreco <- reversetrans + +validrecos = function(freco, rreco){ + +dirSort <- rbind(cbind(dir = "F", freco, valid = NA), cbind(dir = "R", rreco, valid = NA)) +dirSort <- dirSort[order(dirSort$bpst),] + +#if F-R ancestries disagree, get rid of the most telomeric recombination events to get the telomeric ancestries to agree + +telomeric <- TRUE +current.trans <- NA + +for(i in 1:length(dirSort[,1])){ +direction = dirSort[i,1] +comp.dir = ifelse(direction == "F", "R", "F") + +if(telomeric == TRUE){ +if(paste(dirSort[i,4:5], collapse = "-") != paste(dirSort[dirSort$dir == comp.dir,][1,4:5], collapse = "-")){ +dirSort[i,6] <- FALSE}else{ + current.trans <- paste(dirSort[i,4:5], collapse = "-") + telomeric <- FALSE}} else if(is.na(current.trans)){ +current.trans <- paste(dirSort[i,4:5], collapse = "-") } else { +if(paste(dirSort[i,4:5], collapse = "-") == current.trans){ + dirSort[(i-1):i,6] <- TRUE + current.trans <- NA} else { + dirSort[(i-1):i,6] <- FALSE + current.trans <- NA + } + } + } + +dirSort <- dirSort[dirSort$valid == TRUE,] + +output <- NULL +if(length(dirSort[,1]) != 0){ +for(i in 1:(length(dirSort[,1])/2)){ + +outz <- cbind(dirSort[(i*2)-1,2], dirSort[(i*2),2], dirSort[(i*2)-1,3], dirSort[(i*2),3:5], mean(c(dirSort[(i*2)-1,2], dirSort[(i*2),2]))) +colnames(outz) <- c("bp_start", "bp_end", "mark_start", "mark_end", "starthap", "endhap", "bp_mean") +output <- rbind(output, outz) +}} +output +} + +#switch the polarity of the markbp, starthap <-> endhap +reverseSwitch = function(reversetrans, lpath){ + revorder = data.frame(marks = rev(reversetrans[,1]), markbp = rev(reversetrans[,2]), paths = rev(reversetrans[,3]), pathf = rev(reversetrans[,4])) + for(i in 1:length(revorder[,1])){ + revorder[i,2] <- lpath - revorder[i,2] + paths <- revorder[i,3] + pathf <- revorder[i,4] + revorder[i,3] <- pathf + revorder[i,4] <- paths + } + revorder + } + +#### recombinations going forward and reverse #### +recoZ = function(path, position){ + pathl = rle(path)$lengths + pathv = rle(path)$values + + cuml = cumsum(pathl) + cumpos = position[cuml] + nbreaks = (length(cumpos)-1)/2 + starts = cumpos[c(1:nbreaks)*2-1] + starthap = pathv[c(1:nbreaks)*2-1] + endhap = pathv[c(1:nbreaks)*2+1] + markbp = cuml[c(1:nbreaks)*2-1] + data.frame(starts = starts, markbp = markbp, starthap = starthap, endhap = endhap) } + +#### define haplotype after filtering possible recombination events + +cons_haplotype <- function(validrecoZ, npos, Pos){ + +for(i in 1:length(validrecoZ[,1])){ + +if(validrecoZ$mark_start[i] == validrecoZ$mark_end[i]){ +if(Pos[validrecoZ$mark_start[i] + 1] - Pos[validrecoZ$mark_start[i]] > Pos[validrecoZ$mark_start[i]] - Pos[validrecoZ$mark_start[i] - 1]){ + validrecoZ$bp_mean[i] <- validrecoZ$bp_mean[i] - 1 + }else{ + validrecoZ$bp_mean[i] <- validrecoZ$bp_mean[i] + 1 + }}} + +for(i in 1:length(validrecoZ[,1])){ + +for(i in 1:length(validrecoZ[,1])){ +validrecoZ[i,8] <- c(1:(npos-1))[mapply(between.val, low = Pos[c(1:(npos-1))], high = Pos[c(2:npos)], validrecoZ[i,7])] + } +validrecoZ[,9] <- validrecoZ[,8]+1 + } + + +if(length(validrecoZ[,1]) == 1){ +track <- c(rep(validrecoZ[i,5], times = validrecoZ[i,8]), rep(validrecoZ[i,6], times = npos - validrecoZ[i,9] + 1)) +}else{ + validrecoZ <- validrecoZ[order(validrecoZ[,8]),] + starts <- sort(validrecoZ[,8]) + track <- NULL + for(i in 1:length(starts)){ + if(i == 1){ + track <- c(track, rep(validrecoZ[i,5], times = validrecoZ[i,8])) + } + else if(i == length(starts)){ + track <- c(track, rep(validrecoZ[i,5], times = validrecoZ[i,8] - validrecoZ[i-1,9]), rep(validrecoZ[i,6], times = npos - validrecoZ[i,8] + 1)) + } + else{ + track <- c(track, rep(validrecoZ[i,5], times = validrecoZ[i,8] - validrecoZ[i-1,9] + 1)) +}}} +track +} + +between.val <- function(test, low, high){ + (high - test) * (test - low) > 0 + } + +new.HMM_pathFR <- function(phi, d, prob, zero, a, r, Int){ + +hmmpath = HMMeval(phi, d, prob, zero, a, r, Int) + +probR <- prob[rev(c(1:length(d))),] +dR <- rev(d) +dR <- c(NA, dR); dR <- dR[-length(dR)] + +hmmpathr = rev(HMMeval(phi, dR, probR, zero, a, r, Int)) + +cbind(hmmpath, hmmpathr) + +} + +new.HMM_path.consensus <- function(dual.hmm, contig, nstates, Pos){ + +#define breakpoint by transition from haplotype 1 state to memory state going forward and by transition from haplotype 2 state to memory state going in reverse. +npos <- length(Pos) +revPos <- rev(Pos) + +if(sum(c(1:nstates) %in% unique(dual.hmm[,1])) > 1){ + forwardtrans <- recoZ(dual.hmm[,1], Pos) + forwardtrans <- cbind("F", forwardtrans) + } else { forwardtrans <- data.frame("F", NA, NA, sort(unique(dual.hmm[,1]))[1], sort(unique(dual.hmm[,1]))[1])} + +names(forwardtrans) <- c("dir", "bpst", "markst", "starthap", "endhap") + +if(sum(c(1:nstates) %in% unique(dual.hmm[,2])) > 1){ + reverset = recoZ(rev(dual.hmm[,2]), revPos) + reversetrans = reverseSwitch(reverset, npos+1) + reversetrans <- cbind("R", reversetrans) + } else { reversetrans <- data.frame("R", NA, NA, sort(unique(dual.hmm[,2]))[1], sort(unique(dual.hmm[,2]))[1])} + +names(reversetrans) <- c("dir", "bpst", "markst", "starthap", "endhap") +possiblerec <- rbind(forwardtrans, reversetrans) +row.names(possiblerec) <- NULL + +freco <- forwardtrans[,-c(1)] +rreco <- reversetrans[,-c(1)] +c.hap <- list() + +if(is.na(freco[1,1]) | is.na(rreco[1,1])){ +ancestry <- ifelse(is.na(freco[1,1]), freco[1,3], rreco[1,3]) +c.hap[["Ancestry"]] <- rep(ancestry, times = npos) +c.hap[["recoZ"]] <- NULL +} else { + +validrecoZ <- validrecos(freco, rreco) + +if(is.null(validrecoZ)){ +ancestry <- as.numeric(names(table(dual.hmm))[table(dual.hmm) == max(table(dual.hmm))]) +c.hap[["Ancestry"]] <- rep(ancestry, times = npos) +c.hap[["recoZ"]] <- NULL + +}else{ + +c.hap[["Ancestry"]] <- cons_haplotype(validrecoZ, npos, Pos) +c.hap[["recoZ"]] <- validrecoZ +}} +c.hap +} + +#define breakpoint regions +recombs = function(hmmpath, parcount){ + pathl = rle(hmmpath)$lengths + pathv = rle(hmmpath)$values + + cuml = cumsum(pathl) + cumpos = parcount[cuml] + nbreaks = (length(cumpos)-1)/2 + starts = cumpos[c(1:nbreaks)*2-1] + markbp = cuml[c(1:nbreaks)*2-1] + marktrans = cuml[c(1:nbreaks)*2]-cuml[c(1:nbreaks)*2-1] + int = cumpos[c(1:nbreaks)*2]-cumpos[c(1:nbreaks)*2-1] + starthap = pathv[c(1:nbreaks)*2-1] + endhap = pathv[c(1:nbreaks)*2+1] + data.frame(starts = starts, int = int, markbp = markbp, marktrans = marktrans, starthap = starthap, endhap = endhap) +} + + + diff --git a/instructions/.DS_Store b/instructions/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..8cf122dcbd3847e12f72f2d21b665bef36a3d3e2 GIT binary patch literal 6148 zcmeHKJx{|x41F#gidZVKFfiQ6REd8Ol|se9(*A(b4^)ssBf&+2EkBU|RGxjOs#WQV zP_ZNXS zzkr$(G8P$|guFwMluAsgl8+dY(isnVT#>O!n9?En@F96-$tM)C*;zl>a7dBRtux>Z z>@uKdKT22k|J%doe;ee+8E^)UlL2WD#)CeeE$^*+PpNxtpgvMnbX=2glS0R~Vvg0V c_=0M}cuLLh&B~qrr_c@T&}b1LlB6DF6Tf literal 0 HcmV?d00001 From 6c59670d373ace20ebb0e9fd9418d0d463af9750 Mon Sep 17 00:00:00 2001 From: shackett Date: Mon, 17 Sep 2012 12:19:04 -0400 Subject: [PATCH 3/6] ready to merge --- fit-hmm.R | 76 +++++++++++++++++++------------------------------------ 1 file changed, 26 insertions(+), 50 deletions(-) diff --git a/fit-hmm.R b/fit-hmm.R index a2aec70..238d25d 100755 --- a/fit-hmm.R +++ b/fit-hmm.R @@ -3,34 +3,15 @@ ## options(warn=2,error=recover); options(error=quote(q("yes"))) -### manual param filling -setwd("/Users/Sean/Desktop/Andolfatto/9_12_toy") - -dollar0 <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/msg/blah" -dollar0 <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/msg" -toy_indivs <- c("indivA12_AATAAG", "indivE2_CAGCCG", "indivE4_TAGGAG") -direct <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/hmm_data" -outdir <- "/Users/Sean/Desktop/Andolfatto/9_12_toy/hmm_fit" - -#opts <- list(s = "male", i = toy_indivs[1], d = direct, o = outdir, p = 0.1, q = 0.1, r = 0.000001, c = "2L,2R,3L,3R,4", x = "X", y = "all", z = "0,.5,.5", t = 1, H = "both", e = 0.8, j = 30, l = 20) -opts <- list(s = "male", i = toy_indivs[1], d = direct, o = outdir, p = 0.1, q = 0.1, r = 0.000001, c = "2R", x = NULL, y = "all", z = "0,.5,.5", t = 1, H = "ML,viterbi", e = 0.8, j = 30, l = 20) - - - - - - - -#args <- commandArgs() -#dollar0 <- substring(args[grep("^--file=", args)], 8) +args <- commandArgs() +dollar0 <- substring(args[grep("^--file=", args)], 8) source(sprintf("%s/ded.R", dirname(dollar0))) source(sprintf("%s/hmmlib.R", dirname(dollar0))) library("R.methodsS3", lib.loc = dirname(dollar0)) library("R.oo", lib.loc = dirname(dollar0)) - -#opts <- getopts() +opts <- getopts() indivs <- unlist(strsplit(opts$i,split=",")) sex <- opts$s dir <- opts$d @@ -49,7 +30,6 @@ if("viterbi" %in% HMMtype){ HMM_seqpair <- as.numeric(opts$j) HMM_diffthresh <- as.numeric(opts$l) HMM_decay <- as.numeric(opts$e) - vit_bp_calls <- list() } stopifnot(!is.null(indivs), !is.null(dir), !is.null(outdir), length(indivs) == 1) @@ -72,7 +52,6 @@ plotPadding <- 10^(ceiling(log10(aveSpace))-2) alleles <- c("A","C","G","T") -indiv <- indivs for(indiv in indivs) { cat(indiv, "\n") ## if(opts$c == "all") @@ -235,9 +214,6 @@ for(indiv in indivs) { data$est <- apply(prob, 1, which.max) - rand_bp <- round(runif(1, min = 10, max = L-10)) - prob <- rbind(rdirichlet(L - rand_bp, c(0, 0.95, 0.05)), rdirichlet(rand_bp, c(0, 0.05, 0.95))) - ## Posterior probability if("ML" %in% HMMtype){ @@ -251,32 +227,32 @@ for(indiv in indivs) { if(ploidy == 2){zero <- matrix(c(1,6,3,4), ncol =2, byrow = TRUE)} if(ploidy == 1){zero <- NULL} - Int <- 1 - - if(length(prob[,1]) > 1){ - #Run Viterbi-HMM forward and backwards and output the state sequence - dual.hmm <- new.HMM_pathFR(phi, d, prob, zero, HMM_decay, r, Int) - datapos <- data$pos - - #Compare recombination events called by foward and reverse HMM runs and determine a consensus, likely set of recombinations - conc.hap <- new.HMM_path.consensus(dual.hmm, contig, K, data$pos) + Int <- 1 - Pr.z.given.y_vit <- t(sapply(1:L, function(x){ifelse(1:K == conc.hap$Ancestry[x], 1, 0)})) - - if(length(unique(conc.hap$Ancestry)) != 1){ - #gradient based on ambiguous breakpoint boundaries + if(length(prob[,1]) > 1){ + #Run Viterbi-HMM forward and backwards and output the state sequence + dual.hmm <- new.HMM_pathFR(phi, d, prob, zero, HMM_decay, r, Int) + datapos <- data$pos + + #Compare recombination events called by foward and reverse HMM runs and determine a consensus, likely set of recombinations + conc.hap <- new.HMM_path.consensus(dual.hmm, contig, K, data$pos) + + Pr.z.given.y_vit <- t(sapply(1:L, function(x){ifelse(1:K == conc.hap$Ancestry[x], 1, 0)})) + + if(length(unique(conc.hap$Ancestry)) != 1){ + #gradient based on ambiguous breakpoint boundaries - for(recon in 1:length(conc.hap$recoZ[,1])){ - reco <- conc.hap$recoZ[recon,] - interpol_p <- 1 - (datapos[reco$mark_start:reco$mark_end] - datapos[reco$mark_start])/(datapos[reco$mark_end] - datapos[reco$mark_start]) - - Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$starthap] <- interpol_p - Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$endhap] <- 1-interpol_p - }} + for(recon in 1:length(conc.hap$recoZ[,1])){ + reco <- conc.hap$recoZ[recon,] + interpol_p <- 1 - (datapos[reco$mark_start:reco$mark_end] - datapos[reco$mark_start])/(datapos[reco$mark_end] - datapos[reco$mark_start]) + + Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$starthap] <- interpol_p + Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$endhap] <- 1-interpol_p + } + } + } + colnames(Pr.z.given.y_vit) <- paste("viterbi Pr(", ancestries, "|y)") } - colnames(Pr.z.given.y_vit) <- paste("viterbi Pr(", ancestries, "|y)") - vit_bp_calls[[indiv]][[contig]] <- conc.hap$recoZ - } if("ML" %in% HMMtype){data <- cbind(data, Pr.z.given.y_ML)} if("viterbi" %in% HMMtype){data <- cbind(data, Pr.z.given.y_vit)} From 86ee908a48197667c22c7e55a66b6810219f6044 Mon Sep 17 00:00:00 2001 From: shackett Date: Tue, 27 Nov 2012 09:35:14 -0500 Subject: [PATCH 4/6] adding old changes --- fit-hmm.R | 5 ++--- hmmlib.R | 13 ++++++++++++- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/fit-hmm.R b/fit-hmm.R index 238d25d..8151c03 100755 --- a/fit-hmm.R +++ b/fit-hmm.R @@ -257,8 +257,7 @@ for(indiv in indivs) { if("ML" %in% HMMtype){data <- cbind(data, Pr.z.given.y_ML)} if("viterbi" %in% HMMtype){data <- cbind(data, Pr.z.given.y_vit)} - data <- cbind(data, Pr.z.given.y_ML, Pr.z.given.y_vit) - attr(data, "badpos") <- badpos + attr(data, "badpos") <- badpos dataa[[contig]] <- data } @@ -316,7 +315,7 @@ for(indiv in indivs) { ### divvy up homozygous and heterozygous blocks byBlocks <- breakpoint.width(x, y[,par1homo_col], y[,par2homo_col], indiv=indiv, contig=contig, conf1=.05 ,conf2=.95, hmmRunning); - if (is.null(byBlocks[["bps"]])==F) { breakpoints <- rbind(breakpoints,byBlocks[["bps"]]); } + if (nrow(byBlocks[["bps"]])>0) { breakpoints <- rbind(breakpoints,byBlocks[["bps"]]); } ### plot like.par1 <- contig_data[contig_data$read_allele==contig_data$par1ref,]$pos; diff --git a/hmmlib.R b/hmmlib.R index b24692e..a6ae7ac 100755 --- a/hmmlib.R +++ b/hmmlib.R @@ -242,7 +242,18 @@ breakpoint.width <- function(x, y1, y2, indiv=NA, contig=NA, conf1=.05 ,conf2=.9 width1 <- findBreak(indiv,contig,y1,y1>=conf1,y1<=conf2,conf1,conf2,'homozygous_par1') width2 <- findBreak(indiv,contig,y2,y2>=conf1,y2<=conf2,conf1,conf2,'homozygous_par2') - list(blocks=cbind(rbind(width1[["blocks"]],width2[["blocks"]]), hmmtype),bps=cbind(rbind(width1[["bps"]],width2[["bps"]]), hmmtype)) + tmp <- list() + if(length(width1[["blocks"]]) + length(width2[["blocks"]]) != 0){ + tmp[["blocks"]] <- cbind(rbind(width1[["blocks"]],width2[["blocks"]]), hmmtype) + }else{ + tmp[["blocks"]] <- rbind(width1[["blocks"]],width2[["blocks"]]) + } + if(length(width1[["bps"]]) + length(width2[["bps"]]) != 0){ + tmp[["bps"]] <- cbind(rbind(width1[["bps"]],width2[["bps"]]), hmmtype) + }else{ + tmp[["bps"]] <- rbind(width1[["bps"]],width2[["bps"]]) + } + tmp } From 8d50e594d039760fea51dc9b375173cb784a1c5c Mon Sep 17 00:00:00 2001 From: shackett Date: Tue, 27 Nov 2012 09:58:11 -0500 Subject: [PATCH 5/6] changes to run either viterbi or ML, not both, in order to prevent branching into all downstream mapping methods --- fit-hmm.R | 65 +++++++++++++++++++++++++------------------------------ hmmlib.R | 15 ++----------- 2 files changed, 32 insertions(+), 48 deletions(-) diff --git a/fit-hmm.R b/fit-hmm.R index 8151c03..e9febd7 100755 --- a/fit-hmm.R +++ b/fit-hmm.R @@ -21,12 +21,15 @@ deltapar2 <- as.numeric(opts$q) priors <- unlist(strsplit(opts$z,split=",")) theta <- as.numeric(opts$t) -HMMtype <- unlist(strsplit(opts$H, ",")) +HMMtype <- as.character(opts$H) +if(!(HMMtype %in% c("ML", "viterbi"))){ + print(paste("HMMtype must be either ML or viterbi, current option is:", HMMtype)) + } -if("ML" %in% HMMtype){ +if(HMMtype == "ML"){ rfac <- as.numeric(opts$r) } -if("viterbi" %in% HMMtype){ +if(HMMtype == "viterbi"){ HMM_seqpair <- as.numeric(opts$j) HMM_diffthresh <- as.numeric(opts$l) HMM_decay <- as.numeric(opts$e) @@ -146,8 +149,8 @@ for(indiv in indivs) { r <- 1 / contigLengths[1,"length"] ## Arbitrarily use the first contig for unassembled contigs } - if("ML" %in% HMMtype){ - d <- c(NA, diff(data$pos)) + d <- c(NA, diff(data$pos)) + if(HMMtype == "ML"){ p <- 1 - exp(-r*d*rfac) Pi <- array(dim=c(L,K,K), dimnames=list(NULL, ancestries, ancestries)) if(ploidy == 2) { @@ -216,19 +219,18 @@ for(indiv in indivs) { ## Posterior probability - if("ML" %in% HMMtype){ + if(HMMtype == "ML"){ hmm <- forwardback.ded(Pi=Pi, delta=phi, prob=prob) #hmm <- forwardback.ded(Pi=Pi, delta=rep(1/K, K), prob=prob) - Pr.z.given.y_ML <- exp(hmm$logalpha + hmm$logbeta - hmm$LL) - colnames(Pr.z.given.y_ML) <- paste("ML Pr(", ancestries, "|y)") - } - if("viterbi" %in% HMMtype){ + Pr.z.given.y <- exp(hmm$logalpha + hmm$logbeta - hmm$LL) + } + if(HMMtype == "viterbi"){ - if(ploidy == 2){zero <- matrix(c(1,6,3,4), ncol =2, byrow = TRUE)} + if(ploidy == 2){zero <- matrix(c(1,6,3,4), ncol=2, byrow = TRUE)} if(ploidy == 1){zero <- NULL} - + + #Int could be replaced by a vector to add interference, but its vestigial as is. Int <- 1 - if(length(prob[,1]) > 1){ #Run Viterbi-HMM forward and backwards and output the state sequence dual.hmm <- new.HMM_pathFR(phi, d, prob, zero, HMM_decay, r, Int) @@ -237,7 +239,7 @@ for(indiv in indivs) { #Compare recombination events called by foward and reverse HMM runs and determine a consensus, likely set of recombinations conc.hap <- new.HMM_path.consensus(dual.hmm, contig, K, data$pos) - Pr.z.given.y_vit <- t(sapply(1:L, function(x){ifelse(1:K == conc.hap$Ancestry[x], 1, 0)})) + Pr.z.given.y <- t(sapply(1:L, function(x){ifelse(1:K == conc.hap$Ancestry[x], 1, 0)})) if(length(unique(conc.hap$Ancestry)) != 1){ #gradient based on ambiguous breakpoint boundaries @@ -246,17 +248,15 @@ for(indiv in indivs) { reco <- conc.hap$recoZ[recon,] interpol_p <- 1 - (datapos[reco$mark_start:reco$mark_end] - datapos[reco$mark_start])/(datapos[reco$mark_end] - datapos[reco$mark_start]) - Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$starthap] <- interpol_p - Pr.z.given.y_vit[reco$mark_start:reco$mark_end,reco$endhap] <- 1-interpol_p + Pr.z.given.y[reco$mark_start:reco$mark_end,reco$starthap] <- interpol_p + Pr.z.given.y[reco$mark_start:reco$mark_end,reco$endhap] <- 1-interpol_p } } } - colnames(Pr.z.given.y_vit) <- paste("viterbi Pr(", ancestries, "|y)") } - - if("ML" %in% HMMtype){data <- cbind(data, Pr.z.given.y_ML)} - if("viterbi" %in% HMMtype){data <- cbind(data, Pr.z.given.y_vit)} - + + colnames(Pr.z.given.y) <- paste("viterbi Pr(", ancestries, "|y)") + data <- cbind(data, Pr.z.given.y) attr(data, "badpos") <- badpos dataa[[contig]] <- data @@ -266,7 +266,6 @@ for(indiv in indivs) { cat("OK\n") } - contigLengths <- contigLengths[plot.contigs,] ## Track the width of breakpoints @@ -277,20 +276,16 @@ for(indiv in indivs) { plotfile <- file.path(outdir, indiv, paste(indiv, "hmmprob.pdf", sep="-")) if(file.exists(plotfile)) { cat("plot already exists\n") ; next } pdf(file=plotfile, width=7, height=1.5) - par(mar=c(2,2.5,0.5,0.5),bg="transparent",cex.main=.68,cex.lab=.8,font.lab=2,cex.axis=.38,mgp=c(1,.000001,0),xaxs="i") - for(hmmRunning in HMMtype){ - plot(0, 0, xlab="", ylab="", col="transparent", xlim=c(1,sum(as.numeric(contigLengths$length)) + plotPadding*(length(plot.contigs)+1)), ylim=c(-1.01,1.01), axes=F) axis(side=2,at=c(-1,0,1),labels=c("","",""),col="gray38"); mtext(c("par2","par1"),side=2,line=.68,at=c(-1,1),font=2,cex=.8,col=c("blue","red"),las=2); - mtext(paste("HMM:", hmmRunning, collapse = " "), side = 2, line = -0.8, at = -1.4, font=2, cex=0.6, las=2) - box(col="gray68"); + box(col="gray68"); current_start <- plotPadding; - for(contig in plot.contigs) { + for(contig in plot.contigs) { mtext(side=1,at=current_start,contig,font=2,cex=.8,line=1,xpd=T,adj=0) current_end <- current_start + contigLengths[contigLengths$chr == contig,"length"] - 1; @@ -311,31 +306,31 @@ for(indiv in indivs) { if (sum(names(dataa) %in% contig)!=0) { contig_data <- dataa[[contig]]; x <- contig_data$pos - y <- contig_data[,paste(hmmRunning, "Pr(", ancestries, "|y)")] + y <- contig_data[,paste("Pr(", ancestries, "|y)")] ### divvy up homozygous and heterozygous blocks - byBlocks <- breakpoint.width(x, y[,par1homo_col], y[,par2homo_col], indiv=indiv, contig=contig, conf1=.05 ,conf2=.95, hmmRunning); - if (nrow(byBlocks[["bps"]])>0) { breakpoints <- rbind(breakpoints,byBlocks[["bps"]]); } + byBlocks <- breakpoint.width(x, y[,par1homo_col], y[,par2homo_col], indiv=indiv, contig=contig, conf1=.05 ,conf2=.95); + if (is.null(byBlocks[["bps"]])==F) { breakpoints <- rbind(breakpoints,byBlocks[["bps"]]); } ### plot like.par1 <- contig_data[contig_data$read_allele==contig_data$par1ref,]$pos; like.par2 <- contig_data[contig_data$read_allele==contig_data$par2ref,]$pos; plot.posterior(x+current_start, y, ancestries, like.par1+current_start, like.par2+current_start, bounds=c(1,contigLengths[contigLengths$chr==contig,]$length)+current_start-1, subtract=current_start, tickwidth=5*10^7) - + ### report mismatch fraction for homozygous regions if (nrow(byBlocks[["blocks"]])>0) { ## plot fraction of par1/(par1+par2) among informative markers (between -1 and 1) - matchMismatch <- rbind(matchMismatch, cbind(reportCounts(contig_data, as.vector(byBlocks[["blocks"]][,"V1"]), as.vector(byBlocks[["blocks"]][,"V4"]), as.numeric(as.vector(byBlocks[["blocks"]][,"V7"])), as.numeric(as.vector(byBlocks[["blocks"]][,"V8"]))), hmmtype = hmmRunning)) + matchMismatch <- rbind(matchMismatch, reportCounts(contig_data, as.vector(byBlocks[["blocks"]][,"V1"]), as.vector(byBlocks[["blocks"]][,"V4"]), as.numeric(as.vector(byBlocks[["blocks"]][,"V7"])), as.numeric(as.vector(byBlocks[["blocks"]][,"V8"])))) } } current_start <- current_start + contigLengths[contigLengths$chr == contig,"length"] + plotPadding; if (contig != plot.contigs[length(plot.contigs)]) { - abline(v=current_start-(plotPadding/2),col="gray68",lwd=1) + abline(v=current_start-(plotPadding/2),col="gray68",lwd=1) } } - } + etc <- "" main <- sprintf("%s (%s): delta=(%.0e, %.0e)", indiv, sex, deltapar1, deltapar2) dev.off() diff --git a/hmmlib.R b/hmmlib.R index a6ae7ac..7c07ad6 100755 --- a/hmmlib.R +++ b/hmmlib.R @@ -238,22 +238,11 @@ findBreak <- function(indiv,contig,y,z1,z2,conf1,conf2,parentage) { list(blocks=as.data.frame(blocks),bps=as.data.frame(widths)); } -breakpoint.width <- function(x, y1, y2, indiv=NA, contig=NA, conf1=.05 ,conf2=.95, hmmtype) { +breakpoint.width <- function(x, y1, y2, indiv=NA, contig=NA, conf1=.05 ,conf2=.95) { width1 <- findBreak(indiv,contig,y1,y1>=conf1,y1<=conf2,conf1,conf2,'homozygous_par1') width2 <- findBreak(indiv,contig,y2,y2>=conf1,y2<=conf2,conf1,conf2,'homozygous_par2') - tmp <- list() - if(length(width1[["blocks"]]) + length(width2[["blocks"]]) != 0){ - tmp[["blocks"]] <- cbind(rbind(width1[["blocks"]],width2[["blocks"]]), hmmtype) - }else{ - tmp[["blocks"]] <- rbind(width1[["blocks"]],width2[["blocks"]]) - } - if(length(width1[["bps"]]) + length(width2[["bps"]]) != 0){ - tmp[["bps"]] <- cbind(rbind(width1[["bps"]],width2[["bps"]]), hmmtype) - }else{ - tmp[["bps"]] <- rbind(width1[["bps"]],width2[["bps"]]) - } - tmp + list(blocks=rbind(width1[["blocks"]],width2[["blocks"]]),bps=rbind(width1[["bps"]],width2[["bps"]])) } From 5d5867f86e511d00dc663a87fd44fcb12f610e57 Mon Sep 17 00:00:00 2001 From: shackett Date: Tue, 27 Nov 2012 10:12:26 -0500 Subject: [PATCH 6/6] fixed a wrong column name --- fit-hmm.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/fit-hmm.R b/fit-hmm.R index e9febd7..1fdb547 100755 --- a/fit-hmm.R +++ b/fit-hmm.R @@ -255,7 +255,7 @@ for(indiv in indivs) { } } - colnames(Pr.z.given.y) <- paste("viterbi Pr(", ancestries, "|y)") + colnames(Pr.z.given.y) <- paste("Pr(", ancestries, "|y)") data <- cbind(data, Pr.z.given.y) attr(data, "badpos") <- badpos dataa[[contig]] <- data