OSCR

Extensive genetic interactions (epistasis) linked to alcohol use disorder in a high-risk population.

Code ↔ Paper

2 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 2 matches
  1. [1] § Methods › Bi-clustering of SNP x SNP epistasis analysis results ↔ BiClustering.py, lines 1127–1192 · score 0.52 · chromosome pair, bi clustering, smaller, conservative, hypergeometric, intervals
  2. [2] § Methods › Replication in all of us research database ↔ hg19_to_hg38_ranges.py, lines 408–433 · score 0.51 · interacting intervals, ACAF, hg38, boundaries, hg19, mapped

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,800 lines · 80 KB · no license · 1 match

  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Mon Feb 19 02:54:48 2024
  5. @author: staslist
  6. """
  8. import re
  9. import csv
  10. import unittest
  11. import numpy as np
  12. import os
  13. #import time
  14. from scipy.stats import hypergeom
  15. from datetime import datetime
  16. from numba import jit
  18. # Bi-Clustering Algorithm Code
  20. # Assume I have a file that lists all indices in the matrix that are 1.
  22. # GLOBAL VARIABLES
  24. # For tested interaction matrix, assume that every interaction has been tested (exhaustive algorithm)
  26. def merge_list(input_list:list):
  27. output = []
  28. eliminated_indeces = set()
  29. outer_index = 0
  30. while outer_index < len(input_list):
  31. if(outer_index in eliminated_indeces):
  32. outer_index += 1
  33. continue
  35. #to_merge_indeces = []
  36. current_set = input_list[outer_index]
  37. inner_index = 0
  38. while inner_index < len(input_list):
  39. if(inner_index in eliminated_indeces):
  40. inner_index += 1
  41. continue
  43. if(input_list[outer_index] == input_list[inner_index]):
  44. inner_index += 1
  45. continue
  47. if(len(input_list[outer_index].intersection(input_list[inner_index])) > 0):
  48. #to_merge_indeces.append(inner_index)
  49. current_set = current_set.union(input_list[inner_index])
  50. eliminated_indeces.add(inner_index)
  52. inner_index += 1
  54. output.append(current_set)
  56. outer_index += 1
  58. #print(output)
  59. if(output == input_list):
  60. return output
  61. else:
  62. return merge_list(output)
  64. def generate_gene_networks(gene_inter_file:str):
  65. # Convert a collection of gene pairs into a set of gene networks
  66. # The interacting gene pairs are stored in a text file, one gene pair per line
  68. gene_pairs = []
  69. gene_networks = []
  70. with open(gene_inter_file) as csv_file:
  71. csv_reader = csv.reader(csv_file, delimiter=',')
  72. for row in csv_reader:
  73. gene_pairs.append({row[0], row[1]})
  75. gene_networks = merge_list(gene_pairs)
  77. return gene_networks
  79. def check_collection_overlap(collection1, collection2):
  80. for ele in collection1:
  81. if(ele in collection2):
  82. return True
  84. return False
  86. def generate_gene_blocks(fname:str, distance:int):
  87. # Take in a plink range file and aggregate together genes/regulatory elements
  88. # that are located closely together into blocks.
  89. # Write out blocks, one per line, into a new file.
  91. ranges = dict()
  92. pairs = []
  93. blocks = []
  95. with open(fname) as csv_file:
  96. csv_reader = csv.reader(csv_file, delimiter=' ')
  97. for row in csv_reader:
  98. ranges[row[3]] = ((row[0]), int(row[1]), int(row[2]))
  101. for k,v in ranges.items():
  102. chrom1 = v[0]
  103. start1 = v[1]
  104. end1 = v[2]
  105. paired = False
  106. for k2,v2 in ranges.items():
  107. if(v == v2):
  108. continue
  110. #print(k, k2)
  112. chrom2 = v2[0]
  113. start2 = v2[1]
  114. end2 = v2[2]
  116. # Detect one element within another element
  117. if(chrom1 == chrom2 and ((start1 < start2 < end1) or (start1 < end2 < end1)
  118. or (start2 < start1 < end2) or (start2 < end1 < end2))):
  119. pairs.append((k,k2))
  120. paired = True
  121. # Detect neighboring elements
  122. elif(chrom1 == chrom2 and (abs(start2 - end1)<=distance or abs(start1 - end2)<=distance)):
  123. pairs.append((k,k2))
  124. paired = True
  126. if(not paired):
  127. blocks.append({k})
  129. #print("Pairs: ", pairs)
  132. used_pairs = []
  133. # Now conglomerate the pairs into blocks
  134. for pair in pairs:
  135. reverse_pair = (pair[1], pair[0])
  136. if(pair in used_pairs or reverse_pair in used_pairs):
  137. continue
  139. curr_block = {pair[0], pair[1]}
  140. for pair2 in pairs:
  141. if(pair2 in used_pairs):
  142. pass
  143. elif(pair == pair2):
  144. pass
  145. elif(check_collection_overlap(curr_block, pair2)):
  146. curr_block.add(pair2[0])
  147. curr_block.add(pair2[1])
  149. used_pairs.append(pair2)
  152. used_pairs.append(pair)
  153. blocks.append(curr_block)
  155. return blocks
  158. def summation_func(n:int):
  159. # Sum (n-1) + (n-2) + ... + 2 + 1 + 0
  160. assert(n > 1)
  161. result = 0
  162. n = n-1
  163. while n > 0:
  164. result += n
  165. n -= 1
  167. return result
  169. def generate_multicpu_files(anno_files:list, bim_file:str, total_markers:int, num_cpus:int,
  170. out_dir:str, chrom1:str, chrom2:str, pval_cutoff:float):
  171. interval_length = total_markers//num_cpus
  173. i_start = 0
  174. i_end = total_markers
  175. current_i = 0
  176. i_intervals = []
  177. while current_i < total_markers:
  178. if( (current_i + interval_length) < total_markers):
  179. i_intervals.append((current_i, current_i + interval_length))
  180. current_i = current_i + interval_length
  181. else:
  182. i_intervals.append((current_i, total_markers))
  183. current_i = total_markers
  185. counter = 0
  186. for i_interval in i_intervals:
  187. i_s = str(i_interval[0])
  188. i_e = str(i_interval[1])
  189. tot_marks = str(total_markers)
  190. fname_out_py = out_dir + 'BiClustering_Launcher_chr' + chrom1 + '_chr' + chrom2
  191. fname_out_py += '_' + str(counter) + '.py'
  192. with open(fname_out_py, 'w') as writer:
  193. writer.write('from BiClustering import *\n')
  194. writer.write('import sys\n')
  195. writer.write('if __name__ == "__main__":\n')
  196. writer.write('\tassert sys.version_info.major == 3\n')
  197. writer.write('\tassert sys.version_info.minor >= 7\n')
  198. writer.write("\tindir = '/gpfs/group/home/slistopad/BiClustering/'\n")
  199. r = 1
  200. for anno_file in anno_files:
  201. writer.write("\tmarker_file" + str(r) + " = indir + '" + str(anno_file) + "'\n")
  202. r += 1
  203. r = 1
  204. writer.write('\tmarker_files = [')
  205. while r <= len(anno_files):
  206. writer.write('marker_file' + str(r))
  207. if(r != len(anno_files)):
  208. writer.write(',')
  209. r += 1
  210. writer.write(']\n')
  211. writer.write("\tbim_file = indir + '" + str(bim_file) + "'\n")
  212. writer.write("\tchrom1,chrom2 = '" + chrom1 + "','" + chrom2 + "'\n")
  213. writer.write('\tN,n,inter_matrix,up_tri = initialize_matrices2(bim_file,marker_files, chrom1, chrom2, ' + str(pval_cutoff) + ')\n')
  214. writer.write("\tout_dir = '/gpfs/group/home/slistopad/BiClustering/'\n")
  215. writer.write('\ti_start, i_end = '+i_s+', '+i_e+'\n')
  216. writer.write('\tkm_results = compute_k_m_parallel(60, inter_matrix, up_tri, i_start, i_end, N, n)\n')
  217. writer.write('\tcompute_interval_pval_parallel(km_results, N, n, i_start, i_end, 60, out_dir, chrom1=chrom1, chrom2=chrom2)\n')
  219. if(counter == 0):
  220. writer.write('\tpval_results = dict()\n')
  221. writer.write('\ti_intervals = [')
  222. counter_inner = 0
  223. for i_inter in i_intervals:
  224. counter_inner += 1
  225. writer.write('('+str(i_inter[0])+','+str(i_inter[1])+')')
  226. if(counter_inner < len(i_intervals)):
  227. writer.write(',')
  228. else:
  229. writer.write(']\n')
  230. writer.write('\tfor i_pair in i_intervals:\n')
  231. writer.write("\t\tfilename = out_dir + 'chr" + chrom1 + "_chr" + chrom2 + "_pval_results_' + str(i_pair[0]) + '_' + str(i_pair[1]) + '.csv'\n")
  232. writer.write('\t\twith open(filename) as csv_file:\n')
  233. writer.write("\t\t\tcsv_reader = csv.reader(csv_file, delimiter=',')\n")
  234. writer.write("\t\t\tfor row in csv_reader:\n")
  235. writer.write("\t\t\t\tpval_results[(int(row[0]),int(row[1]),int(row[2]),int(row[3]))] = float(row[6])\n")
  237. writer.write("\ts_inter_pairs = sorted(pval_results.items(), key = lambda x: abs(x[1]), reverse = False)\n")
  238. writer.write("\ttrimmed_inter_pairs = trim_intervals(s_inter_pairs)\n")
  239. writer.write("\tfilename = out_dir + 'biclustering_results_chr" + chrom1 + "_chr" + chrom2 + ".txt'\n")
  240. writer.write("\twith open(filename, 'w') as writer:\n")
  241. writer.write("\t\tfor inter_pair in trimmed_inter_pairs:\n")
  242. writer.write("\t\t\twriter.write(str(inter_pair[0][0])+','+str(inter_pair[0][1])+','+str(inter_pair[0][2])+','+str(inter_pair[0][3])+','+str(inter_pair[1]) + '\\n')\n")
  244. counter += 1
  246. fname_out_sh = out_dir + 'BiClustering_Launcher_chr' + chrom1 + '_chr' + chrom2 + '.sh'
  247. with open(fname_out_sh, 'w') as writer:
  248. writer.write("#!/bin/sh\n")
  249. writer.write("#SBATCH --job-name=SL_BiClustering_Launcher_chr" + chrom1 + '_chr' + chrom2 + "\n")
  250. writer.write("#SBATCH --array=1-" + str(len(i_intervals)-1) + "\n")
  251. writer.write("#SBATCH --ntasks=1\n")
  252. writer.write("#SBATCH --mem=48gb\n")
  253. writer.write("#SBATCH --time=240:00:00\n")
  254. writer.write("#SBATCH --partition=shared\n\n")
  255. writer.write("cd $SLURM_SUBMIT_DIR\n")
  256. writer.write("module load use.own\n")
  257. writer.write("module load python/3.8.3\n")
  258. writer.write("export PYTHONPATH=/gpfs/home/slistopad/.local/lib/python3.8/site-packages:$PYTHONPATH\n")
  259. writer.write("python3 " + 'BiClustering_Launcher_chr' + chrom1 + '_chr' + chrom2 + '_${SLURM_ARRAY_TASK_ID}.py' + '\n')
  260. writer.write("python3 " + 'BiClustering_Launcher_chr' + chrom1 + '_chr' + chrom2 + '_0.py' + '\n')
  263. def get_number_of_interacting_snps_in_interval(biclustering_file:str, bim_file:str, gene_location_file:str,
  264. chrom1:str, chrom2:str, interaction_files:list, pval_cutoff:float):
  265. interacting_pairs_per_interval = dict()
  267. chrom1_snps, chrom2_snps = parse_plink_bim_file(bim_file, chrom1, chrom2)
  269. chrom1_ranges, chrom2_ranges = [], []
  271. with open(biclustering_file) as csv_file:
  272. csv_reader = csv.reader(csv_file, delimiter=',')
  273. for row in csv_reader:
  274. range1 = (int(row[0]), int(row[0]) + int(row[2]))
  275. chrom1_ranges.append(range1)
  276. range2 = (int(row[1]), int(row[1]) + int(row[3]))
  277. chrom2_ranges.append(range2)
  280. interactions = parse_remma_interaction_data(interaction_files, pval_cutoff, chrom1, chrom2)
  282. #print('goodbye')
  284. print('chrom1_ranges: ', chrom1_ranges)
  285. print('chrom2_ranges: ', chrom2_ranges)
  287. gene_reg_locations = []
  289. with open(gene_location_file) as csv_file:
  290. csv_reader = csv.reader(csv_file, delimiter=' ')
  291. for row in csv_reader:
  292. gene_reg_locations.append((row[0],row[1],row[2],row[3]))
  294. i = 0
  295. num_interval_pairs = len(chrom1_ranges)
  296. while i < num_interval_pairs:
  297. interval1 = chrom1_snps[chrom1_ranges[i][0]:chrom1_ranges[i][1]]
  298. interval2 = chrom2_snps[chrom2_ranges[i][0]:chrom2_ranges[i][1]]
  299. #print(len(interval1))
  300. #print(len(interval2))
  302. #print(interval1)
  303. #print(interval2)
  305. #print(interactions)
  307. for snp in interval1:
  308. for snp2 in interval2:
  309. if( (snp, snp2) in interactions or (snp2, snp) in interactions ):
  311. chrom1_loc1 = snp.split(':')
  312. chrom1 = chrom1_loc1[0]
  313. loc1 = chrom1_loc1[1]
  314. element1 = ''
  316. for ele_loc in gene_reg_locations:
  317. if(chrom1 == ele_loc[0] and int(loc1) >= int(ele_loc[1]) and int(loc1) <= int(ele_loc[2])):
  318. element1 += ele_loc[3] + ' '
  320. #print(element1)
  322. chrom2_loc2 = snp2.split(':')
  323. chrom2 = chrom2_loc2[0]
  324. loc2 = chrom2_loc2[1]
  325. element2 = ''
  328. for ele_loc in gene_reg_locations:
  329. if(chrom2 == ele_loc[0] and int(loc2) >= int(ele_loc[1]) and int(loc2) <= int(ele_loc[2])):
  330. element2 += ele_loc[3] + ' '
  332. #print(element2)
  334. if (element1, element2) not in interacting_pairs_per_interval:
  335. interacting_pairs_per_interval[(element1, element2)] = 1
  336. else:
  337. interacting_pairs_per_interval[(element1, element2)] += 1
  338. i += 1
  339. print(interacting_pairs_per_interval)
  341. return interacting_pairs_per_interval
  343. def get_msigdb_enrichment(gene_pairs:set, out_dir:str):
  344. #print(pairs)
  345. pairs_map = dict()
  346. fname = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Pathways/c2.all.v2023.2.Hs.symbols.gmt'
  347. with open(fname) as csv_file:
  348. csv_reader = csv.reader(csv_file, delimiter='\t')
  349. for row in csv_reader:
  350. pathway_name = row[0]
  351. pathway_genes = row[2:]
  352. for pair in gene_pairs:
  353. if(pair[0] in pathway_genes and pair[1] in pathway_genes):
  354. try:
  355. pairs_map[pair].append(pathway_name)
  356. except KeyError:
  357. pairs_map[pair] = [pathway_name]
  359. fname_out = out_dir + 'biclustering_results_msigdb_enrichment.txt'
  360. with open(fname_out, 'w') as writer:
  361. for k,v in pairs_map.items():
  362. writer.write(str(k) + ' : ' + str(v) + '\n')
  364. def parse_biclustering_annotation(bicluster_file_annot:str):
  365. element_pairs_dict = dict()
  366. element_pairs_annotated = []
  367. element_pairs_per_interval = []
  368. unique_interacting_snps = []
  369. interacting_snp_pairs_all_intervals = []
  370. regulatory_elements_all = set()
  371. with open(bicluster_file_annot) as csv_file:
  372. csv_reader = csv.reader(csv_file, delimiter='\t')
  373. for row in csv_reader:
  374. # If the row is not empty
  375. if(len(row) == 0):
  376. continue
  377. if(len(row[0]) > 0 and 'Interacting snp pairs:' not in row[0] and
  378. 'Interacting num snp pairs:' not in row[0] and 'Total number of unique' not in row[0]):
  379. # There must be at least one gene/reg element pair
  380. # It must be surrounded by (" ")
  381. interval_pair = row[0]
  382. interval_pair_split = interval_pair.split("),")
  383. pval = row[1]
  384. element_pairs_annotated2 = []
  385. for element_pair in interval_pair_split:
  386. #print(element_pair)
  387. element2_begin = element_pair.find('element2:')
  388. r = 0
  389. record_element = False
  390. element1_set = set()
  391. element1 = ''
  392. while r < (element2_begin-1):
  393. #print(element_pair)
  394. #print(i)
  395. char = element_pair[r]
  396. if(char == "'"):
  397. if(record_element):
  398. if('GH' in element1 and len(element1) == 11):
  399. regulatory_elements_all.add(element1)
  400. element1 = ''
  401. else:
  402. element1_set.add(element1)
  403. element1 = ''
  404. record_element = not record_element
  405. if(record_element and char!= "'"):
  406. element1 += char
  407. r += 1
  409. #print('Element1: ', element1_list)
  411. r = element2_begin
  412. record_element = False
  413. element2_set = set()
  414. element2 = ''
  415. while r < len(element_pair):
  416. #print(interval_pair)
  417. #print(r)
  418. #print(element_pair_end-2)
  419. #print(char)
  420. char = element_pair[r]
  421. #print(char)
  422. if(char == "'"):
  423. if(record_element):
  424. if('GH' in element2 and len(element2) == 11):
  425. regulatory_elements_all.add(element2)
  426. element2 = ''
  427. else:
  428. element2_set.add(element2)
  429. element2 = ''
  430. record_element = not record_element
  431. if(record_element and char!= "'"):
  432. element2 += char
  433. r += 1
  435. #print('Element2: ', element2_list)
  436. element_pair_tuple = (element1_set, element2_set)
  437. if(element_pair_tuple not in element_pairs_annotated):
  438. element_pairs_annotated.append((element1_set, element2_set))
  439. if(element_pair_tuple not in element_pairs_annotated2):
  440. element_pairs_annotated2.append((element1_set, element2_set))
  442. element_pair_str = str((element1_set, element2_set))
  443. element_pair_str_reverse = str((element2_set, element1_set))
  444. if(element_pair_str in element_pairs_dict):
  445. if(pval < element_pairs_dict[element_pair_str]):
  446. element_pairs_dict[element_pair_str] = pval
  447. elif(element_pair_str_reverse in element_pairs_dict):
  448. if(pval < element_pairs_dict[element_pair_str_reverse]):
  449. element_pairs_dict[element_pair_str_reverse] = pval
  450. else:
  451. element_pairs_dict[element_pair_str] = pval
  452. element_pairs_per_interval.append(element_pairs_annotated2)
  454. elif('Total number of unique' in row[0]):
  455. split_row = row[0].split(';')[1:]
  456. #print(split_row)
  457. unique_interacting_snps_in_interval = set()
  458. for unique_snp in split_row:
  459. if(':' in unique_snp):
  460. unique_interacting_snps_in_interval.add(unique_snp)
  461. unique_interacting_snps.append(unique_interacting_snps_in_interval)
  463. elif('Interacting snp pairs:' in row[0]):
  464. current_row = row[0]
  465. pattern = "('\d+:\d+', '\d+:\d+')"
  466. matches = re.findall(pattern, current_row)
  467. interacting_snp_pairs_all = set()
  468. for snp_pair_string in matches:
  469. snp_pair_string = snp_pair_string.replace("'", "")
  470. match_split = snp_pair_string.split(', ')
  471. #print(match_split)
  472. if((match_split[1], match_split[0]) not in interacting_snp_pairs_all):
  473. interacting_snp_pairs_all.add( (match_split[0], match_split[1]) )
  474. interacting_snp_pairs_all_intervals.append(interacting_snp_pairs_all)
  476. return element_pairs_dict, element_pairs_annotated, unique_interacting_snps, interacting_snp_pairs_all_intervals, element_pairs_per_interval, regulatory_elements_all
  478. def parse_biclustering_results(biclustering_file:str, bim_file:str, gene_location_file:str,
  479. gene_hancer_file:str, chrom1:str, chrom2:str, interaction_files:list,
  480. pval_cutoff:float, out_dir:str, single_file_output:bool = False,
  481. snp_list_output:bool = False, strict_annotation:bool = True):
  482. # WARNING, The pval_cutoff should match the p-value cutoff used to conduct the biclustering analysis.
  483. # Otherwise, the mapping of the interacting intervals to interacting genes will not make sense.
  484. # strict annotation = only map interacting snp pairs to genes/reg elements,
  485. # loose annotation = map all snp pairs contained within interacting interval to genes/reg elements
  487. chrom1_snps, chrom2_snps = parse_plink_bim_file(bim_file, chrom1, chrom2)
  489. chrom1_ranges, chrom2_ranges = [], []
  490. p_values = []
  492. with open(biclustering_file) as csv_file:
  493. csv_reader = csv.reader(csv_file, delimiter=',')
  494. for row in csv_reader:
  495. range1 = (int(row[0]), int(row[0]) + int(row[2]))
  496. chrom1_ranges.append(range1)
  497. range2 = (int(row[1]), int(row[1]) + int(row[3]))
  498. chrom2_ranges.append(range2)
  499. p_values.append(row[4])
  502. # modify this section depending on the gene_location_file used
  503. # herein I assume usage of ucsc-hg19
  504. # first read in and map the transcripts to the gene names
  505. transcript_to_gene_name = dict()
  506. kgXref_file = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Gene_Annotations/ucsc-hg19-kgXref.txt'
  507. with open(kgXref_file) as csv_file:
  508. csv_reader = csv.reader(csv_file, delimiter='\t')
  509. for row in csv_reader:
  510. transcript_to_gene_name[row[0]] = row[1]
  512. gene_reg_locations = []
  513. with open(gene_location_file) as csv_file:
  514. csv_reader = csv.reader(csv_file, delimiter='\t')
  515. for row in csv_reader:
  516. gene_reg_locations.append((row[1][3:],row[2],row[3],transcript_to_gene_name[row[0]]))
  518. # We also need to account for all the regulatory elements located within the genes
  519. # that were used in epistasis analysis
  521. with open(gene_hancer_file) as csv_file:
  522. csv_reader = csv.reader(csv_file, delimiter='\t')
  523. first_line = True
  524. for row in csv_reader:
  525. if(first_line):
  526. first_line = False
  527. continue
  528. gene_reg_locations.append((row[0][3:],row[1],row[2],row[3]))
  530. #print("gene_reg_locations: ", gene_reg_locations)
  532. #print('chrom1_ranges: ', chrom1_ranges)
  533. #print('chrom2_ranges: ', chrom2_ranges)
  535. interactions = parse_remma_interaction_data(interaction_files, pval_cutoff, chrom1, chrom2)
  537. #print("Interactions: ", interactions)
  539. if(not single_file_output):
  540. if(not snp_list_output):
  542. if(strict_annotation):
  543. filename = out_dir + 'biclustering_results_chr' + chrom1 + '_chr' + chrom2 + '_annotated.txt'
  544. else:
  545. filename = out_dir + 'biclustering_results_chr' + chrom1 + '_chr' + chrom2 + '_annotated_loosely.txt'
  546. with open(filename, 'w') as writer:
  547. i = 0
  548. num_interval_pairs = len(chrom1_ranges)
  549. while i < num_interval_pairs:
  551. gene_pair_to_num_interacting_snp_pairs = dict()
  552. gene_pair_to_interacting_snp_pairs = dict()
  554. gene_pairs = set()
  555. interval1 = chrom1_snps[chrom1_ranges[i][0]:chrom1_ranges[i][1]]
  556. interval2 = chrom2_snps[chrom2_ranges[i][0]:chrom2_ranges[i][1]]
  557. p_val = p_values[i]
  559. #print("Interval1: ", interval1)
  560. #print("Interval2: ", interval2)
  562. if(strict_annotation):
  563. for snp in interval1:
  564. for snp2 in interval2:
  565. #print(snp, snp2)
  566. if( (snp, snp2) in interactions or (snp2, snp) in interactions ):
  567. # Need to identify number of interacting snp pairs belonging to each gene pair
  569. #print("Found SNP pair in interactions!")
  570. chrom1_loc1 = snp.split(':')
  571. chrom1 = chrom1_loc1[0]
  572. loc1 = chrom1_loc1[1]
  573. element1 = set()
  575. #print('chrom1: ', chrom1)
  576. #print('loc1: ', loc1)
  578. element_name1 = 'intergenic'
  579. for ele_loc in gene_reg_locations:
  580. if(chrom1 == ele_loc[0] and int(loc1) >= int(ele_loc[1]) and int(loc1) <= int(ele_loc[2])):
  581. # There are a small number of overlapping genes, these are the
  582. # poorly defined genes, such as MICB and HLA-C, for whom the exact
  583. # loci are not yet fully established.
  584. element_name1 = ele_loc[3]
  585. element1.add(element_name1)
  587. #print('Element1: ', element1)
  588. if(len(element1) == 0):
  589. element1.add(element_name1)
  591. chrom2_loc2 = snp2.split(':')
  592. chrom2 = chrom2_loc2[0]
  593. loc2 = chrom2_loc2[1]
  594. element2 = set()
  596. #print('chrom2: ', chrom2)
  597. #print('loc2: ', loc2)
  599. element_name2 = 'intergenic'
  600. for ele_loc in gene_reg_locations:
  601. if(chrom2 == ele_loc[0] and int(loc2) >= int(ele_loc[1]) and int(loc2) <= int(ele_loc[2])):
  602. element_name2 = ele_loc[3]
  603. element2.add(element_name2)
  605. #print('Element2: ', element2)
  606. if(len(element2) == 0):
  607. element2.add(element_name2)
  609. for ele_name1 in element1:
  610. for ele_name2 in element2:
  611. if((ele_name1, ele_name2) not in gene_pair_to_num_interacting_snp_pairs and
  612. (ele_name2, ele_name1) not in gene_pair_to_num_interacting_snp_pairs):
  613. gene_pair_to_num_interacting_snp_pairs[(ele_name1, ele_name2)] = 1
  614. gene_pair_to_interacting_snp_pairs[(ele_name1, ele_name2)] = [(snp, snp2)]
  615. elif((ele_name1, ele_name2) in gene_pair_to_num_interacting_snp_pairs):
  616. gene_pair_to_num_interacting_snp_pairs[(ele_name1, ele_name2)] += 1
  617. gene_pair_to_interacting_snp_pairs[(ele_name1, ele_name2)].append((snp, snp2))
  618. elif((ele_name2, ele_name1) in gene_pair_to_num_interacting_snp_pairs):
  619. gene_pair_to_num_interacting_snp_pairs[(ele_name2, ele_name1)] += 1
  620. gene_pair_to_interacting_snp_pairs[(ele_name2, ele_name1)].append((snp2, snp))
  622. #print(gene_pair_to_num_interacting_snp_pairs)
  623. gene_pairs.add( ('element1: '+ str(element1), 'element2: ' + str(element2)) )
  624. writer.write(str(gene_pairs) + '\t')
  625. writer.write(p_values[i] + '\n')
  626. writer.write('Interacting num snp pairs: ' + str(gene_pair_to_num_interacting_snp_pairs) + '\n')
  627. writer.write('Interacting snp pairs: ' + str(gene_pair_to_interacting_snp_pairs) + '\n')
  628. current_union = set()
  629. for k,v in gene_pair_to_interacting_snp_pairs.items():
  630. current_union = current_union | set(v)
  631. writer.write('Total number of unique interacting pairs: ' + str(len(current_union)) + ';')
  632. unique_interacting_snps = set()
  633. for unique_snp_pair in current_union:
  634. unique_interacting_snps.add(unique_snp_pair[0])
  635. unique_interacting_snps.add(unique_snp_pair[1])
  636. for unique_snp in unique_interacting_snps:
  637. writer.write(str(unique_snp) + ';')
  638. writer.write('\n\n')
  639. else:
  640. elements1 = set()
  641. for snp in interval1:
  642. chrom1_loc1 = snp.split(':')
  643. chrom1 = chrom1_loc1[0]
  644. loc1 = chrom1_loc1[1]
  645. element1 = set()
  647. for ele_loc in gene_reg_locations:
  648. if(chrom1 == ele_loc[0] and int(loc1) >= int(ele_loc[1]) and int(loc1) <= int(ele_loc[2])):
  649. element1.add(ele_loc[3])
  650. #print(element1)
  652. elements1.add(str(element1))
  653. elements2 = set()
  654. for snp2 in interval2:
  655. chrom2_loc2 = snp2.split(':')
  656. chrom2 = chrom2_loc2[0]
  657. loc2 = chrom2_loc2[1]
  658. element2 = set()
  659. #print(loc2)
  661. for ele_loc in gene_reg_locations:
  662. if(chrom2 == ele_loc[0] and int(loc2) >= int(ele_loc[1]) and int(loc2) <= int(ele_loc[2])):
  663. element2.add(ele_loc[3])
  664. elements2.add(str(element2))
  666. writer.write(str( (elements1,elements2) ) + '\t')
  667. writer.write(p_values[i] + '\n')
  669. i += 1
  670. else:
  671. #print("HELLO!")
  672. filename = out_dir + 'biclustering_results_chr' + chrom1 + '_chr' + chrom2 + '_snplist_interacting_pairs.txt'
  673. with open(filename, 'w') as writer:
  674. num_interval_pairs = len(chrom1_ranges)
  675. i = 0
  676. while i < num_interval_pairs:
  677. gene_pairs = set()
  678. interval1 = chrom1_snps[chrom1_ranges[i][0]:chrom1_ranges[i][1]]
  679. interval2 = chrom2_snps[chrom2_ranges[i][0]:chrom2_ranges[i][1]]
  681. # format for liftover tool
  682. snp1_info = interval1[0].split(':')
  683. snp1_chrom,snp1_loci = snp1_info[0],snp1_info[1]
  684. snp2_info = interval1[-1].split(':')
  685. snp2_chrom,snp2_loci = snp2_info[0],snp2_info[1]
  687. #writer.write(snp1_chrom + ' ' + snp1_loci + ' ' + snp2_loci + '\n')
  689. snp1_info = interval2[0].split(':')
  690. snp1_chrom,snp1_loci = snp1_info[0],snp1_info[1]
  691. snp2_info = interval2[-1].split(':')
  692. snp2_chrom,snp2_loci = snp2_info[0],snp2_info[1]
  694. #writer.write(snp1_chrom + ' ' + snp1_loci + ' ' + snp2_loci + '\n')
  696. for snp in interval1:
  697. for snp2 in interval2:
  698. #print(snp, snp2)
  699. if( (snp, snp2) in interactions or (snp2, snp) in interactions ):
  700. writer.write(snp + ' : ')
  701. writer.write(snp2 + '\n')
  703. snp1_info = snp.split(':')
  704. snp1_chrom,snp1_loci = snp1_info[0],snp1_info[1]
  706. snp2_info = snp2.split(':')
  707. snp2_chrom,snp2_loci = snp2_info[0],snp2_info[1]
  709. print("plink --bfile /gpfs/group/home/slistopad/REMMA/data/native_american/", end='')
  710. print('NA3_Combined_Fixed_Set_score900_db001_and_exp001_1_and_genehancer', end='')
  711. print('_score25_maf001_region_qc3 --ld ' + snp + ' ' + snp2 + ' --r2 ', end='')
  712. print('--out /gpfs/group/home/slistopad/REMMA/data/native_american/Biclustering_SNPLists/', end='')
  713. print('Combined_Fixed2_Set_Biclustering/biclustering_results_ld_' + snp1_chrom, end='')
  714. print('_' + snp1_loci + '_' + snp2_chrom + '_' + snp2_loci)
  716. i += 1
  718. else:
  719. # Note, what is put into the file, are the locations of the genes, which contain
  720. # interacting SNPs within the interacting intervals.
  722. # This file is made for generation of Circo plots. In it we don't want to see the same gene pairs even if they
  723. # belong to different interval pairs.
  724. filename = out_dir + 'biclustering_results_all.txt'
  725. with open(filename, 'a') as writer:
  726. i = 0
  727. tracker = set()
  728. num_interval_pairs = len(chrom1_ranges)
  729. while i < num_interval_pairs:
  730. gene_pairs = set()
  731. interval1 = chrom1_snps[chrom1_ranges[i][0]:chrom1_ranges[i][1]]
  732. interval2 = chrom2_snps[chrom2_ranges[i][0]:chrom2_ranges[i][1]]
  733. p_val = p_values[i]
  735. #print('interval1: ', interval1)
  736. #print('interval2: ', interval2)
  738. for snp in interval1:
  739. for snp2 in interval2:
  740. if( (snp, snp2) in interactions or (snp2, snp) in interactions ):
  742. chrom1_loc1 = snp.split(':')
  743. chrom1 = chrom1_loc1[0]
  744. loc1 = chrom1_loc1[1]
  745. element1 = ''
  747. for ele_loc in gene_reg_locations:
  748. if(chrom1 == ele_loc[0] and int(loc1) >= int(ele_loc[1]) and int(loc1) <= int(ele_loc[2])):
  749. element1 = 'chr' + ele_loc[0] + ',' + ele_loc[1] + ',' + ele_loc[2]
  751. chrom2_loc2 = snp2.split(':')
  752. chrom2 = chrom2_loc2[0]
  753. loc2 = chrom2_loc2[1]
  754. element2 = ''
  756. for ele_loc in gene_reg_locations:
  757. if(chrom2 == ele_loc[0] and int(loc2) >= int(ele_loc[1]) and int(loc2) <= int(ele_loc[2])):
  758. element2 = 'chr' + ele_loc[0] + ',' + ele_loc[1] + ',' + ele_loc[2]
  760. gene_pairs.add((element1, element2))
  762. #print(gene_pairs)
  763. #print(p_val)
  764. for pair in gene_pairs:
  765. if(pair not in tracker):
  766. writer.write(pair[0] + ',' + pair[1] + ',' + str(p_val) + '\n')
  767. tracker.add(pair)
  768. i += 1
  771. def parse_remma_interaction_data(interaction_files:list, pval_cutoff:float, chrom1:str, chrom2:str):
  772. #print(pval_cutoff)
  773. #print(chrom1)
  774. #print(chrom2)
  776. interactions = dict()
  777. for interaction_file in interaction_files:
  778. with open(interaction_file) as csv_file:
  779. csv_reader = csv.reader(csv_file, delimiter=' ')
  780. header = True
  781. for row in csv_reader:
  782. if(header):
  783. header = False
  784. continue
  785. # Check if entry passes p-value cutoff and matches target chromosomes.
  786. if(float(row[18]) < pval_cutoff and ((row[1] == chrom1 and row[8] == chrom2) or
  787. (row[8] == chrom1 and row[1] == chrom2))):
  789. # If this entry already is recorded in interactions (presumably from another
  790. # interaction file), then overwrite the p-value for the record only if the
  791. # new p-value is smaller.
  792. if((row[2], row[9]) in interactions):
  793. if(float(row[18]) < interactions[(row[2], row[9])]):
  794. interactions[(row[2], row[9])] = float(row[18])
  795. elif((row[9], row[2]) in interactions):
  796. if(float(row[18]) < interactions[(row[9], row[2])]):
  797. interactions[(row[9], row[2])] = float(row[18])
  798. else:
  799. interactions[(row[2], row[9])] = float(row[18])
  800. return interactions
  802. def parse_plink_bim_file(bim_file:str, chrom1:str, chrom2:str):
  803. chrom1_snps = []
  804. chrom2_snps = []
  806. with open(bim_file) as csv_file:
  807. csv_reader = csv.reader(csv_file, delimiter='\t')
  808. for row in csv_reader:
  809. if(chrom1 == row[0]):
  810. chrom1_snps.append(row[1])
  811. if(chrom2 == row[0]):
  812. chrom2_snps.append(row[1])
  814. return chrom1_snps, chrom2_snps
  817. def initialize_matrices2(bim_file:str, interaction_files:list, chrom1:str, chrom2:str,
  818. pval_cutoff:float):
  819. # Assume that testing is always done using all relevant SNPs on one chromosome vs
  820. # all relevant SNPs on another chromosome. A chromosome can be also tested against itself.
  821. # The matrix size is number of SNPs on chromosome 1 vs number of SNPs on chromosome 2.
  822. # Relevant SNPs means SNPs that were included in the epistasis analysis.
  824. # Bim file informs us the number of relevant SNPs on each chromosome, and which SNPs pair
  825. # a given (i,j) value in the matrix represents.
  826. # Interaction file informs us which SNP pairs were found to interact.
  827. # Generally it is expected that the number of interactions is small.
  828. # Whenever the chromosomes are different, all SNP pairs are assumed to have been tested
  829. # for interaction.
  830. # If chromsome is the same, the tested interaction matrix is an upper triangular matrix.
  832. # Important assumption regarding bim file, for each chromosome all SNPs are ordered based
  833. # on their location, from beginning of chromosome to end of chromosome.
  835. # The interaction file is assumed to be a REMMA .anno file.
  837. chrom1_snps, chrom2_snps = parse_plink_bim_file(bim_file, chrom1, chrom2)
  839. interactions = parse_remma_interaction_data(interaction_files, pval_cutoff, chrom1, chrom2)
  841. #print(interactions)
  843. upper_triangular = False
  844. if(chrom1 == chrom2):
  845. upper_triangular = True
  847. N = 0
  848. n = 0
  850. total_markers1 = len(chrom1_snps)
  851. total_markers2 = len(chrom2_snps)
  853. if(upper_triangular):
  854. i = 0
  855. while i < total_markers1:
  856. N += i
  857. i += 1
  858. else:
  859. N = total_markers1 * total_markers2
  861. inter_array = np.zeros((total_markers1, total_markers2), dtype=int)
  862. i = 0
  863. while i < total_markers1:
  864. j = 0
  865. while j < total_markers2:
  866. if((chrom1_snps[i], chrom2_snps[j]) in interactions or
  867. (chrom2_snps[j], chrom1_snps[i]) in interactions):
  868. if(upper_triangular):
  869. if(i >= j):
  870. j += 1
  871. continue
  872. #print(i,j)
  873. inter_array[i][j] = 1
  874. n += 1
  876. j += 1
  878. i += 1
  880. print("Interaction matarix initialized.")
  881. print("N: ", N)
  882. print("n: ", n)
  883. print("Chromosome1: ", chrom1)
  884. print("Chromosome2: ", chrom2)
  885. print("P-Value Cutoff: ", pval_cutoff)
  887. now = datetime.now()
  888. current_time = now.strftime("%H:%M:%S")
  889. print("Current Time =", current_time)
  891. return N, n, inter_array, upper_triangular
  894. def initialize_matrices(total_markers:int, marker_file:str, upper_triangular:bool = True):
  895. # Also, we assume that all values in the interaction matrices are 1s and 0s.
  897. # We assume that all interaction and tested-interaction matrices are square.
  899. # We assume that all tested interaction matrices are either upper triangular or full (of 1s).
  900. # Upper triangular tested interaction matrix represents a set of marker being tested against
  901. # itself for epistasis.
  902. # Full tested interaction matrix represents two completely different sets of markers being
  903. # tested against each other for epistasis.
  905. N = 0
  906. n = 0
  907. inter_matrix = dict()
  909. if(upper_triangular):
  910. i = 0
  911. while i < total_markers:
  912. N += i
  913. i += 1
  914. #print("N = ", N)
  915. else:
  916. N = total_markers**2
  918. i = 0
  919. j = 0
  920. while i < total_markers:
  921. while j < total_markers:
  922. inter_matrix[(i,j)] = 0
  923. j += 1
  925. j = 0
  926. i += 1
  928. with open(marker_file) as csv_file:
  929. csv_reader = csv.reader(csv_file, delimiter=',')
  930. for row in csv_reader:
  931. i,j = int(row[0]), int(row[1])
  932. inter_matrix[(i,j)] = 1
  933. n += 1
  935. return N, n, inter_matrix, upper_triangular
  937. def check_overlap(i:int, j:int, a:int, b:int, i2:int, j2:int, a2:int, b2:int):
  938. # Note if two matrices are identical we return false.
  939. if(i == i2 and j == j2 and a==a2 and b==b2):
  940. return False
  942. assert(i >= 0 and j >= 0 and i2 >= 0 and j2 >= 0)
  943. assert(a >= 1 and b>= 1 and a2 >= 1 and b2 >= 1)
  945. origin_test = ( ((i+a-1) >= i2 >= i) and ((j+b-1) >= j2 >= j) ) or ( ((i2+a2-1) >= i >= i2) and ((j2+b2-1) >= j >= j2) )
  946. if(origin_test):
  947. return True
  948. bottom_left_test = ( (i2 <= (i+a-1) <= (i2+a2-1)) and ((j2+b2-1) >= j >= j2) ) or ( (i <= (i2+a2-1) <= (i+a -1)) and ((j+b-1) >= j2 >= j) )
  949. if(bottom_left_test):
  950. return True
  951. bottom_right_test = ( (i2 <= (i+a-1) <= (i2+a2-1)) and ((j2+b2-1) >= (j+b-1) >= j2) ) or ( (i <= (i2+a2-1) <= (i+a-1)) and ((j+b-1) >= (j2+b2-1) >= j) )
  952. if(bottom_right_test):
  953. return True
  954. upper_right_test = ( ((i+a-1) >= i2 >= i) and ((j2+b2-1) >= (j+b-1) >= j2) ) or ( ((i2+a2-1) >= i >= i2) and ((j+b-1) >= (j2+b2-1) >= j) )
  955. if(upper_right_test):
  956. return True
  958. return False
  960. '''
  961. def comp(i:int, j:int, a:int, b:int, val:str, inter_matrix:dict, tested_inter_matrix:dict):
  962. #print("i =",i,";j = ",j,";a = ", a, ";b = ", b)
  963. if(i < 0):
  964. raise ValueError("Invalid i value.")
  965. elif(j < 0):
  966. raise ValueError("Invalid j value.")
  967. elif(a < 0):
  968. raise ValueError("Invalid a value.")
  969. elif(b < 0):
  970. raise ValueError("Invalid b value.")
  972. if(a == 0 or b == 0):
  973. return 0
  974. elif(a == 1 and b == 1):
  975. if(val == 'k'):
  976. return inter_matrix[(i,j)]
  977. elif(val == 'm'):
  978. return tested_inter_matrix[(i,j)]
  979. else:
  980. raise ValueError("This parameter must be either 'k' or 'm'.")
  981. else:
  982. return comp(i,j,a-1,b-1,val,inter_matrix,tested_inter_matrix) + comp(i+a-1,j,1,b-1,val,inter_matrix,tested_inter_matrix) + comp(i,j+b-1,a-1,1,val,inter_matrix,tested_inter_matrix) + comp(i+a-1,j+b-1,1,1,val,inter_matrix,tested_inter_matrix)
  983. '''
  985. def convert_dict_to_2D_numpy_array(input_dict:dict):
  986. # This function converts a dictionary that maps (row_index, column_index) tuple to integer values.
  987. # It does so using the correct ordering.
  988. s_matrix_values = sorted(input_dict.items(), key = lambda x: x[0], reverse = False)
  989. TwoDList = []
  990. curr_inner = []
  991. prev_row = 0
  992. for ele in s_matrix_values:
  993. if(ele[0][0] > prev_row):
  994. TwoDList.append(curr_inner)
  995. curr_inner = [ele[1]]
  996. prev_row = ele[0][0]
  997. else:
  998. curr_inner.append(ele[1])
  1000. return np.array(TwoDList)
  1002. @jit(nopython=True)
  1003. def comp_alt(i:int, j:int, a:int, b:int, inter_matrix:'np.array', upper_triangular:bool = True,
  1004. compute_k:bool = True):
  1005. # Assume that tested interaction matrix is either upper triangular or full of 1s.
  1006. # In that case m_val can be computed quickly. Assume that one of the two is always true.
  1007. #print("i =",i,";j = ",j,";a = ", a, ";b = ", b)
  1008. if(i < 0):
  1009. raise ValueError("Invalid i value.")
  1010. elif(j < 0):
  1011. raise ValueError("Invalid j value.")
  1012. elif(a < 0):
  1013. raise ValueError("Invalid a value.")
  1014. elif(b < 0):
  1015. raise ValueError("Invalid b value.")
  1017. k_val = 0
  1018. m_val = 0
  1021. if(compute_k):
  1022. if(not upper_triangular):
  1023. m_val = a * b
  1024. else:
  1025. i2 = i
  1026. while i2 < (i+a):
  1027. # At which index do 1s start in this row. That is simpy row_index + 1.
  1028. if(j < (i2 + 1)):
  1029. m_val += max(j + b - (i2 + 1), 0)
  1030. else:
  1031. m_val += b
  1032. i2 += 1
  1034. k_val = np.count_nonzero(inter_matrix[i:(i+a),j:(j+b)])
  1036. '''
  1037. i2 = i
  1038. while i2 < (i+a):
  1039. j2 = j
  1040. while j2 < (j+b):
  1041. k_val += inter_matrix[i2,j2]
  1042. j2 += 1
  1043. i2 += 1
  1044. '''
  1046. return k_val, m_val
  1048. @jit(nopython=True)
  1049. def compute_k_m_parallel(max_inter_length:int, inter_array,
  1050. upper_triangular:bool, i_start:int, i_end:int, N:int, n:int):
  1052. # Our only knowledge about inter_matrix is that it consists of 1s and 0s
  1053. # and that 1s are typically very sparce.
  1055. # We have stronger assumptions about tested_inter_matrix which we deem to be either
  1056. # upper trinagular (with 1s and 0s) or completely filled with 1s.
  1057. # Given this we do not need to actually parse tested_inter_matrix to compute m.
  1059. # Convert inter_matrix and tested_inter_matrix into numpy arrays.
  1060. dimensions = inter_array.shape
  1061. total_markers1 = dimensions[0]
  1062. total_markers2 = dimensions[1]
  1064. i = i_start
  1065. j = 0
  1066. a = max_inter_length
  1067. b = max_inter_length
  1069. km_results = dict()
  1071. # filename = out_dir + 'km_results_' + str(i_start) + '_' + str(i_end) + '.csv'
  1072. # with open(filename, 'w') as writer:
  1073. # writer.write('#i,j,a,b,k,m')
  1075. while i < i_end:
  1076. #skipped_comp_k = 0
  1077. #computed_k = 0
  1078. while j < total_markers2:
  1079. # If we the biggest matrix at (i,j) has k = 0, then all of the matrices within it
  1080. # will also have k = 0.
  1081. k_is_zero = False
  1082. # a_thresh and b_thresh represent biggest submatrix that is all zeros
  1083. a_thresh = 0
  1084. b_thresh = 0
  1085. while a >= 1:
  1086. while b >= 1:
  1087. if((i + a) <= total_markers1 and (j+b) <= total_markers2):
  1088. if(k_is_zero and a <= a_thresh and b <= b_thresh):
  1089. #skipped_comp_k += 1
  1090. k,m = comp_alt(i,j,a,b, inter_array, upper_triangular, compute_k=False)
  1091. else:
  1092. k,m = comp_alt(i,j,a,b, inter_array, upper_triangular)
  1093. #computed_k += 1
  1094. if(k == 0):
  1095. k_is_zero = True
  1096. a_thresh = a
  1097. b_thresh = b
  1098. # Do not store the k,m pair if k <= expected value (m*n/N)
  1099. if(k > (m*n/N)):
  1100. km_results[(i,j,a,b)] = (k,m)
  1101. else:
  1102. pass
  1104. # writer.write(str(i)+','+str(j)+','+str(a)+','+str(b)+','+str(k)+','+str(m))
  1106. b -= 1
  1107. b = max_inter_length
  1108. a -= 1
  1109. a = max_inter_length
  1110. b = max_inter_length
  1111. j += 1
  1112. #print('j = ', j)
  1113. #a = max_inter_length
  1114. #b = max_inter_length
  1115. j = 0
  1116. i += 1
  1117. print('i = ', i)
  1118. #print("Number of skipped k computations: ", skipped_comp_k)
  1119. #print("Number of completed k computations: ", computed_k)
  1121. print("Length of km_results: ", len(km_results))
  1122. #now = datetime.now()
  1123. #current_time = now.strftime("%H:%M:%S")
  1124. #print("Current Time =", current_time)
  1125. return km_results
  1127. #@jit(nopython=True)
  1128. def compute_interval_pval_parallel(km_results:dict, N:int, n:int, i_start:int, i_end:int,
  1129. max_inter_length:int, out_dir:str, chrom1:str = '-1',
  1130. chrom2:str = '-1'):
  1131. # Assume that all k values are bigger than expected k-value of (m*n/N)
  1133. computed_pvals = dict()
  1134. # Conservative Bonferonni Correction
  1135. # In practice the number of tests is much smaller
  1136. correction = 0.05 / (N * (max_inter_length**2))
  1137. # this would be correct if we only computed
  1138. # one chromosome pair, but unfortunately we need to account for all pairs
  1140. # hard coded correction:
  1141. #correction = 0.05 / (60055879978 * 3600)
  1144. if(chrom1 == '-1' or chrom2 == '-1'):
  1145. filename = out_dir + 'pval_results_' + str(i_start) + '_' + str(i_end) + '.csv'
  1146. else:
  1147. filename = out_dir + 'chr' + chrom1 + '_chr' + chrom2 + '_pval_results_'
  1148. filename += str(i_start) + '_' + str(i_end) + '.csv'
  1149. with open(filename, 'w') as writer:
  1151. pval_results = dict()
  1152. # i = 0
  1153. for key,val in km_results.items():
  1154. k,m = val[0], val[1]
  1155. max_k = min(m, n)
  1156. # min_k = max(0, m-(N-n))
  1157. pval_results[key] = 0
  1159. # First check if a p-value for the following k,m pair has already been computed.
  1160. try:
  1161. pval_results[key] = computed_pvals[(val[0],m)]
  1162. except KeyError:
  1163. # If it has not been, then compute the p-value and store it.
  1164. # if(k > (m*n/N)):
  1166. pval_results[key] = hypergeom.sf(k-1,N,n,m)
  1167. # while k <= max_k:
  1168. # pval_results[key] += hypergeom.pmf(k,N,n,m)
  1169. # k += 1
  1171. if(pval_results[key] < correction):
  1172. #print(pval_results[key])
  1173. writer.write(str(key[0])+','+str(key[1])+','+str(key[2])+','+str(key[3])+',')
  1174. writer.write(str(val[0]) + ',' + str(m) + ',' + str(pval_results[key]) + '\n')
  1175. # We do not care about interval pairs that have unusually low
  1176. # number of interacting pairs
  1177. #else:
  1178. # while k >= min_k:
  1179. # pval_results[key] += hypergeom.pmf(k, N, n, m)
  1180. # k -= 1
  1182. computed_pvals[(val[0],m)] = pval_results[key]
  1184. # i += 1
  1185. # if(i%1000 == 0):
  1186. # print(i)
  1187. # print(len(computed_pvals))
  1188. print("Computed p-values.")
  1189. # now = datetime.now()
  1190. # current_time = now.strftime("%H:%M:%S")
  1191. # print("Current Time =", current_time)
  1192. return pval_results
  1194. def trim_intervals(sorted_intervals:list):
  1195. # Now go through interval pairs starting from most significant, and remove all overlapping
  1196. # intervals for the current interval
  1197. outer_count = 0
  1198. indeces_to_delete = set()
  1199. for tup in sorted_intervals:
  1200. # Skip this interval pair if it is already marked for delition
  1201. if(outer_count in indeces_to_delete):
  1202. outer_count += 1
  1203. #print('SKIPPED!')
  1204. continue
  1205. inner_count = 0
  1206. curr_tup = tup
  1207. i,j,a,b = curr_tup[0][0],curr_tup[0][1],curr_tup[0][2],curr_tup[0][3]
  1208. pval_outer = curr_tup[1]
  1209. for tup in sorted_intervals:
  1210. if(inner_count in indeces_to_delete):
  1211. inner_count += 1
  1212. #print('SKIPPED!')
  1213. continue
  1214. i2,j2,a2,b2 = tup[0][0],tup[0][1],tup[0][2],tup[0][3]
  1215. pval_inner = tup[1]
  1216. #print(i,j,a,b)
  1217. #print(i2,j2,a2,b2)
  1218. if(check_overlap(i,j,a,b,i2,j2,a2,b2)):
  1219. if(pval_inner < pval_outer):
  1220. indeces_to_delete.add(outer_count)
  1221. else:
  1222. indeces_to_delete.add(inner_count)
  1223. #print("bye!")
  1225. inner_count += 1
  1227. outer_count += 1
  1229. #print(indeces_to_delete)
  1231. culled_inter_pairs = []
  1232. i = 0
  1233. while i < len(sorted_intervals):
  1234. if(i not in indeces_to_delete):
  1235. culled_inter_pairs.append(sorted_intervals[i])
  1236. i += 1
  1238. return culled_inter_pairs
  1240. class TestBiClusterCodeBase(unittest.TestCase):
  1242. def test_get_msigdb_enrichment(self):
  1243. gene_pairs = {('ABAT', 'ACOX1'), ('CSMD1', 'DLGAP1')}
  1244. out_dir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Test/'
  1245. get_msigdb_enrichment(gene_pairs, out_dir)
  1247. fname = out_dir + 'biclustering_results_msigdb_enrichment.txt'
  1248. with open(fname) as msig_fname:
  1249. line = msig_fname.readline()
  1251. exp_line = "('ABAT', 'ACOX1') : ['AFFAR_YY1_TARGETS_UP', "
  1252. exp_line += "'CARRILLOREIXACH_HEPATOBLASTOMA_VS_NORMAL_DN', 'FLECHNER_BIOPSY_KIDNEY_"
  1253. exp_line += "TRANSPLANT_REJECTED_VS_OK_DN', 'LEE_LIVER_CANCER_MYC_TGFA_DN', 'RODRIGUES_"
  1254. exp_line += "DCC_TARGETS_DN', 'RODRIGUES_THYROID_CARCINOMA_ANAPLASTIC_DN', "
  1255. exp_line += "'RODRIGUES_THYROID_CARCINOMA_POORLY_DIFFERENTIATED_DN']\n"
  1256. self.assertEqual(line, exp_line)
  1258. def test_parse_remma_interaction_data(self):
  1259. indir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Test/'
  1260. ad_anno = indir + 'epiAD_ad_NA3_Combined_Strict_Set_Brief.anno'
  1261. dd_anno = indir + 'epiDD_dd_NA3_Combined_Strict_Set_Brief.anno'
  1262. interactions = parse_remma_interaction_data([ad_anno, dd_anno], 1e-6, '1', '20')
  1264. self.assertEqual(len(interactions), 7)
  1265. self.assertEqual(interactions[('1:7863293', '20:8515243')], 3.0037315721644386e-07)
  1266. self.assertEqual(interactions[('1:7863293', '20:8515956')], 3.0037315721644386e-07)
  1268. self.assertEqual(interactions[('1:7865063', '20:33467717')], 8.801979806253102e-07)
  1269. self.assertEqual(interactions[('1:7865063', '20:33488013')], 7.014568812762179e-07)
  1270. self.assertEqual(interactions[('1:7865063', '20:33514465')], 7.070270435022684e-07)
  1271. self.assertEqual(interactions[('1:7865691', '20:33488013')], 9.85547891823227e-09)
  1272. self.assertEqual(interactions[('1:1725760', '20:9588420')], 7.018367391361476e-07)
  1275. def test_get_number_of_interacting_snps_in_interval(self):
  1276. indir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Test/'
  1277. biclustering_file = indir + 'biclustering_results_chr1_chr5.txt'
  1278. bim_file = indir + 'NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.bim'
  1279. gene_location_file = indir + 'Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_ranges.txt'
  1280. marker_files = [indir + 'epiDD_dd_NA3_Combined_Strict_Set.anno']
  1281. inter_pairs_per_interval = get_number_of_interacting_snps_in_interval(biclustering_file, bim_file, gene_location_file,
  1282. '1', '5', marker_files, 1e-7)
  1284. self.assertEqual(inter_pairs_per_interval[('LHX4 ', 'PPP2R2B ')], 12)
  1285. self.assertEqual(inter_pairs_per_interval[('DISC1 ', 'SLIT3 ')], 5)
  1287. def test_parse_biclustering_results(self):
  1288. indir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/REMMA_Results/'
  1289. marker_file1 = 'epiAA_aa1_approx_parallel_merged_NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.anno'
  1290. marker_file2 = 'epiAA_aa2_approx_parallel_merged_NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.anno'
  1291. marker_file3 = 'epiAA_aa3_approx_parallel_merged_NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.anno'
  1292. marker_file4 = 'epiAD_ad_approx_parallel_merged_NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.anno'
  1293. marker_file5 = 'epiDD_dd_approx_parallel_merged_NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.anno'
  1295. marker_files = [indir + marker_file1, indir + marker_file2, indir + marker_file3, indir + marker_file4, indir + marker_file5]
  1297. indir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Test/'
  1298. biclustering_file = indir + 'biclustering_results_chr1_chr3.txt'
  1299. bim_file = indir + 'NA3_Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_maf001_region_qc3.bim'
  1300. gene_location_file = indir + 'Combined_Strict_Set_score900_db001_or_exp001_1_and_genehancer_score10_ranges.txt'
  1301. out_dir = indir
  1303. parse_biclustering_results(biclustering_file, bim_file, gene_location_file, '1', '3', marker_files, 1e-7, out_dir, False)
  1304. filename = indir + 'biclustering_results_chr1_chr3_annotated.txt'
  1305. rows = []
  1306. with open(filename) as csv_file:
  1307. csv_reader = csv.reader(csv_file, delimiter='\t')
  1308. for row in csv_reader:
  1309. rows.append(row)
  1311. self.assertEqual(rows[0], ["{('DTL ', 'SRGAP3 ')}", '9.605131888649506e-62'])
  1312. self.assertEqual(rows[1], ["{('DTL ', 'SRGAP3 ')}", '6.4212833992207265e-27'])
  1313. self.assertEqual(rows[2], ["{('SRGAP2 ', 'ROBO2 ')}", '4.0199636843080785e-19'])
  1315. parse_biclustering_results(biclustering_file, bim_file, gene_location_file, '1', '3', marker_files, 1e-7, out_dir, True)
  1317. biclustering_file = indir + 'biclustering_results_chr4_chr10.txt'
  1318. parse_biclustering_results(biclustering_file, bim_file, gene_location_file, '4', '10', marker_files, 1e-7, out_dir, True)
  1319. biclustering_file = indir + 'biclustering_results_chr4_chr5.txt'
  1320. parse_biclustering_results(biclustering_file, bim_file, gene_location_file, '4', '5', marker_files, 1e-7, out_dir, True)
  1321. biclustering_file = indir + 'biclustering_results_chr8_chr17.txt'
  1322. parse_biclustering_results(biclustering_file, bim_file, gene_location_file, '8', '17', marker_files, 1e-7, out_dir, True)
  1324. filename = indir + 'biclustering_results_all.txt'
  1325. rows = []
  1326. with open(filename) as csv_file:
  1327. csv_reader = csv.reader(csv_file, delimiter='\t')
  1328. for row in csv_reader:
  1329. rows.append(row)
  1331. self.assertEqual(rows[0], ['chr1,212206918,212280187,chr3,9020275,9293369,9.605131888649506e-62'])
  1332. self.assertEqual(rows[1], ['chr1,206514199,206639783,chr3,75984644,77701114,4.0199636843080785e-19'])
  1333. self.assertEqual(rows[2], ['chr4,6320304,6567327,chr10,52748910,54060110,2.993577604860245e-295'])
  1334. self.assertEqual(rows[3], ['chr4,88080213,88143674,chr5,145967066,146463083,5.967210373347171e-208'])
  1335. self.assertEqual(rows[4], ['chr4,20253234,20622788,chr5,15498304,15941900,4.7347537763144387e-20'])
  1337. self.assertTrue(['chr8,8957214,8964665,chr17,28519336,28564986,1.1381827143817773e-19'] in rows)
  1338. self.assertTrue(['chr8,8957214,8964665,chr17,28443018,28446019,1.1381827143817773e-19'] in rows)
  1340. os.remove(indir + 'biclustering_results_all.txt')
  1342. def test_generate_gene_blocks(self):
  1343. fname = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Biclustering/Test_Input/'
  1344. fname += 'Test_Ranges.txt'
  1345. blocks = generate_gene_blocks(fname, 3000)
  1346. expected_blocks = [{'TEST3'}, {'TEST4'}, {'TEST5'}, {'TEST6'}, {'TEST11'}, {'TEST12'},
  1347. {'TEST22'}, {'TEST2', 'TEST1'}, {'TEST9', 'TEST10', 'TEST8', 'TEST7'},
  1348. {'TEST13', 'TEST14', 'TEST16', 'TEST18'},
  1349. {'TEST15', 'TEST19', 'TEST21', 'TEST20', 'TEST17'}]
  1351. self.assertTrue(len(blocks) == len(expected_blocks))
  1352. for block in expected_blocks:
  1353. self.assertTrue(block in blocks)
  1355. blocks = generate_gene_blocks(fname, 200)
  1356. expected_blocks = [{'TEST1'}, {'TEST2'}, {'TEST3'}, {'TEST4'}, {'TEST5'},
  1357. {'TEST6'}, {'TEST8'}, {'TEST9'}, {'TEST11'}, {'TEST12'},
  1358. {'TEST16'}, {'TEST18'}, {'TEST20'}, {'TEST21'}, {'TEST22'},
  1359. {'TEST10', 'TEST7'}, {'TEST13', 'TEST14'},
  1360. {'TEST19', 'TEST15', 'TEST17'}]
  1362. self.assertTrue(len(blocks) == len(expected_blocks))
  1363. for block in expected_blocks:
  1364. self.assertTrue(block in blocks)
  1366. blocks = generate_gene_blocks(fname, 1000)
  1367. expected_blocks = [{'TEST1'}, {'TEST2'}, {'TEST3'}, {'TEST4'}, {'TEST5'}, {'TEST6'},
  1368. {'TEST11'}, {'TEST12'}, {'TEST22'},
  1369. {'TEST9', 'TEST10', 'TEST8', 'TEST7'},
  1370. {'TEST13', 'TEST14', 'TEST18', 'TEST16'},
  1371. {'TEST15', 'TEST19', 'TEST21', 'TEST20', 'TEST17'}]
  1373. self.assertTrue(len(blocks) == len(expected_blocks))
  1374. for block in expected_blocks:
  1375. self.assertTrue(block in blocks)
  1378. def test_check_overlap(self):
  1379. # Simple Cases
  1380. self.assertTrue(check_overlap(166,621,10,40,167,620,9,42))
  1381. self.assertTrue(check_overlap(0,491,36,11,0,491,35,11))
  1382. self.assertTrue(check_overlap(10,10,5,5,5,12,8,5))
  1383. self.assertTrue(check_overlap(0,0,2,2,1,1,2,2))
  1384. self.assertTrue(check_overlap(2,2,3,6,1,1,3,3))
  1385. self.assertFalse(check_overlap(0,0,2,2,3,2,1,1))
  1386. self.assertFalse(check_overlap(0,0,10,10,40,40,2,2))
  1388. # Challenging/Border Cases
  1389. self.assertFalse(check_overlap(0,0,2,2,0,0,2,2))
  1390. self.assertFalse(check_overlap(0,0,2,2,1,2,1,1))
  1391. self.assertFalse(check_overlap(0,0,2,2,2,2,1,1))
  1392. self.assertFalse(check_overlap(2,2,3,6,1,1,1,3))
  1394. def test_initialize_matrices(self):
  1395. total_markers = 10
  1396. fname = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Biclustering/Test_Input/TenByTen.txt'
  1397. N, n, inter_matrix, upper_tri = initialize_matrices(total_markers, fname)
  1399. inter_matrix_exp = dict()
  1400. i, j = 0, 0
  1401. while i < 10:
  1402. j = 0
  1403. while j < 10:
  1404. inter_matrix_exp[(i,j)] = 0
  1405. j += 1
  1406. i += 1
  1408. inter_matrix_exp[(2,8)] = 1
  1409. inter_matrix_exp[(5,8)] = 1
  1410. inter_matrix_exp[(5,9)] = 1
  1412. self.assertEqual(N, 45)
  1413. self.assertEqual(n, 3)
  1414. self.assertTrue(upper_tri)
  1416. i, j = 0, 0
  1417. while i < 10:
  1418. j = 0
  1419. while j < 10:
  1420. self.assertEqual(inter_matrix[(i,j)], inter_matrix_exp[(i,j)])
  1421. j += 1
  1422. i += 1
  1424. N, n, inter_matrix, upper_tri = initialize_matrices(total_markers, fname, False)
  1426. self.assertEqual(N, 100)
  1427. self.assertEqual(n, 3)
  1428. self.assertFalse(upper_tri)
  1430. i, j = 0, 0
  1431. while i < 10:
  1432. j = 0
  1433. while j < 10:
  1434. self.assertEqual(inter_matrix[(i,j)], inter_matrix_exp[(i,j)])
  1435. j += 1
  1436. i += 1
  1438. def test_initialize_matrices2(self):
  1439. indir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Biclustering/Test_Input/'
  1440. bim_file = indir + 'Test_Data.bim'
  1441. inter_file = indir + 'Test_Interaction.anno'
  1442. chrom1 = '1'
  1443. chrom2 = '1'
  1445. N, n, inter_array, upper_tri = initialize_matrices2(bim_file, [inter_file], chrom1, chrom2, 1e-5)
  1446. self.assertEqual(N, 435)
  1447. self.assertEqual(n, 1)
  1448. self.assertTrue(upper_tri)
  1450. expected_inter_array = np.zeros((inter_array.shape[0],inter_array.shape[1]))
  1451. expected_inter_array[12][19] = 1
  1453. i = 0
  1454. while i < inter_array.shape[0]:
  1455. j = 0
  1456. while j < inter_array.shape[1]:
  1457. self.assertEqual(inter_array[i][j], expected_inter_array[i][j])
  1458. j += 1
  1459. i += 1
  1461. indir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Biclustering/Test_Input/'
  1462. bim_file = indir + 'Test_Data.bim'
  1463. inter_file = indir + 'Test_Interaction.anno'
  1464. chrom1 = '2'
  1465. chrom2 = '3'
  1467. N, n, inter_array, upper_tri = initialize_matrices2(bim_file, [inter_file], chrom1, chrom2, 1e-5)
  1469. self.assertEqual(N, 900)
  1470. self.assertEqual(n, 5)
  1471. self.assertFalse(upper_tri)
  1473. expected_inter_array = np.zeros((inter_array.shape[0],inter_array.shape[1]))
  1474. for k in [(1,0), (8,11), (10,11), (12,14), (14,14)]:
  1475. expected_inter_array[k[0]][k[1]] = 1
  1477. i = 0
  1478. while i < inter_array.shape[0]:
  1479. j = 0
  1480. while j < inter_array.shape[1]:
  1481. self.assertEqual(inter_array[i][j], expected_inter_array[i][j])
  1482. j += 1
  1483. i += 1
  1485. N, n, inter_array, upper_tri = initialize_matrices2(bim_file, [inter_file], chrom1, chrom2, 1e-6)
  1487. self.assertEqual(N, 900)
  1488. self.assertEqual(n, 1)
  1489. self.assertFalse(upper_tri)
  1491. expected_inter_array = np.zeros((inter_array.shape[0],inter_array.shape[1]))
  1492. for k in [(1,0)]:
  1493. expected_inter_array[k[0]][k[1]] = 1
  1495. i = 0
  1496. while i < inter_array.shape[0]:
  1497. j = 0
  1498. while j < inter_array.shape[1]:
  1499. self.assertEqual(inter_array[i][j], expected_inter_array[i][j])
  1500. j += 1
  1501. i += 1
  1503. chrom1 = '3'
  1504. chrom2 = '2'
  1505. N, n, inter_array, upper_tri = initialize_matrices2(bim_file, [inter_file], chrom1, chrom2, 1e-5)
  1507. self.assertEqual(N, 900)
  1508. self.assertEqual(n, 5)
  1509. self.assertFalse(upper_tri)
  1511. expected_inter_array = np.zeros((inter_array.shape[0],inter_array.shape[1]))
  1512. for k in [(0,1), (11,8), (11,10), (14,12), (14,14)]:
  1513. expected_inter_array[k[0]][k[1]] = 1
  1515. i = 0
  1516. while i < inter_array.shape[0]:
  1517. j = 0
  1518. while j < inter_array.shape[1]:
  1519. self.assertEqual(inter_array[i][j], expected_inter_array[i][j])
  1520. j += 1
  1521. i += 1
  1523. chrom1 = '2'
  1524. chrom2 = '3'
  1525. inter_file2 = indir + 'Test_Interaction2.anno'
  1526. N, n, inter_array, upper_tri = initialize_matrices2(bim_file, [inter_file, inter_file2], chrom1, chrom2, 1e-5)
  1527. self.assertEqual(N, 900)
  1528. self.assertEqual(n, 6)
  1529. self.assertFalse(upper_tri)
  1531. expected_inter_array = np.zeros((inter_array.shape[0],inter_array.shape[1]))
  1532. for k in [(1,0), (8,11), (10,11), (12,14), (14,14), (27, 25)]:
  1533. expected_inter_array[k[0]][k[1]] = 1
  1535. i = 0
  1536. while i < inter_array.shape[0]:
  1537. j = 0
  1538. while j < inter_array.shape[1]:
  1539. self.assertEqual(inter_array[i][j], expected_inter_array[i][j])
  1540. j += 1
  1541. i += 1
  1544. N, n, inter_array, upper_tri = initialize_matrices2(bim_file, [inter_file, inter_file2], chrom1, chrom2, 1e-8)
  1545. self.assertEqual(N, 900)
  1546. self.assertEqual(n, 2)
  1547. self.assertFalse(upper_tri)
  1549. expected_inter_array = np.zeros((inter_array.shape[0],inter_array.shape[1]))
  1550. for k in [(8,11), (27, 25)]:
  1551. expected_inter_array[k[0]][k[1]] = 1
  1553. i = 0
  1554. while i < inter_array.shape[0]:
  1555. j = 0
  1556. while j < inter_array.shape[1]:
  1557. self.assertEqual(inter_array[i][j], expected_inter_array[i][j])
  1558. j += 1
  1559. i += 1
  1562. def test_comp_alt(self):
  1563. inter_matrix = dict()
  1564. i, j = 0, 0
  1565. while i < 10:
  1566. j = 0
  1567. while j < 10:
  1568. inter_matrix[(i,j)] = 0
  1569. j += 1
  1570. i += 1
  1572. inter_matrix[(2,8)] = 1
  1573. inter_matrix[(5,8)] = 1
  1574. inter_matrix[(5,9)] = 1
  1576. inter_array = convert_dict_to_2D_numpy_array(inter_matrix)
  1578. #print(inter_array)
  1580. self.assertEqual(comp_alt(0,0,2,2, inter_array, upper_triangular=True), (0,1))
  1581. self.assertEqual(comp_alt(0,0,2,2, inter_array, upper_triangular=False), (0,4))
  1583. self.assertEqual(comp_alt(0,0,2,4, inter_array, upper_triangular=True), (0,5))
  1584. self.assertEqual(comp_alt(0,0,2,4, inter_array, upper_triangular=False), (0,8))
  1586. self.assertEqual(comp_alt(0,0,3,10, inter_array, upper_triangular=True), (1,24))
  1587. self.assertEqual(comp_alt(0,0,3,10, inter_array, upper_triangular=False), (1,30))
  1589. self.assertEqual(comp_alt(5,0,4,4, inter_array, upper_triangular=True), (0,0))
  1590. self.assertEqual(comp_alt(5,3,4,4, inter_array, upper_triangular=True), (0,1))
  1591. self.assertEqual(comp_alt(5,6,2,4, inter_array, upper_triangular=True), (2,7))
  1592. self.assertEqual(comp_alt(5,8,1,1, inter_array, upper_triangular=True), (1,1))
  1593. self.assertEqual(comp_alt(5,8,1,2, inter_array, upper_triangular=True), (2,2))
  1594. self.assertEqual(comp_alt(5,8,2,1, inter_array, upper_triangular=True), (1,2))
  1596. # Performance testing
  1597. '''
  1598. inter_matrix = dict()
  1600. i, j = 0, 0
  1601. while i < 2000:
  1602. j = 0
  1603. while j < 2000:
  1604. inter_matrix[(i,j)] = 0
  1605. j += 1
  1606. i += 1
  1608. inter_matrix[(2,8)] = 1
  1609. inter_matrix[(5,8)] = 1
  1610. inter_matrix[(5,9)] = 1
  1612. inter_array = convert_dict_to_2D_numpy_array(inter_matrix)
  1614. start = time.time()
  1615. self.assertEqual(comp_alt(0,0,1850,1900, inter_array, upper_triangular = False), (3, 3515000))
  1616. print(time.time()-start)
  1617. '''
  1619. def test_compute_k_m_parallel(self):
  1621. inter_array = np.zeros((10,10))
  1622. for k in [(2,8), (5,8), (5,9)]:
  1623. inter_array[k[0]][k[1]] = 1
  1625. # Note, technically, the last parameter should be 3, but we enter 0
  1626. # to force all k,m pairs to be written.
  1627. km_results = compute_k_m_parallel(4, inter_array, True, 0, 9, 45, 0)
  1629. self.assertEqual(km_results[(0,5,4,4)], (1,16))
  1630. self.assertEqual(km_results[(0,8,3,1)], (1,3))
  1631. self.assertEqual(km_results[(1,7,2,2)], (1,4))
  1632. self.assertEqual(km_results[(2,5,4,4)], (2,15))
  1633. self.assertEqual(km_results[(2,8,4,2)], (3,8))
  1634. self.assertEqual(km_results[(4,7,3,3)], (2,9))
  1636. self.assertEqual(km_results[(5,5,4,4)], (1,6))
  1637. self.assertEqual(km_results[(5,5,3,4)], (1,6))
  1638. self.assertEqual(km_results[(5,5,2,4)], (1,5))
  1639. self.assertEqual(km_results[(5,5,1,4)], (1,3))
  1640. self.assertEqual(km_results[(5,6,4,4)], (2,10))
  1641. self.assertEqual(km_results[(5,6,4,3)], (1,6))
  1642. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1643. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1644. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1645. self.assertEqual(km_results[(5,6,3,3)], (1,6))
  1646. self.assertEqual(km_results[(5,6,2,4)], (2,7))
  1648. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1649. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1650. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1651. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1652. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1653. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1654. self.assertEqual(km_results[(5,6,3,4)], (2,9))
  1656. self.assertEqual(km_results[(5,6,2,4)], (2,7))
  1657. self.assertEqual(km_results[(5,8,1,1)], (1,1))
  1658. self.assertEqual(km_results[(5,8,1,2)], (2,2))
  1659. self.assertEqual(km_results[(5,8,2,1)], (1,2))
  1661. inter_array = np.zeros((10,10))
  1662. for k in [(2,8)]:
  1663. inter_array[k[0]][k[1]] = 1
  1665. # Note, technically, the last parameter should be 1, but we enter 0
  1666. # to force all k,m pairs to be written.
  1667. km_results = compute_k_m_parallel(4, inter_array, True, 0, 9, 45, 0)
  1668. self.assertEqual(len(km_results), 63)
  1670. def test_compute_interval_pval_parallel(self):
  1671. inter_array = np.zeros((10,10))
  1672. for k in [(2,8), (5,8), (5,9)]:
  1673. inter_array[k[0]][k[1]] = 1
  1675. # Note, technically, the last parameter should be 3, but we enter 0
  1676. # to force all k,m pairs to be written.
  1677. km_results = compute_k_m_parallel(4, inter_array, True, 0, 9, 45, 0)
  1678. out_dir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Biclustering/Test_Output/'
  1679. pval_results = compute_interval_pval_parallel(km_results, 45, 3, 0, 9, 4, out_dir)
  1680. #print(pval_results)
  1682. self.assertAlmostEqual(pval_results[(0,5,4,4)], 0.7424947145877379)
  1683. self.assertAlmostEqual(pval_results[(0,8,3,1)], 0.19097956307258637)
  1684. self.assertAlmostEqual(pval_results[(1,7,2,2)], 0.24876673713883019)
  1685. self.assertAlmostEqual(pval_results[(2,5,4,4)], 0.25405214940098664)
  1686. self.assertAlmostEqual(pval_results[(2,8,4,2)], 0.0039464411557434895)
  1687. self.assertAlmostEqual(pval_results[(4,7,3,3)], 0.09725158562367894)
  1689. self.assertAlmostEqual(pval_results[(5,5,4,4)], 0.35595489781536405)
  1690. self.assertAlmostEqual(pval_results[(5,5,3,4)], 0.35595489781536405)
  1691. self.assertAlmostEqual(pval_results[(5,5,2,4)], 0.30373502466525815)
  1692. self.assertAlmostEqual(pval_results[(5,5,1,4)], 0.19097956307258632)
  1693. self.assertAlmostEqual(pval_results[(5,6,4,4)], 0.11945031712473582)
  1694. self.assertAlmostEqual(pval_results[(5,6,4,3)], 0.35595489781536405)
  1695. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1696. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1697. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1698. self.assertAlmostEqual(pval_results[(5,6,3,3)], 0.35595489781536405)
  1699. self.assertAlmostEqual(pval_results[(5,6,2,4)], 0.058703312191684544)
  1701. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1702. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1703. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1704. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1705. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1706. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1707. self.assertAlmostEqual(pval_results[(5,6,3,4)], 0.09725158562367894)
  1709. self.assertAlmostEqual(pval_results[(5,6,2,4)], 0.058703312191684544)
  1710. self.assertAlmostEqual(pval_results[(5,8,1,1)], 0.06666666666666664)
  1711. self.assertAlmostEqual(pval_results[(5,8,1,2)], 0.003030303030303028)
  1712. self.assertAlmostEqual(pval_results[(5,8,2,1)], 0.13030303030303028)
  1715. def test_trim_intervals(self):
  1716. in_dir = 'C:/Stas/LabWork/Bioinformatics/Projects/Ch5_AI_AUD/Biclustering/Test_Input/'
  1717. filename = in_dir + 'pval_results_0_1.csv'
  1719. pval_results = dict()
  1720. with open(filename) as csv_file:
  1721. csv_reader = csv.reader(csv_file, delimiter=',')
  1722. for row in csv_reader:
  1723. pval_results[(int(row[0]),int(row[1]),int(row[2]),int(row[3]))] = float(row[6])
  1724. s_inter_pairs = sorted(pval_results.items(), key = lambda x: abs(x[1]), reverse = False)
  1725. #print(len(s_inter_pairs))
  1726. #print(s_inter_pairs)
  1727. trimmed_inter_pairs = trim_intervals(s_inter_pairs)
  1728. #print(trimmed_inter_pairs)
  1730. #print(len(trimmed_inter_pairs))
  1731. #print(trimmed_inter_pairs)
  1732. expected_trimmed_inter_pairs = [((0, 816, 48, 8), 2.5317372178516744e-24),
  1733. ((0, 999, 35, 3), 8.785313396735072e-16),
  1734. ((0, 494, 35, 8), 3.520342278584121e-15)]
  1736. i = 0
  1737. for inter_pair in trimmed_inter_pairs:
  1738. self.assertEqual(inter_pair, expected_trimmed_inter_pairs[i])
  1739. i += 1
  1741. # COMPLETE THIS TEST CASE
  1742. filename = in_dir + 'chr1_chr2_pval_results.csv'
  1743. pval_results = dict()
  1744. with open(filename) as csv_file:
  1745. csv_reader = csv.reader(csv_file, delimiter=',')
  1746. for row in csv_reader:
  1747. pval_results[(int(row[0]),int(row[1]),int(row[2]),int(row[3]))] = float(row[6])
  1748. s_inter_pairs = sorted(pval_results.items(), key = lambda x: abs(x[1]), reverse = False)
  1749. trimmed_inter_pairs = trim_intervals(s_inter_pairs)
  1750. #print(trimmed_inter_pairs)
  1752. expected_trimmed_inter_pairs = [((1, 41, 57, 59), 1.5052159828764328e-29),
  1753. ((65, 53, 35, 47), 5.109419105661788e-22),
  1754. ((45, 1, 55, 33), 1.5051851637750398e-21)]
  1756. i = 0
  1757. for inter_pair in trimmed_inter_pairs:
  1758. self.assertEqual(inter_pair, expected_trimmed_inter_pairs[i])
  1759. i += 1
  1761. def test_convert_dict_to_2D_numpy_array(self):
  1762. tested_inter_matrix = dict()
  1763. i, j = 0, 0
  1764. while i < 10:
  1765. j = 0
  1766. while j < 10:
  1767. if (j > i):
  1768. tested_inter_matrix[(i,j)] = 1
  1769. else:
  1770. tested_inter_matrix[(i,j)] = 0
  1771. j += 1
  1772. i += 1
  1774. array = convert_dict_to_2D_numpy_array(tested_inter_matrix)
  1775. #print(array)
  1776. #for x in np.nditer(array):
  1777. # print(x)
  1779. def test_suite():
  1780. # Unit Tests
  1781. unit_test_suite = unittest.TestSuite()
  1782. unit_test_suite.addTest(TestBiClusterCodeBase('test_initialize_matrices'))
  1783. unit_test_suite.addTest(TestBiClusterCodeBase('test_check_overlap'))
  1784. unit_test_suite.addTest(TestBiClusterCodeBase('test_comp_alt'))
  1785. unit_test_suite.addTest(TestBiClusterCodeBase('test_convert_dict_to_2D_numpy_array'))
  1786. unit_test_suite.addTest(TestBiClusterCodeBase('test_compute_k_m_parallel'))
  1787. unit_test_suite.addTest(TestBiClusterCodeBase('test_compute_interval_pval_parallel'))
  1788. unit_test_suite.addTest(TestBiClusterCodeBase('test_generate_gene_blocks'))
  1789. unit_test_suite.addTest(TestBiClusterCodeBase('test_trim_intervals'))
  1790. unit_test_suite.addTest(TestBiClusterCodeBase('test_initialize_matrices2'))
  1791. unit_test_suite.addTest(TestBiClusterCodeBase('test_parse_biclustering_results'))
  1793. unit_test_suite.addTest(TestBiClusterCodeBase('test_get_number_of_interacting_snps_in_interval'))
  1794. unit_test_suite.addTest(TestBiClusterCodeBase('test_parse_remma_interaction_data'))
  1795. unit_test_suite.addTest(TestBiClusterCodeBase('test_get_msigdb_enrichment'))
  1797. runner = unittest.TextTestRunner()
  1798. runner.run(unit_test_suite)

BiClustering.py at commit effebb8, no license · at the source

Overview

Authors: Stanislav Listopad1, Qian Peng1
ORCID iDs: Stanislav Listopad, Qian Peng
  1. Department of Neuroscience, The Scripps Research Institute, La Jolla, CA USA
Institutions: Scripps Research Institute (United States)
Journal: Communications biology, volume 9, issue 1, article 750
Dates: received 10 March 2025; accepted 21 March 2026; published online 6 April 2026
Type: Research article · Language: English
License: CC BY-NC-ND
Identifiers: DOI 10.1038/s42003-026-09971-7 · PMID 41942605 · PMCID PMC13230857 · OpenAlex W4405907369
Open access: gold, a free copy (OpenAlex)
Preprint: doi:10.1101/2024.12.29.24319759
Status: code verified
Categories: genetics / omics (modality), human (organism), other condition (population)
Methods: Statistics, Smoothing, state filtering, decompositions, Connectivity
Keywords: Diagnostic markers, Addiction, Genome-wide association studies, Quantitative trait
MeSH: Alcoholism*, Epistasis, Genetic*, Genetic Predisposition to Disease*, Female, Genome-Wide Association Study, Humans, Male, Polymorphism, Single Nucleotide (* major topic)
Topic: Genetic Associations and Epidemiology (Genetics, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: NIAAA NIH HHS (T32 AA007456); U.S. Department of Health & Human Services | NIH | National Institute on Alcohol Abuse and Alcoholism (NIAAA) (T32AA007456); U.S. Department of Health &amp; Human Services | NIH | National Institute on Alcohol Abuse and Alcoholism (T32AA007456); U.S. Department of Health & Human Services | NIH | National Institute on Drug Abuse (NIDA) (DP1DA054373); U.S. Department of Health &amp; Human Services | NIH | National Institute on Drug Abuse (DP1DA054373); NIDA NIH HHS (DP1 DA054373)
Citations: cited by 3 papers (Europe PMC); 90 references in the paper

Abstract

The abstract is not reproduced here: the paper's license (CC BY-NC-ND) does not allow it. Read it in the paper, at the publisher or on Europe PMC.

Repositories

Its files are read in the Code ↔ Paper reader above, with 2 matches between paragraphs and lines of code.

staslist/Biclustering_Epistasis

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: effebb849231932f9d91d26d654dac8b09b13c4d, 8 January 2026
Languages: Python (6)
Size: 38 files, 6 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (2 files), Matplotlib (1 file), Numba (1 file), SciPy (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
7 files
At the source: github.com/staslist/Biclustering_Epistasis

njdjyxz/geno-bicluster

License: none: the authors keep all their rights
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: b0fb2f378c52280c258ea09a640f4aa47a86aa7a, 10 October 2025
Languages: Python (6)
Size: 8 files, 6 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, environment (pyproject.toml), tests
Not found: license file, CITATION.cff, continuous integration, documentation
Tools: Numba (4 files), NumPy (4 files), SciPy (4 files), pandas (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
7 files
At the source: github.com/njdjyxz/geno-bicluster

Zenodo 18960797

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 2 files
Software Heritage: not checked
Found in: “Code availability”
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (6 files), Numba (5 files), SciPy (5 files), Matplotlib (1 file), pandas (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
14 files
At the source:

Zenodo 18961359

License: CC-BY-4.0
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Size: 2 files
Software Heritage: not checked
Found in: the references
Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (6 files), Numba (5 files), SciPy (5 files), Matplotlib (1 file), pandas (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers (HTTP 200)
  • 29 September 2026: the link answers (HTTP 200)
14 files
At the source:

Code availability statement

The paper has a code availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

Read it in the paper: doi.org/10.1038/s42003-026-09971-7.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 4 repositories of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 36 scripts, each with its path and the digest of its content;
  • 2 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability statement

The paper has a data availability statement. Its license (CC BY-NC-ND) does not allow reproducing it here; in short, from what the harvester recognized in it:

  • it says that the data are available on request

Read it in the paper: doi.org/10.1038/s42003-026-09971-7.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 2 authors, 4 keywords, 8 MeSH terms, 6 funders, 88 references.

Cite

This paper

Listopad, S., & Peng, Q. (2026). Extensive genetic interactions (epistasis) linked to alcohol use disorder in a high-risk population. Communications biology, 9(1), 750. https://doi.org/10.1038/s42003-026-09971-7

BibTeX

@article{listopad2026extensive,
author = {Listopad, Stanislav and Peng, Qian},
title = {{Extensive genetic interactions (epistasis) linked to alcohol use disorder in a high-risk population}},
journal = {Communications biology},
year = {2026},
month = apr,
volume = {9},
number = {1},
pages = {750},
publisher = {Nature Publishing Group},
issn = {2399-3642},
doi = {10.1038/s42003-026-09971-7},
url = {https://doi.org/10.1038/s42003-026-09971-7},
pmid = {41942605},
pmcid = {PMC13230857}
}

RIS

TY - JOUR
AU - Listopad, Stanislav
AU - Peng, Qian
TI - Extensive genetic interactions (epistasis) linked to alcohol use disorder in a high-risk population
T2 - Communications biology
J2 - Commun Biol
PY - 2026
DA - 2026/04/06
VL - 9
IS - 1
SP - 750
SN - 2399-3642
PB - Nature Publishing Group
DO - 10.1038/s42003-026-09971-7
UR - https://doi.org/10.1038/s42003-026-09971-7
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s42003-026-09971-7",
"type": "article-journal",
"title": "Extensive genetic interactions (epistasis) linked to alcohol use disorder in a high-risk population",
"container-title": "Communications biology",
"author": [
{
"family": "Listopad",
"given": "Stanislav"
},
{
"family": "Peng",
"given": "Qian"
}
],
"container-title-short": "Commun Biol",
"volume": "9",
"issue": "1",
"page": "750",
"DOI": "10.1038/s42003-026-09971-7",
"PMID": "41942605",
"PMCID": "PMC13230857",
"ISSN": "2399-3642",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s42003-026-09971-7",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
6
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1002/gps3.70017 [code]
Shared genetic variants across substance use disorders implicate common neurobiological pathways, a genome-wide mixed methods study.
Journal: General psychiatry
In common: pandas, SciPy, Matplotlib, 1 other tool, genetics / omics, other condition, 3 references
[2] doi:10.1038/s43856-026-01510-z [code]
Mapping genetic convergence across brain structure, mental health, and cardiometabolic disease.
Journal: Communications medicine
In common: Numba, pandas, SciPy, 2 other tools, genetics / omics, other condition, 2 references
[3] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: Numba, pandas, SciPy, 2 other tools, genetics / omics, other condition, 1 reference
[4] doi:10.1038/s41467-026-71803-3 [code]
Charting the transition from in vitro gliogenesis to the in vivo maturation of human glial progenitor cells transplanted into the hypomyelinated mouse brain.
Journal: Nature communications
In common: Numba, pandas, SciPy, 2 other tools, genetics / omics, 1 reference
[5] doi:10.1038/s41467-026-76676-0 [code]
Determinants of functional burden pleiotropy and gene dosage responses across human traits.
Journal: Nature communications
In common: Numba, pandas, SciPy, 2 other tools, genetics / omics, 1 reference
[6] doi:10.1038/s41398-026-04137-9 [code]
An integrative mendelian randomisation and drug mechanism framework for target prioritisation and therapeutic repurposing in major depression.
Journal: Translational psychiatry
In common: Numba, pandas, SciPy, 2 other tools, genetics / omics, 1 reference
[7] doi:10.1038/s41514-026-00391-9 [code]
Region-specific transcriptional signatures of brain aging in the absence of neuropathology at the single-cell level.
Journal: npj aging
In common: Numba, pandas, SciPy, 2 other tools, genetics / omics, 1 reference
[8] doi:10.1038/s41467-026-73996-z [code]
Genetic architecture of white matter microstructure captured by unsupervised deep representation learning of fractional anisotropy maps.
Journal: Nature communications
In common: pandas, SciPy, Matplotlib, 1 other tool, genetics / omics, 2 references
[9] doi:10.1038/s41467-026-71682-8 [code]
GWAS meta-analysis of cerebrospinal fluid Alzheimer's biomarkers reveals loci regulating lipids, brain volume and autophagy.
Journal: Nature communications
In common: pandas, SciPy, NumPy, genetics / omics, 3 references
[10] doi:10.1016/j.xgen.2026.101284 [code]
NERINE reveals rare variant associations in gene networks across phenotypes and implicates an SNCA-PRL-LRRK2 subnetwork in Parkinson's disease.
Journal: Cell genomics
In common: pandas, SciPy, Matplotlib, 1 other tool, other condition, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.