Skip to content

Instantly share code, notes, and snippets.

View nievergeltlab's full-sized avatar

Adam X. Maihofer nievergeltlab

16 followers · 0 following
View GitHub Profile
@nievergeltlab
nievergeltlab / ldsc_results_extractor.sh
Created January 5, 2023 19:40
Extract LDSC results from a list of .log files. Makes it easy to put many results into a table.
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
nievergeltlab / replicate_region.r
Created August 17, 2022 20:37
Determine if a SNP is within a locus by applying a TRUE/FALSE criteria to each required bound. If all bounds are true, the SNP is within a given locus..
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
nievergeltlab / nb_logodds.r
Created June 14, 2022 16:35
Following Sroka 2018, likelihood calculation for a nb glm model with a log odds link This provides odds ratios for count data!
#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
nievergeltlab / get_ns.r
Created June 30, 2021 16:46
Intersect phenotype and covariate files to get total N from each phenotype file.
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
nievergeltlab / 0_forest_plot.r
Created June 29, 2021 19:09
Extract summary data from GWAS to produce a forest plot. Contains code for pivoting alleles to the same strand across datasets.
args <- commandArgs(trailingOnly = TRUE)
snpresults <- args[1]
header_classes <- args[2]
descriptor <-args[3]
outfile <- args[4]
# snpresults="results_filtered/forest_plots/rs10266297"
# header_classes="results_filtered/forest_plots/header_classes.csv"
# descriptor="results_filtered/forest_plots/forestplot_descriptor.csv"
@nievergeltlab
nievergeltlab / optimal_p.r
Created June 7, 2021 16:49
library(MASS)
library(data.table)
require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
library(lmtest)
#if it complains about this package:
#install.packages('zoo',repos="https://cran.r-project.org")
#install.packages('lmtest',repos="https://cran.r-project.org")
@nievergeltlab
nievergeltlab / job_example.sh
Last active May 25, 2021 20:00
IF you want to run batches of jobs on a single node, but need multiple runs (e.g. 6 core processor, but 20 processes to run), this is a mostly efficient way to achieve this. Not as good as a dedicated system, but close as long as jobs are equal in ru
inc1=0
inc2=0
while [ $inc1 -le 100 ]
do
inc1=$(($inc2 + 1))
inc2=$((inc1 + 5))
echo $inc1 $inc2
for i in $(seq $inc1 1 $inc2 )
do
@nievergeltlab
nievergeltlab / 00_launch_plink.sh
Last active May 25, 2021 20:03
@nievergeltlab
nievergeltlab / prs_nl.r
Created May 24, 2021 21:11
library(MASS)
library(data.table)
require(foreign)
require(nnet)
require(ggplot2)
require(reshape2)
library(lmtest)
#if it complains about this package:
#install.packages('zoo',repos="https://cran.r-project.org")
#install.packages('lmtest',repos="https://cran.r-project.org")
@nievergeltlab
nievergeltlab / ukbb_gwas
Created May 24, 2021 21:07
#PLINK2 with dosages
for chr in 2 # {1..22}
do
./plink2 --bgen /mnt/ukbb/adam/bgen/ukb_imp_chr"$chr"_v3.bgen ref-first\
--maf 0.01 \
--geno 0.05 \
--keep UKB_ptsd_eur_unrelated_m0_may13_2021.pheno \
--sample /mnt/ukbb/adam/ptsd/preimputation_genotypes/ukb41209_imp_chr1_v3_s487320.sample \
--pheno UKB_ptsd_eur_unrelated_m0_may13_2021.pheno --pheno-name PCL_sum_imp \