liangyy

# confidence interval of p-value
# ref: https://gist.github.com/hakyim/38431b74c6c0bf90c12f
## copied from ref:
##   the jth order statistic from a
##   uniform(0,1) sample
##   has a beta(j,n-j+1) distribution
##   (Casella & Berger, 2002,
##   2nd edition, pg 230, Duxbury)
##   this portion was posted by anonymous on
##   http://gettinggeneticsdone.blogspot.com/2009/11/qq-plots-of-p-values-in-r-using-ggplot2.html

liangyy

liangyy

qqplot_by_group <- function(pval, group, ...) {
  n <- length(pval)
  pexp <- rank(pval) / n
  df <- data.frame(p.val = pval, grp = group) %>% group_by(grp) %>% mutate(p.exp = rank(p.val) / (n() + 1)) %>% ungroup()
  p <- ggplot(df) + 
    geom_point(aes(x = -log10(p.exp), y = -log10(p.val), color = grp), ...) + 
    geom_hline(yintercept = -log10(0.05 / n)) + 
    geom_abline(slope = 1, intercept = 0, linetype = 2)
  return(p)
}

liangyy

# simple and fast linear regression solver
# y ~ x + covariate
# input:
# y: n x 1 (vector)
# x: n x k (matrix)
# covariate: n x m (optional; matrix)
# 
# output:
# coefficient of x[i] in y ~ x[i] + covariate 
# and corresponding p-value (from t value)

liangyy

load_mixqtl_data_bundle = function(data_dir = 'sample_data/', gene_id = 'ENSG00000198951.6', cis_window = 10000){
  if(!dir.exists(data_dir)) {
    message('data_dir: ', data_dir, ' does not exist! Exit!')
    quit()
  }
  haplotype1 = 'geuvadis_22_hap1.txt.gz'
  haplotype2 = 'geuvadis_22_hap2.txt.gz'
  total_count = 'expression_count.txt.gz'
  allele_specific_count_hap1 = 'geuvadis.asc.h1.tsv.gz'
  allele_specific_count_hap2 = 'geuvadis.asc.h2.tsv.gz'

liangyy

run_gwas = function(X, y) {
  a_y = sd(y)
  b_x = apply(X, 2, sd)
  y_tilde = y / a_y
  # message(1)
  x_tilde = sweep(X, 2, FUN = '/', b_x)
  # message(2)
  xtx = colSums(x_tilde ^ 2)
  xbar = colMeans(x_tilde)
  ybar = mean(y_tilde)

liangyy

# with clang binary downloaded from: http://r.research.att.com/libs/clang-4.0.0-darwin15.6-Release.tar.gz
env CC=/path/to/local/clang4/bin/clang pip install rpy2

liangyy

load_gtex_color_board = function() {
  gtex_color = read.csv('https://bitbucket.org/yanyul/rotation-at-imlab/raw/b49d35b44c0ab936496733268b8578db77eb7f73/data/gtex_tissue_colors.csv', stringsAsFactors = FALSE)
  gtex_color$tissue_site_detail_id[gtex_color$tissue_site_detail_id == 'Cells_Transformed_fibroblasts'] = 'Cells_Cultured_fibroblasts'
  color_board = paste0('#', gtex_color$tissue_color_hex)
  names(color_board) = gtex_color$tissue_site_detail_id
  return(color_board)
}

liangyy

wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/dbsnp_138.b37.vcf.gz.md5
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/dbsnp_138.b37.vcf.gz
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/Mills_and_1000G_gold_standard.indels.b37.vcf.gz.md5
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.indels.b37.vcf.gz
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.indels.b37.vcf.gz.md5
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.snps.high_confidence.b37.vcf.gz
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/1000G_phase1.snps.high_confidence.b37.vcf.gz.md5
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/hapmap_3.3.b37.vcf.gz
wget -c ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/b37/hapmap_3.3.b37.v

liangyy

bcftools query -f query -f chr%CHROM\\t%ID\\t%POS\\t%REF\\t%ALT\\t%QUAL[\\t%DS]\\n test3_vcf.vcf.gz | gzip > test3_to_predixcan_format.gz
bcftools query -l test3_vcf.vcf.gz | tr '_' '\t' > test3_to_predixcan_format.sample