nievergeltlab

#K=pop prevalence
#P=proportion of cases in study
#hsq=Heritability estimate (on observed scale)
#bigT = liability threshold
#tau = density of gaussian
K=0.0659
P=0.0659
h2=0.0365
zv <- dnorm(qnorm(K))

nievergeltlab

for files in $(ls | grep log | grep -v munge | grep ldsc)
do
 h2snp=$(cat $files | grep "Observed scale" | sed 's/(//g' | sed 's/)//g' | awk '{print $5,$6, $5/$6}')
 intercept=$(cat $files | grep "Intercept" | sed 's/(//g' | sed 's/)//g' | awk '{print $2,$3, ($2-1)/$3}' )
 ratio=$(cat $files | grep "Ratio" | sed 's/(//g' | sed 's/)//g' | awk '{print $2,$3, $2/$3}' )
 lambda=$(cat $files | grep "Lambda" | sed 's/(//g' | sed 's/)//g' | awk '{print $3}' )
 chisq=$(cat $files | grep "Mean Chi" | sed 's/(//g' | sed 's/)//g' | awk '{print $3}' )
 outname=$(echo $files | sed 's/.log//g')

nievergeltlab

  BOLT-LMM_v2.3.2/bolt \
    --bed=pca/UKB_tinnitus_eur_unrelated_{1..22}.bed \
    --bim=pca/UKB_tinnitus_eur_unrelated_{1..22}.bim \
    --fam=pca/UKB_tinnitus_eur_unrelated_1.fam  \
    --LDscoresFile=BOLT-LMM_v2.3.2/tables/LDSCORE.1000G_EUR.tab.gz \
    --geneticMapFile=BOLT-LMM_v2.3.2/tables/genetic_map_hg19_withX.txt.gz \
    --lmmForceNonInf \
    --numThreads=6 \
    --bgenFile=/mnt/ukbb/adam/bgen/ukb_imp_chr{1..22}_v3.bgen \
    --bgenMinMAF=1e-3 \

nievergeltlab

d1$replicated <- NA
#All 3 must be true
for (i in 1:dim(d1)[1])
{
 d1[i,]$replicated <- any(d1[i,]$chr == d2$CHR & d1[i,]$pos >= d2$START & d1[i,]$pos <= d2$STOP)
 }
 

nievergeltlab

#Y is an Nx1 vector of the phenotype
#x is a NxN matrix of covariates
#b is beta (supplied by optim)
#d is a dispersion parameter (supplied by optim)
likfun1 <- function(y,x,theta)
{ 
 #the beta parameters will be the first N-1 inputs supplied by optim, dispersion the last
 b=theta[-length(theta)]
 d=theta[length(theta)]
 loglik1 <- sum ( log(gamma(y+d)) -log(gamma(y + 1)) - log(gamma(d)) +

nievergeltlab

for file in $(ls pheno | grep .pheno)
do
 fname=$(echo $file | awk 'BEGIN {FS="_"}{print $2}')
 fname2=$(echo $file | awk 'BEGIN {FS="_"}{print $3}' | sed 's/.pheno//g')
 
# scount=$(awk '{print $3}' pheno/$file | grep -v NA | wc -l | awk '{print $1}')
 scount=$(LC_ALL=C join <(awk '{print $1"_"$2,$3}' pheno/$file  | LC_ALL=C sort -k1b,1 ) <( awk '{print $1"_"$2}' pheno/p2_"$fname"_eur_"$fname2"_agesex.cov | LC_ALL=C sort -k1b,1) | awk '{print $2}' | grep -v NA | wc -l | awk '{print $1}')
 
 echo $fname $fname2 $scount 
 head -n1 pheno/$file

nievergeltlab

ls ../bgn | sed 's/.bed//g' | sed 's/.bim//g' | sed 's/.fam//g' | sort -u | awk '{print "../bgn/"$1".bed","../bgn/"$1".bim","../bgn/"$1".fam"}' > mergelist.txt

../plink --merge-list mergelist.txt --maf 0.00000001 --geno 0.02 --make-bed --out pgbd_nomono_goodgeno

nievergeltlab

nAA=76
nAa=1501
naa=7823
n=nAA+nAa+naa
pA=(2*nAA+nAa)/(2*n)
pa=1-pA
obs= c(nAA,nAa,naa)
exp=c(n*pA^2,2*n*pA*pa,n*pa^2)
val=sum(((obs-exp)^2)/exp)
1-pchisq(val,df=1)

nievergeltlab

LC_ALL=C sort -k1b,1 file1.txt > file1.txt.sorted
LC_ALL=C sort -k1b,1 file2.txt > file2.txt.sorted
LC_ALL=C join file1.txt.sorted file2.txt.sorted > joined.txt

nievergeltlab

expand.grid(study= c("pgc","hou"), imputation = c("unimputed","imputed"), subjs = c("pgbd","pgbd_va"), pheno = c("enter_M_v2","nonresponder","eventnonresponder"))