Added AB Matrix: Antibody model

models/protein/matrices.go

@@ -1185,7 +1185,6 @@ func WAGMats() (dmat *mat.Dense, pi []float64) {
 // PLoS ONE. 2007 Jun 6;2:e503.
 // [thanks to Sergei L. Kosakovsky]
 // Translated from HYPHY to Phyml format by Federico Abascal.
-//
 func HIVBMats() (dmat *mat.Dense, pi []float64) {
 	var i, j, naa int
 	naa = 20
@@ -1413,3 +1412,237 @@ func HIVBMats() (dmat *mat.Dense, pi []float64) {
 	dmat = mat.NewDense(naa, naa, m)
 	return
 }
+
+// Antibody-Specific Model of Amino Acid Substitution
+// From AB Model publication (supplementary material online)
+// Alexander Mirsky, Linda Kazandjian, Maria Anisimova,
+// Antibody-Specific Model of Amino Acid Substitution for Immunological Inferences from Alignments of Antibody Sequences,
+// Molecular Biology and Evolution, Volume 32, Issue 3, March 2015, Pages 806–819,
+// https://doi.org/10.1093/molbev/msu340
+func ABMats() (dmat *mat.Dense, pi []float64) {
+	var i, j, naa int
+	naa = 20
+	m := make([]float64, naa*naa)
+	pi = make([]float64, naa)
+
+	m[1*20+0] = 1.784266e-01
+	m[2*20+0] = 9.291290e-02
+	m[2*20+1] = 7.829130e-01
+	m[3*20+0] = 1.241095e+00
+	m[3*20+1] = 5.795374e-02
+	m[3*20+2] = 7.185182e+00
+	m[4*20+0] = 8.929181e-03
+	m[4*20+1] = 1.821885e-01
+	m[4*20+2] = 1.374268e-06
+	m[4*20+3] = 2.340019e-02
+	m[5*20+0] = 1.992269e-01
+	m[5*20+1] = 1.923901e+00
+	m[5*20+2] = 8.705989e-02
+	m[5*20+3] = 1.843856e-01
+	m[5*20+4] = 1.046446e-08
+	m[6*20+0] = 9.521821e-01
+	m[6*20+1] = 6.273863e-02
+	m[6*20+2] = 5.038373e-01
+	m[6*20+3] = 7.426619e+00
+	m[6*20+4] = 7.519215e-11
+	m[6*20+5] = 3.691671e+00
+	m[7*20+0] = 1.851951e+00
+	m[7*20+1] = 1.089400e+00
+	m[7*20+2] = 4.868901e-01
+	m[7*20+3] = 2.112400e+00
+	m[7*20+4] = 5.891123e-02
+	m[7*20+5] = 5.516340e-02
+	m[7*20+6] = 1.389370e+00
+	m[8*20+0] = 5.241316e+00
+	m[8*20+1] = 1.049550e+01
+	m[8*20+2] = 1.405444e+01
+	m[8*20+3] = 1.126995e+01
+	m[8*20+4] = 3.963388e+00
+	m[8*20+5] = 8.908434e+00
+	m[8*20+6] = 7.298080e+00
+	m[8*20+7] = 9.139518e+00
+	m[9*20+0] = 1.140412e-01
+	m[9*20+1] = 3.245175e-01
+	m[9*20+2] = 1.762721e+00
+	m[9*20+3] = 3.916999e-02
+	m[9*20+4] = 6.594967e-04
+	m[9*20+5] = 6.712736e-06
+	m[9*20+6] = 1.029959e-04
+	m[9*20+7] = 3.560482e-02
+	m[9*20+8] = 4.706586e+00
+	m[10*20+0] = 6.969101e-02
+	m[10*20+1] = 3.932002e-01
+	m[10*20+2] = 2.769442e-02
+	m[10*20+3] = 3.020502e-02
+	m[10*20+4] = 6.079219e-03
+	m[10*20+5] = 6.802781e-01
+	m[10*20+6] = 1.283121e-03
+	m[10*20+7] = 2.157936e-02
+	m[10*20+8] = 5.879103e+00
+	m[10*20+9] = 1.601123e+00
+	m[11*20+0] = 7.388355e-02
+	m[11*20+1] = 7.549240e+00
+	m[11*20+2] = 6.190318e+00
+	m[11*20+3] = 6.622772e-02
+	m[11*20+4] = 3.722878e-16
+	m[11*20+5] = 3.030805e+00
+	m[11*20+6] = 3.608816e+00
+	m[11*20+7] = 5.504400e-02
+	m[11*20+8] = 1.455741e+00
+	m[11*20+9] = 5.059793e-01
+	m[11*20+10] = 2.158451e-02
+	m[12*20+0] = 6.299271e-02
+	m[12*20+1] = 3.362326e-01
+	m[12*20+2] = 3.972173e-02
+	m[12*20+3] = 3.357577e-02
+	m[12*20+4] = 7.213178e-03
+	m[12*20+5] = 1.233336e-03
+	m[12*20+6] = 7.659566e-02
+	m[12*20+7] = 2.187264e-02
+	m[12*20+8] = 2.298295e+00
+	m[12*20+9] = 1.096748e+01
+	m[12*20+10] = 5.647985e+00
+	m[12*20+11] = 1.238634e+00
+	m[13*20+0] = 1.130146e-01
+	m[13*20+1] = 8.208677e-02
+	m[13*20+2] = 1.955446e-01
+	m[13*20+3] = 1.031734e-01
+	m[13*20+4] = 1.993818e-01
+	m[13*20+5] = 1.496610e-03
+	m[13*20+6] = 5.288625e-02
+	m[13*20+7] = 1.984772e-01
+	m[13*20+8] = 5.642309e+00
+	m[13*20+9] = 2.714705e+00
+	m[13*20+10] = 3.390618e+00
+	m[13*20+11] = 4.649035e-03
+	m[13*20+12] = 3.947940e+00
+	m[14*20+0] = 1.800713e+00
+	m[14*20+1] = 3.498713e-01
+	m[14*20+2] = 7.342554e-03
+	m[14*20+3] = 1.509482e-01
+	m[14*20+4] = 4.878395e-03
+	m[14*20+5] = 7.426909e-01
+	m[14*20+6] = 2.889815e-02
+	m[14*20+7] = 7.915056e-02
+	m[14*20+8] = 1.049496e+01
+	m[14*20+9] = 5.016568e-02
+	m[14*20+10] = 1.149931e+00
+	m[14*20+11] = 9.948994e-03
+	m[14*20+12] = 7.417279e-02
+	m[14*20+13] = 3.556198e-01
+	m[15*20+0] = 9.988358e-01
+	m[15*20+1] = 1.926435e+00
+	m[15*20+2] = 7.348346e+00
+	m[15*20+3] = 5.822988e-01
+	m[15*20+4] = 2.639482e-01
+	m[15*20+5] = 5.906405e-04
+	m[15*20+6] = 6.776709e-02
+	m[15*20+7] = 9.984215e-01
+	m[15*20+8] = 5.439116e+00
+	m[15*20+9] = 6.007607e-01
+	m[15*20+10] = 1.580539e-01
+	m[15*20+11] = 8.688405e-02
+	m[15*20+12] = 1.861354e-02
+	m[15*20+13] = 9.813064e-01
+	m[15*20+14] = 1.284651e+00
+	m[16*20+0] = 2.912317e+00
+	m[16*20+1] = 1.135258e+00
+	m[16*20+2] = 2.147175e+00
+	m[16*20+3] = 1.516881e-01
+	m[16*20+4] = 3.225214e-06
+	m[16*20+5] = 1.202094e-01
+	m[16*20+6] = 6.016624e-02
+	m[16*20+7] = 7.862767e-02
+	m[16*20+8] = 3.443285e+00
+	m[16*20+9] = 3.087152e+00
+	m[16*20+10] = 5.702792e-01
+	m[16*20+11] = 1.039298e+00
+	m[16*20+12] = 1.415612e+00
+	m[16*20+13] = 3.674486e-02
+	m[16*20+14] = 9.057112e-01
+	m[16*20+15] = 3.058575e+00
+	m[17*20+0] = 7.939549e-02
+	m[17*20+1] = 5.724286e-01
+	m[17*20+2] = 7.310937e-04
+	m[17*20+3] = 1.423897e-02
+	m[17*20+4] = 4.440833e-01
+	m[17*20+5] = 4.332983e-05
+	m[17*20+6] = 2.252612e-02
+	m[17*20+7] = 1.386853e-01
+	m[17*20+8] = 7.013890e+00
+	m[17*20+9] = 6.318748e-02
+	m[17*20+10] = 3.378544e-01
+	m[17*20+11] = 8.024263e-03
+	m[17*20+12] = 1.011149e-01
+	m[17*20+13] = 2.199856e-01
+	m[17*20+14] = 5.516074e-03
+	m[17*20+15] = 1.385142e-01
+	m[17*20+16] = 1.412361e-02
+	m[18*20+0] = 1.433528e-01
+	m[18*20+1] = 1.711315e-01
+	m[18*20+2] = 2.622763e+00
+	m[18*20+3] = 9.078338e-01
+	m[18*20+4] = 7.741612e-01
+	m[18*20+5] = 2.737091e-02
+	m[18*20+6] = 1.240642e-01
+	m[18*20+7] = 2.295842e-01
+	m[18*20+8] = 2.055414e+01
+	m[18*20+9] = 2.903165e-01
+	m[18*20+10] = 1.521320e-01
+	m[18*20+11] = 7.109973e-02
+	m[18*20+12] = 2.246759e-03
+	m[18*20+13] = 7.074464e+00
+	m[18*20+14] = 1.992133e-01
+	m[18*20+15] = 8.104751e-01
+	m[18*20+16] = 9.984255e-02
+	m[18*20+17] = 6.121284e-01
+	m[19*20+0] = 3.774477e+00
+	m[19*20+1] = 1.366145e-01
+	m[19*20+2] = 4.931206e-02
+	m[19*20+3] = 4.076074e-01
+	m[19*20+4] = 2.243512e-02
+	m[19*20+5] = 9.047737e-03
+	m[19*20+6] = 5.795409e-01
+	m[19*20+7] = 4.228200e-01
+	m[19*20+8] = 6.890244e+00
+	m[19*20+9] = 7.926675e+00
+	m[19*20+10] = 3.595310e+00
+	m[19*20+11] = 3.493440e-02
+	m[19*20+12] = 4.396720e+00
+	m[19*20+13] = 1.643946e+00
+	m[19*20+14] = 2.217442e-01
+	m[19*20+15] = 7.477041e-02
+	m[19*20+16] = 2.166054e-01
+	m[19*20+17] = 9.663569e-02
+	m[19*20+18] = 5.010635e-01
+
+	for i = 0; i < naa; i++ {
+		for j = 0; j < i; j++ {
+			m[j*naa+i] = m[i*naa+j]
+		}
+	}
+
+	pi[0] = 6.541704e-02
+	pi[1] = 4.708366e-02
+	pi[2] = 3.168984e-02
+	pi[3] = 4.688141e-02
+	pi[4] = 2.150693e-02
+	pi[5] = 4.240711e-02
+	pi[6] = 2.842211e-02
+	pi[7] = 1.005278e-01
+	pi[8] = 9.812606e-03
+	pi[9] = 3.424424e-02
+	pi[10] = 6.222565e-02
+	pi[11] = 4.844488e-02
+	pi[12] = 1.760370e-02
+	pi[13] = 3.478555e-02
+	pi[14] = 3.962469e-02
+	pi[15] = 1.280566e-01
+	pi[16] = 8.199314e-02
+	pi[17] = 3.393045e-02
+	pi[18] = 7.586119e-02
+	pi[19] = 4.948141e-02
+
+	dmat = mat.NewDense(naa, naa, m)
+	return
+}

models/protein/model.go

@@ -14,6 +14,7 @@ const (
 	MODEL_LG
 	MODEL_WAG
 	MODEL_HIVB
+	MODEL_AB
 
 	BL_MIN  = 1.e-08
 	BL_MAX  = 100.0
@@ -52,8 +53,11 @@ func NewProtModel(model int, usegamma bool, alpha float64) (*ProtModel, error) {
 		m, pi = WAGMats()
 	case MODEL_HIVB:
 		m, pi = HIVBMats()
+	case MODEL_AB:
+		m, pi = ABMats()
+
 	default:
-		return nil, fmt.Errorf("This protein model is not implemented")
+		return nil, fmt.Errorf("this protein model is not implemented")
 	}
 	return &ProtModel{
 		pi,
@@ -84,6 +88,8 @@ func ModelStringToInt(model string) int {
 		return MODEL_WAG
 	case "hivb":
 		return MODEL_HIVB
+	case "ab":
+		return MODEL_AB
 	default:
 		return -1
 	}
@@ -100,7 +106,7 @@ func (model *ProtModel) InitModel(aafreqs []float64) error {
 	ns := 20
 
 	if model.mat == nil || model.pi == nil {
-		return fmt.Errorf("Matrices have not been initialized")
+		return fmt.Errorf("matrices have not been initialized")
 	}
 
 	if aafreqs != nil && len(aafreqs) != ns {
@@ -129,7 +135,7 @@ func (model *ProtModel) InitModel(aafreqs []float64) error {
 
 	model.eigen = &mat.Eigen{}
 	if ok = model.eigen.Factorize(model.mat, mat.EigenRight); !ok {
-		return fmt.Errorf("Problem during matrix decomposition")
+		return fmt.Errorf("problem during matrix decomposition")
 	}
 	model.reigenvect = mat.NewDense(ns, ns, nil)
 	model.leigenvect = mat.NewDense(20, 20, nil)