#########################################################################################
#											#
#                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


#########################################################################################
#											#
#                Merlin									#
#											#
#########################################################################################


# make a zygosity file for the twins MZ=1 DZ=2
awk ' { print  $1, $2, "10 11", $3, $4}  ' twzyg.dat > zyg.ped
# add in the parents
awk ' { print  $1 }' zyg.ped | sort | uniq > temp
awk ' { print  $1,  "10 0 0 1 0"}  ' temp >> zyg.ped
awk ' { print  $1,  "11 0 0 2 0"}  ' temp >> zyg.ped

# make a dat file for zyg.ped
echo 'Z Zygosity' > zyg.dat

# make a dat file for f.ped
echo 'T Trait' > f.dat
awk '{ print "M", $2}' AP_7m_snp.map >>f.dat

# make a map file in merlin format - CHR SNP cM (BP/1M)
awk '{print "0", $2, $4/1000000}' AP_7m_snp.map > merlin.map

# run association - score test
merlin -p zyg.ped,f.ped -d zyg.dat,f.dat -m merlin.map  --singlepoint --trim --fastAssoc --tabulate --prefix fastAssoc

# run association - Likelihood-ratio test
merlin -p zyg.ped,f.ped -d zyg.dat,f.dat -m merlin.map  --singlepoint --trim --assoc --tabulate --prefix assoc