# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python [conda env:.conda-ExoCarto] *
#     language: python
#     name: conda-env-.conda-ExoCarto-py
# ---

# # Imports & Modules

# +
# General imports
import os
import sys
import pandas as pd
import scipy
pd.options.mode.chained_assignment = None  # default='warn'
import numpy as np
from tqdm import tqdm
sys.path.append('../')
# Other imports

tqdm.pandas()
import yaml
import json

# Figures imports
import matplotlib
from matplotlib.lines import Line2D   
import matplotlib.patches as mpatches
import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.ticker as mtick
from statannot import add_stat_annotation
import matplotlib.font_manager as font_manager
import string


# Font settings
font_dirs = ['/home/weber/Fonts', ]
font_files = font_manager.findSystemFonts(fontpaths=font_dirs)
font_list = font_manager.createFontList(font_files)
font_manager.fontManager.ttflist.extend(font_list)

pvalues_cutoff = json.load(open('config/config_files.json'))['pvalues_cutoff']
pvalues_cutoff = [[e, str(e)] for e in pvalues_cutoff]

import matplotlib.font_manager as font_manager

from matplotlib import rcParams
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Arial']
rcParams['font.weight'] = 'light'
rcParams['axes.facecolor'] = 'white'
rcParams['font.size'] = 18
# -

# # Load files

# ## Homo sapiens

# ### Update with count & write file

# +
# H sapiens

hs_genes = pd.read_parquet("/gstock/GeneIso/V2/Genes.parquet")
hs_mrna = pd.read_parquet("/gstock/GeneIso/V2/mRNA.parquet")
hs_exons = pd.read_parquet("/gstock/GeneIso/V2/Exons.parquet")
hs_exons['ranges'] = hs_exons['Exon region start (bp)'].astype(str) + '-' + hs_exons['Exon region end (bp)'].astype(str)

hs_cds = pd.read_parquet("/gstock/GeneIso/V2/CDS.parquet")
hs_cds['ranges'] = hs_cds['Exon region start (bp)'].astype(str) + '-' + hs_cds['Exon region end (bp)'].astype(str)

# TODO > REMOVE
## Used to compute nb of completely different exon boundaries / gene => differ from exon events where overlap is possible
hs_cds = hs_cds.groupby('GeneID').progress_apply(overlap_apply)
hs_cds.to_parquet("/gstock/GeneIso/V2/CDS_update.parquet")
hs_exons = hs_exons.groupby('GeneID').progress_apply(overlap_apply)
hs_exons.to_parquet("/gstock/GeneIso/V2/Exons_update.parquet")

hs_five_UTR = pd.read_parquet("/gstock/GeneIso/V2/5_UTR.parquet")
hs_three_UTR = pd.read_parquet("/gstock/GeneIso/V2/3_UTR.parquet")

# TODO > REMOVE
hs_five_UTR['ranges']  = hs_five_UTR["5' UTR start"].astype(int).astype(str) + '-' +  hs_five_UTR["5' UTR end"].astype(int).astype(str)
hs_three_UTR['ranges'] = hs_three_UTR["3' UTR start"].astype(int).astype(str) + '-' + hs_three_UTR["3' UTR end"].astype(int).astype(str)
hs_five_UTR = hs_five_UTR.groupby('GeneID').progress_apply(overlap_apply)
hs_five_UTR.to_parquet("/gstock/GeneIso/V2/5_UTR_update.parquet")

hs_three_UTR = hs_three_UTR.groupby('GeneID').progress_apply(overlap_apply)
hs_three_UTR.to_parquet("/gstock/GeneIso/V2/3_UTR_update.parquet")

hs_introns = pd.read_parquet("/gstock/GeneIso/V2/Introns.parquet")
hs_introns.loc[hs_introns['Strand'] == 1, 'Ordinal_nb_inverted'] = hs_introns.loc[hs_introns['Strand'] == 1, 'Ordinal_nb_inverted'] + 1
# -

# ## Mus Musculus

# +
genes = pd.read_parquet("/gstock/GeneIso/V2/Genes_MM.parquet")
# mrna = pd.read_parquet("/gstock/GeneIso/V2/mRNA.parquet")
exons = pd.read_parquet("/gstock/GeneIso/V2/Exons_MM.parquet")
exons['ranges'] = exons['Exon region start (bp)'].astype(str) + '-' + exons['Exon region end (bp)'].astype(str)


cds = pd.read_parquet("/gstock/GeneIso/V2/CDS_MM.parquet")
cds['ranges'] = cds['Exon region start (bp)'].astype(str) + '-' + cds['Exon region end (bp)'].astype(str)

# TODO > REMOVE
cds = cds.groupby('GeneID').progress_apply(overlap_apply)
exons = exons.groupby('GeneID').progress_apply(overlap_apply)

cds.to_parquet("/gstock/GeneIso/V2/CDS_update_MM.parquet")
exons.to_parquet("/gstock/GeneIso/V2/Exons_update_MM.parquet")


five_UTR = pd.read_parquet("/gstock/GeneIso/V2/5_UTR_MM.parquet")
three_UTR = pd.read_parquet("/gstock/GeneIso/V2/3_UTR_MM.parquet")
introns = pd.read_parquet("/gstock/GeneIso/V2/Introns_MM.parquet")

# TODO > REMOVE
## Used to compute nb of completely different exon boundaries / gene => differ from exon events where overlap is possible

five_UTR['ranges'] = five_UTR['UTR_start'].astype(int).astype(str) + '-' + five_UTR['UTR_end'].astype(int).astype(str)
three_UTR['ranges'] = three_UTR['UTR_start'].astype(int).astype(str) + '-' + three_UTR['UTR_end'].astype(int).astype(str)
five_UTR = five_UTR.groupby('GeneID').progress_apply(overlap_apply)
three_UTR = three_UTR.groupby('GeneID').progress_apply(overlap_apply)

five_UTR.to_parquet("/gstock/GeneIso/V2/5_UTR_update_MM.parquet")
three_UTR.to_parquet("/gstock/GeneIso/V2/3_UTR_update_MM.parquet")




#TODO
# introns.loc[introns['Strand'] == 1, 'Ordinal_nb_inverted'] = introns.loc[introns['Strand'] == 1, 'Ordinal_nb_inverted'] + 1
# -

# # Descritive stats

# ## Genes

genes_new = genes[['Miso_siso', 'GeneID', 'Gene name', 'ensembl_gene_id', 'chromosome_name', 'start_position', 'end_position']].drop_duplicates()
genes_new['Gene_length'] = genes_new['end_position'] - genes_new['start_position']

describe = pd.concat(
    [
        hs_genes.groupby('Miso_siso')['Gene_length'].describe().round(0).astype(int),
        genes_new.groupby('Miso_siso')['Gene_length'].describe().round(0).astype(int),
    ],
    axis=0
)
describe.index = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe

# ## mRNA

mrna = genes[['Miso_siso', 'GeneID', 'Gene name', 'ensembl_gene_id', 'ensembl_transcript_id', 'chromosome_name', 'transcript_start', 'transcript_end', 'transcript_length']].drop_duplicates()

mrna_describe = pd.concat(
    [
        hs_mrna.groupby('Miso_siso')['Transcript length (including UTRs and CDS)'].describe().round(0).astype(int).T,
        mrna.groupby('Miso_siso')['transcript_length'].describe().round(0).astype(int).T,
    ],
    axis=1
)
mrna_describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
mrna_describe

# ## Exons

# ### Length

# +
hs_exons = hs_exons.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)', 'Length'])
hs_exons = hs_exons.loc[hs_exons['Length'] > 0]

exons = exons.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)', 'Length'])
exons = exons.loc[exons['Length'] > 0]

describe = pd.concat(
    [
        hs_exons.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
        exons.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
    ],
    axis=1
)
describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe
# -

# ### Nb

# +
hs_exons = hs_exons.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)', 'Length'])
hs_exons = hs_exons.loc[hs_exons['Length'] > 0]

exons = exons.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)', 'Length'])
exons = exons.loc[exons['Length'] > 0]

describe = pd.concat(
    [
        hs_exons.groupby(['Miso_siso', 'GeneID'])['Exon region start (bp)'].nunique().groupby('Miso_siso').describe().round(0).astype(int).T,
        exons.groupby(['Miso_siso', 'GeneID'])['Exon region start (bp)'].nunique().groupby('Miso_siso').describe().round(0).astype(int).T,
    ],
    axis=1
)
describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe
# -

# ## CDS

# +
hs_cds = hs_cds.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)', 'Length'])
hs_cds = hs_cds.loc[hs_cds['Length'] > 0]

cds = cds.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)', 'Length'])
cds = cds.loc[cds['Length'] > 0]

describe = pd.concat(
    [
        hs_cds.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
        cds.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
    ],
    axis=1
)
describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe
# -

# ## Introns length

# +
hs_introns = hs_introns.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Introns'])
hs_introns = hs_introns.loc[hs_introns['Length'] > 0]

introns = introns.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Introns'])
introns = introns.loc[introns['Length'] > 0]

describe = pd.concat(
    [
        hs_introns.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
        introns.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
    ],
    axis=1
)
describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe
# -

# ## Introns length according ordinal

# +
hs_introns = pd.read_parquet("/gstock/GeneIso/V2/Introns.parquet")
hs_introns.loc[hs_introns['Strand'] == 1, 'Ordinal_nb_inverted'] = hs_introns.loc[hs_introns['Strand'] == 1, 'Ordinal_nb_inverted'] + 1
hs_introns = hs_introns.loc[hs_introns['Ordinal_nb'] <= 5] 

introns = pd.read_parquet("/gstock/GeneIso/V2/Introns_MM.parquet")
introns.loc[introns['Strand'] == 1, 'Ordinal_nb_inverted'] = introns.loc[introns['Strand'] == 1, 'Ordinal_nb_inverted'] + 1
introns = introns.loc[introns['Ordinal_nb'] <= 5] 

pd.concat(
    [
        hs_introns.groupby(['Miso_siso', 'Ordinal_nb'])['Length'].describe().round(0).astype(int),
        introns.groupby(['Miso_siso', 'Ordinal_nb'])['Length'].describe().round(0).astype(int)

    ], axis=0
).T
# -

# ## 5' UTR length

# +
hs_five_UTR = hs_five_UTR.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)'])
hs_five_UTR = hs_five_UTR.loc[hs_five_UTR['Length'] > 0]

five_UTR = five_UTR.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)'])
five_UTR = five_UTR.loc[five_UTR['Length'] > 0]

describe = pd.concat(
    [
        hs_five_UTR.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
        five_UTR.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
    ],
    axis=1
)
describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe
# -

# ## 3' UTR length

# +
hs_three_UTR = hs_three_UTR.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)'])
hs_three_UTR = hs_three_UTR.loc[hs_three_UTR['Length'] > 0]

three_UTR = three_UTR.drop_duplicates(subset=['Miso_siso', 'GeneID', 'Exon region start (bp)', 'Exon region end (bp)'])
three_UTR = three_UTR.loc[three_UTR['Length'] > 0]

describe = pd.concat(
    [
        three_UTR.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
        three_UTR.groupby('Miso_siso')['Length'].describe().round(0).astype(int).T,
    ],
    axis=1
)
describe.columns = pd.MultiIndex.from_product([['HS', 'MM'], ['MISOG', 'SISOG']])
describe
# -

# # MISOG SISOG category comparison for HS & MM

# +
# Count nb of transcripts per MM gene
new_count = genes.groupby('ensembl_gene_id')['ensembl_transcript_id'].nunique().rename('transcript_count_mm').reset_index()
new_count.loc[new_count['transcript_count_mm'] > 1, 'Miso_siso_mm'] = 'Miso'
new_count.loc[new_count['transcript_count_mm'] == 1, 'Miso_siso_mm'] = 'Siso'

# Merge with previous DF and retrieve Miso Siso HS status
new_count = pd.merge(genes, new_count,on='ensembl_gene_id')[['GeneID', 'ensembl_gene_id', 'Miso_siso', 'Miso_siso_mm']].drop_duplicates()
# new_count.groupby('Miso_siso')['Miso_siso_mm'].value_counts()
new_count['Miso_siso_full'] = 'HS:' + new_count['Miso_siso'] + '-MM:' + new_count['Miso_siso_mm']
new_count
# -

# # MM ENCODE Expression

# +
encode_mm = pd.read_csv('/gstock/biolo_datasets/ENCODE/ENCODE_SRS_TPM_summary.tsv.gz', compression='gzip', sep='\t') 

# Ratio = % of 535 experiments of short read seq where ENST is expressed > 0.1 TPM
## Like GTEx, cutoff = expression > 0.1 TPM in >= 20% exp
encode_mm = encode_mm.loc[encode_mm['TPM_ratio'] >= 0.2]
encode_mm
# -

# ## Merge 

# +
# Merge MM genes with MM genes expressed 

new_genes_mm = pd.merge(
    genes,
    encode_mm.rename({'ENST': 'ensembl_transcript_id'}, axis=1),
    on='ensembl_transcript_id'
)
# -

# ## New count expressed

# +
# Same as before but to retieve expressed genes

new_count = new_genes_mm.groupby('ensembl_gene_id')['ensembl_transcript_id'].nunique().rename('transcript_count_mm').reset_index()
new_count.loc[new_count['transcript_count_mm'] > 1, 'Miso_siso_mm'] = 'Miso'
new_count.loc[new_count['transcript_count_mm'] == 1, 'Miso_siso_mm'] = 'Siso'
new_count = pd.merge(genes, new_count,on='ensembl_gene_id')[['GeneID', 'ensembl_gene_id', 'Miso_siso', 'Miso_siso_mm']].drop_duplicates()
# new_count.groupby('Miso_siso')['Miso_siso_mm'].value_counts()
new_count['Miso_siso_full'] = 'HS:' + new_count['Miso_siso'] + '-MM:' + new_count['Miso_siso_mm']
new_count_expr = new_count.copy()
# -

# Before ENCODE expr
new_count

# After ENCODE expr
new_count_expr

# +
# Merge to compare MM genes that change of status

merge_merge = pd.merge(
    new_count[['GeneID', 'ensembl_gene_id', 'Miso_siso_mm']],
    new_count_expr[['GeneID', 'ensembl_gene_id', 'Miso_siso_mm']].rename({'Miso_siso_mm' : 'Miso_siso_mm_expr'}, axis=1),
    on=['GeneID', 'ensembl_gene_id']
)

# Export to reload later
merge_merge.to_csv('/gstock/GeneIso/V2_with_MM_expression/misog_sisog_mm_expr.tsv.gz', compression='gzip', sep='\t', index=False)
merge_merge
# -

# Identify MM genes that change of status before and after ENCODE expr
changing_status_mm_genes = pd.merge(
    merge_merge.loc[(merge_merge['Miso_siso_mm'] == "Miso") & (merge_merge['Miso_siso_mm_expr'] == "Siso")],
    genes[['GeneID', 'Gene name']],
    on='GeneID'
).drop_duplicates()
changing_status_mm_genes

# ## Count of genes that change after expression in MM

# MM genes before & after expr
merge_merge.groupby(['Miso_siso_mm'])['Miso_siso_mm_expr'].value_counts()

# Final distribution for Miso Siso HS & Miso Siso MM after expr
new_count.groupby('Miso_siso')['Miso_siso_mm'].value_counts()

# +
# Prepare for export to process with stringdb command of laetitia

new_count_to_export = pd.merge(new_count, genes[['GeneID', 'Gene name', 'mmusculus_homolog_associated_gene_name']].drop_duplicates(), on='GeneID')

for condition in new_count_to_export.Miso_siso_full.unique().tolist():
    print(condition)
    print('/gstock/GeneIso/V2_with_MM_expression/Gene_lists_STRING/{}_HS_gene_names.txt'.format(condition.replace(':', '-')))
    new_count_to_export.loc[new_count_to_export['Miso_siso_full'] == condition, 'Gene name'].to_csv('/gstock/GeneIso/V2_with_MM_expression/Gene_lists_STRING/{}_HS_gene_names.txt'.format(condition.replace(':', '-')), index=False, header=False)
#     new_count_to_export.loc[new_count_to_export['Miso_siso_full'] == condition, 'mmusculus_homolog_associated_gene_name'].to_csv('/gstock/GeneIso/V2/Gene_lists_STRING/{}_MM_gene_names.txt'.format(condition.replace(':', '-')), index=False, header=False)

StarcoderdataJupyterScriptsDedupFiltered