Skip to content
This repository has been archived by the owner on Dec 29, 2023. It is now read-only.

Commit

Permalink
Add antibiotics effects (#31)
Browse files Browse the repository at this point in the history
  • Loading branch information
Fumire committed Mar 20, 2023
1 parent 98ffc73 commit d1d3d33
Show file tree
Hide file tree
Showing 3 changed files with 99 additions and 33 deletions.
27 changes: 12 additions & 15 deletions jwlee230/Program/Makefile
Original file line number Diff line number Diff line change
Expand Up @@ -943,6 +943,9 @@ Output/Step44/everything.%.pdf: Python/step44.py Output/Step44/everything.%.DAT.
Output/Step44/singleton.%.pdf: Python/step44.py Output/Step44/singleton.%.DAT.tsv | build.log
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) 1> $@.stdout 2> $@.stderr

Output/Step44/singleton.%.Antibiotics.Mouth.pdf: Python/step44-6.py Output/Step44/singleton.%.Antibiotics.Mouth.DAT.tsv | build.log
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) 1> $@.stdout 2> $@.stderr

Output/Step44/singleton.%.box.pdf: Python/step44-4.py Output/Step44/singleton.%.tsv Output/Step44/singleton.%.DAT.tsv Output/Metadata/everything.tsv | build.log
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) 1> $@.stdout 2> $@.stderr

Expand Down Expand Up @@ -1589,23 +1592,17 @@ step64: step64-clustering step64-dendrogram
Output/Step65 Output/Step65_Updown:
mkdir -p $@

Output/Step65/everything.%.Mouth.Up.pdf: Python/step65.py Output/Step44/everything.%.Mouth.DAT.tsv Output/Step44/everything.%.EL.Mouth.DAT.tsv Output/Step44/everything.%.EF.Mouth.DAT.tsv Output/Step44/everything.%.LF.Mouth.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" --up 1> $@.stdout 2> $@.stderr

Output/Step65/everything.%.Mouth.Down.pdf: Python/step65.py Output/Step44/everything.%.Mouth.DAT.tsv Output/Step44/everything.%.EL.Mouth.DAT.tsv Output/Step44/everything.%.EF.Mouth.DAT.tsv Output/Step44/everything.%.LF.Mouth.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" --down 1> $@.stdout 2> $@.stderr

Output/Step65/everything.%.Cervix.Up.pdf: Python/step65.py Output/Step44/everything.%.Cervix.DAT.tsv Output/Step44/everything.%.EL.Cervix.DAT.tsv Output/Step44/everything.%.EF.Cervix.DAT.tsv Output/Step44/everything.%.LF.Cervix.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" --up 1> $@.stdout 2> $@.stderr
Output/Step65/everything.%.Mouth.pdf: Python/step65.py Output/Step44/everything.%.Mouth.DAT.tsv Output/Step44/everything.%.EL.Mouth.DAT.tsv Output/Step44/everything.%.EF.Mouth.DAT.tsv Output/Step44/everything.%.LF.Mouth.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" 1> $@.stdout 2> $@.stderr

Output/Step65/everything.%.Cervix.Down.pdf: Python/step65.py Output/Step44/everything.%.Cervix.DAT.tsv Output/Step44/everything.%.EL.Cervix.DAT.tsv Output/Step44/everything.%.EF.Cervix.DAT.tsv Output/Step44/everything.%.LF.Cervix.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" --down 1> $@.stdout 2> $@.stderr
Output/Step65/everything.%.Cervix.pdf: Python/step65.py Output/Step44/everything.%.Cervix.DAT.tsv Output/Step44/everything.%.EL.Cervix.DAT.tsv Output/Step44/everything.%.EF.Cervix.DAT.tsv Output/Step44/everything.%.LF.Cervix.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" 1> $@.stdout 2> $@.stderr

Output/Step65/everything.%.Vagina.Up.pdf: Python/step65.py Output/Step44/everything.%.Vagina.DAT.tsv Output/Step44/everything.%.EL.Vagina.DAT.tsv Output/Step44/everything.%.EF.Vagina.DAT.tsv Output/Step44/everything.%.LF.Vagina.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" --up 1> $@.stdout 2> $@.stderr
Output/Step65/everything.%.Vagina.pdf: Python/step65.py Output/Step44/everything.%.Vagina.DAT.tsv Output/Step44/everything.%.EL.Vagina.DAT.tsv Output/Step44/everything.%.EF.Vagina.DAT.tsv Output/Step44/everything.%.LF.Vagina.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" 1> $@.stdout 2> $@.stderr

Output/Step65/everything.%.Vagina.Down.pdf: Python/step65.py Output/Step44/everything.%.Vagina.DAT.tsv Output/Step44/everything.%.EL.Vagina.DAT.tsv Output/Step44/everything.%.EF.Vagina.DAT.tsv Output/Step44/everything.%.LF.Vagina.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" --down 1> $@.stdout 2> $@.stderr
Output/Step65/singleton.%.Antibiotics.Mouth.pdf: Python/step65.py Output/Step44/singleton.%.Mouth.DAT.tsv Output/Step44/singleton.%.Antibiotics.Mouth.DAT.tsv | build.log Output/Step65
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "PTB" "Antibiotics" 1> $@.stdout 2> $@.stderr

Output/Step65_Updown/everything.%.Mouth.tsv: Python/step65-1.py Output/Step44/everything.%.Mouth.DAT.tsv Output/Step44/everything.%.EL.Mouth.DAT.tsv Output/Step44/everything.%.EF.Mouth.DAT.tsv Output/Step44/everything.%.LF.Mouth.DAT.tsv | build.log Output/Step65_Updown
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Early PTB vs. Late PTB+Normal" "Early PTB vs. Late PTB" "Earlry PTB vs. Normal" "Late PTB vs. Normal" 1> $@.stdout 2> $@.stderr
Expand All @@ -1619,7 +1616,7 @@ Output/Step65_Updown/everything.%.Vagina.tsv: Python/step65-1.py Output/Step44/e
Output/Step65_Updown/singleton.%.Mouth.tsv: Python/step65-1.py Output/Step44/everything.%.Mouth.DAT.tsv Output/Step44/singleton.%.Mouth.DAT.tsv | build.log Output/Step65_Updown
$(DOCKER) $(PYTHON) $(addprefix /,$^ $@) --annotation "Twin" "Singleton" 1> $@.stdout 2> $@.stderr

step65-multiple: Output/Step65/everything.DADA2.homd.Mouth.Up.pdf Output/Step65/everything.DADA2.homd.Mouth.Down.pdf Output/Step65/everything.DADA2.homd.uncorrected.Mouth.Up.pdf Output/Step65/everything.DADA2.homd.uncorrected.Mouth.Down.pdf Output/Step65/everything.DADA2.homd.Cervix.Up.pdf Output/Step65/everything.DADA2.homd.Cervix.Down.pdf Output/Step65/everything.DADA2.homd.uncorrected.Cervix.Up.pdf Output/Step65/everything.DADA2.homd.uncorrected.Cervix.Down.pdf Output/Step65/everything.DADA2.homd.Vagina.Up.pdf Output/Step65/everything.DADA2.homd.Vagina.Down.pdf Output/Step65/everything.DADA2.homd.uncorrected.Vagina.Up.pdf Output/Step65/everything.DADA2.homd.uncorrected.Vagina.Down.pdf
step65-multiple: Output/Step65/everything.DADA2.homd.Mouth.pdf Output/Step65/everything.DADA2.homd.uncorrected.Mouth.pdf Output/Step65/everything.DADA2.homd.Cervix.pdf Output/Step65/everything.DADA2.homd.uncorrected.Cervix.pdf Output/Step65/everything.DADA2.homd.Vagina.pdf Output/Step65/everything.DADA2.homd.uncorrected.Vagina.pdf
step65-updown: Output/Step65_Updown/everything.DADA2.homd.Mouth.tsv Output/Step65_Updown/everything.DADA2.homd.uncorrected.Mouth.tsv Output/Step65_Updown/everything.DADA2.homd.Cervix.tsv Output/Step65_Updown/everything.DADA2.homd.uncorrected.Cervix.tsv Output/Step65_Updown/everything.DADA2.homd.Vagina.tsv Output/Step65_Updown/everything.DADA2.homd.uncorrected.Vagina.tsv
step65: step65-multiple step65-updown
.PHONY: step65-multiple step65-updown step65
Expand Down
79 changes: 79 additions & 0 deletions jwlee230/Program/Python/step44-6.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,79 @@
"""
step44-6.py: Volcano plot with differetially abundant taxa with blue/red
"""
import argparse
import warnings
import adjustText
import matplotlib
import matplotlib.pyplot
import numpy
import pandas
import step00

ratio_threshold = 2
p_threshold = 0.05


if __name__ == "__main__":
parser = argparse.ArgumentParser()

parser.add_argument("input", help="Input TSV file", type=str)
parser.add_argument("output", help="Output PDF file", type=str)

args = parser.parse_args()

if not args.input.endswith(".tsv"):
raise ValueError("INPUT must end with .TSV!!")
elif not args.output.endswith(".pdf"):
raise ValueError("Output file must end with .PDF!!")

input_data = pandas.read_csv(args.input, sep="\t", index_col=0)
input_data["-log10(p)"] = -1 * numpy.log10(input_data["padj"])
input_data["simple_name"] = list(map(step00.consistency_taxonomy, list(input_data.index)))
print(input_data)

ceil = numpy.ceil(max(numpy.absolute(input_data["log2FoldChange"])))
input_data.to_csv(args.output.replace(".pdf", ".list.tsv"), sep="\t")

matplotlib.use("Agg")
matplotlib.rcParams.update(step00.matplotlib_parameters)
warnings.filterwarnings("error")

fig, ax = matplotlib.pyplot.subplots(figsize=(24, 24))
texts = list()

down_results = input_data.loc[((input_data["log2FoldChange"] < numpy.log2(1 / ratio_threshold)) & (input_data["-log10(p)"] > (-1 * numpy.log10(p_threshold)))), :].sort_values("-log10(p)", ascending=False)
up_results = input_data.loc[((input_data["log2FoldChange"] > numpy.log2(ratio_threshold)) & (input_data["-log10(p)"] > (-1 * numpy.log10(p_threshold)))), :].sort_values("-log10(p)", ascending=False)
ns_results = input_data.loc[(((input_data["log2FoldChange"] < numpy.log2(ratio_threshold)) & (input_data["log2FoldChange"] > numpy.log2(1 / ratio_threshold))) | (input_data["-log10(p)"] < (-1 * numpy.log10(p_threshold)))), :]

print(up_results.sort_values("simple_name"))
print(down_results.sort_values("simple_name"))

matplotlib.pyplot.scatter(ns_results["log2FoldChange"], ns_results["-log10(p)"], s=400, c="gray", marker="o", edgecolors=None)
matplotlib.pyplot.scatter(up_results["log2FoldChange"], up_results["-log10(p)"], s=400, c="tab:red", marker="o", edgecolors=None)
matplotlib.pyplot.scatter(down_results["log2FoldChange"], down_results["-log10(p)"], s=400, c="tab:blue", marker="o", edgecolors=None)

for index, row in down_results.iterrows():
texts.append(matplotlib.pyplot.text(row["log2FoldChange"], row["-log10(p)"], step00.simplified_taxonomy(index), color="tab:blue", fontsize="xx-small"))

for index, row in up_results.iterrows():
texts.append(matplotlib.pyplot.text(row["log2FoldChange"], row["-log10(p)"], step00.simplified_taxonomy(index), color="tab:red", fontsize="xx-small"))

matplotlib.pyplot.xlabel("log2(PTB/Normal)")
matplotlib.pyplot.ylabel("-log10(adj. p)")
matplotlib.pyplot.axvline(numpy.log2(1 / ratio_threshold), color="k", linestyle="--", linewidth=5)
matplotlib.pyplot.axvline(numpy.log2(ratio_threshold), color="k", linestyle="--", linewidth=5)
matplotlib.pyplot.axhline(-1 * numpy.log10(p_threshold), color="k", linestyle="--", linewidth=5)
matplotlib.pyplot.xlim((-ceil, ceil))
matplotlib.pyplot.grid(True)
matplotlib.pyplot.title(f"Up: {up_results.shape[0]}; Down: {down_results.shape[0]}")
matplotlib.pyplot.tight_layout()

adjustText.adjust_text(texts, arrowprops=dict(arrowstyle="-", color="black", alpha=0.3), lim=step00.big, ax=ax)
matplotlib.pyplot.tight_layout()

fig.savefig(args.output)
matplotlib.pyplot.close(fig)

output_results = input_data.loc[((input_data["log2FoldChange"] < numpy.log2(1 / ratio_threshold)) | (input_data["log2FoldChange"] > numpy.log2(ratio_threshold))) & (input_data["-log10(p)"] > (-1 * numpy.log10(p_threshold)))].sort_values("log2FoldChange")
output_results.to_csv(args.output.replace(".pdf", ".selected.tsv"), sep="\t")
26 changes: 8 additions & 18 deletions jwlee230/Program/Python/step65.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,12 +2,13 @@
step65.py: Venn diagram for differentially abundant taxa of DESeq2
"""
import argparse
import typing
import matplotlib
import matplotlib.pyplot
import numpy
import pandas
import tqdm
import venn
import upsetplot
import step00

ratio_threshold = 2
Expand All @@ -21,10 +22,6 @@
parser.add_argument("output", help="Output PDF file", type=str)
parser.add_argument("--annotation", help="Annotated name(s)", type=str, nargs="+")

updown = parser.add_mutually_exclusive_group(required=True)
updown.add_argument("--up", action="store_true", default=False)
updown.add_argument("--down", action="store_true", default=False)

args = parser.parse_args()

if list(filter(lambda x: not x.endswith(".tsv"), args.input)):
Expand All @@ -34,28 +31,21 @@
elif len(args.input) != len(args.annotation):
raise ValueError("Length of annotation must match length of input!!")

input_dict = dict()
input_dict: typing.Dict[str, typing.Set[str]] = dict()
for file, annot in tqdm.tqdm(zip(args.input, args.annotation)):
input_data = pandas.read_csv(file, sep="\t", index_col=0)
if args.up:
input_data = input_data.loc[(input_data["log2FoldChange"] > numpy.log2(ratio_threshold)) & (input_data["padj"] < p_threshold)]
elif args.down:
input_data = input_data.loc[(input_data["log2FoldChange"] < numpy.log2(1 / ratio_threshold)) & (input_data["padj"] < p_threshold)]
else:
raise Exception("Something went wrong!!")
input_dict[annot] = set(input_data.index)
input_dict[f"{annot}-Up"] = set(input_data.loc[(input_data["log2FoldChange"] > numpy.log2(ratio_threshold)) & (input_data["padj"] < p_threshold)].index)
input_dict[f"{annot}-Down"] = set(input_data.loc[(input_data["log2FoldChange"] < numpy.log2(1 / ratio_threshold)) & (input_data["padj"] < p_threshold)].index)

print("Union:", len(set.union(*input_dict.values())))

matplotlib.use("Agg")
matplotlib.rcParams.update(step00.matplotlib_parameters)

fig, ax = matplotlib.pyplot.subplots(figsize=(10, 10))
fig = matplotlib.pyplot.figure(figsize=(56, 24))

if set.union(*input_dict.values()):
venn.venn(input_dict, fmt=step00.venn_format, ax=ax)

matplotlib.pyplot.tight_layout()
upsetplot.plot(upsetplot.from_contents(input_dict), fig=fig, show_counts="%d", show_percentages=True, element_size=None)

fig.savefig(args.output)
fig.savefig(args.output, bbox_inches="tight")
matplotlib.pyplot.close(fig)

0 comments on commit d1d3d33

Please sign in to comment.