<jupyter_start><jupyter_text>Plotting ContactmapsUsing these functions, contact maps can be constructed from the pairs files to visualize your transcripts/isoforms of interest. First, we visualize the differences among the following methods (KARR-seq, PARIS and RIC-seq) by looking at the Telomerase RNA component, TERC, transcript. By proving a dot-bracket notation file, we can overlap secondary structures onto the contact map. Here, we use the known secondary structure obtained from [Rfam](https://rfam.org/family/RF00024) to illustrate.<jupyter_code>data_KARRseq = ["data/inner/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HEK293T_M13_B01.dedup.pairs.gz"]
data_PARIS = ["data/inner/hrefseq/MAPQ1_SPAN0/PARIS-HEK293T_B97_B01.nodedup.pairs.gz"]
data_RIC = ["data/inner/hrefseq/MAPQ1_SPAN0/RIC-HeLa-rRNA-_M94_B01.dedup.pairs.gz"]

fig, ax = plt.subplots(1, 3, figsize=(18, 8))

# KARR-seq (Left)
plot_contactmap(data_KARRseq, p_chrom1="NR_001566", p_resolution=10, p_max=3, p_figsize=(6, 6),
                p_structure=False, p_genome="hg19", p_display="full", p_balance=False,
                p_dot_size=0.5, f_dot="data/sec_struct/TERC.dot",
                f_size="data/reference/hg19/hg19_refseq_wo_version.sizes", p_ax=(fig, ax[0]), p_verbose=False)

# PARIS (Center)
plot_contactmap(data_PARIS, p_chrom1="NR_001566", p_resolution=10, p_max=3, p_figsize=(6, 6),
                p_structure=False, p_genome="hg19", p_display="full", p_balance=False,
                p_dot_size=0.5, f_dot="/mnt/chenga/projects/rnarna/analysis/TechComparison/TERC.dot",
                f_size="data/reference/hg19/hg19_refseq_wo_version.sizes", p_ax=(fig, ax[1]), p_verbose=False)

# RIC-seq (Right)
plot_contactmap(data_RIC, p_chrom1="NR_001566", p_resolution=10, p_max=3, p_figsize=(6, 6),
                p_structure=False, p_genome="hg19", p_display="full", p_balance=False,
                p_dot_size=0.5, f_dot="/mnt/chenga/projects/rnarna/analysis/TechComparison/TERC.dot",
                f_size="data/reference/hg19/hg19_refseq_wo_version.sizes", p_ax=(fig, ax[2]), p_verbose=False)<jupyter_output><empty_output><jupyter_text>Clustering chimeric reads and classifying them into loops, left stripes and right stripes To cluster the chimeric groups, we first define a similarity function by counting the amount of overlapd between two chimeric pairs. We can then construct the graph by building the similarity matrix, and then cluster the chimeric reads into chimeric groups by using a greedy modularity algorithm. We can then classify the chimeric groups into loops, left stripes and right stripes.<jupyter_code>from analysis.clustering.cluster_pairs import cluster_chimeric_pairs
data_txt = "data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2-refold_M13_R01.txt.gz"

dict_kwargs = {
    "p_iid": "NM_001404", "p_genome": None, 
    "p_long": 200, "p_span": 10, "p_overlap": 10, "p_min_reads": 3, 
    "p_min_interactions": 200, "p_max_interactions": 20000, 
    "f_sizes": None, "f_tsne": None, "f_pairs": data_txt, "f_output": None, "f_bed": None
}
results = cluster_chimeric_pairs(**dict_kwargs)<jupyter_output>NM_001404 None 200 10 10 3 200 20000 None None /mnt/chenga/projects/rnarna/data/star/pairs/se/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2-refold_M13_R01.txt.gz None None
13
[['NM_001404', 471, 540, 'NM_001404', 823, 890], ['NM_001404', 816, 898, 'NM_001404', 986, 1074], ['NM_001404', 275, 336, 'NM_001404', 478, 566], ['NM_001404', 1198, 1287, 'NM_001404', 1365, 1467], ['NM_001404', 458, 541, 'NM_001404', 1367, 1444], ['NM_001404', 926, 998, 'NM_001404', 1401, 1470], ['NM_001404', 847, 906, 'NM_001404', 1260, 1407], ['NM_001404', 409, 495, 'NM_001404', 1093, 1195], ['NM_001404', 171, 222, 'NM_001404', 942, 1071], ['NM_001404', 614, 691, 'NM_001404', 1187, 1236], ['NM_001404', 937, 976, 'NM_001404', 1176, 1296], ['NM_001404', 1017, 1195, 'NM_001404', 1185, 1219], ['NM_001404', 194, 287, 'NM_001404', 1373, 1417]]<jupyter_text>Plotting Contactmaps with ArcBandsTo visualize the chimeric groups, we can plot them as arcbands and align it above the heatmap. Here, we will compare an orthologue pair from mouse (Eef1g, mESC) and human (EEF1G, HepG2). Using a BLAST analysis, even though they have a 91.4% precentage identity, their general topologies are similar, suggesting that structure are conserved and plausibly functional. If the secondary structure is not available, you can use RNAfold to predict the ensemble base-pairing matrix by setting`p_structure=True`.Here we plot the mouse Eef1g ...<jupyter_code>data_mouse = ["data/pairs/mrefseq/MAPQ1_SPAN0/G1_kethoxal-mESC-refold_M13_R01.dedup.pairs.gz", 
              "data/pairs/mrefseq/MAPQ1_SPAN0/G1_kethoxal-mESC-refold_M13_R02.dedup.pairs.gz"]

f_cluster_mouse = "data/clustering/output/mrefseq/G1_kethoxal-mESC-refold_M13_C01.dedup.classification.txt"

plot_contactmap_w_arcbands(data_mouse, f_cluster_mouse, p_chrom1="NM_026007.4", p_resolution=25, p_max=3, p_figsize=(10, 10),
                           p_structure=True, p_genome=None, p_display="upper", p_balance=False,
                           f_fasta = "data/reference/mm10/mm10_refseq.fa",
                           f_size = "data/reference/mm10/mm10_refseq.sizes",
                           p_dot_size=1.0, p_multi=False, p_verbose=False)<jupyter_output><empty_output><jupyter_text>... and then for human EEF1G<jupyter_code>data_human = ["data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2-refold_M13_R01.dedup.pairs.gz", 
              "data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2-refold_M13_R02.dedup.pairs.gz"]
f_cluster_human = "data/clustering/output/hrefseq/G1_kethoxal-HepG2-refold_M13_C01.dedup.classification.txt"

plot_contactmap_w_arcbands(data_human, f_cluster_human, p_chrom1="NM_001404", p_resolution=25, p_max=3, p_figsize=(10, 10),
                           p_structure=True, p_genome="hg19", p_display="upper", p_balance=False,
                           f_fasta = "data/reference/hg19/hg19_refseq_wo_version.fasta",
                           f_size = "data/reference/hg19/hg19_refseq_wo_version.sizes",
                           p_dot_size=1.0, p_multi=False, p_verbose=False)<jupyter_output><empty_output><jupyter_text>Differential MapsTo identify the features separating the in-vitro and in-vivo contact maps, we can use the plot_diffmap function.<jupyter_code>data_vitro = ["data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2-refold_M13_R01.dedup.pairs.gz"] 
data_vivo = ["data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2_M13_R01.dedup.pairs.gz"]

fig, ax = plt.subplots(1, 2, figsize=(18, 8))

plot_contactmap(data_vitro, p_chrom1="NM_001085", p_resolution=25, p_max=5, p_figsize=(6, 6),
             p_structure=True, p_genome="hg19", p_display="full", p_balance=False,
             p_dot_size=1.0, 
             f_size="data/reference/hg19/hg19_refseq_wo_version.sizes", p_ax=(fig, ax[0]), p_verbose=False)

plot_contactmap(data_vivo, p_chrom1="NM_001085", p_resolution=25, p_max=5, p_figsize=(6, 6),
             p_structure=True, p_genome="hg19", p_display="full", p_balance=False,
             p_dot_size=1.0, 
             f_size="data/reference/hg19/hg19_refseq_wo_version.sizes", p_ax=(fig, ax[1]), p_verbose=False)
data_pairs = ["data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2-refold_M13_R01.dedup.pairs.gz", 
              "data/pairs/hrefseq/MAPQ1_SPAN0/G1_kethoxal-HepG2_M13_R01.dedup.pairs.gz"]

plot_diffmap(data_pairs, p_chrom1="NM_001085", p_resolution=25, p_max=5, p_figsize=(6, 6),
             p_structure=True, p_genome="hg19", p_display="upper", p_balance=False,
             p_dot_size=1.0, p_diff="norm_diff",
             f_size="data/reference/hg19/hg19_refseq_wo_version.sizes", p_verbose=True)<jupyter_output>['NM_001085:0-1628|NM_001085:0-1628']
Dim: (66, 66) | Total Counts: 2165
Dim: (66, 66) | Total Counts: 1635
1 1
(66, 66)
(66, 66)
{'fraction_non_alphanumeric': 0.12669567443050933, 'fraction_numerical': 0.09508574353724085, 'mean_word_length': 4.658942795076032, 'pattern_counts': {'":': 13, '<': 19, '<?xml version=': 0, '>': 19, 'https://': 1, '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': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '458227', 'license': 'MIT', 'n_tokens_mistral': 3570, 'n_tokens_neox': 2987, 'n_words': 555, 'path': 'example/NB_visualize.ipynb', 'repo_name': 'ouyang-lab/KARR-seq', 'size': '7814'}