##################################################################################################################
##################################################################################################################
## PLINK script: recoding the data (for extraction later) and do GxE with dichotomous Exposure variable

plink --bfile Transferrin --chr 3 --from-kb 134840 --to-kb 135052 --recodeA --out  TFsnps

### Analysis using ALL SNPs

plink --bfile Transferrin --gxe --covar covD.dat --pheno Tr.pheno --out fullGxEanalysis

### Analysis using Only SNPs in TF gene

plink --bfile Transferrin --chr 3 --from-kb 134840 --to-kb 135052 --gxe --covar covD.dat --pheno Tr.pheno --out TFGxEanalysisD

### Try to run interaction analysis with quantitative Exposure

plink --bfile Transferrin --chr 3 --from-kb 134840 --to-kb 135052 --gxe --covar covC.dat --pheno Tr.pheno --out TFGxEanalysisC



##################################################################################################################
##################################################################################################################
### R Scripts : Reading data into R and getting rid of missing phenotypes

snps = read.table("TFsnps.raw", header = T) 
phenotypes = read.table("Tr.pheno", header = F)
envC = read.table("covC.dat", header = F)

missing.phenotypes = which(is.na(phenotypes[,3]))
Z = as.matrix(snps[-missing.phenotypes,7:30])
y = as.vector(phenotypes[-missing.phenotypes,3])
E = as.vector(envC[-missing.phenotypes,3])

sum(is.na(Z)) ## this gives us the number of SNPs that are missing

w.missing = which(apply(is.na(Z),1,sum) > 0) ## This figures out which subjects have ANY missing values

Z = Z[-w.missing,]  ## now we get rid of the subjects with any missing SNPs
y = y[-w.missing]   ## BE CAREFUL: because we are over-writing the variable, don't run this line more than once by accident.
E = E[-w.missing]

sum(is.na(Z))
sum(is.na(y)) ## Now we see that there is no missingness in either the Z or the y.


##################################################################################################################
##################################################################################################################
### R Scripts : Marginal Association Test for GxE 

pvs = rep(NA, ncol(Z))

for (j in 1:ncol(Z)){
    mod = lm(y~Z[,j]*E)
    pvs[j] = summary(mod)$coef[4,4] ### Need to make sure plucking off the correct value
}

pvs


##################################################################################################################
##################################################################################################################
### R Scripts : 2-df test

pvs.2df = rep(NA, ncol(Z))

for (j in 1:ncol(Z)) {
    mod0 = lm(y~E)
    mod1 = lm(y~Z[,j]*E)
    pvs.2df[j] = anova(mod1,mod0)$P[2]
}
pvs.2df


##################################################################################################################
##################################################################################################################
### R Scripts : Simulate outcome

set.seed(150)
n = length(y)
y.new = 5*(E^2) + rnorm(n, 0, 1)


pvs = rep(NA, ncol(Z))

for (j in 1:ncol(Z)){
    mod = lm(y.new~Z[,j]*E)
    pvs[j] = summary(mod)$coef[4,4] ### Need to make sure plucking off the correct value
}

pvs

##################################################################################################################
##################################################################################################################
## PLINK script: testing for gene-gene interactions

### Analysis using ALL SNPs

plink --bfile Transferrin --epistasis --pheno Tr_dich.pheno --out fullGxGanalysis


### Analysis using Only SNPs in TF gene

plink --bfile Transferrin --epistasis --chr 3 --from-kb 134840 --to-kb 135052 --pheno Tr_dich.pheno --out TFGxGanalysis


### Case Only analysis

plink --bfile Transferrin --fast-epistasis --case-only --chr 3 --from-kb 134840 --to-kb 135052 --pheno Tr_dich.pheno --out TFGxG_CO_analysis 

plink --bfile Transferrin --fast-epistasis --case-only --gap 1  --chr 3 --from-kb 134840 --to-kb 135052 --pheno Tr_dich.pheno --out TFGxG_CO_analysis