<jupyter_start><jupyter_text>Estimating fixation index in trans-eQTLs- toc: true - badges: true- comments: true  BackgroundTrans-eQTLs can be confounded by population substructure. To get an idea, I am looking at the fixation index of the trans-eQTLs. I calculated $F_{ST}$ from the 1000 Genome Project using 668 EUR (including CEU, FIN, GBR, IBS, TSI) subjects and 806 AFR (including ESN, GWD, LWK, MSL, YRI) subjects. I used the mapping between GTEx variant IDs and dbGap 151 hg38 mappings for rsids. I used the [hg38 mapped version of 1KG](ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/GRCh38_positions/).  Path to all data files<jupyter_code>#collapse-hide
teqtl_fst_file = '/scratch/sbanerj/trans-eqtl/fst_from_1kg/trans_eqtls/trans_eqtls_EUR_AFR_fst.txt'
gwas_fst_file = '/scratch/sbanerj/trans-eqtl/fst_from_1kg/gwas/gwas_SNPs_EUR_AFR_fst.txt'
random_fst_file = '/scratch/sbanerj/trans-eqtl/fst_from_1kg/trans_eqtls/random_snps_EUR_AFR_fst.txt'
chr1_fst_file = '/scratch/sbanerj/trans-eqtl/fst_from_1kg/fst/EUR-AFR_chr1.weir.fst'
teqtl_results = '/usr/users/sbanerj/trans-eQTL/gtex_v8_tejaas_trans_eQTLs_20200309/protein_coding_lncRNA_fixed_gamma_alt-haa-spl-pan-wb_combined_cut5e-8'
chr1_rr_file = '/scratch/sbanerj/trans-eqtl/dev-pipeline/gtex_v8_202003/as/tejaas/raw_std/permnull_sb0.1_knn30/chr1/rr.txt'
tissue_file = 'gtex/tissue_table.txt'
json_file = 'gtex/gtex_v8_metadata.json'
nsample_file = "gtex/tissue_nsamples.txt"
popsample_file = 'gtex/tissue_nsamples_with_pop.txt'<jupyter_output><empty_output><jupyter_text>Python libraries<jupyter_code>#collapse-hide
import os
import numpy as np
import json
import matplotlib.pyplot as plt
from utils import mpl_stylesheet
mpl_stylesheet.banskt_presentation(fontfamily = 'latex-clearsans', fontsize = 18, colors = 'banskt', dpi = 72)
from utils import utils<jupyter_output><empty_output><jupyter_text>Read metadata for our trans-eQTLs<jupyter_code>#collapse-hide
tshorts, tfulls, tstrings = utils.read_tissues(tissue_file)
with open(json_file) as instream:
    gtex_meta = json.load(instream)
tissue_colors = dict()
tissue_names = dict()
tissue_nsamples = dict()

for tshort, tfull, tstring in zip(tshorts, tfulls, tstrings):
    if tshort in tshorts:
        tissue_names[tshort] = tstring
        tissue_colors[tshort] = "#" + gtex_meta[tfull]["colorHex"]
        
tissue_nsamples = dict()       
tissue_noneur_samples = dict()
tissue_afr_samples = dict()
tissue_frac_noneur_to_eur = dict()
with open(popsample_file, 'r') as instream:
    header = instream.readline().split()
    totidx = header.index('TOT')
    euridx = header.index('EUR')
    afridx = header.index('AFR')
    for line in instream:
        linesplit = line.strip().split()
        tshort = linesplit[0].strip()
        ntot = int(linesplit[totidx])
        neur = int(linesplit[euridx])
        tissue_nsamples[tshort] = ntot
        tissue_noneur_samples[tshort] = ntot - neur
        tissue_afr_samples[tshort] = int(linesplit[afridx])
        tissue_frac_noneur_to_eur[tshort] = (ntot - neur) / neur
        
brain_tissues = ['bam', 'ban', 'bca', 'bceh', 'bce', 'bco', 'bfr', 'bhi', 'bhy', 'bnu', 'bpu', 'bsp', 'bsu']
altsb_tissues = ['haa', 'pan', 'spl', 'wb']<jupyter_output><empty_output><jupyter_text>Define some functions to read results<jupyter_code>#collapse-hide
def get_fstmap(fst_file, headerline = 0):
    fstmap = dict()
    for c in range(1,23):
        fstmap[f'chr{c}'] = {}

    with open(fst_file, 'r') as instream:
        for line in instream:
            linesplit = line.strip().split()
            chrm = int(linesplit[0])
            bppos = int(linesplit[1])
            fst = float(linesplit[2])
            if np.isnan(fst): fst = 0
            fstmap[f'chr{chrm}'][bppos] = fst
    return fstmap

def read_transeqtls(filename):
    idcol = 0
    chrcol = 1
    bpcol = 2
    pvcol = 7
    rsidlist = list()
    chrmlist = list()
    bplist   = list()
    pvallist = list()
    with open(filename, 'r') as instream:
        next(instream)
        for line in instream:
            linesplit = line.strip().split()
            rsid = linesplit[idcol]
            chrm = int(linesplit[chrcol])
            bppos = int(linesplit[bpcol])
            pval = float(linesplit[pvcol])
            rsidlist.append(rsid)
            chrmlist.append(chrm)
            bplist.append(bppos)
            pvallist.append(pval)
    return rsidlist, chrmlist, pvallist, bplist<jupyter_output><empty_output><jupyter_text>Read the results for trans-eQTLs<jupyter_code>#collapse-hide
teqtl_fstmap = get_fstmap(teqtl_fst_file)
random_fstmap = get_fstmap(random_fst_file)
gwas_fstmap = get_fstmap(gwas_fst_file)
teqtl_fst = [x for c in range(1,23) for k, x in teqtl_fstmap[f'chr{c}'].items()]
random_fst = [x for c in range(1,23) for k, x in random_fstmap[f'chr{c}'].items()]
gwas_fst = [x for c in range(1,23) for k, x in gwas_fstmap[f'chr{c}'].items()]

mean_fst_per_tissue = dict()
tissue_nteqtls_with_fst = dict()
for tshort in tshorts:
    _fst = list()
    resfile = os.path.join(teqtl_results, tshort, 'trans_eqtls_ldpruned.txt')
    gwvarids, gwchrms, gwpvals, gwbppos = read_transeqtls(resfile)
    for i, varid in enumerate(gwvarids):
        chrm = gwchrms[i]
        bppos = gwbppos[i]
        if bppos in teqtl_fstmap[f'chr{chrm}']:
            _fst.append(teqtl_fstmap[f'chr{chrm}'][bppos])
    if _fst:
        tissue_nteqtls_with_fst[tshort] = len(_fst)
        mean_fst_per_tissue[tshort] = np.mean(np.array(_fst))<jupyter_output><empty_output><jupyter_text>Read the results for p-values<jupyter_code>#collapse-hide
chr1_fstmap = get_fstmap(chr1_fst_file)

pvalue_lowfst = list()
pvalue_highfst = list()

gwvarids, gwchrms, gwpvals, gwbppos = read_transeqtls(chr1_rr_file)
for i, varid in enumerate(gwvarids):
    chrm = gwchrms[i]
    bppos = gwbppos[i]
    if bppos in chr1_fstmap[f'chr{chrm}']:
        this_fst = chr1_fstmap[f'chr{chrm}'][bppos]
        if this_fst <= 0.5:
            pvalue_lowfst.append(gwpvals[i])
        else:
            pvalue_highfst.append(gwpvals[i])<jupyter_output><empty_output><jupyter_text>Plot * Compare $F_{ST}$ distributions for all our trans-eQTLs with random SNPs. * Check if the mean of $F_{ST}$ depends on the number of trans-eQTLs discovered. Right panel: color according to GTEx tissues, size according to number of samples<jupyter_code>#collapse-hide
fig = plt.figure(figsize = (12, 6))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)

ax1.hist(gwas_fst, alpha = 0.7, bins = 50, lw = 0, density = True, label = 'GWAS SNPs')
ax1.hist(random_fst, alpha = 0.5, bins = 50, lw = 0, density = True, label = 'Random SNPs')
ax1.hist(teqtl_fst, alpha = 0.7, bins = 50, lw = 0, density = True, label = 'Trans-eQTLs')

ax1.legend()
ax1.set_xlabel('$F_{ST}$')
ax1.set_ylabel('Density')

whichtissues = [x for x in tshorts if x in mean_fst_per_tissue]
yvals = [mean_fst_per_tissue[x] for x in whichtissues]
xvals = [tissue_frac_noneur_to_eur[x] for x in whichtissues]
# xvals = [np.log2(tissue_nteqtls_with_fst[x]) for x in whichtissues]
sizes = [tissue_nsamples[x] / 2 for x in whichtissues]
# sizes = [tissue_sample_frac_noneur_to_eur[x] * 1000 for x in whichtissues]
mcolors = [tissue_colors[x] for x in whichtissues]
ax2.scatter(xvals, yvals, color = mcolors, s = sizes, edgecolors = 'black', zorder = 10)
ax2.set_xlabel('non-European / European')
ax2.set_ylabel('Mean $F_{ST}$')
# xticks = [2, 4, 6, 8, 10, 12]
# ax2.set_xticks(xticks)
# ax2.set_xticklabels([int(2**x) for x in xticks])

plt.tight_layout()
plt.show()<jupyter_output><empty_output><jupyter_text>* Compare Tejaas P-values of SNPs with high $F_{ST}$ and low $F_{ST}$<jupyter_code>#collapse-hide
fig = plt.figure(figsize = (12, 6))
ax1 = fig.add_subplot(111)

ax1.hist(pvalue_highfst, bins = 100, lw = 0, density = True, alpha = 0.7, label = "High $F_{ST}$")
ax1.hist(pvalue_lowfst, bins = 100, lw = 0, density = True, alpha = 0.7, label = "Low $F_{ST}$")
ax1.set_xlabel("P-value")
ax1.set_ylabel("Density")
ax1.legend()

plt.tight_layout()
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1278204810314902, 'fraction_numerical': 0.021696007934540046, 'mean_word_length': 4.296782665791201, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 33, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '92259', 'license': 'Apache-2.0', 'n_tokens_mistral': 3403, 'n_tokens_neox': 3269, 'n_words': 647, 'path': '_notebooks/2020-06-08-estimating-fixation-index-in-trans-eqtls.ipynb', 'repo_name': 'banskt/trans-eqtl-ideas', 'size': '8066'}