De novo structural variants in autism spectrum disorder disrupt distal regulatory interactions of neuronal genes.
The 7 matches
- [1] § Methods › Variant scoring with SuPreMo-Akita ↔ scripts/SuPreMo.py, lines 312–387 · score 0.79 · reverse complement, map corresponds, contact frequency maps, alternate allele, Akita scores, augmented
- [2] § Methods › Variant exclusions ↔ scripts/SuPreMo.py, lines 134–188 · score 0.76 · generate mutated sequences, prediction models, genome folding, SuPreMo, pipeline, kb
- [3] § Results › ASD dnSVs disrupt neighboring neuronal regulatory element contacts ↔ walkthroughs/SuPreMo-Akita_weighted_scores.ipynb, lines 155–193 · score 0.63 · ROI bins, weight track, weighted scoring, disruption track, maps
- [4] § Results › ASD dnSVs disrupt neighboring neuronal regulatory element contacts ↔ scripts/get_Akita_scores_utils.py, lines 305–337 · score 0.63 · ROI bins, weight track, weighted scoring, disruption track, maps
- [5] § Methods › Akita model ↔ scripts/SuPreMo.py, lines 6–64 · score 0.60 · H1hESC, GM12878, HCT116, HFF, Spearman, MSE
- [6] § Methods › CREint weighted scoring ↔ walkthroughs/SuPreMo-Akita_weighted_scores.ipynb, lines 155–193 · score 0.60 · weight track, weighted score, disruption track, BED, ROI, bin
- [7] § Methods › Criteria for variant prioritization ↔ scripts/SuPreMo.py, lines 312–387 · score 0.53 · reverse complement, predicted contacts, augmenting, shifting, sequence, Akita
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 · 876 lines · 40 KB · no license · 4 matches
- #!/usr/bin/env python
- # SuPreMo: Pipeline for generating mutated sequences for input into predictive models and for scoring variants for disruption to genome folding.
- # Written in Python v 3.7.11
- '''
- usage: SuPreMo [-h] [--sequences SEQUENCES] [--fa FASTA] [--roi ROI] [--roi_scales ROI_SCALES [ROI_SCALES ...]] [--genome {hg19,hg38}]
- [--Akita_cell_types {HFF,H1hESC,GM12878,IMR90,HCT116} [{HFF,H1hESC,GM12878,IMR90,HCT116} ...]] [--scores {mse,corr,ssi,scc,ins,di,dec,tri,pca} [{mse,corr,ssi,scc,ins,di,dec,tri,pca} ...]]
- [--shift_by SHIFT_WINDOW [SHIFT_WINDOW ...]] [--shifts_file SHIFTS_FILE] [--file OUT_FILE] [--dir OUT_DIR] [--limit SVLEN_LIMIT] [--seq_len SEQ_LEN]
- [--revcomp {no_revcomp,add_revcomp,only_revcomp}] [--augment] [--get_seq] [--get_tracks] [--get_maps] [--get_Akita_scores] [--nrows NROWS]
- Input file
- Pipeline for generating mutated sequences for input into predictive models and for scoring variants for disruption to genome folding.
- positional arguments:
- Input file Input file with variants. Accepted formats are:
- VCF file,
- TSV file (SVannot output),
- BED file,
- TXT file.
- Can be gzipped. Coordinates should be 1-based and left open (,], except for SNPs (follow vcf 4.1/4.2 specifications). See custom_perturbations.ipynb for BED or TXT file input specifications.
- options:
- -h, --help show this help message and exit
- --sequences SEQUENCES
- Input fasta file with sequences. This file is one outputted by SuPreMo.
- (default: None)
- --fa FASTA Path to reference genome fasta file. Default: data/{genome}.fa, where genome is hg38 or hg19.
- (default: None)
- --roi ROI Path to text file with regions of intererst (roi) to use for scoring. Tab-delimited text file with required columns: chrom, start, end or chrom, start1, end1, start2, end2 for paired regions. Alternatively, use 'genes' to weight 2 kb window around protein-coding gene transcription start sites.
- (default: None)
- --roi_scales ROI_SCALES [ROI_SCALES ...]
- Scale by which regions of interest (roi) are upweighted compared to the rest of the map. Default is 10 meaning regions of interest have weights of 10 whereas the rest of the bins have weights of 1.
- (default: [10])
- --genome {hg19,hg38} Genome to be used: hg19 or hg38. (default: ['hg38'])
- --Akita_cell_types {HFF,H1hESC,GM12878,IMR90,HCT116} [{HFF,H1hESC,GM12878,IMR90,HCT116} ...]
- Cell types to make Akita predictions for: HFF, H1hESC, GM12878, IMR90, and/or HCT116.
- (default: ['HFF'])
- --scores {mse,corr,ssi,scc,ins,di,dec,tri,pca} [{mse,corr,ssi,scc,ins,di,dec,tri,pca} ...]
- Method(s) used to calculate disruption scores. Use abbreviations as follows:
- mse: Mean squared error
- corr: Spearman correlation
- ssi: Structural similarity index measure
- scc: Stratum adjusted correlation coefficient
- ins: Insulation
- di: Directionality index
- dec: Contact decay
- tri: Triangle method
- pca: Principal component method.
- (default: ['mse', 'corr'])
- --shift_by SHIFT_WINDOW [SHIFT_WINDOW ...]
- Values for shifting prediciton windows inputted as space-separated integers (e.g. -1 0 1). Values outside of range -450000 ≤ x ≤ 450000 will be ignored. Prediction windows at the edge of chromosome arms will only be shifted in the direction that is possible (ex. for window at chrom start, a -1 shift will be treated as a 1 shift since it is not possible to shift left).
- (default: [0])
- --shifts_file SHIFTS_FILE
- Path to file with values to shift each variant by. Required when not all variants are shifted by the same amount. File should be a text file with the same number of rows as the input file and a column which contains a shift value (negative or positive integer) by which to shift the variants in the corresponding row of the input file. (default: None)
- --file OUT_FILE Prefix for output files. Saved files will overwrite any existing files.
- (default: SuPreMo)
- --dir OUT_DIR Output directory. If directory already exists, files will be saved in existing directory. If the same files already exists in that directory, new files will overwrite them. (default: SuPreMo_output)
- --limit SVLEN_LIMIT Maximum length of variants to be scored. Filtering out variants that are too big can save time and memory. If not specified, will be set to 2/3 of seq_len.
- (default: None)
- --seq_len SEQ_LEN Length for sequences to generate. Default value is based on Akita requirement. If non-default value is set, get_Akita_scores must be false.
- (default: 1048576)
- --revcomp {no_revcomp,add_revcomp,only_revcomp}
- Option to use the reverse complement of the sequence:
- no_revcomp: no, only use the standard sequence;
- add_revcomp: yes, use both the standard sequence and its reverse complement;
- only_revcomp: yes, only use the reverse complement of the sequence.
- The reverse complement of the sequence is only taken with 0 shift.
- (default: ['no_revcomp'])
- --augment
- Only applicable if --get_Akita_scores is specified. Get the mean and median scores from sequences with specified shifts and reverse complement. If augment is used but shift and revcomp are not specified, the following four sequences will be used:
- 1) no augmentation: 0 shift and no reverse complement,
- 2) +1bp shift and no reverse complement,
- 3) -1bp shift and no reverse complement,
- 4) 0 shift and take reverse complement.
- (default: False)
- --get_seq Save sequences for the reference and alternate alleles in fa file format. If --get_seq is not specified, must specify --get_Akita_scores.
- Sequence name format: {var_index}_{shift}_{revcomp_annot}_{seq_index}_{var_rel_pos}
- var_index: input row number, followed by _0, _1, etc for each allele of variants with multiple alternate alleles;
- shift: integer that window is shifted by;
- revcomp_annot: present only if reverse complement of sequence was taken;
- seq_index: index for sequences generated for that variant: 0-1 for non-BND reference and alternate sequences and 0-2 for BND left and right reference sequence and alternate sequence;
- var_rel_pos: relative position of variant in sequence: list of two for non-BND variant positions in reference and alternate sequence and an integer for BND breakend position in reference and alternate sequences.
- There are 2-3 entries per prediction (2 for non-BND variants and 3 for BND variants).
- To read fasta file: pysam.Fastafile(filename).fetch(seqname, start, end).upper().
- To get sequence names in fasta file: pysam.Fastafile(filename).references.
- (default: False)
- --get_tracks Save disruption score tracks (448 bins) in npy file format. Only possible for mse and corr scores.
- Dictionary item name format: {var_index}_{track}_{shift}_{revcomp_annot}
- var_index: input row number, followed by _0, _1, etc for each allele of variants with multiple alternate alleles;
- track: disruption score track specified;
- shift: integer that window is shifted by;
- revcomp_annot: present only if reverse complement of sequence was taken.
- There is 1 entry per prediction: a 448x1 array.
- To read into a dictionary in python: np.load(filename, allow_pickle="TRUE").item()
- (default: False)
- --get_maps Save predicted contact frequency maps in npy file format.
- Dictionary item name format: {var_index}_{shift}_{revcomp_annot}
- var_index: input row number, followed by _0, _1, etc for each allele of variants with multiple alternate alleles;
- shift: integer that window is shifted by;
- revcomp_annot: present only if reverse complement of sequence was taken.
- There is 1 entry per prediction. Each entry contains the following: 2 (3 for chromosomal rearrangements) arrays that correspond to the upper right triangle of the predicted contact frequency maps, the relative variant position in the map, and the first coordinate of the sequence that the map corresponds to.
- To read into a dictionary in python: np.load(filename, allow_pickle="TRUE").item()
- (default: False)
- --get_Akita_scores Get disruption scores. If --get_Akita_scores is not specified, must specify --get_seq. Scores saved in a dataframe with the same number of rows as the input. For multiple alternate alleles, the scores are separated by a comma. To convert the scores from strings to integers, use float(x), after separating rows with multiple alternate alleles. Scores go up to 20 decimal points.
- (default: False)
- --nrows NROWS Number of rows (perturbations) to read at a time from input. When dealing with large inputs, selecting a subset of rows to read at a time allows scores to be saved in increments and uses less memory. Files with scores and filtered out variants will be temporarily saved in output direcotry. The file names will have a suffix corresponding to the set of nrows (0-based), for example for an input with 2700 rows and with nrows = 1000, there will be 3 sets. At the end of the run, these files will be concatenated into a comprehensive file and the temporary files will be removed.
- (default: 1000)
- '''
- # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
- # Parse through arguments
- import argparse
- class CustomFormatter(argparse.RawTextHelpFormatter, argparse.ArgumentDefaultsHelpFormatter):
- pass
- parser = argparse.ArgumentParser(
- prog = 'SuPreMo',
- formatter_class = CustomFormatter,
- description='''Pipeline for generating mutated sequences for input into predictive models and for scoring variants for disruption to genome folding.''')
- parser.add_argument('in_file',
- metavar = 'Input file',
- help = ''' Input file with variants. Accepted formats are:
- VCF file,
- TSV file (SVannot output),
- BED file,
- TXT file.
- Can be gzipped. Coordinates should be 1-based and left open (,], except for SNPs (follow vcf 4.1/4.2 specifications). See custom_perturbations.ipynb for BED or TXT file input specifications.
- ''',
- type = str)
- parser.add_argument('--sequences',
- dest = 'sequences',
- help = ''' Input fasta file with sequences. This file is one outputted by SuPreMo.
- ''',
- type = str,
- required = False)
- parser.add_argument('--fa',
- dest = 'fasta',
- help = '''Path to reference genome fasta file. Default: data/{genome}.fa, where genome is hg38 or hg19.
- ''',
- type = str,
- required = False)
- parser.add_argument('--roi',
- dest = 'roi',
- help = '''Path to text file with regions of intererst (roi) to use for scoring. Tab-delimited text file with required columns: chrom, start, end or chrom, start1, end1, start2, end2 for paired regions. Alternatively, use 'genes' to weight 2 kb window around protein-coding gene transcription start sites.
- ''',
- type = str,
- required = False)
- parser.add_argument('--roi_scales',
- dest = 'roi_scales',
- nargs = '+',
- help = '''Scale by which regions of interest (roi) are upweighted compared to the rest of the map. Default is 10 meaning regions of interest have weights of 10 whereas the rest of the bins have weights of 1.
- ''',
- type = int,
- default = [10],
- required = False)
- parser.add_argument('--genome',
- dest = 'genome',
- nargs = 1,
- help = '''Genome to be used: hg19 or hg38.''',
- type = str,
- choices = ['hg19', 'hg38'],
- default = ['hg38'],
- required = False)
- parser.add_argument('--Akita_cell_types',
- dest = 'Akita_cell_types',
- nargs = '+',
- help = '''Cell types to make Akita predictions for: HFF, H1hESC, GM12878, IMR90, and/or HCT116.
- ''',
- type = str,
- choices = ['HFF', 'H1hESC', 'GM12878', 'IMR90', 'HCT116'],
- default = ['HFF'],
- required = False)
- parser.add_argument('--scores',
- dest = 'scores',
- nargs = '+',
- help = '''
- Method(s) used to calculate disruption scores. Use abbreviations as follows:
- mse: Mean squared error
- corr: Spearman correlation
- ssi: Structural similarity index measure
- scc: Stratum adjusted correlation coefficient
- ins: Insulation
- di: Directionality index
- dec: Contact decay
- tri: Triangle method
- pca: Principal component method.
- ''',
- type = str,
- choices = ['mse', 'corr', 'ssi', 'scc', 'ins', 'di', 'dec', 'tri', 'pca'],
- default = ['mse', 'corr'],
- required = False)
- parser.add_argument('--shift_by',
- dest = 'shift_window',
- nargs = '+',
- help = '''Values for shifting prediciton windows inputted as space-separated integers (e.g. -1 0 1). Values outside of range -450000 ≤ x ≤ 450000 will be ignored. Prediction windows at the edge of chromosome arms will only be shifted in the direction that is possible (ex. for window at chrom start, a -1 shift will be treated as a 1 shift since it is not possible to shift left).
- ''',
- type = int,
- default = [0],
- required = False)
- parser.add_argument('--shifts_file',
- dest = 'shifts_file',
- help = 'Path to file with values to shift each variant by. Required when not all variants are shifted by the same amount. File should be a text file with the same number of rows as the input file and a column which contains a shift value (negative or positive integer) by which to shift the variants in the corresponding row of the input file.',
- type = str,
- required = False)
- parser.add_argument('--file',
- dest = 'out_file',
- help = '''Prefix for output files. Saved files will overwrite any existing files.
- ''',
- type = str,
- default = 'SuPreMo',
- required = False)
- parser.add_argument('--dir',
- dest = 'out_dir',
- help = 'Output directory. If directory already exists, files will be saved in existing directory. If the same files already exists in that directory, new files will overwrite them.',
- type = str,
- default = 'SuPreMo_output',
- required = False)
- parser.add_argument('--limit',
- dest = 'svlen_limit',
- help = '''Maximum length of variants to be scored. Filtering out variants that are too big can save time and memory. If not specified, will be set to 2/3 of seq_len.
- ''',
- type = int,
- required = False)
- parser.add_argument('--seq_len',
- dest = 'seq_len',
- help = '''Length for sequences to generate. Default value is based on Akita requirement. If non-default value is set, get_Akita_scores must be false.
- ''',
- type = int,
- default = 1048576,
- required = False)
- parser.add_argument('--revcomp',
- dest = 'revcomp',
- nargs = 1,
- help = '''
- Option to use the reverse complement of the sequence:
- no_revcomp: no, only use the standard sequence;
- add_revcomp: yes, use both the standard sequence and its reverse complement;
- only_revcomp: yes, only use the reverse complement of the sequence.
- The reverse complement of the sequence is only taken with 0 shift.
- ''',
- type = str,
- choices = ['no_revcomp', 'add_revcomp', 'only_revcomp'],
- default = ['no_revcomp'],
- required = False)
- parser.add_argument('--augment',
- dest = 'augment',
- help = '''
- Only applicable if --get_Akita_scores is specified. Get the mean and median scores from sequences with specified shifts and reverse complement. If augment is used but shift and revcomp are not specified, the following four sequences will be used:
- 1) no augmentation: 0 shift and no reverse complement,
- 2) +1bp shift and no reverse complement,
- 3) -1bp shift and no reverse complement,
- 4) 0 shift and take reverse complement.
- ''',
- action='store_true',
- required = False)
- parser.add_argument('--get_seq',
- dest = 'get_seq',
- help = '''Save sequences for the reference and alternate alleles in fa file format. If --get_seq is not specified, must specify --get_Akita_scores.
- Sequence name format: {var_index}_{shift}_{revcomp_annot}_{seq_index}_{var_rel_pos}
- var_index: input row number, followed by _0, _1, etc for each allele of variants with multiple alternate alleles;
- shift: integer that window is shifted by;
- revcomp_annot: present only if reverse complement of sequence was taken;
- seq_index: index for sequences generated for that variant: 0-1 for non-BND reference and alternate sequences and 0-2 for BND left and right reference sequence and alternate sequence;
- var_rel_pos: relative position of variant in sequence: list of two for non-BND variant positions in reference and alternate sequence and an integer for BND breakend position in reference and alternate sequences.
- There are 2-3 entries per prediction (2 for non-BND variants and 3 for BND variants).
- To read fasta file: pysam.Fastafile(filename).fetch(seqname, start, end).upper().
- To get sequence names in fasta file: pysam.Fastafile(filename).references.
- ''',
- action = 'store_true',
- required = False)
- parser.add_argument('--get_tracks',
- dest = 'get_tracks',
- help = '''Save disruption score tracks (448 bins) in npy file format. Only possible for mse and corr scores.
- Dictionary item name format: {var_index}_{track}_{shift}_{revcomp_annot}
- var_index: input row number, followed by _0, _1, etc for each allele of variants with multiple alternate alleles;
- track: disruption score track specified;
- shift: integer that window is shifted by;
- revcomp_annot: present only if reverse complement of sequence was taken.
- There is 1 entry per prediction: a 448x1 array.
- To read into a dictionary in python: np.load(filename, allow_pickle="TRUE").item()
- ''',
- action='store_true',
- required = False)
- parser.add_argument('--get_maps',
- dest = 'get_maps',
- help = '''Save predicted contact frequency maps in npy file format.
- Dictionary item name format: {var_index}_{shift}_{revcomp_annot}
- var_index: input row number, followed by _0, _1, etc for each allele of variants with multiple alternate alleles;
- shift: integer that window is shifted by;
- revcomp_annot: present only if reverse complement of sequence was taken.
- There is 1 entry per prediction. Each entry contains the following: 2 (3 for chromosomal rearrangements) arrays that correspond to the upper right triangle of the predicted contact frequency maps, the relative variant position in the map, and the first coordinate of the sequence that the map corresponds to.
- To read into a dictionary in python: np.load(filename, allow_pickle="TRUE").item()
- ''',
- action='store_true',
- required = False)
- parser.add_argument('--get_Akita_scores',
- dest = 'get_Akita_scores',
- help = '''Get disruption scores. If --get_Akita_scores is not specified, must specify --get_seq. Scores saved in a dataframe with the same number of rows as the input. For multiple alternate alleles, the scores are separated by a comma. To convert the scores from strings to integers, use float(x), after separating rows with multiple alternate alleles. Scores go up to 20 decimal points.
- ''',
- action='store_true',
- required = False)
- parser.add_argument('--nrows',
- dest = 'nrows',
- help = '''Number of rows (perturbations) to read at a time from input. When dealing with large inputs, selecting a subset of rows to read at a time allows scores to be saved in increments and uses less memory. Files with scores and filtered out variants will be temporarily saved in output direcotry. The file names will have a suffix corresponding to the set of nrows (0-based), for example for an input with 2700 rows and with nrows = 1000, there will be 3 sets. At the end of the run, these files will be concatenated into a comprehensive file and the temporary files will be removed.
- ''',
- type = int,
- default = 1000,
- required = False)
- args = parser.parse_args()
- in_file = args.in_file
- input_sequences = args.sequences
- fasta_path = args.fasta
- roi = args.roi
- roi_scales = args.roi_scales
- genome = args.genome[0]
- Akita_cell_types = args.Akita_cell_types
- scores_to_use = args.scores
- shift_by = args.shift_window
- shifts_file = args.shifts_file
- out_file = args.out_file
- out_dir = args.out_dir
- svlen_limit = args.svlen_limit
- seq_len = args.seq_len
- revcomp = args.revcomp[0]
- augment = args.augment
- get_seq = args.get_seq
- get_tracks = args.get_tracks
- get_maps = args.get_maps
- get_Akita_scores = args.get_Akita_scores
- var_set_size = args.nrows
- __version__ = '1.0'
- # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
- # Adjust inputs from arguments
- # Handle paths
- import os
- from pathlib import Path
- # This file path and repo path
- repo_path = Path(__file__).parents[1]
- # Output directory and file
- if not os.path.exists(out_dir):
- os.mkdir(out_dir)
- out_file = os.path.join(out_dir, out_file)
- # Data path
- chrom_lengths_path = f'{repo_path}/data/chrom_lengths_{genome}'
- centromere_coords_path = f'{repo_path}/data/centromere_coords_{genome}'
- # Handle argument dependencies
- use_roi = False
- if fasta_path is None:
- fasta_path = f'{repo_path}/data/{genome}.fa'
- if seq_len != 1048576:
- get_Akita_scores = False
- if svlen_limit is None:
- svlen_limit = 2/3*seq_len
- elif svlen_limit > 2/3*seq_len:
- raise ValueError("Maximum SV length limit should not be >2/3 of sequence length.")
- if not get_seq and not get_Akita_scores:
- raise ValueError('Either get_seq and/or get_Akita_scores must be specified.')
- if roi is not None and not get_Akita_scores:
- raise ValueError('get_Akita_scores must be specified to use roi.')
- # Adjust shift input: Remove shifts that are outside of allowed range
- max_shift = 0.4*seq_len
- shift_by = [x for x in shift_by if x > -max_shift and x < max_shift]
- # Adjust input for taking the reverse complement
- if revcomp == 'no_revcomp':
- revcomp_decision = [False]
- elif revcomp == 'add_revcomp':
- revcomp_decision = [False, True]
- elif revcomp == 'only_revcomp':
- revcomp_decision = [True]
- if augment and shift_by == [0] and revcomp == 'no_revcomp':
- shift_by = [-1,0,1]
- revcomp_decision = [False, True]
- revcomp_decision_i = revcomp_decision
- import pysam
- if input_sequences is not None:
- seq_names = pysam.Fastafile(input_sequences).references
- # Create dictionaries to save sequences, maps, and disruption score tracks, if specified
- if get_seq:
- sequences = {}
- if get_maps:
- variant_maps = {}
- if get_Akita_scores == False:
- get_Akita_scores = True
- print('Must get scores to get maps. --get_Akita_scores was not specified but will be applied.')
- if get_tracks:
- variant_tracks = {}
- if get_Akita_scores == False:
- get_Akita_scores = True
- print('Must get scores to get tracks. --get_Akita_scores was not specified but will be applied.')
- # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
- # Read in (and adjust) data
- import pandas as pd
- chrom_lengths = pd.read_table(chrom_lengths_path, header = None, names = ['CHROM', 'chrom_max'])
- centromere_coords = pd.read_table(centromere_coords_path, sep = '\t')
- fasta_open = pysam.Fastafile(fasta_path)
- # Assign necessary values to variables across module
- # Reading utilities
- import sys
- sys.path.insert(0, 'scripts/')
- import reading_utils
- reading_utils.var_set_size = var_set_size
- # SuPreMo utilities
- import get_seq_utils
- get_seq_utils.fasta_open = fasta_open
- get_seq_utils.chrom_lengths = chrom_lengths
- get_seq_utils.centromere_coords = centromere_coords
- get_seq_utils.svlen_limit = svlen_limit
- get_seq_utils.seq_length = seq_len
- get_seq_utils.half_patch_size = round(seq_len/2)
- # SuPreMo-Akita utilities
- if get_Akita_scores:
- import get_Akita_scores_utils
- get_Akita_scores_utils.chrom_lengths = chrom_lengths
- get_Akita_scores_utils.centromere_coords = centromere_coords
- if roi is None or roi_scales == [0]:
- use_roi = False
- else:
- use_roi = True
- import get_roi_utils
- get_roi_utils.pred_len = get_Akita_scores_utils.target_length_cropped * get_Akita_scores_utils.bin_size
- get_roi_utils.d = get_Akita_scores_utils.bin_size * 5
- get_Akita_scores_utils.roi_coords_BED = get_roi_utils.get_roi(roi, genome)
- from pybedtools import BedTool
- get_Akita_scores_utils.BedTool = BedTool
- get_Akita_scores_utils.roi_scales = roi_scales
- import sys
- import numpy as np
- nt = ['A', 'T', 'C', 'G']
- var_set = 0
- var_set_list = []
- print(f'Log file being saved here: {out_file}_log')
- while True:
- # Read in variants
- variants = reading_utils.read_input(in_file, var_set)
- if len(variants) == 0:
- break
- if shifts_file is not None:
- variants = pd.concat([variants,
- pd.read_csv(shifts_file, sep = '\t', low_memory=False, names = ['shift_by'],
- skiprows = 1 + var_set*var_set_size, nrows = var_set_size)],
- axis = 1)
- # Index input based on row number and create output with same indexes
- variants['var_index'] = list(range(var_set*var_set_size, var_set*var_set_size + len(variants)))
- variants['var_index'] = variants['var_index'].astype(str)
- # If there are multiple alternate alleles, split those into new rows and indexes
- if any([',' in x for x in variants.ALT]):
- variants = (variants
- .set_index(['CHROM', 'POS', 'REF', 'var_index'])
- .apply(lambda x: x.str.split(',').explode())
- .reset_index())
- g = variants.groupby(['var_index'])
- variants.loc[g['var_index'].transform('size').gt(1),
- 'var_index'] += '-'+g.cumcount().astype(str)
- variant_scores = pd.DataFrame({'var_index':variants.var_index})
- if use_roi:
- variant_roi_ids = pd.DataFrame({'var_index':variants.var_index})
- # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
- # Filter out variants that cannot be scored
- # Get indexes for variants to exclude
- # Exclude mitochondrial variants
- chrM_var = pd.DataFrame({'var_index' : list(variants[variants.CHROM == 'chrM'].var_index),
- 'reason' : ' Mitochondrial chromosome.'})
- # Exclude variants larger than limit
- if 'SVLEN' in variants.columns:
- too_long_var = pd.DataFrame({'var_index' : [y for x,y in zip(variants.SVLEN, variants.var_index)
- if not pd.isnull(x) and abs(int(x)) > svlen_limit],
- 'reason' : f' SV longer than {svlen_limit}.'})
- unsuitable_var = pd.DataFrame({'var_index' : [y for x,y,z in zip(variants.SVTYPE, variants.var_index, variants.ALT)
- if not pd.isnull(x) and
- x not in ["DEL", "DUP", "INV", "BND"] and
- all([g not in nt for g in z])],
- 'reason' : ' SV type not compatible.'})
- else:
- too_long_var = pd.DataFrame()
- unsuitable_var = pd.DataFrame()
- filtered_out = pd.concat([chrM_var, too_long_var], axis = 0)
- filtered_out = pd.concat([filtered_out, unsuitable_var], axis = 0)
- filtered_out.var_index = filtered_out.var_index.astype('str')
- # Save filtered out variants into file
- filtered_out.to_csv(f'{out_file}_filtered_out_{var_set}', sep = ':', index = False, header = False)
- # Exclude
- variants = variants[[x not in filtered_out.var_index.values for x in variants.var_index]]
- # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
- # Run: Make Akita predictions and calculate disruption scores
- # Create log file to save standard output with error messages
- std_output = sys.stdout
- log_file = open(f'{out_file}_log_{var_set}','w')
- sys.stdout = log_file
- # Loop through each row (not index) and get disruption scores
- for i in range(len(variants)):
- variant = variants.iloc[i]
- var_index = variant.var_index
- CHR = variant.CHROM
- POS = variant.POS
- REF = variant.REF
- ALT = variant.ALT
- if 'SVTYPE' in variants.columns:
- END = variant.END
- SVTYPE = variant.SVTYPE
- SVLEN = variant.SVLEN
- else:
- END = np.nan
- SVTYPE = np.nan
- SVLEN = 0
- if shifts_file is not None:
- shift_by = [int(variant.shift_by)]
- for shift in shift_by:
- # Take reverse complement only with 0 shift
- if shift != 0 & True in revcomp_decision:
- revcomp_decision_i = [False]
- else:
- revcomp_decision_i = revcomp_decision
- for revcomp in revcomp_decision_i:
- try:
- if revcomp:
- revcomp_annot = '_revcomp'
- else:
- revcomp_annot = ''
- if input_sequences is not None:
- # Generate sequences_i from sequence input
- if revcomp_annot == '':
- sequence_names = [x for x in seq_names if x.startswith(f'{var_index}_{shift}') and
- 'revcomp' not in x]
- elif revcomp_annot == '_revcomp':
- sequence_names = [x for x in seq_names if x.startswith(f'{var_index}_{shift}{revcomp_annot}')]
- sequences_i = []
- for sequence_name in sequence_names:
- sequences_i.append(pysam.Fastafile(input_sequences).fetch(sequence_name, 0, seq_len).upper())
- sequences_i.append([int(x) for x in sequence_name.split('[')[1].split(']')[0].split('_')])
- else:
- # Create sequences_i from variant input
- sequences_i = get_seq_utils.get_sequences_SV(CHR, POS, REF, ALT, END, SVTYPE, shift, revcomp)
- if get_seq:
- # Get relative position of variant in sequence
- var_rel_pos = str(sequences_i[-1]).replace(', ', '_')
- for ii in range(len(sequences_i[:-1][:3])):
- sequences[f'{var_index}_{shift}{revcomp_annot}_{ii}_{var_rel_pos}'] = sequences_i[:-1][ii]
- if get_Akita_scores:
- scores = get_Akita_scores_utils.get_scores(CHR, POS, SVTYPE, SVLEN,
- sequences_i, scores_to_use,
- shift, revcomp,
- get_tracks, get_maps, use_roi, Akita_cell_types)
- if get_tracks:
- for track in [x for x in scores.keys() if 'track' in x]:
- variant_tracks[f'{var_index}_{track}_{shift}{revcomp_annot}'] = scores[track]
- del scores[track]
- if get_maps:
- for map in [x for x in scores.keys() if 'map' in x]:
- variant_maps[f'{var_index}_{map}_{shift}{revcomp_annot}'] = scores[map]
- del scores[map]
- if shifts_file is not None:
- shift = 'shifted'
- if use_roi:
- for roi_ids in [x for x in scores.keys() if 'roi_id' in x]:
- variant_roi_ids.loc[variant_scores.var_index == var_index, f'roi_ids_{shift}'] = scores[roi_ids]
- del scores[roi_ids]
- for score in scores:
- variant_scores.loc[variant_scores.var_index == var_index,
- f'{score}_{shift}{revcomp_annot}'] = scores[score]
- print(str(var_index) + ' (' + str(shift) + f' shift{revcomp_annot})')
- except Exception as e:
- print(str(var_index) + ' (' + str(shift) + f' shift{revcomp_annot})' + ': Error:', e)
- pass
- # Write standard output with error messages and warnings to log file
- sys.stdout = std_output
- log_file.close()
- # Combine results from all sets
- # Write sequences to fasta file
- if get_seq:
- if var_set == 0:
- sequences_all = sequences.copy()
- else:
- sequences_all.update(sequences)
- # Write scores to data frame
- if get_Akita_scores:
- # Take average of augmented sequences
- if augment:
- for score in scores:
- cols = [x for x in variant_scores.columns if score in x]
- variant_scores[f'{score}_mean'] = variant_scores[cols].mean(axis = 1)
- variant_scores[f'{score}_median'] = variant_scores[cols].median(axis = 1)
- variant_scores.drop(cols, axis = 1, inplace = True)
- # Convert scores from float to string so you can merge scores for variants with multiple alleles
- for col in variant_scores.iloc[:,1:].columns:
- variant_scores[col] = [format(x, '.20f') for x in variant_scores[col]]
- # Join scores for alternate alleles, separated by a comma
- if any(['-' in x for x in variant_scores.var_index]):
- variant_scores['var_index'] = variant_scores.var_index.str.split('-').str[0]
- variant_scores = (variant_scores
- .set_index(['var_index'], drop = False)
- .rename(columns = {'var_index':'var_index2'})
- .groupby('var_index2')
- .transform(','.join)
- .reset_index()
- .drop_duplicates())
- if var_set == 0:
- variant_scores.to_csv(f'{out_file}_scores_{var_set}', sep = '\t', index = False)
- if use_roi:
- variant_roi_ids.to_csv(f'{out_file}_roi_ids_{var_set}', sep = '\t', index = False)
- else:
- variant_scores.to_csv(f'{out_file}_scores_{var_set}', sep = '\t', index = False, header = False)
- if use_roi:
- variant_roi_ids.to_csv(f'{out_file}_roi_ids_{var_set}', sep = '\t', index = False, header = False)
- if get_tracks:
- if var_set == 0:
- variant_tracks_all = variant_tracks.copy()
- else:
- variant_tracks_all.update(variant_tracks)
- if get_maps:
- if var_set == 0:
- variant_maps_all = variant_maps.copy()
- else:
- variant_maps_all.update(variant_maps)
- var_set_list.append(var_set)
- var_set += 1
- # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
- # Save results
- # Write sequences to fasta file
- if get_seq:
- sequences_fasta = open(f'{out_file}_sequences.fa','w')
- for seq_name, sequence in sequences_all.items():
- seq_name_line = ">" + seq_name + "\n"
- sequences_fasta.write(seq_name_line)
- sequence_line = sequence + "\n"
- sequences_fasta.write(sequence_line)
- sequences_fasta.close()
- # Write disruption tracks and/or predictions to array
- if get_tracks:
- np.save(f'{out_file}_tracks.npy', variant_tracks_all)
- if get_maps:
- np.save(f'{out_file}_maps.npy', variant_maps_all)
- # Combine subset files into one
- os.system(f'rm -f {out_file}_filtered_out; \
- for file in {out_file}_filtered_out_*; \
- do cat "$file" >> {out_file}_filtered_out && rm "$file"; \
- done')
- os.system(f'rm -f {out_file}_log; \
- for file in {out_file}_log_*; \
- do cat "$file" >> {out_file}_log && rm "$file"; \
- done')
- if get_Akita_scores:
- os.system(f'rm -f {out_file}_scores; \
- for file in {out_file}_scores_*; \
- do cat "$file" >> {out_file}_scores && rm "$file"; \
- done')
- if use_roi:
- os.system(f'rm -f {out_file}_roi_ids; \
- for file in {out_file}_roi_ids_*; \
- do cat "$file" >> {out_file}_roi_ids && rm "$file"; \
- done')
- # Adjust log file to only have 1 row per variant
- if os.path.exists(f'{out_file}_log'):
- log_file = pd.read_csv(f'{out_file}_log', names = ['output'], sep = '\t')
- # Move warnings (printed 1 line before variant) to variant line
- indexes = np.array([[index, index+1] for (index, item) in enumerate(log_file.output) if item.startswith('Warning')])
- if len(indexes) != 0:
- log_file.loc[indexes[:,1],'output'] = [x+': '+y for x,y in zip(list(log_file.loc[indexes[:,1],'output']),
- list(log_file.loc[indexes[:,0],'output']))]
- log_file.drop(indexes[:,0], axis = 0, inplace = True)
- log_file.to_csv(f'{out_file}_log', sep = '\t', header = None, index = False)
SuPreMo.py at commit 9799281, no license · at the source
Overview
- Gladstone Institute of Data Science and Biotechnology, San Francisco, California 94158, USA
- Department of Epidemiology and Biostatistics, University of California San Francisco, California 94158, USA
- Institute for Human Genetics, University of California San Francisco, San Francisco, California 94143, USA
- Department of Neurology, University of California San Francisco, San Francisco, California 94143, USA
- Weill Institute for Neurosciences, University of California San Francisco, San Francisco, California 94158, USA
- Bakar Computational Health Sciences Institute, University of California, San Francisco, California 94143, USA
- Chan Zuckerberg Biohub, San Francisco, California 94158, USA
Abstract
Three-dimensional genome organization plays a critical role in gene regulation, and disruptions can lead to developmental disorders by altering the contact between genes and their distal regulatory elements. Structural variants (SVs) can disturb local genome organization, such as the merging of topologically associating domains upon boundary deletion. Testing large numbers of SVs experimentally for their effects on chromatin structure and gene expression is time and cost prohibitive. To address this, we propose a computational approach to predict SV impacts on genome folding, which can help prioritize causal hypotheses for functional testing. We develop a weighted scoring method that measures chromatin contact changes specifically affecting regions of interest, such as regulatory elements or promoters, and implement it in the SuPreMo-Akita software. With this tool, we rank hundreds of de novo SVs (dnSVs) from autism spectrum disorder (ASD) individuals and their unaffected siblings based on predicted disruptions to nearby neuronal regulatory interactions. This reveals that putative cis-regulatory element interactions (CREints) are more disrupted by dnSVs from ASD probands versus unaffected siblings. We prioritize candidate variants that disrupt ASD CREints and validate our top-ranked locus using isogenic excitatory neurons with and without the dnSV, confirming accurate predictions of disrupted chromatin contacts. This study suggests that disrupted genome folding is a potential genetic mechanism in a subset of ASD cases and provides a general strategy for prioritizing variants predicted to disrupt regulatory interactions across tissues.
Reproduced under the paper's license (CC BY-NC), from the paper cited above.
Repository
Its files are read in the Code ↔ Paper reader above, with 7 matches between paragraphs and lines of code.
ketringjoni/SuPreMo
9799281ec9b4ea5ea03702931e57950903f75424, 14 November 2024Availability: 1 check, the latest on 26 September 2026: the link answers
- 26 September 2026: the link answers
16 files
- scripts/
SuPreMo.py , Python, 876 lines, 4 matches - scripts/
get_Akita_scores_utils.p , Python, 718 lines, 1 matchy - scripts/
get_roi_utils.py , Python, 222 lines - scripts/
get_seq_utils.py , Python, 1,033 lines - scripts/
hicrep.py , Python, 84 lines - scripts/
plotting_utils.py , Python, 691 lines - scripts/
reading_utils.py , Python, 184 lines - scripts/
scoring.py , Python, 746 lines - scripts/
test_install_SuPreMo-Aki , Python, 41 linesta.py - scripts/
test_install_SuPreMo.py , Python, 16 lines - walkthroughs/
SuPreMo-Akita_walkthroug , Jupyter, 343 linesh.ipynb - walkthroughs/
SuPreMo-Akita_weighted_s , Jupyter, 354 lines, 2 matchescores.ipynb - walkthroughs/
SuPreMo_walkthrough.ipyn , Jupyter, 358 linesb - walkthroughs/
custom_perturbations.ipy , Jupyter, 231 linesnb - walkthroughs/
example_application.ipyn , Jupyter, 396 linesb - README.md, Text, 290 lines
Code availability
SuPreMo V2 code is available as Supplemental Code and at GitHub (https://
Reproduced under the paper's license (CC BY-NC), from the paper cited above.
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:
- 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
- 15 scripts, each with its path and the digest of its content;
- 7 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
Datasets cited
- geo:GSE281283, at NCBI GEO; found in the text, “Data access”
Other data links
- ncbi.nlm.nih.gov/
geo , NCBI; found in the text, “Data access”
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 3, 28 September 2026
- Authors: added Xingjie Ren (0000-0002-2348-4162); Amanda Everitt (0000-0001-9720-1922); Yin Shen (0000-0001-9901-5613); removed Xingjie Ren; Amanda Everitt; Yin Shen
Version 1, 27 September 2026: the first record
Recorded: type, language, journal, volume, issue, pages, dates, 5 authors, 9 MeSH terms, 10 funders, 55 references.
Cite
This paper
Gjoni, K., Ren, X., Everitt, A., Shen, Y., & Pollard, K. S. (2026). De novo structural variants in autism spectrum disorder disrupt distal regulatory interactions of neuronal genes. Genome research, 36(9), 1825-1835. https://
BibTeX
@article{gjoni2026de,
author = {Gjoni, Ketrin and Ren, Xingjie and Everitt, Amanda and Shen, Yin and Pollard, Katherine S},
title = {{De novo structural variants in autism spectrum disorder disrupt distal regulatory interactions of neuronal genes}},
journal = {Genome research},
year = {2026},
month = sep,
volume = {36},
number = {9},
pages = {1825--1835},
publisher = {Cold Spring Harbor Laboratory Press},
issn = {1088-9051},
doi = {10.1101/
url = {https://
pmid = {42680562},
pmcid = {PMC13534236}
}
RIS
TY - JOUR
AU - Gjoni, Ketrin
AU - Ren, Xingjie
AU - Everitt, Amanda
AU - Shen, Yin
AU - Pollard, Katherine S
TI - De novo structural variants in autism spectrum disorder disrupt distal regulatory interactions of neuronal genes
T2 - Genome research
J2 - Genome Res
PY - 2026
DA - 2026/
VL - 36
IS - 9
SP - 1825
EP - 1835
SN - 1088-9051
PB - Cold Spring Harbor Laboratory Press
DO - 10.1101/
UR - https://
LA - en
ER -
CSL-JSON
{
"id": "10.1101/
"type": "article-journal",
"title": "De novo structural variants in autism spectrum disorder disrupt distal regulatory interactions of neuronal genes",
"container-title": "Genome research",
"author": [
{
"family": "Gjoni",
"given": "Ketrin"
},
{
"family": "Ren",
"given": "Xingjie"
},
{
"family": "Everitt",
"given": "Amanda"
},
{
"family": "Shen",
"given": "Yin"
},
{
"family": "Pollard",
"given": "Katherine S"
}
],
"container-title-short":
"volume": "36",
"issue": "9",
"page": "1825-1835",
"DOI": "10.1101/
"PMID": "42680562",
"PMCID": "PMC13534236",
"ISSN": "1088-9051",
"publisher": "Cold Spring Harbor Laboratory Press",
"URL": "https://
"language": "en",
"issued": {
"date-parts": [
[
2026,
9,
1
]
]
}
}
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.1016/j.xhgg.2026.100652 [code]
- CRISPR-engineered deletion of POGZ alters transcription factor binding at promoters of genes involved in synaptic signaling.Journal: HGG advancesIn common: pysam, BEDTools, pandas, 3 other tools, autism, cellular / molecular, 6 references
- [2] doi:10.1038/s41467-026-76675-1 [code]
- Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.Journal: Nature communicationsIn common: pysam, Biopython, BEDTools, 5 other tools, autism, 3 references
- [3] doi:10.1038/s41467-026-71877-z [code]
- Hi-Compass: a depth-aware deep learning framework for predicting cell-type-specific 3D genome organization from single-cell to spatial resolution.Journal: Nature communicationsIn common: pysam, BEDTools, scikit-image, 3 other tools, 4 references
- [4] doi:10.21203/rs.3.rs-9927928/v1 [code]
- Genome-wide and allele-resolved maps of the radial architecture of the mouse genomeJournal: Research Square (preprint)In common: pysam, BEDTools, pandas, 3 other tools, 5 references
- [5] doi:10.1038/s41467-026-75700-7 [code]
- Gene regulatory innovations from transposable elements in primate cerebellum development.Journal: Nature communicationsIn common: pysam, Biopython, BEDTools, 6 other tools, cellular / molecular
- [6] doi:10.1126/sciadv.aed2952 [code]
- Activation of transposable elements is linked to a region- and cell type-specific interferon response in Parkinson's disease.Journal: Science advancesIn common: pysam, Biopython, BEDTools, 4 other tools, cellular / molecular, 3 references
- [7] doi:10.1038/s41592-026-03057-2 [code]
- CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.Journal: Nature methodsIn common: pysam, Biopython, BEDTools, 6 other tools
- [8] doi:10.1186/s13059-026-04177-w [code]
- Genomic sequence evolution underlying human neocortical interareal diversification.Journal: Genome biologyIn common: pysam, BEDTools, TensorFlow, 5 other tools, cellular / molecular, 1 reference
- [9] doi:10.1038/s41586-026-10391-0 [code]
- Cell-type-targeted mitochondrial transplantation rescues cell degeneration.Journal: NatureIn common: Biopython, TensorFlow, scikit-learn, 4 other tools, cellular / molecular, 3 references
- [10] doi:10.1038/s41586-026-10512-9 [code]
- Astrocyte glucocorticoid receptor signalling restricts neuronal plasticity.Journal: NatureIn common: pysam, Biopython, BEDTools, 4 other tools, cellular / molecular, 1 reference
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.
Claim this paper
Correct its record
Say what each link of this record is, remove the ones that are not the paper's, add the ones that are missing. The correction becomes a new version of the record, in its Versions section.
Validate its tracing map
You validate the map as this page shows it: 1 repository of the authors' code, each at its verified commit and with its license, 15 scripts, and 7 matches between paragraphs and code (see the Code and Map sections). It then receives a DOI on Zenodo, with you (your ORCID iD) and OSCR as its creators; the code itself is not deposited.
The map's fingerprint: sha256:43f6515919d296ca…
Add the badge to its README
The badge links the code to this page. Copy one of these into the README of the paper's code: only you decide where it goes, and nothing is changed for you.
Markdown
[, paste the snippet at the top, then “Commit changes…” and, to review it first, “Create a new branch and start a pull request”. You open the pull request; OSCR asks for no permission.
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.
