Modified files for strain-specific residual variance option #269

jb621-star · 2023-03-20T14:12:36Z

These changes include the modifications to the math inside mvlmm, and attempts to add an option for the user to input their own strain-specific residual variance column within the files gemma_io.cpp and param.cpp. Thanks for your help!

… option similar to cvt

…eplaced in the math with the null Ve

Small formatting change to -residvar

Updated resid_var math to CalcQi

Updated most functions except MphEM and remaining P-val functions

Continued to make appropriate math changes to MphCalcSigma, UpdateV, MphCalcLogL, MphEM, MphCalcP, MphCalcBeta, and other functions.

Finished theoretical changes in other functions and replaced all "eps_eval" entries with "eval_vec, sigmasq".

Replaced any lingering mentions of "eps_eval" with "eval_vec, sigmasq".

Changed mentions of "eps_eval" with "sigmasq". Still need to figure out how to get eigenvectors from the same place "eval" is derived.

Changed mention of eps_eval as a vector to "gsl_matrix *sigmasq"

Changed the one mention of "eps_eval" to "eval_vec" and "sigmasq". May actually need to return to this file in order for the code to actually run.

eps_eval --> eval_vec, sigmasq

Changed eps_eval to sigmasq and corrected any identified regions where dimensions did not match.

Replaced input of CopyResid from eps_eval to sigmasq

aartaka

Here's what I've noticed. Most of the errors you're getting are due to mismatch between the arguments function requires and the arguments you provide, so ensure these are in sync. The rest I've tried to locate and comment on. Let's have a call (coordinate via Matrix) in case there's something I missed.

aartaka · 2024-10-20T04:25:48Z

src/param.h

@@ -136,6 +136,7 @@ class PARAM {
  string file_anno; // Optional.
  string file_gxe;  // Optional.
  string file_cvt;  // Optional.
+  string file_residvar; // Optional


Notice you called it file_residvar here, but are using it as file_resid elsewhere. A typo?

src/gemma.cpp

aartaka · 2024-10-20T04:27:33Z

src/gemma.cpp

@@ -2795,22 +2816,22 @@ void GEMMA::BatchRun(PARAM &cPar) {
          if (!cPar.file_bfile.empty()) {
            // PLINK analysis
            if (cPar.file_gxe.empty()) {
-              cLmm.AnalyzePlink(U, eval, UtW, &UtY_col.vector, W,
+              cLmm.AnalyzePlink(U, eval, eval_vec, sigmasq, UtW, &UtY_col.vector, W,


You didn't define eval_vec anywhere earlier, that's why compiler is shouting. What is the intended use for eval_vec?

Also note that sigmasq is not defined either.

And then there's also required/provided argument mismatch. AnalyzePlink requires (U, eval, UtW, Uty, W, y, gwasnps), while you provide (U, eval, eval_vec, sigmasq, UtW, Uty, W, y, gwasnps). I guess AnalyzePlink needs to be modified to accept new arguments first.

aartaka · 2024-10-20T04:33:02Z

src/mvlmm.h

@@ -35,6 +35,7 @@ class MVLMM {

  string file_bfile;
  string file_geno;
+  //string file_residvar maybe?


Yes, if you need to use it (although you seems to only use residuals as sigmasq, so maybe no need for it.) And don't forget to add a line assigning it to CopyFromParam in case you add this field.

Changed "file_residvar" back to "file_resid"

Transferred all instances of "eval_vec" to "U".

Changed all instances of "_residvar" to "resid" to keep things consistent.

Changed all instances of "residvar" to "resid" to keep things consistent

Changed inputs to "CalcLmmVgVeBeta", "CalcLambda", and "AnalyzePlink" according to their use elsewhere (i.e. no more "eval_vec").

Changed residvar to resid for consistency.

Changed all mentions of residvar to resid for consistency

Changed all residvar to resid for consistency.

Changed "residvar" to "resid" and removed unnecessary mentions of "eval_vec" from other functions.

Removed repeat parameter mentions of eigenvectors U in Analyze functions

Fixed examples of mishandled data types from when defining "scalar" data type to: double scalar = gsl_matrix_get(sigmasq, 0, 0) / ve; gsl_matrix_scale(Sigma, scalar); double epsilon = gsl_matrix_get(Sigma, k, k);

Fixed remaining cases of how the "scalar" variable is defined, and how many other functions are attributed.

Cleaned up various parameter specifications within function definitions.

Corrected CalcLambda and CalcLmmVgVeBeta mentions that were actually from lmm no mvlmm

sigmasq needs to be in the AnalyzeBimbam and AnalyzePlink function definitions.

const gsl_matrix *sigmasq, needs to be in these functions due to MphEM calls within other functions in mvlmm.cpp.

Removed duplicate mention of U within AnalyzePlinkGXE

Updated AnalyzePlink and AnalyzeBimbam function calls to include "sigmasq".

Added variable declaration for "sigmasq" matrix.

Added sigmasq variable calls for other places.

Converted "resid" from vector<double to vector<vector<double>>

Removed the (I think) unnecessary input of U within the "CalcLmmVgVeBeta" function call.

Updated residual variance file reading call from ReadFile_column() to ReadFile_resid().

Bug fixes after session with Artyom

Big fixes after session with Artyom

Bug fixes after session with Artyom

Major bug fixes after session with Artyom

Changed ReadFile_resid's parameter expectations.

Redefined ReadFile_resid() to expect a gsl_matrix instead of <vector<vector<double>> for resid. This seems to be a consistent sticking point.

Changed ReadFile_resid so that gsl_matrix operations are being used for resid instead of vector<vector<double>> operations.

Threw out the following additional definition of ReadFile_resid: void ReadFile_resid(const string &file_resid, bool &error, gsl_matrix *sigmasq) { debug_msg("entered"); igzstream infile(file_resid.c_str(), igzstream::in); if (!infile) { cout << "error! fail to open the residual variance file: " << file_resid << endl; error = true; return; } size_t n_row = sigmasq->size1, n_col = sigmasq->size2, i_row = 0, i_col = 0; gsl_matrix_set_zero(sigmasq); //change this so that sigmasq = V_e somehow string line; char *ch_ptr; double d; while (getline(infile, line)) { if (i_row == n_row) { cout << "error! number of rows in the residual variance file is larger " << "than expected." << endl; error = true; } i_col = 0; ch_ptr = strtok((char *)line.c_str(), " ,\t"); while (ch_ptr != NULL) { if (i_col == n_col) { cout << "error! number of columns in the residual variance file " << "is larger than expected, for row = " << i_row << endl; error = true; } d = atof(ch_ptr); gsl_matrix_set(sigmasq, i_row, i_col, d); i_col++; ch_ptr = strtok(NULL, " ,\t"); } i_row++; } infile.close(); infile.clear(); return; }

Added additional ReadFile_resid call back, changing mentions of sigmasq back to resid.

Edited user documentation for "-resid" option.

Changed mention of "resid" to "sigmasq" within the function "CheckResid"

Created new version of "CheckResid" function that hopefully helps in debugging the empty "sigmasq" issue and commented out seemingly redundant "CheckResid" function initiation.

Restored ReadFile_resid to its original, because the other caused compile errors.

Jameson Blount and others added 20 commits March 20, 2023 09:52

Modified files for strain-specific residual variance option

7c70ec5

Fixed 'cvt' variable and missing 0 < geno < 2 checks

ca93fb2

Fixed errors except for MphInitial errors

763a0b5

Fixed parameter names, attempted to integrate residual variance as an…

18e15e3

… option similar to cvt

More bug fixes

72f5a24

Added PARAM 'CheckResidvar' so that if no vector is supplied, it is r…

f76326b

…eplaced in the math with the null Ve

Update gemma.cpp

58e370f

Small formatting change to -residvar

Update mvlmm.cpp

924daa5

Updated resid_var math to CalcQi

Update mvlmm.cpp

db72327

Updated most functions except MphEM and remaining P-val functions

Update mvlmm.cpp

363ec00

Continued to make appropriate math changes to MphCalcSigma, UpdateV, MphCalcLogL, MphEM, MphCalcP, MphCalcBeta, and other functions.

Update mvlmm.cpp

c40699b

Finished theoretical changes in other functions and replaced all "eps_eval" entries with "eval_vec, sigmasq".

Update gemma.cpp

96d5f17

Replaced any lingering mentions of "eps_eval" with "eval_vec, sigmasq".

Update gemma.cpp

112d499

Update gemma_io.cpp

6f1d2cd

Changed mentions of "eps_eval" with "sigmasq". Still need to figure out how to get eigenvectors from the same place "eval" is derived.

Update gemma_io.h

057d27a

Changed mention of eps_eval as a vector to "gsl_matrix *sigmasq"

Update lmm.cpp

a8352e8

Changed the one mention of "eps_eval" to "eval_vec" and "sigmasq". May actually need to return to this file in order for the code to actually run.

Update mvlmm.h

f858b30

eps_eval --> eval_vec, sigmasq

Update param.cpp

587c181

Changed eps_eval to sigmasq and corrected any identified regions where dimensions did not match.

Update param.h

e85449d

Update param.h

ad54957

Replaced input of CopyResid from eps_eval to sigmasq

aartaka reviewed Oct 20, 2024

View reviewed changes

jb621-star added 9 commits October 22, 2024 12:13

Update param.h

b0da2e4

Changed "file_residvar" back to "file_resid"

Update gemma.cpp

9580e68

Transferred all instances of "eval_vec" to "U".

Update gemma_io.cpp

1b9c986

Changed all instances of "_residvar" to "resid" to keep things consistent.

Update gemma_io.h

e193600

Changed all instances of "residvar" to "resid" to keep things consistent

Update lmm.cpp

8409e49

Changed inputs to "CalcLmmVgVeBeta", "CalcLambda", and "AnalyzePlink" according to their use elsewhere (i.e. no more "eval_vec").

Update param.h

06b5f62

Changed residvar to resid for consistency.

Update param.cpp

be69a55

Changed all mentions of residvar to resid for consistency

Update gemma.cpp

2b25992

Changed all residvar to resid for consistency.

Update mvlmm.h

e1d92d0

Changed "residvar" to "resid" and removed unnecessary mentions of "eval_vec" from other functions.

jb621-star added 30 commits October 28, 2024 16:09

Update mvlmm.cpp

7347f95

Removed repeat parameter mentions of eigenvectors U in Analyze functions

Update mvlmm.cpp

6c284c4

Fixed examples of mishandled data types from when defining "scalar" data type to: double scalar = gsl_matrix_get(sigmasq, 0, 0) / ve; gsl_matrix_scale(Sigma, scalar); double epsilon = gsl_matrix_get(Sigma, k, k);

Update mvlmm.cpp

ed2db02

Fixed remaining cases of how the "scalar" variable is defined, and how many other functions are attributed.

Update mvlmm.cpp

8b43a6b

Update mvlmm.cpp

95dace7

Cleaned up various parameter specifications within function definitions.

Update param.cpp

d397275

Update mvlmm.cpp

e3768cb

Corrected CalcLambda and CalcLmmVgVeBeta mentions that were actually from lmm no mvlmm

Update mvlmm.cpp

32f4188

sigmasq needs to be in the AnalyzeBimbam and AnalyzePlink function definitions.

Update mvlmm.h

baff409

const gsl_matrix *sigmasq, needs to be in these functions due to MphEM calls within other functions in mvlmm.cpp.

Update mvlmm.cpp

ad18a8d

Removed duplicate mention of U within AnalyzePlinkGXE

Update gemma.cpp

b4a4641

Updated AnalyzePlink and AnalyzeBimbam function calls to include "sigmasq".

Update gemma.cpp

b198ae6

Added variable declaration for "sigmasq" matrix.

Update gemma.cpp

56fc6c5

Added sigmasq variable calls for other places.

Update param.h

f8764ae

Converted "resid" from vector<double to vector<vector<double>>

Update mvlmm.cpp

dd20afe

Removed the (I think) unnecessary input of U within the "CalcLmmVgVeBeta" function call.

Update param.cpp

c75b7d2

Updated residual variance file reading call from ReadFile_column() to ReadFile_resid().

Update gemma.cpp

882c9cd

Bug fixes after session with Artyom

Update mvlmm.cpp

f6e9299

Big fixes after session with Artyom

Update mvlmm.h

73412d2

Bug fixes after session with Artyom

Update param.h

4bfc619

Bug fixes after session with Artyom

Update param.cpp

2516635

Major bug fixes after session with Artyom

Update gemma_io.h

0876bf0

Changed ReadFile_resid's parameter expectations.

Update gemma_io.h

66a5a37

Redefined ReadFile_resid() to expect a gsl_matrix instead of <vector<vector<double>> for resid. This seems to be a consistent sticking point.

Update gemma_io.cpp

b6bfb79

Changed ReadFile_resid so that gsl_matrix operations are being used for resid instead of vector<vector<double>> operations.

Update gemma_io.cpp

ea678d7

Added additional ReadFile_resid call back, changing mentions of sigmasq back to resid.

Update gemma.cpp

40b31bc

Edited user documentation for "-resid" option.

Update param.cpp

bbb458b

Changed mention of "resid" to "sigmasq" within the function "CheckResid"

Update gemma_io.cpp

84cfda9

Created new version of "CheckResid" function that hopefully helps in debugging the empty "sigmasq" issue and commented out seemingly redundant "CheckResid" function initiation.

Update gemma_io.cpp

65dab2e

Restored ReadFile_resid to its original, because the other caused compile errors.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Modified files for strain-specific residual variance option #269

Modified files for strain-specific residual variance option #269

jb621-star commented Mar 20, 2023

aartaka left a comment

aartaka Oct 20, 2024

aartaka Oct 20, 2024

aartaka Oct 20, 2024

aartaka Oct 20, 2024

Modified files for strain-specific residual variance option #269

Are you sure you want to change the base?

Modified files for strain-specific residual variance option #269

Conversation

jb621-star commented Mar 20, 2023

aartaka left a comment

Choose a reason for hiding this comment

aartaka Oct 20, 2024

Choose a reason for hiding this comment

aartaka Oct 20, 2024

Choose a reason for hiding this comment

aartaka Oct 20, 2024

Choose a reason for hiding this comment

aartaka Oct 20, 2024

Choose a reason for hiding this comment