add in new HMM plots

Author: xji3

.DS_Store

@@ -237,3 +237,34 @@ colnames(summary.mat) <- c("lnL", "max lnL", "tract length",
 "Pr(S_0)", "Pr(S_1)",
 "df", "d^2f", "c.i. tract length", "c.i tract length")
 summary.mat
+source('~/.active-rstudio-document', echo=TRUE)
+# plot one paralog
+paralog = "EDN_ECP"
+lnL.ratio <- as.vector(read.table(paste("./summary/", paralog, "_MG94_nonclock_HMM_log_posterior_ratio.txt", sep = "")))
+Viterbi.path <- as.vector(read.table(paste("./summary/", paralog, "_MG94_nonclock_HMM_Viterbi_path.txt", sep = "")))
+lnL.surface <- as.vector(read.table(paste("./summary/", paralog, "_MG94_nonclock_HMM_lnL_surface.txt", sep = "")))
+IGC.sw.lnL <- as.vector(read.table(paste("./summary/", paralog, "_MG94_nonclock_sw_lnL.txt", sep = "")))
+Force.sw.lnL <- as.vector(read.table(paste("./summary/Force_", paralog, "_MG94_nonclock_sw_lnL.txt", sep = "")))
+plot(lnL.ratio[, 1], xlab = "codon number", ylab = "log values",
+type = "l", col = "black", lty = 1,
+main = paste(paralog, " HMM result"),
+ylim = c(min(-0.5, min(lnL.ratio)), max(lnL.ratio)))
+lines(1:dim(Viterbi.path)[1], Viterbi.path[, 1], type = "S", lty = 2, col = "red")
+#lines(1:dim(IGC.sw.lnL)[1], IGC.sw.lnL[, 2] - Force.sw.lnL[, 2], type = "l", lty = 3, col = "red")
+legend(1, max(lnL.ratio), legend = c("log posterior ratio", "Viterbi path"),
+lty = c(1, 2), col = c( "black", "red"))
+plot(-lnL.surface[, 1], xlab = "tract length in nucleotide", ylab= "lnL", type = "l", col = "black", lty = 1,
+main = paste(paralog, " lnL surface"))
+summary.mat <- read.table("./HMM_tract_MG94_nonclock_summary.txt")
+# Now calculate standard deviation of lnP
+lnP <- log(3.0 / summary.mat[, 3])
+sd.lnP <- 1.0 / sqrt(-summary.mat[, 7])
+low.cf <- exp(lnP - 1.96 * sd.lnP)
+up.cf  <- exp(lnP + 1.96 * sd.lnP)
+up.cf[up.cf > 1] <- 1.0
+summary.mat <- cbind(summary.mat, 3.0 / up.cf, 3.0 / low.cf)
+rownames(summary.mat) <- c("EDN_ECP")
+colnames(summary.mat) <- c("lnL", "max lnL", "tract length",
+"Pr(S_0)", "Pr(S_1)",
+"df", "d^2f", "c.i. tract length", "c.i tract length")
+summary.mat

.Rhistory

@@ -68,23 +68,29 @@ def main(args):
     seq_index_file = './' + '_'.join(paralog) + '_seq_index.txt'
     state_list = ['No IGC event (Si = 0)','At least one IGC event (Si > 0)']
 
-    x_3 = np.array([-0.69727878,
-                  -0.53710801,
-                  -0.72400474,
-                  0.72385788,
-                  -0.18983735,
-                  -2.5199719,
-                  -2.01452935,
-                  -1.0337633,
-                  -3.29369029,
-                  -1.75318807,
-                  -3.25869777,
-                  -2.27341043,
-                  -4.20160402,
-                  -4.110472,
-                  -0.5])
+    save_name = './save/EDN_ECP_Ind_MG94_IGC_save.txt'
+    summary_file = './summary/EDN_ECP_Ind_MG94_HMM_summary.txt'
+    summary_file_1D = './summary/EDN_ECP_Ind_MG94_HMM_1D_summary.txt'
+    Ind_MG94_IGC = IndCodonGeneconv( newicktree, alignment_file, paralog, Model = 'MG94', Force = Force, clock = None, save_path = './save/',
+                             save_name = save_name)
+    Ind_MG94_IGC.get_mle(True, True, 0, 'BFGS')
+    x = np.concatenate((Ind_MG94_IGC.x, [0.0]))
+
+    IGC_sitewise_lnL_file = summary_path + '_'.join(paralog) + '_Ind_MG94_nonclock_sw_lnL.txt'
+    NOIGC_sitewise_lnL_file = summary_path + 'NOIGC_' + '_'.join(paralog) + '_Ind_MG94_nonclock_sw_lnL.txt'
+        
+    Ind_MG94_IGC.get_sitewise_loglikelihood_summary(IGC_sitewise_lnL_file, False)
+    Ind_MG94_IGC.get_sitewise_loglikelihood_summary(NOIGC_sitewise_lnL_file, True)
+
+    
+    test = HMMJSGeneconv(save_file, newicktree, alignment_file, paralog, summary_path, x, save_path, IGC_sitewise_lnL_file, NOIGC_sitewise_lnL_file,                         state_list, seq_index_file)
+    
+
+    log_p_list = np.log(3.0/np.array(range(3, 1001)))
+    plot_file = './plot/HMM_' + '_'.join(paralog) + '_lnL_1D_surface.txt'
+    test.plot_tract_p(log_p_list, plot_file)
+
+    test.get_mle(display = True, two_step = True, One_Dimension = True)
+    test.get_summary(summary_file_1D)
+    
 
-    test = HMMJSGeneconv(save_file, newicktree, alignment_file, paralog, summary_path, x_3, save_path, IGC_save_name, Force_save_name,
-                         IGC_sitewise_lnL_file, Force_sitewise_lnL_file,
-                         state_list, seq_index_file)
-    test.get_mle(display = True, two_step = False, OneDimension = True)

Run_IS_HMM.py

@@ -0,0 +1,19 @@
+1.694970388693524455e+03
+2.909982296947770108e-01
+2.434371217926649678e-01
+2.068434784667338988e-01
+2.587211700458240116e-01
+2.070543014300672890e+00
+8.467952031245288280e-01
+5.280029360273003070e-01
+2.446307008266619676e-02
+4.633131449367799276e-02
+1.964016947431149462e-01
+3.266001629134150375e-01
+3.265333578314633806e-02
+1.562930819382047642e-01
+3.460162330033541428e-02
+9.009211334423779249e-02
+1.432976146309197611e-02
+1.603597178254540304e-02
+# ll pi_a pi_c pi_g pi_t kappa omega tau eta tract_p (N0,N1) (N0,Tamarin) (N1,N2) (N1,Macaque) (N2,N3) (N2,Orangutan) (N3,Chimpanzee) (N3,Gorilla)