UW-GAC / GENESIS

GENetic EStimation and Inference in Structured samples (GENESIS): Statistical methods for analyzing genetic data from samples with population structure and/or relatedness
https://bioconductor.org/packages/GENESIS
34 stars 13 forks source link

pcrelate issue for large sample size #70

Closed GraceSheng closed 3 years ago

GraceSheng commented 3 years ago

I was running 76K samples for 33,765 SNPs, using PLINK bed/bim/fam files. PCrelate failed with the following error message:

Error in rbindlist(l, use.names, fill, idcol) : 
  Total rows in the list is 2165981000 which is larger than the maximum number of rows, currently 2147483647
Calls: pcrelate ... pcrelate -> .local -> .pcrelate -> rbind -> rbind -> rbindlist
In addition: Warning message:
In .pcrelate(gdsobj, pcs = pcs, scale = scale, ibd.probs = ibd.probs,  :
  small.samp.correct can only be used when all samples are analyzed in one block and `scale != none`
Execution halted
run_PCair_PCrelate_1st_iter_V1_largemem.R failed

It appears to hit the maximum row size (2^31-1) for a R data.table. Does this mean we cannot run pcrelate for more than ~65K samples? Any advice or insights? Thank you!

smgogarten commented 3 years ago

We have run pcrelate on larger numbers of samples by calling its sub-functions separately, which also allows for running sample blocks in parallel on a compute cluster; see #38 for a description of how to do it. This allows filtering of the results for each block to only keep pairs above a certain kinship threshold before combining them.

We should really add this threshold option to the pcrelate function so it can be done more easily, but it might be a while before we have time to do that.

GraceSheng commented 3 years ago

Thank you Stephanie for suggesting the workaround. As I intend to keep all the pairwise kc rows, I tried an alternative approach: basically I use a list of data table to hold kinBtwn and write kinBtwn out to a file after the pcrelate call.

I browsed the post processing codes for kinBtwn (skipped the ones that are applicable to small sample size), both "correct k2 for HW departure" and "correctK0" appear to be doing row based corrections which does not need info from all the other pairs in the table. So I just run each of my kinBtwn data table through these corrections.

I tested this on a smaller dataset with both standard pcrelate and my custom pcrelate, their kinBtwn table looked identical after resorting by ID1 and ID2. Though my testing might not be sufficient.

Since my knowledge of pcrelate is very limited, is there any pitfall that I might have overlooked? I plan to only use mine when sample size is > ~65K.

Thank you for your time!

mconomos commented 3 years ago

The correctK0 function is certainly fine to be called by block as it operates only on each row. The correctK2 function is also fine, as long as you also have the estimated kinSelf data loaded for each sample in the pair -- the f estimates from kinSelf are used as part of the correction.

GraceSheng commented 3 years ago

Thank you for confirming. yes, the full kinSelf data table is always passed to correctK2 which should cover each sample in the pair for each of my kinBtwn data tables.