Skip to content

Commit

Permalink
--r[2]-[un]phased, phased-r2 log-likelihood bugfix, plink 1.9 cm-wind…
Browse files Browse the repository at this point in the history 
…ow bugfix
  • Loading branch information
@chrchang
chrchang committed Dec 12, 2023
1 parent db19b76 commit 6678559
Show file tree
Hide file tree
Showing 6 changed files with 863 additions and 263 deletions.
6 changes: 3 additions & 3 deletions 1.9/plink.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,7 +93,7 @@

static const char ver_str[] =
#ifdef STABLE_BUILD
"PLINK v1.90b7.1"
"PLINK v1.90b7.2"
#else
"PLINK v1.90p"
#endif
Expand All @@ -105,10 +105,10 @@ static const char ver_str[] =
#else
" 32-bit"
#endif
" (5 Nov 2023)";
" (11 Dec 2023)";
static const char ver_str2[] =
// include leading space if day < 10, so character length stays the same
" "
""
#ifdef STABLE_BUILD
" " // (don't want this when version number has two trailing digits)
#else
Expand Down
3 changes: 2 additions & 1 deletion 1.9/plink_ld.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5516,7 +5516,8 @@ int32_t ld_report_dprime(pthread_t* threads, Ld_info* ldip, FILE* bedfile, uintp
}
} else {
cur_marker_cm = marker_cms[marker_uidx1_tmp] + window_cm;
while ((marker_pos[marker_uidx2_fwd2] <= cur_marker_pos) && (marker_cms[marker_uidx2_fwd2] <= window_cm)) {
// bugfix (11 Dec 2023)
while ((marker_pos[marker_uidx2_fwd2] <= cur_marker_pos) && (marker_cms[marker_uidx2_fwd2] <= cur_marker_cm)) {
marker_uidx2_fwd = marker_uidx2_fwd2;
window_lead_ct++;
if (marker_uidx2_fwd == chrom_last) {
Expand Down
4 changes: 2 additions & 2 deletions 2.0/plink2.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ static const char ver_str[] = "PLINK v2.00a6"
#elif defined(USE_AOCL)
" AMD"
#endif
" (23 Nov 2023)";
" (11 Dec 2023)";
static const char ver_str2[] =
// include leading space if day < 10, so character length stays the same
""
Expand Down Expand Up @@ -7239,7 +7239,7 @@ int main(int argc, char** argv) {
if (unlikely(EnforceParamCtRange(argvk[arg_idx], param_ct, 1, 1))) {
goto main_ret_INVALID_CMDLINE_2A;
}
reterr = AllocFname(argvk[arg_idx], flagname_p, 0, &pc.vcor_info.ld_snp_list_fname);
reterr = AllocFname(argvk[arg_idx + 1], flagname_p, 0, &pc.vcor_info.ld_snp_list_fname);
if (unlikely(reterr)) {
goto main_ret_1;
}
Expand Down
65 changes: 9 additions & 56 deletions 2.0/plink2_common.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -304,14 +304,17 @@ void PopulateRescaledDosageF(const uintptr_t* genoarr, const uintptr_t* dosage_p

void DosagePhaseinfoPatch(const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dosage_vec, const uintptr_t* dphase_present, uint32_t sample_ct, SDosage* dphase_delta) {
const uint32_t sample_ctl = BitCtToWordCt(sample_ct);
uintptr_t dosage_present_word = 0;
for (uint32_t widx = 0; widx != sample_ctl; ++widx) {
uintptr_t phasepresent_nodphase_word = phasepresent[widx];
if (dphase_present) {
phasepresent_nodphase_word &= ~dphase_present[widx];
}
if (phasepresent_nodphase_word) {
const uintptr_t phaseinfo_word = phaseinfo[widx];
const uintptr_t dosage_present_word = dosage_present[widx];
if (dosage_present) {
dosage_present_word = dosage_present[widx];
}
const uint32_t sample_idx_offset = widx * kBitsPerWord;
do {
const uint32_t sample_idx_lowbits = ctzw(phasepresent_nodphase_word);
Expand All @@ -333,7 +336,7 @@ void DosagePhaseinfoPatch(const uintptr_t* phasepresent, const uintptr_t* phasei
}
}

void PopulateDenseDphase(const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dense_dosage_vec, const uintptr_t* dphase_present, const SDosage* dphase_delta, uint32_t sample_ct, uint32_t phasepresent_ct, uint32_t dphase_ct, SDosage* dense_dphase_delta) {
void PopulateDenseDphase(const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dense_dosage_vec, const uintptr_t* dphase_present, const SDosage* dphase_delta, uint32_t sample_ct, uint32_t phasepresent_ct, uint32_t dosage_ct, uint32_t dphase_ct, SDosage* dense_dphase_delta) {
const uint32_t dosagev_ct = DivUp(sample_ct, kDosagePerVec);
ZeroDphaseArr(dosagev_ct * kDosagePerVec, dense_dphase_delta);
if (dphase_ct) {
Expand All @@ -348,60 +351,7 @@ void PopulateDenseDphase(const uintptr_t* phasepresent, const uintptr_t* phasein
dphase_present = nullptr;
}
if (phasepresent_ct) {
DosagePhaseinfoPatch(phasepresent, phaseinfo, dosage_present, dense_dosage_vec, dphase_present, sample_ct, dense_dphase_delta);
}
}

void DosagePhaseinfoPatchSubset(const uintptr_t* sample_include, const uint32_t* sample_include_cumulative_popcounts, const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dosage_vec, const uintptr_t* dphase_present, uint32_t raw_sample_ct, SDosage* dphase_delta) {
const uint32_t raw_sample_ctl = BitCtToWordCt(raw_sample_ct);
for (uint32_t widx = 0; widx != raw_sample_ctl; ++widx) {
const uintptr_t sample_include_word = sample_include[widx];
uintptr_t phasepresent_nodphase_word = phasepresent[widx] & sample_include_word;
if (dphase_present) {
phasepresent_nodphase_word &= ~dphase_present[widx];
}
if (phasepresent_nodphase_word) {
const uintptr_t phaseinfo_word = phaseinfo[widx];
const uintptr_t dosage_present_word = dosage_present[widx];
const uint32_t sample_idx_offset = sample_include_cumulative_popcounts[widx];
do {
const uintptr_t shifted_bit = phasepresent_nodphase_word & (-phasepresent_nodphase_word);
const uint32_t sample_idx = sample_idx_offset + PopcountWord(sample_include_word & (shifted_bit - 1));
int32_t cur_diff = kDosageMid;
if (dosage_present_word & shifted_bit) {
cur_diff = DosageHomdist(dosage_vec[sample_idx]);
}
// todo: verify the compiler optimizes this well
if (!(phaseinfo_word & shifted_bit)) {
cur_diff = -cur_diff;
}
dphase_delta[sample_idx] = cur_diff;
phasepresent_nodphase_word ^= shifted_bit;
} while (phasepresent_nodphase_word);
}
}
}

void PopulateDenseDphaseSubset(const uintptr_t* sample_include, const uint32_t* sample_include_cumulative_popcounts, const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dense_dosage_vec, const uintptr_t* dphase_present, const SDosage* dphase_delta, uint32_t raw_sample_ct, uint32_t sample_ct, uint32_t phasepresent_ct, uint32_t dphase_ct, SDosage* dense_dphase_delta) {
const uint32_t dosagev_ct = DivUp(sample_ct, kDosagePerVec);
ZeroDphaseArr(dosagev_ct * kDosagePerVec, dense_dphase_delta);
if (dphase_ct) {
uintptr_t widx = 0;
uintptr_t cur_bits = dphase_present[0];
for (uint32_t dphase_idx = 0; dphase_idx != dphase_ct; ++dphase_idx) {
const uintptr_t shifted_bit = BitIter1y(dphase_present, &widx, &cur_bits);
const uintptr_t sample_include_word = sample_include[widx];
if (sample_include_word & shifted_bit) {
const uint32_t sample_idx = sample_include_cumulative_popcounts[widx] + PopcountWord(sample_include_word & (shifted_bit - 1));
const int32_t dphase_delta_val = dphase_delta[dphase_idx];
dense_dphase_delta[sample_idx] = dphase_delta_val;
}
}
} else {
dphase_present = nullptr;
}
if (phasepresent_ct) {
DosagePhaseinfoPatchSubset(sample_include, sample_include_cumulative_popcounts, phasepresent, phaseinfo, dosage_present, dense_dosage_vec, dphase_present, raw_sample_ct, dense_dphase_delta);
DosagePhaseinfoPatch(phasepresent, phaseinfo, dosage_ct? dosage_present : nullptr, dense_dosage_vec, dphase_present, sample_ct, dense_dphase_delta);
}
}

Expand Down Expand Up @@ -4252,6 +4202,9 @@ void GetBpWindow(const uintptr_t* variant_include, const uint32_t* variant_bps,

PgenGlobalFlags GflagsVfilter(const uintptr_t* variant_include, const unsigned char* vrtypes, uint32_t raw_variant_ct, PgenGlobalFlags input_gflags) {
PgenGlobalFlags read_gflags = kfPgenGlobal0;
if (!vrtypes) {
return input_gflags;
}
const uintptr_t* vrtypes_alias = R_CAST(const uintptr_t*, vrtypes);
const uint32_t raw_variant_ctl = BitCtToWordCt(raw_variant_ct);
uint32_t mask_multiply = ((input_gflags & kfPgenGlobalHardcallPhasePresent)? 0x10 : 0) + ((input_gflags & kfPgenGlobalDosagePresent)? 0x60 : 0) + ((input_gflags & kfPgenGlobalDosagePhasePresent)? 0x80 : 0);
Expand Down
4 changes: 1 addition & 3 deletions 2.0/plink2_common.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -429,9 +429,7 @@ void PopulateRescaledDosage(const uintptr_t* genoarr, const uintptr_t* dosage_pr

void PopulateRescaledDosageF(const uintptr_t* genoarr, const uintptr_t* dosage_present, const Dosage* dosage_main, float slope, float intercept, float missing_val, uint32_t sample_ct, uint32_t dosage_ct, float* expanded_dosages);

void PopulateDenseDphase(const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dense_dosage_vec, const uintptr_t* dphase_present, const SDosage* dphase_delta, uint32_t sample_ct, uint32_t phasepresent_ct, uint32_t dphase_ct, SDosage* dense_dphase_delta);

void PopulateDenseDphaseSubset(const uintptr_t* sample_include, const uint32_t* sample_include_cumulative_popcounts, const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dense_dosage_vec, const uintptr_t* dphase_present, const SDosage* dphase_delta, uint32_t raw_sample_ct, uint32_t sample_ct, uint32_t phasepresent_ct, uint32_t dphase_ct, SDosage* dense_dphase_delta);
void PopulateDenseDphase(const uintptr_t* phasepresent, const uintptr_t* phaseinfo, const uintptr_t* dosage_present, const Dosage* dense_dosage_vec, const uintptr_t* dphase_present, const SDosage* dphase_delta, uint32_t sample_ct, uint32_t phasepresent_ct, uint32_t dosage_ct, uint32_t dphase_ct, SDosage* dense_dphase_delta);

// assumes trailing bits of genoarr are zeroed out
HEADER_INLINE uint32_t AtLeastOneHetUnsafe(const uintptr_t* genoarr, uint32_t sample_ct) {
Expand Down
Loading

0 comments on commit 6678559

Please sign in to comment.