gibran hemani explodecomputer

RSID lookup with tabix

Strategy:

Extract the chr/pos/rsid from dbsnp into a .bed file. Here it makes a separate .bed file per chromosome for parallelisation of lookups
Create a tabix index for the .bed files
Extract the chr/pos from the target VCF file
Use tabix to query the target chr/pos against the tabix indexed .bed files. Parellelised across chromosomes using GNU parallel.
Update the VCF file with the extracted RSIDs

	library(dplyr)
	library(data.table)
	library(seqminer)
	library(stringr)

	dat <- fread("/mnt/project/analysis/bmi-alcohol-ancestry/phen_extract.csv")
	str(dat)

	anc <- fread("/mnt/project/data/ancestry/ukb_InferredAncestry.txt")
	str(anc)

	#!/bin/bash

	set -e

	MY_GITHUB_URL=$1

	if [ -z "${MY_GITHUB_URL}" ]; then
	echo "Please MY_GITHUB_URL as the argument to this script"
	exit 1;
	fi

	#!/bin/bash

	#!/bin/bash

	for chr in {1..22}
	do
	# Get the list of variants that pass the threshold
	bcftools filter -e 'INFO/R2<0.8' /mnt/project/Bulk/Imputation/Imputation\ from\ genotype\ \(TOPmed\)/helper_files/ukb21007_c${chr}_b0_v1.sites.vcf.gz \| \
	bcftools filter -e 'INFO/AF<0.01' \| \
	bcftools filter -e 'INFO/AF>0.99' \| \

	library(readxl)
	library(dplyr)

	download.file("https://www.biorxiv.org/content/biorxiv/early/2022/06/18/2022.06.17.496443/DC2/embed/media-2.xlsx?download=true", here("data", "media-2.xlsx"))
	a <- read_xlsx(here("data", "media-2.xlsx"), sheet="ST6", skip=2) %>%
	dplyr::select(
	exposure="Assay Target",
	id.exposure="Assay Target",
	SNP="rsID",
	beta.exposure="BETA (discovery, wrt. A1)",

	library(dplyr)

	standardise <- function(d, ea_col="ea", oa_col="oa", beta_col="beta", eaf_col="eaf", chr_col="chr", pos_col="pos", vid_col="vid")
	{
	toflip <- d[[ea_col]] > d[[oa_col]]
	d[[eaf_col]][toflip] <- 1 - d[[eaf_col]][toflip]
	d[[beta_col]][toflip] <- d[[beta_col]][toflip] * -1
	temp <- d[[oa_col]][toflip]
	d[[oa_col]][toflip] <- d[[ea_col]][toflip]
	d[[ea_col]][toflip] <- temp

	library(ieugwasr)
	library(dplyr)

	ukb <- tophits("ukb-d-I9_IHD")
	bbj <- tophits("bbj-a-159", pop="EAS", r2=0.0011)
	table(duplicated(c(ukb$rsid, bbj$rsid)))

	cl <- bind_rows(ukb, bbj) %>% select(rsid, p) %>% ld_clump()

	a <- associations(cl$rsid, c("bbj-a-159", "ukb-d-I9_IHD"))

	library(dplyr)

	#' Genarte genotypes for nuclear families
	#'
	#' @param af Vector of allele frequencies for SNPs (length = number of SNPs to simulate)
	#' @param nfam Number of families to generate
	#'
	#' @return List of genotype matrices for each family member type (mum, dad, sib1, sib2)
	make_families <- function(af, nfam)
	{

	library(lubridate)
	library(dplyr)
	library(knitr)

	#' Create agenda
	#'
	#' Alternating weeks of Talks and Discussions.
	#' Presenter for a talk leads discussion in the following week.
	#'
	#' @param start Start date in format "yyyy/mm/dd"

	# MR-Base (i.e. using the cloud)

	library(TwoSampleMR)
	metab_ids <- c("met-e-1", "met-e-2")

	# Extract data - from API
	exp_dat <- extract_instruments("alz_id")
	out_dat <- extract_outcome_data(exp_dat$SNP, metab_ids)

	# Harmonise and run