Skip to content
Snippets Groups Projects
Verified Commit 05d7f527 authored by nfontrod's avatar nfontrod
Browse files

src/07_ribosome_density.py: changes to display only ribosome density for one treatment at a time

parent 6db81916
No related branches found
No related tags found
No related merge requests found
...@@ -23,8 +23,10 @@ class Config_file: ...@@ -23,8 +23,10 @@ class Config_file:
htseq_rnaseq = htseq_fold / "htseq_count_rna" htseq_rnaseq = htseq_fold / "htseq_count_rna"
htseq_ribo = htseq_fold / "htseq_count_ribo" htseq_ribo = htseq_fold / "htseq_count_ribo"
result = Path(__file__).parents[1] / "results" result = Path(__file__).parents[1] / "results"
tdd_braf = result / "TDD_analysis" / "correlation" / "TDD_BRAF.txt" tdd_braf_itdd = result / "TDD_analysis" / "correlation" / \
tdd_dmso = result / "TDD_analysis" / "correlation" / "TDD_DMSO.txt" "TDD_BRAF_BRAF_DOWN.txt"
tdd_dmso_dtdd = result / "TDD_analysis" / "correlation" / \
"TDD_DMSO_BRAF_UP.txt"
output = result / "TDD_analysis" / "density_figures" output = result / "TDD_analysis" / "density_figures"
...@@ -35,7 +37,7 @@ def load_htseqfiles(mfolder: Path) -> List[Path]: ...@@ -35,7 +37,7 @@ def load_htseqfiles(mfolder: Path) -> List[Path]:
:param mfolder: A folder containing htseq count file :param mfolder: A folder containing htseq count file
:return: The list of files in `mfolder` :return: The list of files in `mfolder`
>>> res = load_htseqfiles(Config_file.htseq_ribo) >>> res = lrgoad_htseqfiles(Config_file.htseq_ribo)
>>> len(res) >>> len(res)
6 6
>>> [x.stem for x in res[0:2]] >>> [x.stem for x in res[0:2]]
...@@ -203,10 +205,9 @@ def avg_replicate(df_density: pd.DataFrame) -> pd.DataFrame: ...@@ -203,10 +205,9 @@ def avg_replicate(df_density: pd.DataFrame) -> pd.DataFrame:
3 D DMSO CTRL 2.125000 3 D DMSO CTRL 2.125000
4 E DMSO CTRL 2.600000 4 E DMSO CTRL 2.600000
""" """
df = df_density.drop("replicate", axis=1) \ return df_density.drop(
.groupby(["gene", "treatment", "group"]) \ "replicate", axis=1).groupby(
.mean().reset_index() ["gene", "treatment", "group"]).mean().reset_index()
return df
def statistical_analysis(df_density: pd.DataFrame, output: Path, def statistical_analysis(df_density: pd.DataFrame, output: Path,
...@@ -225,7 +226,7 @@ def statistical_analysis(df_density: pd.DataFrame, output: Path, ...@@ -225,7 +226,7 @@ def statistical_analysis(df_density: pd.DataFrame, output: Path,
treat = df_density["treatment"].unique()[0] treat = df_density["treatment"].unique()[0]
outfig = output / f"diag_{treat}_{list_type}.pdf" outfig = output / f"diag_{treat}_{list_type}.pdf"
outpval = output / f"diag_pval_{treat}_{list_type}.txt" outpval = output / f"diag_pval_{treat}_{list_type}.txt"
stat_s = r(f""" stat_s = r("""
require("DHARMa") require("DHARMa")
require("glmmTMB") require("glmmTMB")
function(data, outfig, outpval){{ function(data, outfig, outpval){{
...@@ -257,32 +258,35 @@ def create_density_figure(df_density: pd.DataFrame, gene_file: Path, ...@@ -257,32 +258,35 @@ def create_density_figure(df_density: pd.DataFrame, gene_file: Path,
df_density dataframe df_density dataframe
:param output_file: File where the figure will be created :param output_file: File where the figure will be created
""" """
sns.set(context="talk") sns.set_theme(context="talk", font_scale=2, style="whitegrid")
df_mean = avg_replicate(df_density) df_mean = avg_replicate(df_density)
df_mean["log10_density"] = np.log(df_mean["density"] + 1) df_mean["log10_density"] = np.log(df_mean["density"] + 1)
g = sns.catplot(x="treatment", y="log10_density", hue="group", ctreatment = "DMSO" if "UP" in gene_file.stem else "BRAF"
data=df_mean, kind="violin", height=10, aspect=1.7, df_mean = df_mean[df_mean["treatment"] == ctreatment]
palette={"TDD_BRAF": "#a600ff", "CTRL": "dimgrey", pval = statistical_analysis(
"TDD_DMSO": "#ff8400"}, cut=True) df_density[df_density["treatment"] == ctreatment].copy(),
pdmso = statistical_analysis(
df_density[df_density["treatment"] == "DMSO"].copy(),
output_file.parent, gene_file.stem)
pbraf = statistical_analysis(
df_density[df_density["treatment"] == "BRAF"].copy(),
output_file.parent, gene_file.stem) output_file.parent, gene_file.stem)
pdmso = "%.2e" % pdmso g = sns.catplot(x="group", y="log10_density",
pbraf = "%.2e" % pbraf data=df_mean, kind="violin", height=10, aspect=1.7,
g.fig.suptitle(f"Ribosome density " palette={"TDD_BRAF_BRAF_DOWN": "grey", "CTRL": "white",
f"between {gene_file.stem} and CTRL" "TDD_DMSO_BRAF_UP": "grey"}, cut=True)
f"genes in BRAF and DMSO condition\n" pval = "%.2e" % pval
f"(p_treatment_BRAF = {pbraf} & p_treatment_DMSO = {pdmso}") title = f"Ribosome density " \
g.set(ylabel="log10(density +1)") f"between {gene_file.stem} and CTRL" \
f"genes in {ctreatment} condition\n" \
f"(pval = {pval})"
g.fig.subplots_adjust(top=0.98)
g.fig.suptitle(title, fontsize=10)
g.set(ylabel="log10(density +1)", xlabel="")
g.savefig(output_file) g.savefig(output_file)
def density_figure_maker(): def density_figure_maker():
"""
Create the violin indicating the ribosome density figure
"""
df = create_density_df() df = create_density_df()
list_files = [Config_file.tdd_braf, Config_file.tdd_dmso] list_files = [Config_file.tdd_braf_itdd, Config_file.tdd_dmso_dtdd]
for mfile in list_files: for mfile in list_files:
dfc = add_group_column(df.copy(), mfile) dfc = add_group_column(df.copy(), mfile)
Config_file.output.mkdir(exist_ok=True) Config_file.output.mkdir(exist_ok=True)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment