<jupyter_start><jupyter_text>Library design for CTP-05, Intronic RNAby <NAME>This library design is based on target regions designed by Stephen  0. Imports<jupyter_code>%run "E:\Users\puzheng\Documents\Startup_py3.py"
sys.path.append(r"E:\Users\puzheng\Documents")

import ImageAnalysis3
from ImageAnalysis3 import get_img_info, visual_tools, corrections, library_tools

from ImageAnalysis3.library_tools import LibraryDesigner as ld
from ImageAnalysis3.library_tools import LibraryTools as lt

%matplotlib notebook
print(os.getpid())<jupyter_output>95264<jupyter_text>1. Load probes and convert to dic<jupyter_code># folder for genomic info
genome_folder = r'E:\Genomes\hg38'

# Library directories
pool_folder = r'X:\Libraries\CTP-07'
# folder for sub-pool
library_folder = os.path.join(pool_folder, 'chr21_intronic')
# folder for fasta sequences
sequence_folder = os.path.join(library_folder, 'sequences')

# folder to save result probes
save_folder = os.path.join(library_folder, 'reports')
import csv
source_folder = r'E:\Users\puzheng\Documents\Libraries\CTP-05\Candidate_Probes'
probe_file = 'Chr21genes_96Intron.csv'
probes = []
with open(os.path.join(source_folder, probe_file),'r') as csvfile:
    reader = csv.reader(csvfile, delimiter=',', quotechar='|')
    header = next(reader,None)
    for row in reader:
        _source_probe = {h.lower():info for h,info in zip(header, row)}
        _probe = {}
        for _key,_value in _source_probe.items():
            if 'position' in _key and 'position' not in _probe:
                _probe['position'] = int(_value)
            if 'gene' in _key and 'gene' not in _probe:
                if '_' in _value:
                    _probe['gene'] = _value.split('_')[0]
                else:
                    _probe['gene'] = _value
            if 'sequence' in _key and 'sequence' not in _probe:
                _probe['sequence'] = _value
        probes.append(_probe)
# merge into dictionary
probe_dic = {}
for _probe in probes:
    if _probe['gene'] not in probe_dic:
        probe_dic[_probe['gene']] = [_probe]
        
    else:
        probe_dic[_probe['gene']].append(_probe)

        
position_file = "chr21_intron_positions.csv"
position_dic = {}
with open(os.path.join(source_folder, position_file),'r') as csvfile:
    reader = csv.reader(csvfile, delimiter=',', quotechar='|')
    for row in reader:
        position_dic[row[0]] = int(row[1])
        
# sort dictionary and append probe id, position
for _gene, _pb_list in probe_dic.items():
    _sorted_pb_list = sorted(_pb_list, key=lambda v:int(v['position']))
    # extract position
    _gene_start = position_dic[_gene]
    _gene_end = _sorted_pb_list[-1]['position'] + len(_sorted_pb_list[-1]['sequence'])+_gene_start
    for _i, _probe in enumerate(_sorted_pb_list):
        _sorted_pb_list[_i]['id'] = _i
        _sorted_pb_list[_i]['region'] = f'chr21:{_gene_start}-{_gene_end}'
        
print(len(probe_dic))
#>chr21:26842909-26847909_gene_ADAMTS1_pb_0_pos_0_readouts_[NDB_1147_u,NDB_1147_u,NDB_1147_u]

probe_dic['APP']<jupyter_output><empty_output><jupyter_text>3. Screening genes and probes<jupyter_code># keep the most 5' probes and exclude genes without enough probes
min_num_probe = 12 # this number is from MERFISH experience as well as Long Cai 2018 paper
max_num_probe = 80 # this is just an arbitrary number larger than 2x min_num_probe
# initialize
kept_probe_dic = {}
for _gene, _pb_list in sorted(probe_dic.items(), key=lambda v:position_dic[v[0]]):
    if len(_pb_list) < min_num_probe:
        print(f"Gene: {_gene} has probes less than {min_num_probe}, skip")
    else:
        kept_probe_dic[_gene] = sorted(_pb_list, key=lambda p:int(p["position"]))[:min(len(_pb_list), max_num_probe)]
print("Number of genes kept:", len(kept_probe_dic))
# save 
pickle.dump(kept_probe_dic, open(os.path.join(library_folder, 'filtered_candidate_dict.pkl'), 'wb'))<jupyter_output><empty_output><jupyter_text>save to fasta file<jupyter_code># biopython for SeqRecord
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# blast
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio.Blast import NCBIXML<jupyter_output><empty_output><jupyter_text>4. Append barcode info<jupyter_code># design gene_readout_dict
gene_readout_dict = {_gene:['u'+str(i)]*3 for i,_gene in enumerate(kept_probe_dic)}
# load primers and readouts
reload(library_tools.probes)
primers = library_tools.probes.load_primers([4,3])
unique_readouts = library_tools.probes.load_readouts(len(gene_readout_dict), 'NDB_new', _num_colors=3, _start_id=29)
readout_dict = {'u':unique_readouts}
kept_probe_dic['APP']
# Assemble probes
reload(library_tools.probes)
cand_probes, readout_summary = library_tools.probes.Assemble_probes(library_folder, kept_probe_dic, gene_readout_dict,
                                                                    readout_dict, primers, rc_targets=True, save=True)<jupyter_output>- Assemble probes by given target sequences, readouts and primers.
-- included readout types: ['u']
--- assemblying 80 probes in region: U2AF1L5
--- assemblying 80 probes in region: KCNE1B
--- assemblying 76 probes in region: CYP4F29P
--- assemblying 80 probes in region: ANKRD20A11P
--- assemblying 64 probes in region: HSPA13
--- assemblying 80 probes in region: NRIP1
--- assemblying 80 probes in region: USP25
--- assemblying 80 probes in region: CXADR
--- assemblying 80 probes in region: BTG3
--- assemblying 80 probes in region: C21orf91
--- assemblying 80 probes in region: NCAM2
--- assemblying 80 probes in region: MRPL39
--- assemblying 80 probes in region: JAM2
--- assemblying 80 probes in region: GABPA
--- assemblying 51 probes in region: ATP5PF
--- assemblying 80 probes in region: APP
--- assemblying 73 probes in region: ADAMTS1
--- assemblying 80 probes in region: ADAMTS5
--- assemblying 80 probes in region: N6AMT1
--- assemblying 80 probes in region: LTN1
--- assemblying 80 pro[...]<jupyter_text>5. Check quality<jupyter_code># biopython for SeqRecord
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# blast
from Bio.Blast.Applications import NcbiblastnCommandline
from Bio.Blast import NCBIXML 
# folder for this library
# candidate full-length probe filename
candidate_full_name = 'candidate_probes.fasta'
# load full probes
full_records = []
with open(os.path.join(library_folder, candidate_full_name), 'r') as handle:
    for record in SeqIO.parse(handle, "fasta"):
        full_records.append(record)

print(f"Total probe loaded: {len(full_records)}")
import ImageAnalysis3.library_tools.quality_check as qc
reload(qc)

primer_check = qc._check_primer_usage(full_records, primers[0], primers[1])
print(f"primer_check:{primer_check}")
reg_size_dic, len_check = qc._check_region_size(full_records, min_size=15)
print(f"len_check:{len_check}")
reg_readout_dic, reg2readout_check = qc._check_region_to_readouts(full_records, readout_dict)
print(f"reg2readout_check:{reg2readout_check}")
readout_reg_dic, readout2reg_check = qc._check_readout_to_region(reg_readout_dic, full_records, readout_dict, target_len=30)
print(f"readout2reg_check:{readout2reg_check}")
int_map = qc._construct_internal_map(full_records, library_folder)
readout_count_dic, readout_count_check = qc._check_readout_in_probes(readout_reg_dic, reg_size_dic, int_map, readout_dict)
print(f"readout_count_check:{readout_count_check}")
kept_records, removed_count = qc._check_between_probes(full_records, int_map)
# save kept records
with open(os.path.join(library_folder, 'filtered_full_probes.fasta'), 'w') as output_handle:
    SeqIO.write(kept_records, output_handle, "fasta")<jupyter_output>-- Checking primer usage, total probes: 7430
primer_check:True
gene: U2AF1L5 -> 80
gene: KCNE1B -> 80
gene: CYP4F29P -> 76
gene: ANKRD20A11P -> 80
gene: HSPA13 -> 64
gene: NRIP1 -> 80
gene: USP25 -> 80
gene: CXADR -> 80
gene: BTG3 -> 80
gene: C21orf91 -> 80
gene: NCAM2 -> 80
gene: MRPL39 -> 80
gene: JAM2 -> 80
gene: GABPA -> 80
gene: ATP5PF -> 51
gene: APP -> 80
gene: ADAMTS1 -> 73
gene: ADAMTS5 -> 80
gene: N6AMT1 -> 80
gene: LTN1 -> 80
gene: USP16 -> 80
gene: CCT8 -> 77
gene: BACH1 -> 80
gene: TIAM1 -> 80
gene: SOD1 -> 80
gene: SCAF4 -> 80
gene: HUNK -> 80
gene: MIS18A -> 80
gene: URB1 -> 80
gene: EVA1C -> 80
gene: SYNJ1 -> 80
gene: PAXBP1 -> 80
gene: IFNAR2 -> 80
gene: IL10RB -> 80
gene: IFNAR1 -> 80
gene: IFNGR2 -> 80
gene: SON -> 80
gene: DONSON -> 75
gene: ITSN1 -> 80
gene: CRYZL1 -> 80
gene: MRPS6 -> 80
gene: RCAN1 -> 50
gene: RUNX1 -> 80
gene: SETD4 -> 54
gene: CBR3 -> 76
gene: MORC3 -> 80
gene: CHAF1B -> 80
gene: SIM2 -> 80
gene: HLCS -> 80
gene: TTC3 -> 80
gene: DSCR9 -> 74
ge[...]<jupyter_text>blast<jupyter_code>qc.Blast_probes(kept_records, library_folder)
kept_pbs, blast_keep_dic, hard_count_list, soft_count_list = qc.Screening_Probes_by_Blast(library_folder, 60,
                                                                                          hard_thres=29)
kept_pbs
final_probe_folder = os.path.join(library_folder, 'final_probes')
if not os.path.exists(final_probe_folder):
    os.makedirs(final_probe_folder)
reload(qc)

primer_check = qc._check_primer_usage(kept_pbs, primers[0], primers[1])
print(f"primer_check:{primer_check}")
reg_size_dic, len_check = qc._check_region_size(kept_pbs, min_size=15)
print(f"len_check:{len_check}")
reg_readout_dic, reg2readout_check = qc._check_region_to_readouts(kept_pbs, readout_dict)
print(f"reg2readout_check:{reg2readout_check}")
readout_reg_dic, readout2reg_check = qc._check_readout_to_region(reg_readout_dic, kept_pbs, readout_dict, target_len=30)
print(f"readout2reg_check:{readout2reg_check}")
int_map = qc._construct_internal_map(kept_pbs, library_folder)
readout_count_dic, readout_count_check = qc._check_readout_in_probes(readout_reg_dic, reg_size_dic, int_map, readout_dict)
print(f"readout_count_check:{readout_count_check}")
kept_records, removed_count = qc._check_between_probes(kept_pbs, int_map)
# save kept records
with open(os.path.join(final_probe_folder, 'extra_filtered_full_probes.fasta'), 'w') as output_handle:
    SeqIO.write(kept_records, output_handle, "fasta")
len(kept_records)
str(kept_records[0].seq)[40:70]<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11122994652406418, 'fraction_numerical': 0.0351968886728245, 'mean_word_length': 4.817873303167421, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 79, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '513065', 'license': 'MIT', 'n_tokens_mistral': 4083, 'n_tokens_neox': 3681, 'n_words': 850, 'path': 'Library_design/CTP-07/CTP-07_Chr21_intronic_RNA_newNDB.ipynb', 'repo_name': 'shiwei23/Chromatin_Analysis_Scripts', 'size': '10287'}