import numpy as np
from snipar.utilities import make_id_dict
from snipar.pedigree import find_individuals_with_sibs
from snipar.gtarray import gtarray
[docs]def get_indices_given_ped(ped, gts_ids, imp_fams=None, ids=None, sib=False, verbose=False):
"""
Used in get_gts_matrix_given_ped to get the ids of individuals with observed/imputed parental genotypes and, if sib=True, at least one genotyped sibling.
It returns those ids along with the indices of the relevant individuals and their first degree relatives in the observed genotypes (observed indices),
and the indices of the imputed parental genotypes for those individuals.
"""
# Made dictionary for observed genotypes
gts_id_dict = make_id_dict(gts_ids)
# If IDs not provided, use all individuals with observed genotypes
if ids is None:
ids = gts_ids
# Find individuals with genotyped siblings
if sib:
# Look in full genotype sample in case some genotyped sibs are not in ids
ids = gts_ids
ids = find_individuals_with_sibs(ids, ped, gts_ids, return_ids_only=True)
if verbose:
print('Found ' + str(ids.shape[0]) + ' individuals with genotyped siblings')
### Find parental status
if verbose:
print('Checking for observed/imputed parental genotypes')
par_status, gt_indices, fam_labels = find_par_gts(ids, ped, gts_id_dict, imp_fams=imp_fams)
# Find which individuals can be used
none_missing = np.min(gt_indices, axis=1) >= 0
N = np.sum(none_missing)
if N == 0:
raise ValueError(
'No individuals with phenotype observations and complete observed/imputed genotype observations')
# Take those that can be used
gt_indices = gt_indices[none_missing, :]
par_status = par_status[none_missing, :]
ids = ids[none_missing]
parcount = np.sum(par_status==0,axis=1)
if verbose:
print(str(N) + ' individuals with phenotype observations and complete observed/imputed genotype observations')
print(str(np.sum(parcount==0))+' individuals with imputed but no observed parental genotypes')
print(str(np.sum(parcount==1))+' individuals with one observed and one imputed parent')
print(str(np.sum(parcount==2))+' individuals with both parents observed')
# Find indices of individuals and their parents in observed genotypes
observed_indices = np.sort(np.unique(np.hstack((gt_indices[:, 0],
gt_indices[par_status[:, 0] == 0, 1],
gt_indices[par_status[:, 1] == 0, 2]))))
# Get indices of imputed parents
imp_indices = np.sort(np.unique(np.hstack((gt_indices[par_status[:, 0] == 1, 1],
gt_indices[par_status[:, 1] == 1, 2]))))
# Return ids with imputed/observed parents
return ids, observed_indices, imp_indices, parcount
[docs]def find_par_gts(pheno_ids, ped, gts_id_dict, imp_fams=None):
"""
Used in get_gts_matrix to find whether individuals have imputed or observed parental genotypes, and to
find the indices of the observed/imputed parents in the observed/imputed genotype arrays.
'par_status' codes whether an individual has parents that are observed or imputed or neither.
'gt_indices' records the relevant index of the parent in the observed/imputed genotype arrays
'fam_labels' records the family of the individual based on the pedigree
"""
# Whether mother and father have observed/imputed genotypes
par_status = np.zeros((pheno_ids.shape[0],2),dtype=int)
par_status[:] = -1
# Indices of obsered/imputed genotypes in relevant arrays
gt_indices = np.zeros((pheno_ids.shape[0],3),dtype=int)
gt_indices[:] = -1
## Build dictionaries
# Where each individual is in the pedigree
ped_dict = make_id_dict(ped,1)
# Where the imputed data is for each family
if imp_fams is not None:
fam_dict = make_id_dict(imp_fams)
# Store family ID of each individual
fam_labels = np.zeros((pheno_ids.shape[0]),dtype=ped.dtype)
# Find status and find indices
for i in range(0,pheno_ids.shape[0]):
# Find index in genotypes
if pheno_ids[i] in gts_id_dict:
gt_indices[i,0] = gts_id_dict[pheno_ids[i]]
# Find index in pedigree
if pheno_ids[i] in ped_dict:
ped_i = ped[ped_dict[pheno_ids[i]], :]
fam_labels[i] = ped_i[0]
# Check for observed father
if ped_i[2] in gts_id_dict:
gt_indices[i,1] = gts_id_dict[ped_i[2]]
par_status[i,0] = 0
# Check for observed mother
if ped_i[3] in gts_id_dict:
gt_indices[i, 2] = gts_id_dict[ped_i[3]]
par_status[i,1] = 0
# If parent not observed, look for imputation
if imp_fams is not None:
if ped_i[0] in fam_dict:
imp_index = fam_dict[ped_i[0]]
# Check if this is imputation of father, or mother, or both
if ped_i[4] == 'False' and not par_status[i,0] == 0:
gt_indices[i, 1] = imp_index
par_status[i, 0] = 1
if ped_i[5] == 'False' and not par_status[i,1] == 0:
gt_indices[i, 2] = imp_index
par_status[i, 1] = 1
return par_status, gt_indices, fam_labels
[docs]def make_gts_matrix(gts, par_status, gt_indices, imp_gts=None, parsum = False):
"""
Used in get_gts_matrix to construct the family based genotype matrix given
observed/imputed genotypes. 'gt_indices' has the indices in the observed/imputed genotype arrays;
and par_status codes whether the parents are observed (0) or imputed (1).
"""
if np.min(gt_indices)<0:
raise(ValueError('Missing genotype index'))
if imp_gts is None:
if np.max(par_status)==1:
raise(ValueError('No imputed parental genotypes provided'))
N = gt_indices.shape[0]
if parsum:
gdim = 2
else:
gdim = 3
G = np.zeros((N,gdim,gts.shape[1]),np.float32)
# Proband genotypes
G[:,0,:] = gts[gt_indices[:,0],:]
# Paternal genotypes
G[par_status[:, 0] == 0, 1 ,:] = gts[gt_indices[par_status[:, 0] == 0, 1], :]
if imp_gts is not None:
G[par_status[:, 0] == 1, 1, :] = imp_gts[gt_indices[par_status[:, 0] == 1, 1], :]
# Maternal genotypes
if parsum:
G[par_status[:, 1] == 0, 1, :] += gts[gt_indices[par_status[:, 1] == 0, 2], :]
if imp_gts is not None:
G[par_status[:, 1] == 1, 1, :] += imp_gts[gt_indices[par_status[:, 1] == 1, 2], :]
else:
G[par_status[:, 1] == 0, 2, :] = gts[gt_indices[par_status[:, 1] == 0, 2], :]
if imp_gts is not None:
G[par_status[:, 1] == 1, 2, :] = imp_gts[gt_indices[par_status[:, 1] == 1, 2], :]
return G
[docs]def get_fam_means(ids,ped,gts,gts_ids,remove_proband = True, return_famsizes = False):
"""
Used in get_gts_matrix to find the mean genotype in each sibship (family) for each SNP or for a PGS.
The gtarray that is returned is indexed based on the subset of ids provided from sibships of size 2 or greater.
If remove_proband=True, then the genotype/PGS of the index individual is removed from the fam_mean given for that individual.
"""
ids, ids_fams, gts_fams = find_individuals_with_sibs(ids, ped, gts_ids)
fams = np.unique(ids_fams)
fams_dict = make_id_dict(fams)
# Compute sums of genotypes in each family
fam_sums = np.zeros((fams.shape[0],gts.shape[1]),dtype=gts.dtype)
fam_counts = np.zeros((fams.shape[0]),dtype=int)
for i in range(0,fams.shape[0]):
fam_indices = np.where(gts_fams==fams[i])[0]
fam_sums[i,:] = np.sum(gts[fam_indices,:],axis=0)
fam_counts[i] = fam_indices.shape[0]
# Place in vector corresponding to IDs
if remove_proband:
gts_id_dict = make_id_dict(gts_ids)
G_sib = np.zeros((ids.shape[0],gts.shape[1]),dtype = np.float32)
for i in range(0,ids.shape[0]):
fam_index = fams_dict[ids_fams[i]]
G_sib[i,:] = fam_sums[fam_index,:]
n_i = fam_counts[fam_index]
if remove_proband:
G_sib[i,:] = G_sib[i,:] - gts[gts_id_dict[ids[i]],:]
n_i = n_i-1
G_sib[i,:] = G_sib[i,:]/float(n_i)
if return_famsizes:
return [gtarray(G_sib, ids),fam_counts,fam_sums]
else:
return gtarray(G_sib,ids)