#########################################################################################
#											#
#                Plink									#
#											#
#########################################################################################


# sibs need to be 'nonfounders' for QFAM to work, thus putting in dummy "10" and "11" pat/mat codes in PED
awk ' { print $1, $2,"10 11", $5, $6, $7, $8, $9, $10, $11, $12, $13, $14, $15, $16 }  ' AP_7m_snp.ped > f.ped
cp AP_7m_snp.map f.map

# create 'cluster' file for pairs:
awk ' { print $1,$2,$1 } ' f.ped > f.clst

# Basic linear regression
plink --file f --linear --out reg1

# Huber-White sandwich variance estimator (based on twin clusters)
plink --file f --linear --within f.clst --out reg2

# QFAM: within family
plink --file f --qfam --mperm 5000 --out qfam1

# QFAM: between family
plink --file f --qfam-between --mperm 5000 --out qfam2

# QFAM: total test
plink --file f --qfam-total --mperm 5000 --out qfam3

# conclusions

# p-values get slightly less interesting after adjusting
cat reg1.assoc.linear
cat reg2.assoc.linear

# the qfam-total is v. similar to the HW-adjusted estimates, that makes sense 
cat qfam1.qfam.within.perm
cat qfam2.qfam.between.perm
cat qfam3.qfam.total.perm
cat reg2.assoc.linear