alglib::ap_error

[phil-grayson]

Hi again,

Once I worked through the previous issues with my test data, I ramped up to run a full analysis on 266 samples (140 case and 126 controls). Everything was running smoothly and it produced all expected output (including significant hits this time!), but I was printing stdout/stderr to a file and discovered the following in that file:

terminate called after throwing an instance of 'alglib::ap_error'
./runHawk: line 53: 33265 Aborted                 (core dumped) $hawkDir/log_reg_case.out -t $noThread $sex_confounder_arg -p $noPC > pvals_case_top.txt
terminate called after throwing an instance of 'alglib::ap_error'
./runHawk: line 53: 35031 Aborted                 (core dumped) $hawkDir/log_reg_control.out -t $noThread $sex_confounder_arg -p $noPC > pvals_control_top.txt

As I said, all expected output files were generated regardless of this error, but I'm concerned that those results might not be trustworthy if something failed during the run. Can you provide any insight on this?

I'm running on a Linux machine with 64 CPUs and 252 Gb and running the exact same code that previously worked (without error) for my ~20 sample test dataset. I haven't modified anything for threading and according to htop it looks like HAWK is running on 16 threads by default on this system.

Thanks!

[atifrahman]

It appears that our new cpp implementation has a bug that we didn't encounter while analyzing our test datasets. I'm guessing at some point the external library alglib had an underflow or something and crashed. So, there might been more significant hits to be found.

Would it be possible to share the following files? We can take a look pcs.evec gwas_info.txt gwas_eigenstratX.ind case_total_kmers.txt control_total_kmers.txt case_out_w_bonf_top.kmerDiff control_out_w_bonf_top.kmerDiff

In the meantime, it may be worth running the old R implementation below (the relevant portions from runHawk_old). You'll need to edit the R scripts to specify the confounders (sex, no. PCs). It'll take more time than the cpp version but shouldn't be too bad. Let me know if you need help modifying the R scripts... the parameter to makeCluster is the number of cores.

hawkDir=/home/atif/exp_hawk_tool/HAWK eigenstratDir=/home/atif/exp_hawk_tool/EIG6.0.1-Hawk isDiploid=1

noInd=$(cat sorted_files.txt | wc -l);

noPC=2 noThread=32 useSexCounfounder=False covFile=""

caseCount=$(cat case_sorted_files.txt | wc -l); controlCount=$(cat control_sorted_files.txt | wc -l);

Rscript $hawkDir/log_reg_case.R Rscript $hawkDir/log_reg_control.R

paste pvals_case_top.txt case_out_w_bonf_top.kmerDiff > pvals_case_top_merged.txt sort -g -k 1 -t $'\t' pvals_case_top_merged.txt > pvals_case_top_merged_sorted.txt

paste pvals_control_top.txt control_out_w_bonf_top.kmerDiff > pvals_control_top_merged.txt sort -g -k 1 -t $'\t' pvals_control_top_merged.txt > pvals_control_top_merged_sorted.txt

$hawkDir/convertToFasta_bf_correction.out

[phil-grayson]

Thanks for the follow-up! I submitted it to re-run when I posted my issue and it has completed with the same error again. Another thing that is worth mentioning given that you didn't ask for these files specifically is that the .fasta files that come out of this run have the pvalue as the labels, but the count instead of the ACTGs for the sequence. They look like this:

>0.000000e+00
6842
>0.000000e+00
5727
>0.000000e+00
5060

For the files, pcs.evec needed to be changed to a .txt to be uploaded. case_out and control_out are a few gb each (even gzipped), so it won't let me upload them here. Is there somewhere else I could put them? Thanks for looking into this for me!

case_total_kmers.txt control_total_kmers.txt pcs.evec.txt gwas_info.txt

[atifrahman]

Oh... interesting! Can you upload the large files somewhere (Google Drive or something like that) and share with me atif.bd@gmail.com?

Without your help, we wouldn't be able to debug our method. So, thank you!

[phil-grayson]

Thanks! It's on it's way!

[phil-grayson]

Just wanted to give this a bump. Based on the files I've sent, is it still reasonable to run the R scripts or is there a larger bug at play? Thanks!

[atifrahman]

Sorry, we had a busy week! We have found an error - we're getting an underflow from time to time but haven't figured out how to solve it with the C++ implementation.

Can you please let us know the parameters you want to use to run the C++ code (number of PCs to use as confounders, whether to use sex as a confounder and if there are any other confounder you want to add)? We'll run the R scripts ourselves. We need to do that to debug the C++ version anyway.

[phil-grayson]

Thank you very much for following up! I've attached my runHawk (had to change it to a .txt to upload). I was just using the default parameters (so it looks like 2 PCs?). As far as I know, all I've changed in that executable is hawkDir and eigenstratDir. This analysis is searching for sex differences, so we don't want sex as a confounder. Let me know if there's anything else you need for your testing! runHawk.txt