jianyangqt / gcta

GCTA software
GNU General Public License v3.0
87 stars 26 forks source link

Abnormal memory requirements while making GRMs #91

Open cjfei18 opened 3 months ago

cjfei18 commented 3 months ago

Dear GCTA team, I've encountered an issue when making GRMs during the GREML analysis. I was using the 1.94.0 beta version of GTCA and retried with 1.94.1 version but the problem persisted. My code is as follows:

for i in {1..600}; do
gcta64 \
--mbfile plink_file.list \
--make-grm-alg 1 \
--make-grm-part 600 ${i} \
--threads 128 \
--extract maf4_ld1_snplist.txt \
--out step3/maf4_ld1 
done

It worked well for i from 1 to 598, but a segmentation fault occurred when i=599. The log file is as follows:

Reading PLINK BIM file from [chr21.bim]...
37339960 SNPs to be included from BIM file(s).
Reading PLINK BIM file from [chr22.bim]...
37922580 SNPs to be included from BIM file(s).
Get 12166065 SNPs from list [maf4_ld1_snplist.txt].
After extracting SNP, 12166065 SNPs remain.
Error: can't allocate enough memory to store the (parted) GRM: 10915231.959338GB required.
An error occurs, please check the options or data

It requested 10915231.959338GB of memory to store the GRM, which seems unreasonable. The same error occurred for i=600.
A similar segmentation error was reported in issue #26 and I'm unsure if this is the same issue. Could you please help me resolve this problem? Thank you in advance for your time and assistance!

longmanz commented 3 months ago

Hi, Thank you for reporting this. This indeed seems to be a bug. I would suggest a few things to see if we can workaround it:

  1. Use a smaller set of SNPs to compute the GRM. Currently you are using m = 12166065 which is a lot. You may consider using HapMap3 SNPs only (~1m SNPs). Usually using common HapMap3 SNPs should be good enough to reconstruct the relationship matrix.
  2. You may consider parallelize your jobs instead of doing it in a loop.
  3. Using a smaller number for the total part and smaller umber of threads. Even if you are computing GRM for UK Biobank, you should not need a total of 600 parts and 128 threads. Usually ~10 threads per part and 200 parts should be enough (given you are parallelizing your jobs).
cjfei18 commented 3 months ago

Hi, Thank you for your quick response and suggestions! Because we are doing a WGS analysis and want rare variants to be covered, we used a large set of SNPs, which is similar in magnitude to the number of SNPs used in Wainschtein et al. (PMID: 35256806). Due to the large number of SNPs, we split the analysis into many parts and set a higher number of threads. I apologize for the confusion caused by the for i in {1..600}; do line in my code; I actually ran the jobs in parallel rather than using just one loop. I tried setting 10 threads as you suggested, but I still encountered the same error. However, when I divided the jobs into 300 parts with 10 threads each, the 300th job ran successfully. I am wondering if it would be acceptable to use the 300th job from the 300 parts instead of the 599th and 600th jobs from the 600 parts? Thank you again for your guidance!

longmanz commented 3 months ago

Hi, Thank you for your information and clarification. And thank you for checking the scenario of using 300 parts. Unfortunately, you will not be able to merge the 300th part instead of the 599 + 600 parts. You will need to re-run the 300 parts altogether. Sorry about that. But the good news is that this strategy can workaround the bug you previously encountered. So I would suggest you re-running the job with 300 parts.