<jupyter_start><jupyter_text>1- Merging genotyping datasets (*Immuno* and *NeuroX*)<jupyter_code>from scipy.stats import norm
from matplotlib_venn import venn3

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import subprocess
import math
import sys
import os<jupyter_output><empty_output><jupyter_text>Introduction During this step of the preprocessing phase, we will merge *Immuno* and *NeuroX* genotyping data in a single dataset.Genetic variant positions in *Immuno* datasets were mapped on hg18 genome assembly, while *NeuroX* SNPs were mapped on hg19 assembly. This could cause issues and inconsistencies when merging together the two datasets. To solve this issue we have to lift *Immuno* variant coordinates from hg18 to the corresponding hg19 genomic positions. After having lifted *Immuno* SNP positions from hg18 to hg19 coordinates, we will merge the two datasets. The merged datasets will be used during the following steps.  Exploring *Immuno* and *NeuroX* Each dataset has 3 files:- **BIM**: contains SNPs, their genomic coordinates (chromosome and position) and phase (heterozygote or homozygous)    - **FAM**: contains subjects information (e.g. ID, gender, etc.)    - **BED**: binary file storing a matrix linking subjects and SNPs<jupyter_code>immuno_path = "../../data/genotyping/IMMUNO/"
neurox_path = "../../data/genotyping/NEUROX/"<jupyter_output><empty_output><jupyter_text>We begin by exploring raw genotyping data stored in *Immuno* and *NeuroX* BIM files.Let's start with *Immuno*.<jupyter_code>immuno_bim_fn = os.path.join(immuno_path, "IMMUNO.bim")
neurox_bim_fn = os.path.join(neurox_path, "NEUROX.bim")

immuno_bim = pd.read_csv(immuno_bim_fn, sep="\t", header=None)
neurox_bim = pd.read_csv(neurox_bim_fn, sep="\t", header=None)
immuno_bim.head(n=10)<jupyter_output><empty_output><jupyter_text>Now, let's look at *NeuroX*.<jupyter_code>neurox_bim.head(n=10)<jupyter_output><empty_output><jupyter_text>Lifting *Immuno* SNP positions to hg19 genome assembly As mentioned above, *Immuno* variant positions are mapped on hg18, while those of *NeuroX* are mapped on hg19 genome assembly. To solve this problem we lift *Immuno* SNP positions to the corresponding hg19 genomic locations. To lift SNP positions we use the web-based tool liftOver at the UCSC web-server (https://genome.ucsc.edu/cgi-bin/hgLiftOver).Since liftOver accepts as input UCSC BED files and not BIMs, we first need to convert *Immuno* BIM in a corresponding UCSC BED file.<jupyter_code># minimal UCSC BED file requires
#
# CHROM    START    STOP    FEATURE_ID 

tmp_chr = immuno_bim.iloc[:,0]  # chromosomes
tmp_pos = immuno_bim.iloc[:,3]  # SNP positions
tmp_name = immuno_bim.iloc[:,1] # SNP IDs

# append 'chr' in front of chromosome numbers (required by liftOver)
tmp_chr_str = [''.join(["chr", str(c)]) for c in tmp_chr.values]

# replace chr23, chr24, chr25 and chr26 with chrX, chrY, chrX and chrMT
tmp_chr_str2 = [c.replace('23', 'X') for c in tmp_chr_str]
tmp_chr_str3 = [c.replace('24', 'Y') for c in tmp_chr_str2]
tmp_chr_str4 = [c.replace('25', 'X') for c in tmp_chr_str3]
tmp_chr_str5 = [c.replace('26', 'MT') for c in tmp_chr_str4]

# write the resulting UCSC BED file
bed = pd.concat(
    [pd.DataFrame(tmp_chr_str5), tmp_pos, tmp_pos+1, tmp_name],
    axis=1
)
print(bed.head(n=5))
print(bed.tail(n=5))
bed.to_csv(
    os.path.join(immuno_path, "IMMUNO_tolift.bed"), sep="\t",
    index=False, header=False
)<jupyter_output>0        3        3             1
0  chr1   898835   898836  imm_1_898835
1  chr1  1108138  1108139  vh_1_1108138
2  chr1  1110294  1110295  vh_1_1110294
3  chr1  1125105  1125106     rs9729550
4  chr1  1130298  1130299     rs1815606
            0          3          3           1
196519   chrX  154889906  154889907   rs3093525
196520   chrX  154889917  154889918   rs3093526
196521   chrX  154892630  154892631   rs2742301
196522   chrX  154907376  154907377   rs2981835
196523  chrMT       3721       3722  MitoA3721G<jupyter_text>Now we can map *Immuno* positions to hg19 genome assembly using liftOver (https://genome.ucsc.edu/cgi-bin/hgLiftOver)LiftOver mapped succesfully 196,520 of IMMUNO SNPs of the original 196,524 genetic variants.The four unmapped SNPs were:- imm_3_50875337- rs1574660- rs1131012- MitoA3721GBefore proceeding, we rename the downloaded UCSC BED file with the lifted SNPs as ```IMMUNO_lifted.bed```.<jupyter_code>bed_lifted = pd.read_csv(os.path.join(immuno_path, "IMMUNO_lifted.bed"), sep="\t", header=None)

# we keep SNPs whose conversion to hg19 positions was succesful
good_snps = bed_lifted.iloc[:,3]
good_snps.to_csv(
    os.path.join(immuno_path, "good_snps.txt"),
    index=False, header=False
)<jupyter_output><empty_output><jupyter_text>Now, we keep the mapped SNPs and create the corresponding BIM, FAM and BED files (named ```IMMUNO_hg19```).<jupyter_code>!plink --bfile {os.path.join(immuno_path, "IMMUNO")} --update-map {os.path.join(immuno_path, "IMMUNO_lifted.bed")} 2 4 --make-bed --extract {os.path.join(immuno_path, "good_snps.txt")} --out {os.path.join(immuno_path, "IMMUNO_hg19")}<jupyter_output>PLINK v1.90b6.12 64-bit (28 Oct 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 <NAME>, <NAME>   GNU General Public License v3
Logging to ../../data/genotyping/IMMUNO/IMMUNO_hg19.log.
Options in effect:
  --bfile ../../data/genotyping/IMMUNO/IMMUNO
  --extract ../../data/genotyping/IMMUNO/good_snps.txt
  --make-bed
  --out ../../data/genotyping/IMMUNO/IMMUNO_hg19
  --update-map ../../data/genotyping/IMMUNO/IMMUNO_lifted.bed 2 4

16384 MB RAM detected; reserving 8192 MB for main workspace.
196524 variants loaded from .bim file.
523 people (343 males, 180 females) loaded from .fam.
--update-map: 196520 values updated.
Warning: Base-pair positions are now unsorted!
--extract: 196520 variants remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 523 founders and 0 nonfounders present.
Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596[...]<jupyter_text>Merging *Immuno* and *NeuroX* Since both *Immuno* and *NeuroX* SNP positions are now mapped on the same genome assembly (hg19), we can merge the two datasets.<jupyter_code>immuno_bim = pd.read_csv(os.path.join(immuno_path, "IMMUNO_hg19.bim"), sep="\t", header=None)  # BIM with SNPs mapped on hg19
immuno_bim.head(n=10)
neurox_bim.head(n=10)  # look again at NeuroX<jupyter_output><empty_output><jupyter_text>Before merging *Immuno* and *NeuroX*, we keep only those subjects with genotyping data in both the datasets.<jupyter_code>immuno_fam_fn = os.path.join(immuno_path, "IMMUNO_hg19.fam")
neurox_fam_fn = os.path.join(neurox_path, "NEUROX.fam")

immuno_fam = pd.read_csv(immuno_fam_fn, sep=" ", header=None)
neurox_fam = pd.read_csv(neurox_fam_fn, sep=" ", header=None)<jupyter_output><empty_output><jupyter_text>Let's explore the content of the two FAM files.Let's begin with *Immuno*.<jupyter_code># structure 
# FID IID FATHERID MOTHERID GENDER PHENOTYPE
immuno_fam.head(n=5)<jupyter_output><empty_output><jupyter_text>And *NeuroX*.<jupyter_code># structure 
# FID IID FATHERID MOTHERID GENDER PHENOTYPE
neurox_fam.head(n=5)<jupyter_output><empty_output><jupyter_text>We can now retrieve FAMIDs and IIDs of those subjects in both *Immuno* and *NeuroX*.<jupyter_code>immuno_subj = immuno_fam.iloc[:,1].tolist()
neurox_subj = neurox_fam.iloc[:,1].tolist()

# get the subjects appearing in both IMMUNO and NEUROX data
common_subj = set(immuno_subj).intersection(set(neurox_subj))

# write the common subjects to a file
common_subj_fn = "../../data/genotyping/common_subj.txt"
pd.DataFrame(list(zip(common_subj, common_subj)), columns=["FID", "IID"]).to_csv(
    common_subj_fn,
    header=False,
    index=False,
    sep=" "
)<jupyter_output><empty_output><jupyter_text>In *Immuno* there are 523 subjects, while in *NeuroX* there are 619 subjects. The subjects in both the datasets are 520. Before merging *Immuno* and *NeuroX*, we keep data belonging only to the subjects shared by both the datsets.<jupyter_code>!plink --bfile {os.path.join(immuno_path, "IMMUNO_hg19")} --keep {common_subj_fn} --make-bed --out {os.path.join(immuno_path, "IMMUNO_common")}
!plink --bfile {os.path.join(neurox_path, "NEUROX")} --keep {common_subj_fn} --make-bed --out {os.path.join(neurox_path, "NEUROX_common")}<jupyter_output>PLINK v1.90b6.12 64-bit (28 Oct 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 <NAME>, <NAME>   GNU General Public License v3
Logging to ../../data/genotyping/NEUROX/NEUROX_common.log.
Options in effect:
  --bfile ../../data/genotyping/NEUROX/NEUROX
  --keep ../../data/genotyping/common_subj.txt
  --make-bed
  --out ../../data/genotyping/NEUROX/NEUROX_common

16384 MB RAM detected; reserving 8192 MB for main workspace.
267607 variants loaded from .bim file.
619 people (409 males, 210 females) loaded from .fam.
--keep: 520 people remaining.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 520 founders and 0 nonfounders present.
Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989 done.
Warning: 4672 het. haploid genotypes present (see
../../data/genotyping/NEURO[...]<jupyter_text>PPMI documentation tells us that the genotyping was made using custom arrays. Therefore, SNP identifiers often are not named with standard rsIDs. Moreover, we want to keep as much SNPs as possible from the two datasets for upcoming analyses.Thus, to get those SNPs uniquely present in *Immuno*, we need to work on their genomic coordinates.<jupyter_code>#function to rename snps as chrom:pos_WT_AT
def rename_SNP(snpname, chrom, pos, a1, a2):
    if chrom == 0 or pos == 0:
        return snpname  # manage unmapped SNPs on neurox
    return ":".join([str(chrom), "_".join([str(pos), str(a1), str(a2)])])

# the BIM are those containing the common subjects data
immuno_bim = pd.read_csv(
    os.path.join(immuno_path, "IMMUNO_common.bim"), sep="\t",
    header=None
)
neurox_bim = pd.read_csv(
    os.path.join(neurox_path, "NEUROX_common.bim"), sep="\t",
    header=None
)

renamedsnps_immuno = immuno_bim.apply(lambda row : rename_SNP(row[1], row[0], row[3], row[4], row[5]), axis=1)
immuno_coords = immuno_bim
immuno_coords[6] = renamedsnps_immuno

renamedsnps_neurox = neurox_bim.apply(lambda row : rename_SNP(row[1], row[0], row[3], row[4], row[5]), axis=1)
neurox_coords = neurox_bim
neurox_coords[6] = renamedsnps_neurox
immuno_coords.head(n=10)
neurox_coords.head(n=10)<jupyter_output><empty_output><jupyter_text>Let's now rename the SNP IDs as chr:position_WT_AT.<jupyter_code># write name map file for IMMUNO dataset
pd.DataFrame(
    zip(immuno_coords.iloc[:,1].tolist(), immuno_coords.iloc[:,6].tolist())
).to_csv(
    os.path.join(immuno_path, "namemap.txt"),
    sep=" ",
    index=False
)

# write name map file for NEUROX dataset
pd.DataFrame(
    zip(neurox_coords.iloc[:,1].tolist(), neurox_coords.iloc[:,6].tolist())
).to_csv(
    os.path.join(neurox_path, "namemap.txt"),
    sep=" ",
    index=False
)
!plink --bfile {os.path.join(immuno_path, "IMMUNO_common")} --update-name {os.path.join(immuno_path, "namemap.txt")} --make-bed --out {os.path.join(immuno_path, "IMMUNO_renamed")}
!plink --bfile {os.path.join(neurox_path, "NEUROX_common")} --update-name {os.path.join(neurox_path, "namemap.txt")} --make-bed --out {os.path.join(neurox_path, "NEUROX_renamed")}<jupyter_output>PLINK v1.90b6.12 64-bit (28 Oct 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 <NAME>, <NAME>   GNU General Public License v3
Logging to ../../data/genotyping/IMMUNO/IMMUNO_renamed.log.
Options in effect:
  --bfile ../../data/genotyping/IMMUNO/IMMUNO_common
  --make-bed
  --out ../../data/genotyping/IMMUNO/IMMUNO_renamed
  --update-name ../../data/genotyping/IMMUNO/namemap.txt

16384 MB RAM detected; reserving 8192 MB for main workspace.
196520 variants loaded from .bim file.
520 people (341 males, 179 females) loaded from .fam.
--update-name: 196520 values updated, 1 variant ID not present.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 520 founders and 0 nonfounders present.
Calculating allele frequencies... 10111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989 done.
Warning: 66842 het. hapl[...]<jupyter_text>Then, let's continue by taking the SNPs that are only in *IMMUNO* dataset.<jupyter_code># get the names of SNPs only available in IMMUNO and write 
# the SNPs appearing uniquely in IMMUNO to a TXT file
immuno_bim = pd.read_csv(
    os.path.join(immuno_path, "IMMUNO_renamed.bim"),
    sep="\t",
    header=None
)
neurox_bim = pd.read_csv(
    os.path.join(neurox_path, "NEUROX_renamed.bim"),
    sep="\t",
    header=None
)
only_immuno_snps_fn = "../../data/genotyping/onlyIMMUNOsnps.txt"
neurox_snpids_set = set(neurox_bim.iloc[:,1].to_list())
immuno_bim[~immuno_bim[1].isin(neurox_snpids_set)].iloc[:,1].to_csv(
    only_immuno_snps_fn, header=False, index=False
)<jupyter_output><empty_output><jupyter_text>Then, in *Immuno* we keep only those SNPs which are not available in *NeuroX* dataset.<jupyter_code>!plink --bfile {os.path.join(immuno_path, "IMMUNO_renamed")} --extract {only_immuno_snps_fn} --make-bed --out {os.path.join(immuno_path, "IMMUNO_uniquesnps")}<jupyter_output>PLINK v1.90b6.12 64-bit (28 Oct 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 <NAME>, <NAME>   GNU General Public License v3
Logging to ../../data/genotyping/IMMUNO/IMMUNO_uniquesnps.log.
Options in effect:
  --bfile ../../data/genotyping/IMMUNO/IMMUNO_renamed
  --extract ../../data/genotyping/onlyIMMUNOsnps.txt
  --make-bed
  --out ../../data/genotyping/IMMUNO/IMMUNO_uniquesnps

16384 MB RAM detected; reserving 8192 MB for main workspace.
196520 variants loaded from .bim file.
520 people (341 males, 179 females) loaded from .fam.
--extract: 189564 variants remaining.
Warning: At least 402 duplicate IDs in --extract file.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 520 founders and 0 nonfounders present.
Calculating allele frequencies... 1011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798[...]<jupyter_text>We can now merge *Immuno* and *NeuroX* in a single datset.Note that are created also the corresponding BIM, FAM and BED, representing the merged data.<jupyter_code>merged_ppmi_fn = "../../data/genotyping/PPMI_merge"

!plink --bfile {os.path.join(neurox_path, "NEUROX_renamed")} --bmerge {os.path.join(immuno_path, "IMMUNO_uniquesnps")} --make-bed --out {merged_ppmi_fn}
!cat ../../data/genotyping/PPMI_merge.log | grep 'Warning' | grep 'Variants' | grep -o "'.*'" | cut -d' ' -f1,3 | awk '{print $1"\n"$2;}' | tr -d "'" | sort -u > ../../data/genotyping/merge_rm_snps.txt
!plink --bfile {os.path.join(neurox_path, "NEUROX_renamed")} --bmerge {os.path.join(immuno_path, "IMMUNO_uniquesnps")} --exclude ../../data/genotyping/merge_rm_snps.txt --make-bed --out {merged_ppmi_fn}<jupyter_output>PLINK v1.90b6.12 64-bit (28 Oct 2019)          www.cog-genomics.org/plink/1.9/
(C) 2005-2019 <NAME>, <NAME>   GNU General Public License v3
Logging to ../../data/genotyping/PPMI_merge.log.
Options in effect:
  --bfile ../../data/genotyping/NEUROX/NEUROX_renamed
  --bmerge ../../data/genotyping/IMMUNO/IMMUNO_uniquesnps
  --exclude ../../data/genotyping/merge_rm_snps.txt
  --make-bed
  --out ../../data/genotyping/PPMI_merge

16384 MB RAM detected; reserving 8192 MB for main workspace.
520 people loaded from ../../data/genotyping/NEUROX/NEUROX_renamed.fam.
520 people to be merged from
../../data/genotyping/IMMUNO/IMMUNO_uniquesnps.fam.
Of these, 0 are new, while 520 are present in the base dataset.
266454 markers loaded from ../../data/genotyping/NEUROX/NEUROX_renamed.bim.
189564 markers to be merged from
../../data/genotyping/IMMUNO/IMMUNO_uniquesnps.bim.
Of these, 189162 are new, while 402 are present in the base dataset.
Warning: Variants 'NeuroX_PINK1_23bp_del_ex7' a[...]
{'fraction_non_alphanumeric': 0.12101055714895954, 'fraction_numerical': 0.08561664734240557, 'mean_word_length': 4.748158540862855, 'pattern_counts': {'":': 1, '<': 77, '<?xml version=': 0, '>': 78, 'https://': 2, 'lorem ipsum': 0, 'www.': 5, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '296074', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7194, 'n_tokens_neox': 6045, 'n_words': 1527, 'path': 'src/1-Preprocessing/1-Merging_and_Lifting.ipynb', 'repo_name': 'InfOmics/MUVIG', 'size': '16477'}