Skip to content

Instantly share code, notes, and snippets.

View nievergeltlab's full-sized avatar

Adam X. Maihofer nievergeltlab

11 followers · 0 following
View GitHub Profile
@nievergeltlab
nievergeltlab / forest_plot
Last active February 14, 2017 20:15
Forest plot for meta-analysis in R
library(metafor)
y <- scan(what=numeric())
0.273936
0.278002
0.117789
0.170343
se <- scan(what=numeric())
0.0864617
@nievergeltlab
nievergeltlab / name_chromosome_chunk.sh
Created March 13, 2017 16:18
Get chromosome number and chunk number from file name. For ricopili.
chunknum=$(echo $chunk | sed 's/.*\(chr[0-9]*_[0-9]*_[0-9]*\).*/\1/')
chrnum=$(echo $chunk | sed 's/.*\(chr[0-9]*\).*/\1/')
@nievergeltlab
nievergeltlab / ibd_compare
Created March 22, 2017 22:08
2 file IBD, if strand mismatches and allele mismatches
bimfile=pts_meg2_mix_am-qc.bim
grep -P "A\tT" $bimfile > ambiguous_snps.txt
grep -P "T\tA" $bimfile >> ambiguous_snps.txt
grep -P "C\tG" $bimfile >> ambiguous_snps.txt
grep -P "G\tC" $bimfile >> ambiguous_snps.txt
#filer out ambiguous
$plink_location --bfile pts_meg2_mix_am-qc --exclude ambiguous_snps.txt --make-bed --out pts_meg2_mix_am-qc-ambig
@nievergeltlab
nievergeltlab / extract_snp_from_dosages
Created March 24, 2017 18:23
Algorithm: Find which file the SNP is in, then zgrep the SNP out. This script will output a .dosage file for each snp in the list. to combine them into a file, it's of course just a matter of using "cat"
#Requires a list of SNPs written into snplist.txt
dosage_directory=write_path_to_dosage_files_here
for snp in $(cat snplist.txt)
do
filename=$(grep -w -m1 -l $snp "$dosage_directory"/*.map | sed 's/.map/.gz/g' )
zgrep -w -m1 $snp $filename > "$snp".dosage
done
@nievergeltlab
nievergeltlab / analyze_dosages.pbs
Last active April 4, 2017 19:59
Analyze dosages data from Ricopili
#PBS -lnodes=1
#PBS -lwalltime=0:05:00
#!/bin/bash
while getopts a:b:o:p:d: option
do
case "${option}"
in
a) phenotype=${OPTARG};;
@nievergeltlab
nievergeltlab / make_dos.pbs
Last active April 4, 2017 20:10
This makes imputed probabilities from ricopili into a dosage file format Working with job system (TORQUE). Have a list of files. Want to run same operation on all of them. Can run multiple at the same time on a single node The files to be processed are listed in a single file. The number of processes to be run on a single node are noted. The job…
#PBS -lnodes=1
#PBS -lwalltime=0:05:00
#!/bin/bash
while getopts n:o:p:d:s: option
do
case "${option}"
in
o) outputfile=${OPTARG};;
@nievergeltlab
nievergeltlab / tar_wildcards.sh
Created April 8, 2017 23:17
Extract files from tar with wildcards
tar xvzf dnhs_qc_v2_mar7_2017.tgz --wildcards '*-qc.fam' --wildcards '*.header'
@nievergeltlab
nievergeltlab / convert_dosage_to_gmmat_format.sh
Created April 10, 2017 16:00
Logistic Mixed Models with GMMAT
#Script to convert dosage data to GMMAT format. Requires also a list of usuable SNPs as input, but this list does not need to be filtered.
#Create a valid SNPlist
info=0.9
maf=0.01
zcat /home/maihofer/vets/qc/imputation/distribution/vets_eur_analysis_run2/daner_vets_eur_analysis_run2.gz | awk -v info=$info -v maf=$maf '{if (($8 > info) && ($6 > maf) && ($7 > maf) && ($6 < 1-maf) && ($7 < 1-maf) && ($11 != "NA")) print $2}' > european_valid_snps.snplist
echo "SNP" > snp.txt
cat snp.txt european_valid_snps.snplist | LC_ALL=C sort -k 1b,1 > european_valid_snps.sorted.snplist
@nievergeltlab
nievergeltlab / merge_discordance.r
Created April 19, 2017 20:55
Descriptive statistics for genotype merge discordances
#Merge data in PLINK, using merge mode 6 or 7
./plink --bfile YEHUDA --bmerge YEHUDA.bed YEHUDA.bim YEHUDA.fam --merge-mode 7 --out yehude-merge
R
setwd('F:/rutgers_2')
dat <- read.table('yehude-merge.diff', header=T,nr=800000,stringsAsFactors=F)
library(plyr)
#Determine general amount of disagreement for each SNP. Ones that have especially high disagreement may be badly genotyped
quantile(table(dat$SNP))
#Plot it
@nievergeltlab
nievergeltlab / manhattan_plot.r
Last active May 1, 2017 18:51
Manhattan Plots V2 of the plotter will color SNPs within a set kb window of any hits
args <- commandArgs(trailingOnly = TRUE)
scriptloc <- args[1]
results <- args[2]
outfile <- args[3]
snpcol <- args[4]
chrcol <- args[5]
bpcol <- args[6]
pcol <- args[7]