<jupyter_start><jupyter_text>AnnoQ Gene IDs Comparison The following code was designed to be run using input from [AnnoQ](http://annoq.org/), which compiles results from [ANNOVAR](https://annovar.openbioinformatics.org/en/latest/), [SnpEff](http://pcingola.github.io/SnpEff/), and [VEP](https://uswest.ensembl.org/index.html). The goal of this project is to compare Gene IDs assigned by the different annotation tools to the same SNP in order to inform the end users of any potential missing information. **Disclosure**: The following code is broken into pieces. The code ran in the USC HPC was enclosed in a *main()* function in order to automatically run all the Chromosome SNPs  Libraries<jupyter_code>import pandas as pd
import numpy as np
import re<jupyter_output><empty_output><jupyter_text>Functions<jupyter_code># Eliminating repeats
def no_repeats(nparray):
    place = []
    for i in range(0,size):
        temp = []
        temp = np.unique(nparray[i])
        place.append(temp)
    return(place)
        
# Mining gene ID
def extract(lists):
    place =[]
    for i in range(0,size):
        temp = []
        temp = re.findall("ENSG...........", lists[i])
        place.append(temp)
    return(place)

# Looking for differences and writting down the GeneID
def difference(array1,array2):
    place =[]
    for i in range(0,size):
        temp = []
        temp = np.setdiff1d(array1[i], array2[i], assume_unique = True)
        place.append(temp)
    return(place)

# Counting
def counter(arrays):
    place = 0
    for i in range(0,size):
        if(arrays[i].size > 0):
            place = place + 1
    return place<jupyter_output><empty_output><jupyter_text>Data Processing  File Input & Data Mining<jupyter_code># Reading file 
filename = "Test_data.txt"
cdata = pd.read_csv(filename, delimiter="\t", dtype= object)

# Selecting data
AN_ID_1 = cdata["ANNOVAR_ensembl_Gene_ID"]  #Gene ID should be here
AN_ID_2 = cdata["ANNOVAR_ensembl_Closest_gene(intergenic_only)"]    #alternative place for Gene ID

SN_ID_1 = cdata["SnpEff_ensembl_Gene_ID"]   #Gene ID for SnpEff

VP_ID_1 = cdata["VEP_ensembl_Gene_ID"]      #Gene ID for VEP

chr = cdata["chr"] # chromosome number
pos = cdata["pos"] # SNP position
ref = cdata["ref"] # reference
alt = cdata["alt"] # mutation
rs = cdata["rs_dbSNP151"] # rs ID

# Finding length of dataset
size = len(AN_ID_1)<jupyter_output><empty_output><jupyter_text>Gene ID Extraction & Concatination<jupyter_code># Finding the Gene ID in ANNOVAR
AN_ID_tog = []
for i in range(0,size):
    if(AN_ID_1[i] == "."):
        AN_ID_tog.append(AN_ID_2[i])
    else:
        AN_ID_tog.append(AN_ID_1[i])
        
# Extracting Gene ID in ANNOVAR
AN_ID = []
AN_ID = extract(AN_ID_tog)

# Extracting Gene ID in SnpEff
SN_ID = []
SN_ID = extract(SN_ID_1)

# Extracting Gene ID in VEP
VP_ID = []
VP_ID = extract(VP_ID_1)

# Converting lists to arrays
AN_ID = np.array(AN_ID, dtype=object)
SN_ID = np.array(SN_ID, dtype=object)
VP_ID = np.array(VP_ID, dtype=object)

# Getting rid of repeated Gene IDs
AN_ID_clean = no_repeats(AN_ID)
SN_ID_clean = no_repeats(SN_ID)
VP_ID_clean = no_repeats(VP_ID)

# Merging all Gene IDs 
united = []
for i in range(0,size):
  temp = []
  temp = np.concatenate((AN_ID[i], SN_ID[i], VP_ID[i]), axis=None)
  united.append(temp)

# Getting rid of repeats for each SNP
unique_clean = no_repeats(united)<jupyter_output><empty_output><jupyter_text>Comparing unique Gene IDs<jupyter_code># Creating empty arrays
AN_ID_check = []
SN_ID_check = []
VP_ID_check = []

# Parsing through dataset
for i in range(0,size):
  zize = len(unique_clean[i])   # number of Gene IDs
  if zize == 1:                 # if only 1 Gene ID 
    if unique_clean[i] in AN_ID_clean[i]:   ## Check if present
      AN_ID_check.append([zize]) 
    else:
        AN_ID_check.append([])
    if unique_clean[i] in SN_ID_clean[i]:
      SN_ID_check.append([zize]) 
    else:
        SN_ID_check.append([])
    if unique_clean[i] in VP_ID_clean[i]:
      VP_ID_check.append([zize]) 
    else: 
      VP_ID_check.append([])
  elif zize > 1:                 # if multiple Gene IDs
    place_1 = []
    place_2 = []
    place_3 = []
    for x in range(0,zize):      ## Check which ones are present 
      if unique_clean[i][x] in AN_ID_clean[i]: place_1.append(x+1)
      if unique_clean[i][x] in SN_ID_clean[i]: place_2.append(x+1)
      if unique_clean[i][x] in VP_ID_clean[i]: place_3.append(x+1)
    AN_ID_check.append(place_1)
    SN_ID_check.append(place_2)
    VP_ID_check.append(place_3)<jupyter_output>/usr/local/lib/python3.7/dist-packages/ipykernel_launcher.py:29: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison<jupyter_text>Comparison Results Array Table<jupyter_code># Creating results table
d = {"Chr":chr, 
     "Position":pos,
     "Ref":ref,
     "Alt":alt,
     "rs ID":rs,
     "Gene_IDs":unique_clean,
     "ANNOVAR":AN_ID_check,
     "SnpEff":SN_ID_check,
     "VEP":VP_ID_check}

results = pd.DataFrame(d)
results.head()
all_agree = 0 #with GeneID union
two_agree = 0 #with GeneID union
one_agree = 0 #with GeneID union
none_agree= 0 #with GeneID union
AN_missing = 0 #
SN_missing = 0 #
VP_missing = 0

for i in range(0,size):
  zize = len(unique_clean[i])
  temp = 0
  if zize == len(AN_ID_check[i]):
    temp += 1
  if zize == len(SN_ID_check[i]):
    temp += 1
  if zize == len(VP_ID_check[i]):
    temp += 1
  if temp == 3:             # this means all should agree 
    all_agree += 1
  elif temp == 2:           # exactly two agree
    two_agree += 1
  elif temp == 1:           # exactly one agrees
    one_agree += 1
  elif temp == 0:           # no tool has all GeneIDs
    none_agree+= 1          # can hardly see it happening
all_agree
size
two_agree
one_agree
none_agree<jupyter_output><empty_output><jupyter_text>Results File Creation<jupyter_code># Naming output file assuming .gz compression
filename_1 = filename.split("/")
filename_1 = filename_1[-1]
filename_1 = filename_1.split(".")
result_filename = filename_1[0]+"_processed."+filename_1[1] #adding the file extension

# Creating a .txt file with results inside
#results.to_csv(path_or_buf=('/content/'+result_filename), sep="\t", index=False)<jupyter_output><empty_output><jupyter_text>Quantifying Differences<jupyter_code># Comparing SN vs AN
dif_SN_AN = []
dif_SN_AN = difference(SN_ID_clean, AN_ID_clean)
dif_SN_AN_sum = counter(dif_SN_AN)

# Comparing SN vs VP
dif_SN_VP = []
dif_SN_VP = difference(SN_ID_clean, VP_ID_clean)
dif_SN_VP_sum = counter(dif_SN_VP)

# Comparing AN vs SN
dif_AN_SN = []
dif_AN_SN = difference(AN_ID_clean, SN_ID_clean)
dif_AN_SN_sum = counter(dif_AN_SN)

# Comparing AN vs VP
dif_AN_VP = []
dif_AN_VP = difference(AN_ID_clean, VP_ID_clean)
dif_AN_VP_sum = counter(dif_AN_VP)

# Comparing VP vs SN
dif_VP_SN = []
dif_VP_SN = difference(VP_ID_clean, SN_ID_clean)
dif_VP_SN_sum = counter(dif_VP_SN)

# Comparing VP vs AN
dif_VP_AN = []
dif_VP_AN = difference(VP_ID_clean, AN_ID_clean)
dif_VP_AN_sum = counter(dif_VP_AN)<jupyter_output><empty_output><jupyter_text>Summary File Creation<jupyter_code># Summary results
l_1 = "For Chromosome: " + chr[1] + "\n"
l_8 = "There are: " + str(size) + " SNPs\n"
l_9 = "All tools agree with all GeneIDs: " + str(all_agree) + "\n"
l_10= "Two tools agree with all GeneIDs: " + str(two_agree) + "\n"
l_11= "One tool agrees with all GeneIDs: " + str(one_agree) + "\n"
l_12= "No tool agrees with all GeneIDs:  " + str(none_agree) + "\n"
l_2 = str(dif_AN_SN_sum) + "\t in ANNOVAR \t but not in \t SnpEff \n"
l_3 = str(dif_AN_VP_sum) + "\t in ANNOVAR \t but not in \t VEP \n"
l_4 = str(dif_SN_AN_sum) + "\t in SnpEff \t but not in \t ANNOVAR \n"
l_5 = str(dif_SN_VP_sum) + "\t in SnpEff \t but not in \t VEP \n"
l_6 = str(dif_VP_AN_sum) + "\t in VEP    \t but not in \t ANNOVAR \n"
l_7 = str(dif_VP_SN_sum) + "\t in VEP    \t but not in \t SnpEff \n"


# Creating a .txt file with the summary results
chromo = filename_1[0]+"_summary."+filename_1[1]
f = open('/content/'+ chromo,"w+")
f.writelines([l_1, l_8, "\n", l_9, l_10, l_11, l_12, "\n", l_2, l_3, l_4, l_5, l_6, l_7])
f.close()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12252141982864137, 'fraction_numerical': 0.01383108935128519, 'mean_word_length': 3.9671732522796352, 'pattern_counts': {'":': 9, '<': 36, '<?xml version=': 0, '>': 38, 'https://': 2, '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': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '624642', 'license': 'Apache-2.0', 'n_tokens_mistral': 3269, 'n_tokens_neox': 3118, 'n_words': 785, 'path': 'AnnoQ_comparison_IDs_v2.ipynb', 'repo_name': 'quemeb/AnnoQ_Gene-ID_comparison', 'size': '8170'}