<jupyter_start><jupyter_text>Generate Random Genome SequencesThis notebook generates `N` random of the length `K` with the minimal distance `MIN_DIST` to any Ensembl gene.Note that some sequences may be missing in the human genome (= consists of 'N' letters only).  SetupInstallation for colab environment.<jupyter_code>!pip install biopython pyensembl
from google.colab import drive
drive.mount('/content/drive')
!pyensembl install --release 97 --species human
import pandas as pd
import numpy as np
import gzip
from tqdm.notebook import tqdm

from Bio import SeqIO   # for reading fasta files
from pyensembl import EnsemblRelease   # to get the gene list

ENSEMBL_RELEASE = 97
DNA_TOPLEVEL_FASTA_PATH = "/content/drive/My Drive/data/ensembl/Homo_sapiens.GRCh38.dna.toplevel.fa.gz"

# to generate random sequences
N = 50_000    # how many
K = 200       # how long
MIN_DIST = 50 # how far should they be from the gene
OUTPUT_FILE = '/content/drive/My Drive/data/random/random_seqs3.csv'   # where to save them
# OUTPUT_FILE = 'random_seqs.csv'

CHRS = [str(chr) for chr in range(1,23)] + ['X', 'Y', 'MT']<jupyter_output><empty_output><jupyter_text>Get length of each chromosome<jupyter_code>def get_chr_lengths(fasta_path=DNA_TOPLEVEL_FASTA_PATH):
    chr_lengths = {}
    
    with gzip.open(fasta_path, "rt") as handle:
        for record in tqdm(SeqIO.parse(handle, "fasta"), total=24):
            chr_lengths[record.id] = len(record.seq)
            if record.id == "MT": 
                # stop, do not read small contigs
                break
    
    # check that we have all chromosomes
    assert set(chr_lengths.keys()) == set(CHRS) 
    
    return chr_lengths

contig_lengths = get_chr_lengths()
contig_lengths
total_length = pd.Series(contig_lengths).sum()
total_length<jupyter_output><empty_output><jupyter_text>Get gene list<jupyter_code># release 97 uses human reference genome GRCh38
data = EnsemblRelease(ENSEMBL_RELEASE)
human_genes = data.genes()
len(human_genes)
human_genes[0]
human_genes_tuples = [(x.gene_id, x.gene_name, x.biotype, x.contig, x.start, x.end, x.strand) for x in human_genes]
human_genes_table = pd.DataFrame.from_records(human_genes_tuples, columns=["id", "symbol", "biotype", "chr", "start", "end", "strand"])
assert all(human_genes_table.start <= human_genes_table.end)

human_genes_table.head()
excluded_regions = human_genes_table[["chr", "start", "end"]][human_genes_table.chr.isin(CHRS)].copy()
excluded_regions.chr.value_counts()
excluded_regions.head()
total_gene_length = ((excluded_regions.end - excluded_regions.start) + 1).sum()
total_gene_length, total_gene_length/total_length
excluded_regions['start_dist'] = excluded_regions['start'] - MIN_DIST - K
excluded_regions['end_dist'] = excluded_regions['end'] + MIN_DIST<jupyter_output><empty_output><jupyter_text>Check on chr1 how many possible seqs are really there<jupyter_code>c = '1'
c_len = contig_lengths[c]
possible_seq_starts = np.ones(c_len, dtype=bool)

for i in tqdm(range(excluded_regions.shape[0])):
    if excluded_regions.chr.values[i] == c:
        start = max(excluded_regions.start_dist.values[i], 1)
        end = min(excluded_regions.end_dist.values[i], c_len)
        possible_seq_starts[(start-1):end] = 0
possible_seq_starts.sum(), possible_seq_starts.sum()/c_len
biotypes = pd.Series([g.biotype for g in human_genes])
biotypes.value_counts()<jupyter_output><empty_output><jupyter_text>Random sequence generation<jupyter_code>def get_random_chr(chr_lengths: pd.Series):
    chr_probs = chr_lengths / chr_lengths.sum()
    return CHRS[np.argwhere(np.random.multinomial(1, chr_probs))[0][0]]

def is_intersecting(c, pos, df_forbidden):
    intersecting = (df_forbidden.chr.values == c) & (df_forbidden.start_dist.values <= pos) & (df_forbidden.end_dist.values >= pos)
    return intersecting.any()

def get_random_pos(df_forbidden: pd.DataFrame, chr_lengths: pd.Series):
    c = get_random_chr(chr_lengths)
    c_len = chr_lengths[c]
    pos = np.random.randint(c_len) + 1
    
    while is_intersecting(c, pos, df_forbidden):
        pos = np.random.randint(c_len) + 1
    
    return c, pos
chr_lengths = pd.Series(contig_lengths)
get_random_chr(chr_lengths) 
is_intersecting('1', 11_000_000, excluded_regions)
get_random_pos(excluded_regions, chr_lengths)
random_seqs = [get_random_pos(excluded_regions, chr_lengths) for i in tqdm(range(N))]<jupyter_output><empty_output><jupyter_text>Get actual genomic sequences<jupyter_code>seqs = pd.DataFrame.from_records(random_seqs, columns=["chr", "start"])
seqs['end'] = seqs['start'] + K - 1
seqs['seq'] = ''
seqs.head()
def which(self):
    try:
        self = list(iter(self))
    except TypeError as e:
        raise Exception("""'which' method can only be applied to iterables.
        {}""".format(str(e)))
    indices = [i for i, x in enumerate(self) if bool(x) == True]
    return(indices)

with gzip.open(DNA_TOPLEVEL_FASTA_PATH, "rt") as handle:
    for record in tqdm(SeqIO.parse(handle, "fasta"), total=24):
        sel_seqs = which(seqs.chr == record.id)
        for i in sel_seqs:
            seqs.loc[i, "seq"] = str(record.seq[(seqs.start[i]-1):seqs.end[i]])
        
        if record.id == "MT": 
            # stop, do not read small contigs
            break
seqs.head()
len(seqs.seq.values[0])
seqs.seq.values[1]<jupyter_output><empty_output><jupyter_text>Save generated sequences to file<jupyter_code>seqs.to_csv(OUTPUT_FILE, index=False)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12183950164895566, 'fraction_numerical': 0.009893733968486624, 'mean_word_length': 4.8074468085106385, 'pattern_counts': {'":': 2, '<': 31, '<?xml version=': 0, '>': 30, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '301905', 'license': 'Apache-2.0', 'n_tokens_mistral': 2046, 'n_tokens_neox': 1943, 'n_words': 415, 'path': 'data/Random_Genome_Seqs.ipynb', 'repo_name': 'simecek/PseudoDNA_Generator', 'size': '5458'}