--- title: "20-supplementary-files" author: "Kathleen Durkin" date: "2026-06-15" output: github_document: toc: true number_sections: true bookdown::html_document2: theme: cosmo toc: true toc_float: true number_sections: true code_folding: show code_download: true html_document: theme: cosmo toc: true toc_float: true number_sections: true code_folding: show code_download: true editor_options: markdown: wrap: 72 --- ```{r setup, include=FALSE} library(tidyverse) library(readr) library(stringr) library(purrr) ``` Code for any compiling/formating of supplementary figures/tables for the DDE manuscript. # Machinery.fasta reference table ```{r build-machinery-table} protein_names <- c( # DNA methylation "DNMT1" = "DNA (cytosine-5-)-methyltransferase 1", "DNMT3A" = "DNA (cytosine-5-)-methyltransferase 3 alpha", "DNMT3B" = "DNA (cytosine-5-)-methyltransferase 3 beta", "DNMT3L" = "DNA (cytosine-5-)-methyltransferase 3 like", # TET dioxygenases "TET1" = "Methylcytosine dioxygenase TET1", "TET2" = "Methylcytosine dioxygenase TET2", "TET3" = "Methylcytosine dioxygenase TET3", # Methyl-CpG binding domain "MBD1" = "Methyl-CpG-binding domain protein 1", "MBD2" = "Methyl-CpG-binding domain protein 2", "MBD3" = "Methyl-CpG-binding domain protein 3", "MBD4" = "Methyl-CpG-binding domain protein 4", "MBD5" = "Methyl-CpG-binding domain protein 5", "MBD6" = "Methyl-CpG-binding domain protein 6", # UHRF "UHRF1" = "E3 ubiquitin-protein ligase UHRF1", "UHRF2" = "E3 ubiquitin-protein ligase UHRF2", # ZBTB readers "ZBTB33" = "Zinc finger and BTB domain-containing protein 33 (Kaiso)", "ZBTB4" = "Zinc finger and BTB domain-containing protein 4", "ZBTB38" = "Zinc finger and BTB domain-containing protein 38", "ZFP57" = "Zinc finger protein 57", "KLF4" = "Krueppel-like factor 4", "EGR1" = "Early growth response protein 1", "WT1" = "Wilms tumor protein", "CTCF" = "CCCTC-binding factor", # HDACs "HDAC1" = "Histone deacetylase 1", "HDAC2" = "Histone deacetylase 2", "HDAC3" = "Histone deacetylase 3", "HDAC4" = "Histone deacetylase 4", "HDAC5" = "Histone deacetylase 5", "HDAC6" = "Histone deacetylase 6", "HDAC7" = "Histone deacetylase 7", "HDAC8" = "Histone deacetylase 8", "HDAC9" = "Histone deacetylase 9", "HDAC10" = "Histone deacetylase 10", "HDAC11" = "Histone deacetylase 11", # Sirtuins "SIRT1" = "NAD-dependent protein deacetylase sirtuin-1", "SIRT2" = "NAD-dependent protein deacetylase sirtuin-2", "SIRT3" = "NAD-dependent protein deacetylase sirtuin-3", "SIRT4" = "NAD-dependent protein lysine acylase sirtuin-4", "SIRT5" = "NAD-dependent protein lysine desuccinylase sirtuin-5", "SIRT6" = "NAD-dependent protein deacylase sirtuin-6", "SIRT7" = "NAD-dependent protein desuccinylase sirtuin-7", # KATs "KAT1-HAT1" = "Histone acetyltransferase type B catalytic subunit (HAT1/KAT1)", "KAT2A" = "Histone acetyltransferase KAT2A (GCN5)", "KAT2B" = "Histone acetyltransferase KAT2B (PCAF)", "KAT3A-CREBBP" = "CREB-binding protein (CREBBP/KAT3A)", "KAT3B-EP300" = "Histone acetyltransferase p300 (EP300/KAT3B)", "KAT4-TAF1" = "Transcription initiation factor TFIID subunit 1 (TAF1/KAT4)", "KAT5" = "Histone acetyltransferase KAT5 (Tip60)", "KAT6A" = "Histone acetyltransferase KAT6A (MOZ)", "KAT6B" = "Histone acetyltransferase KAT6B (MORF)", "KAT7" = "Histone acetyltransferase KAT7 (HBO1)", "KAT8" = "Histone acetyltransferase KAT8 (MOF)", "KAT9-ELP3" = "Elongator complex protein 3 (ELP3/KAT9)", "KAT12-GTF3C4" = "Transcription factor IIIC subunit 4 (GTF3C4/KAT12)", "KAT13A-NCOA1" = "Nuclear receptor coactivator 1 (NCOA1/KAT13A)", "KAT13B-NCOA3" = "Nuclear receptor coactivator 3 (NCOA3/KAT13B)", "KAT13C-NCOA2" = "Nuclear receptor coactivator 2 (NCOA2/KAT13C)", "KAT13D-CLOCK" = "Circadian locomoter output cycles protein kaput (CLOCK/KAT13D)", "KAT14" = "Histone acetyltransferase KAT14 (CSRP2BP)", # KMTs "KMT1C-EHMT2" = "Histone-lysine N-methyltransferase EHMT2 (G9a/KMT1C)", "KMT1D-EHMT1" = "Histone-lysine N-methyltransferase EHMT1 (GLP/KMT1D)", "KMT1E-SETDB1" = "Histone-lysine N-methyltransferase SETDB1 (KMT1E)", "KMT1F-SETDB2" = "Histone-lysine N-methyltransferase SETDB2 (KMT1F)", "SUV39H1-KMT1A" = "Histone-lysine N-methyltransferase SUV39H1 (KMT1A)", "SUV39H2-KMT1B" = "Histone-lysine N-methyltransferase SUV39H2 (KMT1B)", "KMT2A" = "Histone-lysine N-methyltransferase 2A (MLL1)", "KMT2B" = "Histone-lysine N-methyltransferase 2B (MLL4)", "KMT2C" = "Histone-lysine N-methyltransferase 2C (MLL3)", "KMT2D" = "Histone-lysine N-methyltransferase 2D (MLL2)", "KMT2E" = "Histone-lysine N-methyltransferase 2E (MLL5)", "KMT2F-SETD1A" = "Histone-lysine N-methyltransferase SETD1A (KMT2F)", "KMT2G-SETD1B" = "Histone-lysine N-methyltransferase SETD1B (KMT2G)", "KMT2H-ASH1L" = "Histone-lysine N-methyltransferase ASH1L (KMT2H)", "KMT3A-SETD2" = "Histone-lysine N-methyltransferase SETD2 (KMT3A)", "KMT3B-NSD1" = "Histone-lysine N-methyltransferase NSD1 (KMT3B)", "KMT3C-SMYD2" = "SET and MYND domain-containing protein 2 (SMYD2/KMT3C)", "KMT4-DOT1L" = "Histone-lysine N-methyltransferase DOT1L (KMT4)", "KMT5A\u2014SETD8" = "Histone-lysine N-methyltransferase SETD8 (KMT5A)", "KMT5B" = "Histone-lysine N-methyltransferase KMT5B (SUV420H1)", "KMT5C" = "Histone-lysine N-methyltransferase KMT5C (SUV420H2)", "KMT6-EZH2" = "Histone-lysine N-methyltransferase EZH2 (KMT6)", "KMT7-SETD7/9" = "Histone-lysine N-methyltransferase SETD7 (KMT7)", "KMT8-PRDM2" = "PR domain zinc finger protein 2 (PRDM2/KMT8)", "ASH2L" = "Set1/Ash2 histone methyltransferase complex subunit ASH2", "WHSC1-NSD2" = "Histone-lysine N-methyltransferase NSD2 (WHSC1/MMSET)", "WHSC1L-NSD3" = "Histone-lysine N-methyltransferase NSD3 (WHSC1L1)", "EZH1" = "Histone-lysine N-methyltransferase EZH1", "SETD3" = "Histone methyltransferase SETD3", "SETD4" = "SET domain-containing protein 4", "SETD5" = "SET domain-containing protein 5", "SETD6" = "N-lysine methyltransferase SETD6", "SETMAR" = "Histone-lysine N-methyltransferase SETMAR (Metnase)", "SMYD1" = "SET and MYND domain-containing protein 1", "SMYD3" = "SET and MYND domain-containing protein 3", "SMYD4" = "SET and MYND domain-containing protein 4", "SMYD5" = "SET and MYND domain-containing protein 5", # KDMs "KDM1A" = "Lysine-specific histone demethylase 1A (LSD1)", "KDM1B" = "Lysine-specific histone demethylase 1B (LSD2)", "KDM2A" = "Lysine-specific demethylase 2A", "KDM2B" = "Lysine-specific demethylase 2B", "KDM3A" = "Lysine-specific demethylase 3A", "KDM3B" = "Lysine-specific demethylase 3B", "KDM3C-JMJD1C" = "Lysine-specific demethylase JMJD1C (KDM3C)", "KDM4A" = "Lysine-specific demethylase 4A", "KDM4B" = "Lysine-specific demethylase 4B", "KDM4C" = "Lysine-specific demethylase 4C", "KDM4D" = "Lysine-specific demethylase 4D", "KDM5A" = "Lysine-specific demethylase 5A (JARID1A)", "KDM5B" = "Lysine-specific demethylase 5B (JARID1B)", "KDM5C" = "Lysine-specific demethylase 5C (JARID1C)", "KDM5D" = "Lysine-specific demethylase 5D (JARID1D)", "KDM6A" = "Lysine-specific demethylase 6A (UTX)", "KDM6B" = "Lysine-specific demethylase 6B (JMJD3)", "KDM7A" = "Lysine-specific demethylase 7A (JHDM1D/PHF8-related)", "KDM8" = "Lysine-specific demethylase 8 (JMJD5)", "UTY" = "Inactive histone demethylase UTY", "PHF2" = "PHD finger protein 2 (histone demethylase)", "PHF8" = "Histone lysine demethylase PHF8", "JMJD6" = "Bifunctional arginine demethylase and lysyl-hydroxylase JMJD6", "JMJD7" = "Lysine demethylase JMJD7", "JMJD8" = "JMJD8 protein", "JARID2" = "Protein Jumonji (JARID2)", "HR" = "Lysine demethylase and nuclear receptor corepressor HR", "RIOX1" = "Ribosomal oxygenase 1 (NO66)", "RIOX2" = "Ribosomal oxygenase 2 (MINA53)", "TYW5" = "tRNA wybutosine-synthesizing protein 5", # PRDM "PRDM1" = "PR domain zinc finger protein 1 (Blimp1)", "MECOM-PRDM3" = "EVI1/MDS1 and EVI1 complex locus protein (MECOM/PRDM3)", "PRDM4" = "PR domain zinc finger protein 4", "PRDM5" = "PR domain zinc finger protein 5", "PRDM6" = "PR domain zinc finger protein 6", "PRDM7" = "PR domain zinc finger protein 7", "PRDM8" = "PR domain zinc finger protein 8", "PRDM9" = "PR domain zinc finger protein 9", "PRDM10" = "PR domain zinc finger protein 10", "PRDM11" = "PR domain zinc finger protein 11", "PRDM12" = "PR domain zinc finger protein 12", "PRDM13" = "PR domain zinc finger protein 13", "PRDM14" = "PR domain zinc finger protein 14", "PRDM15" = "PR domain zinc finger protein 15", "PRDM16" = "PR domain zinc finger protein 16", # PRMTs "PRMT2" = "Protein arginine N-methyltransferase 2", "PRMT3" = "Protein arginine N-methyltransferase 3", "PRMT5" = "Protein arginine N-methyltransferase 5", "PRMT6" = "Protein arginine N-methyltransferase 6", "PRMT7" = "Protein arginine N-methyltransferase 7", "PRMT8" = "Protein arginine N-methyltransferase 8", "PRMT9" = "Protein arginine N-methyltransferase 9", # PARPs "PARP2" = "Poly [ADP-ribose] polymerase 2", "PARP3" = "Poly [ADP-ribose] polymerase 3", "PARP4" = "Poly [ADP-ribose] polymerase 4", "PARP5A-TNKS" = "Tankyrase-1 (TNKS/PARP5A)", "PARP5B-TNKS2" = "Tankyrase-2 (TNKS2/PARP5B)", "PARP6" = "Poly [ADP-ribose] polymerase 6", "PARP7-TIPARP" = "TCDD-inducible poly [ADP-ribose] polymerase (TIPARP/PARP7)", "PARP8" = "Poly [ADP-ribose] polymerase 8", "PARP9" = "Poly [ADP-ribose] polymerase 9 (BAL1)", "PARP10" = "Poly [ADP-ribose] polymerase 10", "PARP11" = "Poly [ADP-ribose] polymerase 11", "PARP12" = "Poly [ADP-ribose] polymerase 12", "PARP13-ZC3HAV1" = "Zinc finger CCCH-type antiviral protein 1 (ZAP/PARP13)", "PARP14" = "Poly [ADP-ribose] polymerase 14 (BAL2)", "PARP15" = "Poly [ADP-ribose] polymerase 15 (BAL3)", "PARP16" = "Poly [ADP-ribose] polymerase 16", "PARG" = "Poly(ADP-ribose) glycohydrolase", "ARH1-ADPRH" = "ADP-ribosylhydrolase ARH1", "ARH3-ADPRS" = "ADP-ribosylhydrolase ARH3", "OGT" = "UDP-N-acetylglucosamine--peptide N-acetylglucosaminyltransferase (OGT)", # SUMO "SUMO2" = "Small ubiquitin-related modifier 2", "SUMO3" = "Small ubiquitin-related modifier 3", "SENP1" = "Sentrin-specific protease 1", "SENP2" = "Sentrin-specific protease 2", "SENP3" = "Sentrin-specific protease 3", "SENP5" = "Sentrin-specific protease 5", "SENP6" = "Sentrin-specific protease 6", "SENP7" = "Sentrin-specific protease 7", "SENP8" = "Sentrin-specific protease 8 (NEDD8 deSUMOylation)", # Ubiquitin E1/E2/E3 "UBA1" = "Ubiquitin-activating enzyme E1", "UBA2" = "SUMO-activating enzyme subunit 2", "UBA3" = "NEDD8-activating enzyme E1 catalytic subunit", "UBA5" = "Ubiquitin-like modifier-activating enzyme 5", "UBA6" = "Ubiquitin-activating enzyme E1-like protein", "UBA7" = "Ubiquitin-activating enzyme E1-like protein 2 (UBE1L)", "UBE2A" = "Ubiquitin-conjugating enzyme E2 A", "UBE2D1" = "Ubiquitin-conjugating enzyme E2 D1", "UBE2D2" = "Ubiquitin-conjugating enzyme E2 D2", "UBE2D3" = "Ubiquitin-conjugating enzyme E2 D3", "UBE2V1" = "Ubiquitin-conjugating enzyme E2 variant 1", "UBE2Z" = "Ubiquitin-conjugating enzyme E2 Z", "UBE3A" = "Ubiquitin-protein ligase E3A", "UBE4A" = "Ubiquitin conjugation factor E4 A", "UBR2" = "E3 ubiquitin-protein ligase UBR2", "BRCA1" = "Breast cancer type 1 susceptibility protein (BRCA1)", "HUWE1" = "E3 ubiquitin-protein ligase HUWE1", "RNF8" = "E3 ubiquitin-protein ligase RNF8", "RNF20" = "E3 ubiquitin-protein ligase BRE1A (RNF20)", "RNF40" = "E3 ubiquitin-protein ligase BRE1B (RNF40)", "RNF168" = "E3 ubiquitin-protein ligase RNF168", "TRIM37" = "E3 ubiquitin-protein ligase TRIM37", "DTX3L" = "Deltex E3 ubiquitin ligase 3L (DTX3L)", "CUL4A" = "Cullin-4A", "DDB1" = "DNA damage-binding protein 1", "TOM1" = "Target of Myb protein 1 (TOM1)", # USPs "USP1" = "Ubiquitin carboxyl-terminal hydrolase 1", "USP2" = "Ubiquitin carboxyl-terminal hydrolase 2", "USP3" = "Ubiquitin carboxyl-terminal hydrolase 3", "USP4" = "Ubiquitin carboxyl-terminal hydrolase 4", "USP5" = "Ubiquitin carboxyl-terminal hydrolase 5", "USP6" = "Ubiquitin carboxyl-terminal hydrolase 6", "USP7" = "Ubiquitin carboxyl-terminal hydrolase 7 (HAUSP)", "USP8" = "Ubiquitin carboxyl-terminal hydrolase 8", "USP9Y" = "Ubiquitin carboxyl-terminal hydrolase 9, Y-linked", "USP10" = "Ubiquitin carboxyl-terminal hydrolase 10", "USP11" = "Ubiquitin carboxyl-terminal hydrolase 11", "USP12" = "Ubiquitin carboxyl-terminal hydrolase 12", "USP13" = "Ubiquitin carboxyl-terminal hydrolase 13", "USP14" = "Ubiquitin carboxyl-terminal hydrolase 14", "USP15" = "Ubiquitin carboxyl-terminal hydrolase 15", "USP16" = "Ubiquitin carboxyl-terminal hydrolase 16", "USP18" = "Ubiquitin carboxyl-terminal hydrolase 18", "USP19" = "Ubiquitin carboxyl-terminal hydrolase 19", "USP20" = "Ubiquitin carboxyl-terminal hydrolase 20", "USP21" = "Ubiquitin carboxyl-terminal hydrolase 21", "USP22" = "Ubiquitin carboxyl-terminal hydrolase 22", "USP24" = "Ubiquitin carboxyl-terminal hydrolase 24", "USP25" = "Ubiquitin carboxyl-terminal hydrolase 25", "USP26" = "Ubiquitin carboxyl-terminal hydrolase 26", "USP27X" = "Ubiquitin carboxyl-terminal hydrolase 27, X-linked", "USP28" = "Ubiquitin carboxyl-terminal hydrolase 28", "USP29" = "Ubiquitin carboxyl-terminal hydrolase 29", "USP30" = "Ubiquitin carboxyl-terminal hydrolase 30", "USP31" = "Ubiquitin carboxyl-terminal hydrolase 31", "USP32" = "Ubiquitin carboxyl-terminal hydrolase 32", "USP33" = "Ubiquitin carboxyl-terminal hydrolase 33", "USP34" = "Ubiquitin carboxyl-terminal hydrolase 34", "USP35" = "Ubiquitin carboxyl-terminal hydrolase 35", "USP36" = "Ubiquitin carboxyl-terminal hydrolase 36", "USP37" = "Ubiquitin carboxyl-terminal hydrolase 37", "USP38" = "Ubiquitin carboxyl-terminal hydrolase 38", "USP39" = "Ubiquitin carboxyl-terminal hydrolase 39", "USP40" = "Ubiquitin carboxyl-terminal hydrolase 40", "USP41" = "Ubiquitin carboxyl-terminal hydrolase 41", "USP42" = "Ubiquitin carboxyl-terminal hydrolase 42", "USP43" = "Ubiquitin carboxyl-terminal hydrolase 43", "USP44" = "Ubiquitin carboxyl-terminal hydrolase 44", "USP45" = "Ubiquitin carboxyl-terminal hydrolase 45", "USP46" = "Ubiquitin carboxyl-terminal hydrolase 46", "USP47" = "Ubiquitin carboxyl-terminal hydrolase 47", "USP48" = "Ubiquitin carboxyl-terminal hydrolase 48", "USP49" = "Ubiquitin carboxyl-terminal hydrolase 49", "USP50" = "Ubiquitin carboxyl-terminal hydrolase 50", "USP51" = "Ubiquitin carboxyl-terminal hydrolase 51", "USP53" = "Ubiquitin carboxyl-terminal hydrolase 53", "USP54" = "Ubiquitin carboxyl-terminal hydrolase 54", "PAN2" = "Poly(A)-specific ribonuclease subunit PAN2 (USP52)", # UCH / OTU / other DUBs "UCHL1" = "Ubiquitin carboxyl-terminal hydrolase isozyme L1", "UCHL3" = "Ubiquitin carboxyl-terminal hydrolase isozyme L3", "UCHL4" = "Ubiquitin carboxyl-terminal hydrolase isozyme L4", "UCHL5" = "Ubiquitin carboxyl-terminal hydrolase isozyme L5 (UCH37)", "OTUD1" = "OTU domain-containing protein 1", "OTUD3" = "OTU domain-containing protein 3", "OTUD4" = "OTU domain-containing protein 4", "OTUD5" = "OTU domain-containing protein 5 (DUBA)", "OTUD6A" = "OTU domain-containing protein 6A", "OTUD6B" = "OTU domain-containing protein 6B", "OTUD7A" = "OTU domain-containing protein 7A (CEZANNE2)", "OTUD7B" = "OTU domain-containing protein 7B (CEZANNE)", "BAP1" = "Ubiquitin carboxyl-terminal hydrolase BAP1", "BRCC3" = "Lys-63-specific deubiquitinase BRCC3", "MYSM1" = "Histone H2A deubiquitinase MYSM1", "ATXN3" = "Ataxin-3 (Josephin domain deubiquitinase)", "ATXN3L" = "Ataxin-3-like protein", "JOSD1" = "Josephin domain-containing protein 1", "JOSD2" = "Josephin domain-containing protein 2", "PSMD14" = "26S proteasome non-ATPase regulatory subunit 14 (POH1/Rpn11)", "PSMD7" = "26S proteasome non-ATPase regulatory subunit 7 (Rpn8)", "COPS5" = "COP9 signalosome complex subunit 5 (JAB1/CSN5)", "COPS6" = "COP9 signalosome complex subunit 6 (CSN6)", "STAMBPL1"= "Stam-binding protein-like 1 (AMSH-LP)", "ZRANB1" = "Zinc finger RANBP2-type and RING finger-containing protein 1 (TRABID)", "VCPIP1" = "Deubiquitinating protein VCPIP1", "EIF3H" = "Eukaryotic translation initiation factor 3 subunit H", "EIF3F" = "Eukaryotic translation initiation factor 3 subunit F", "MPND" = "MPN domain-containing protein", "PRPF8" = "Pre-mRNA-processing-splicing factor 8 (PRPF8/SNRNP220)", "ARID1A" = "AT-rich interactive domain-containing protein 1A", "ARID1B" = "AT-rich interactive domain-containing protein 1B", # ADAR "ADAR" = "Double-stranded RNA-specific adenosine deaminase (ADAR1)", "ADARB1" = "Double-stranded RNA-specific editase B1 (ADAR2)", "ADARB2" = "Double-stranded RNA-specific editase B2 (ADAR3)", # m6A writers / readers / erasers "METTL14" = "N6-adenosine-methyltransferase non-catalytic subunit (METTL14)", "METTL16" = "RNA N6-methyladenosine methyltransferase METTL16", "METTL1" = "tRNA (guanine(46)-N(7))-methyltransferase METTL1", "METTL2A" = "tRNA (cytosine(34)-N(4))-methyltransferase METTL2A", "METTL2B" = "tRNA (cytosine(34)-N(4))-methyltransferase METTL2B", "METTL4" = "N6-adenosine-methyltransferase METTL4", "METTL5" = "rRNA adenine N(6)-methyltransferase (METTL5)", "METTL6" = "Methyltransferase-like protein 6", "METTL8" = "Methyltransferase-like protein 8", "WTAP" = "pre-mRNA-splicing regulator WTAP (m6A writer complex)", "VIRMA" = "Vir-like m6A methyltransferase associated (KIAA1429)", "CBLL1" = "Hakai (CBLL1); E3 ubiquitin ligase, m6A writer complex", "ZC3H13" = "Zinc finger CCCH domain-containing protein 13 (m6A writer complex)", "RBM15" = "RNA-binding protein 15", "RBM15B" = "RNA-binding protein 15B", "PCIF1" = "mRNA cap-specific m6Am methyltransferase PCIF1", "FTO" = "Alpha-ketoglutarate-dependent dioxygenase FTO (m6A eraser)", "ALKBH5" = "mRNA demethylase ALKBH5 (m6A eraser)", "ALKBH1" = "tRNA N(1)-methyl adenine demethylase ALKBH1", "ALKBH3" = "tRNA demethylase ALKBH3", "YTHDC1" = "YTH domain-containing protein 1 (m6A nuclear reader)", "YTHDC2" = "YTH domain-containing protein 2 (m6A reader/helicase)", "YTHDF1" = "YTH domain-containing family protein 1 (m6A cytoplasmic reader)", "YTHDF2" = "YTH domain-containing family protein 2 (m6A decay reader)", "YTHDF3" = "YTH domain-containing family protein 3 (m6A reader)", "HNRNPC" = "Heterogeneous nuclear ribonucleoprotein C", "HNRNPA2B1" = "Heterogeneous nuclear ribonucleoproteins A2/B1", "RBMX" = "RNA-binding motif protein, X-linked (HNRNPG)", "IGF2BP1" = "Insulin-like growth factor 2 mRNA-binding protein 1", "IGF2BP2" = "Insulin-like growth factor 2 mRNA-binding protein 2", "IGF2BP3" = "Insulin-like growth factor 2 mRNA-binding protein 3", "FMR1" = "Fragile X mental retardation protein 1", "PRRC2A" = "Proline-rich coiled-coil-containing protein 2A (m6A reader)", "ELF3" = "ETS-related transcription factor ELF3", # RNA methylation (non-m6A) "NSUN2" = "tRNA (cytosine(34)-C(5))-methyltransferase NSUN2", "NSUN3" = "Methyltransferase-like protein NSUN3 (mitochondrial tRNA)", "NSUN4" = "Bifunctional methyltransferase/rRNA maturation factor NSUN4", "NSUN5" = "rRNA methyltransferase NSUN5", "NSUN6" = "tRNA methyltransferase NSUN6", "TRDMT1" = "tRNA (cytosine-5-)-methyltransferase TRDMT1 (DNMT2)", "NOP2" = "Ribosome biogenesis protein NOP2 (NSUN1)", "DKC1" = "H/ACA ribonucleoprotein complex subunit 4 (Dyskerin/DKC1)", "TRMT1" = "tRNA (guanine(26)-N(2))-methyltransferase TRMT1", "TRMT61A" = "tRNA (adenine(58)-N(1))-methyltransferase TRMT61A", "TRMT61B" = "rRNA (adenine(58)-N(1))-methyltransferase TRMT61B", "TRMT10A" = "tRNA (guanine(9)-N(1))-methyltransferase TRMT10A", "TRMT10C" = "Mitochondrial tRNA methyltransferase TRMT10C", "TRMT112" = "tRNA methyltransferase subunit TRMT112", "ALYREF" = "THO complex subunit ALYREF (Aly/REF export factor)", "ZCCHC4" = "Zinc finger CCHC domain-containing protein 4 (rRNA m6A writer)", "NAT10" = "N-acetyltransferase 10 (RNA cytidine acetyltransferase)", "DUS1L" = "tRNA-dihydrouridine synthase 1-like", # Pseudouridine synthases "PUSL1" = "Pseudouridine synthase 1-like", "PUS1" = "tRNA pseudouridine synthase A (PUS1)", "PUS3" = "tRNA pseudouridine synthase C (PUS3)", "PUS7" = "tRNA pseudouridine synthase G (PUS7)", "PUS7L" = "tRNA pseudouridine synthase 7-like", "PUS10" = "tRNA pseudouridine synthase 10 (RPUSD2-related)", "RPUSD1" = "RNA pseudouridylate synthase domain-containing 1", "RPUSD2" = "RNA pseudouridylate synthase domain-containing 2", "RPUSD3" = "RNA pseudouridylate synthase domain-containing 3", "RPUSD4" = "RNA pseudouridylate synthase domain-containing 4", "TRUB1" = "Probable tRNA pseudouridine synthase 1 (TRUB1)", "TRUB2" = "Probable tRNA pseudouridine synthase 2 (TRUB2)", # NATs / acetyltransferases "NAA10" = "N-alpha-acetyltransferase 10 (NatA catalytic subunit)", "NAA50" = "N-alpha-acetyltransferase 50 (NatE catalytic subunit)", "NAA60" = "N-alpha-acetyltransferase 60 (NatF catalytic subunit)", "NAT8" = "N-acetyltransferase 8 (camello-like 1)", "NAT8B" = "N-acetyltransferase 8B (camello-like 2)", "NAT8L" = "N-acetyltransferase 8-like", "ATAT1" = "Alpha-tubulin acetyltransferase 1", "GTF3C2" = "Transcription factor IIIC 110 kDa subunit (TFIIIC110/GTF3C2)", "GTF3C1" = "Transcription factor IIIC 220 kDa subunit (TFIIIC220/GTF3C1)", "DLAT" = "Dihydrolipoyllysine-residue acetyltransferase (PDC-E2)", "ATF2" = "Cyclic AMP-dependent transcription factor ATF-2", "MCM3AP" = "Germinal-center associated nuclear protein (GANP/MCM3AP)", "MAPT" = "Microtubule-associated protein tau (MAPT)", "BRD4" = "Bromodomain-containing protein 4", "ESCO1" = "N-acetyltransferase ESCO1 (establishment of cohesion)", "ESCO2" = "N-acetyltransferase ESCO2", "BLOC1S1" = "Biogenesis of lysosome-related organelles complex 1 subunit 1", "ACAT1" = "Acetyl-CoA acetyltransferase 1 (mitochondrial)", "ACAT2" = "Acetyl-CoA acetyltransferase 2 (cytoplasmic)", # Histone variants / linker histones "H1-0" = "Histone H1.0", "H1-1" = "Histone H1.1 (HIST1H1A)", "H1-2" = "Histone H1.2 (HIST1H1C)", "H1-3" = "Histone H1.3 (HIST1H1D)", "H1-4" = "Histone H1.4 (HIST1H1E)", "H1-5" = "Histone H1.5 (HIST1H1B)", "H1-6" = "Histone H1t", "H1-7" = "Histone H1oo (oocyte-specific)", "H1-8" = "Histone H1.8 (H1foo)", "H1-10" = "Histone H1.10", "H2AX" = "Histone H2A.x", "H2AZ1" = "Histone H2A.z.1", "H2AZ2" = "Histone H2A.z.2", "H2AW" = "Histone H2A.W", "H2AP" = "Histone H2A.P", "H2AJ" = "Histone H2A.J", "H2AB1" = "Histone H2A-Bbd type 1 (H2A.Bbd)", "H2AB2" = "Histone H2A-Bbd type 2", "H2AB3" = "Histone H2A-Bbd type 3", "MACROH2A1" = "Core histone macro-H2A.1 (macroH2A1/H2AFY)", # Chromatin regulators "BAZ1B" = "Bromodomain adjacent to zinc finger domain protein 1B (WSTF)", "BAZ2B" = "Bromodomain adjacent to zinc finger domain protein 2B", "HIRA" = "Histone cell cycle regulator (HIRA)", "WDR5" = "WD repeat-containing protein 5 (H3K4 methyltransferase complex subunit)", # EYA phosphatases "EYA1" = "Eyes absent homolog 1 (H2AX phosphatase)", "EYA2" = "Eyes absent homolog 2", "EYA3" = "Eyes absent homolog 3", "EYA4" = "Eyes absent homolog 4", # Protein phosphatases "PPM1D" = "Protein phosphatase 1D (WIP1)", "PPP6C" = "Serine/threonine-protein phosphatase 6 catalytic subunit", "PPP4C" = "Serine/threonine-protein phosphatase 4 catalytic subunit", "CDCA2" = "Cell division cycle associated 2 (Repo-Man; PP1 targeting)", "PPP1CA" = "Serine/threonine-protein phosphatase PP1-alpha catalytic subunit", "PPP1CB" = "Serine/threonine-protein phosphatase PP1-beta catalytic subunit", "PPP1CC" = "Serine/threonine-protein phosphatase PP1-gamma catalytic subunit", "PPP1R1B" = "Protein phosphatase 1 regulatory subunit 1B (DARPP-32)", "PPP1R1C" = "Protein phosphatase 1 regulatory subunit 1C (I-3)", "PPP1R2" = "Protein phosphatase inhibitor 2", "PPP1R3A" = "Protein phosphatase 1 regulatory subunit 3A (GM)", "PPP1R3B" = "Protein phosphatase 1 regulatory subunit 3B (GL)", "PPP1R3C" = "Protein phosphatase 1 regulatory subunit 3C (PTG)", "PPP1R3D" = "Protein phosphatase 1 regulatory subunit 3D (R6)", "PPP1R3E" = "Protein phosphatase 1 regulatory subunit 3E", "PPP1R3F" = "Protein phosphatase 1 regulatory subunit 3F", "PPP1R3G" = "Protein phosphatase 1 regulatory subunit 3G", "PPP1R7" = "Protein phosphatase 1 regulatory subunit 7", "PPP1R8" = "Protein phosphatase 1 regulatory subunit 8 (NIPP1)", "PPP1R9A" = "Protein phosphatase 1 regulatory subunit 9A (neurabin-1)", "PPP1R9B" = "Protein phosphatase 1 regulatory subunit 9B (neurabin-2/spinophilin)", "PPP1R10" = "Protein phosphatase 1 regulatory subunit 10 (PNUTS)", "PPP1R11" = "Protein phosphatase 1 regulatory subunit 11 (PNUTS-like)", "PPP1R12A" = "Protein phosphatase 1 regulatory subunit 12A (MYPT1)", "PPP1R12B" = "Protein phosphatase 1 regulatory subunit 12B (MYPT2)", "PPP1R12C" = "Protein phosphatase 1 regulatory subunit 12C (MBS85)", "PPP1R13B" = "Protein phosphatase 1 regulatory subunit 13B (ASPP1)", "PPP1R14A" = "Protein phosphatase 1 regulatory subunit 14A (CPI-17)", "PPP1R14B" = "Protein phosphatase 1 regulatory subunit 14B (PHI-1)", "PPP1R14C" = "Protein phosphatase 1 regulatory subunit 14C", "PPP1R14D" = "Protein phosphatase 1 regulatory subunit 14D", "PPP1R15B" = "Protein phosphatase 1 regulatory subunit 15B (GADD34-related)", "PPP1R16A" = "Protein phosphatase 1 regulatory subunit 16A (MYPT3)", "PPP1R16B" = "Protein phosphatase 1 regulatory subunit 16B (TA-MYPT)", # Kinases "AURKB" = "Aurora kinase B", "AURKA" = "Aurora kinase A", "AURKC" = "Aurora kinase C", "JAK2" = "Tyrosine-protein kinase JAK2", "HASPIN" = "Serine/threonine-protein kinase haspin (GSG2)", "BUB1" = "Mitotic checkpoint serine/threonine-protein kinase BUB1", "WEE1" = "Wee1-like protein kinase", "SGO2" = "Shugoshin-like 2 (SGO2)", "SGO1" = "Shugoshin-like 1 (SGO1)", "DLK1" = "Protein delta homolog 1 (Preadipocyte factor 1/Pref-1)", "DAPK3" = "Death-associated protein kinase 3 (ZIPK)", "PKN1" = "Serine/threonine-protein kinase N1 (PKN1/PRK1)", "PKM" = "Pyruvate kinase PKM", "STK4" = "Serine/threonine-protein kinase 4 (MST1/STK4)", "MST1" = "Hepatocyte growth factor-like protein (MST1)", "SLK" = "STE20-like serine/threonine-protein kinase SLK", "PRKAA1" = "5'-AMP-activated protein kinase catalytic subunit alpha-1 (AMPKa1)", "PRKAA2" = "5'-AMP-activated protein kinase catalytic subunit alpha-2 (AMPKa2)", "DBF4" = "Protein DBF4 homolog A (Cdc7 kinase regulatory subunit)", "DYRK1A" = "Dual specificity tyrosine-phosphorylation-regulated kinase 1A", "RPS6KA4" = "Ribosomal protein S6 kinase alpha-4 (MSK2)", "RPS6KA5" = "Ribosomal protein S6 kinase alpha-5 (MSK1)", "PRKCB" = "Protein kinase C beta type (PKCb)", "JKAMP" = "JMJD8-associated regulatory complex component (JKAMP)", "MAPK8" = "Mitogen-activated protein kinase 8 (JNK1)", "MAPK9" = "Mitogen-activated protein kinase 9 (JNK2)", "MDC1" = "Mediator of DNA damage checkpoint protein 1", "MCPH1" = "Microcephalin (MCPH1; DNA damage sensor)", "PAK2" = "Serine/threonine-protein kinase PAK2", "ATM" = "Serine-protein kinase ATM (DNA damage checkpoint)", "ATR" = "Serine/threonine-protein kinase ATR", "CSNK2B" = "Casein kinase II subunit beta", "PRKCD" = "Protein kinase C delta type", "MAP2K1" = "Dual specificity mitogen-activated protein kinase kinase 1 (MEK1)", "MAP3K12" = "Mitogen-activated protein kinase kinase kinase 12 (DLK/MUK)", "MAP3K22" = "Mitogen-activated protein kinase kinase kinase 22", "MAPK1" = "Mitogen-activated protein kinase 1 (ERK2)", "MAPK3" = "Mitogen-activated protein kinase 3 (ERK1)", "MAPK8B" = "Mitogen-activated protein kinase 8b (JNK1b, zebrafish)", "MAPK11" = "Mitogen-activated protein kinase 11 (p38b)", "MAPK12" = "Mitogen-activated protein kinase 12 (p38g/ERK6)", "MAPK13" = "Mitogen-activated protein kinase 13 (p38d/SAPK4)", "MAPK14" = "Mitogen-activated protein kinase 14 (p38a)", "PRKDC" = "DNA-dependent protein kinase catalytic subunit (DNA-PKcs)", "PIM1" = "Proto-oncogene serine/threonine-protein kinase Pim-1", "CDK8" = "Cyclin-dependent kinase 8 (Mediator kinase)", "CHUK" = "Inhibitor of nuclear factor kappa-B kinase subunit alpha (IKKa)", "CDK1" = "Cyclin-dependent kinase 1 (CDC2)", "CHEK1" = "Serine/threonine-protein kinase Chk1", "APBB1" = "Amyloid precursor protein-binding family B member 1 (Fe65)", "YWHAQ" = "14-3-3 protein theta", # ZC3H12 ribonucleases "ZC3H12A" = "Endoribonuclease ZC3H12A (MCPIP1/Regnase-1)", "ZC3H12B" = "Endoribonuclease ZC3H12B (Regnase-2)", "ZC3H12C" = "Endoribonuclease ZC3H12C (Regnase-3)", "ZC3H12D" = "Endoribonuclease ZC3H12D (Regnase-4)", # Other "DCP2" = "mRNA-decapping enzyme subunit DCP2", "YBX1" = "Nuclease-sensitive element-binding protein 1 (YBX1/YB-1)", "RAG1" = "V(D)J recombination-activating protein 1 (RAG1)", "ALG13" = "UDP-N-acetylglucosamine transferase subunit ALG13", "HSPBAP1"= "HSPA binding protein 1 (HSPBAP1)" ) strip_isoform <- function(label) { # Remove trailing isoform number, e.g. Dnmt1-201 -> Dnmt1 sub("-\\d+$", "", label) } species_from_id <- function(eid) { dplyr::case_when( startsWith(eid, "ENSMUSP") ~ "Mus musculus", startsWith(eid, "ENSDARP") ~ "Danio rerio", TRUE ~ "Homo sapiens" ) } lookup_protein_name <- function(gene_name) { key <- toupper(gene_name) if (!is.na(protein_names[key])) return(unname(protein_names[key])) # Strip alias suffixes and try progressively shorter hyphen-joined prefixes simplified <- gsub("[\u2014\u2013/].*", "", key) parts <- strsplit(simplified, "-")[[1]] for (i in seq(length(parts), 1)) { candidate <- paste(parts[seq_len(i)], collapse = "-") if (!is.na(protein_names[candidate])) return(unname(protein_names[candidate])) } "" } # Parse fasta and build table lines <- readLines("../../data/Machinery.fasta") headers <- lines[startsWith(lines, ">")] parsed <- regmatches(headers, regexec("^>(\\S+)\\s+peptide:\\s+(\\S+)", headers)) label <- sapply(parsed, `[`, 2) # e.g. "Dnmt1-201" eid <- sapply(parsed, `[`, 3) # e.g. "ENSMUSP00000004202" gene_name <- strip_isoform(label) protein_name <- sapply(gene_name, lookup_protein_name, USE.NAMES = FALSE) machinery_ref <- data.frame( accession = eid, protein_name = protein_name, gene_name = gene_name, species = species_from_id(eid), stringsAsFactors = FALSE ) write.csv(machinery_ref, "../output/20-supplementary-files/Machinery_reference_table.csv", row.names = FALSE) ``` # Epimachinery category labels for Machinery.fasta table ```{r} machinery <- read_csv("../output/20-supplementary-files/Machinery_reference_table.csv", show_col_types = FALSE) classify_gene <- function(gene) { # Uppercase to handle mouse mixed-case gene names gene <- toupper(gene) cats <- list( "ADP-ribosylation" = list( prefix = c("PARP", "TNKS", "TARG", "MACROD", "OARD", "ADPRH", "ARH", "TIPARP"), exact = c("PARG", "ARTD1") ), "DNA methylation & reading" = list( prefix = c("DNMT", "TET", "MBD", "UHRF", "MECP", "ZBTB"), exact = c("PRDM14", "ZFP57") ), "RNA modification" = list( prefix = c("METTL", "YTHDF", "YTHDC", "ALKBH", "WTAP", "VIRMA", "ZC3H13", "IGF2BP", "NSUN", "DKC", "TRMT", "PUS", "NAT10", "ADAR", "TRUB", "RPUSD"), exact = c("FTO", "RBM15", "RBM15B", "PCIF1", "BUD23", "WDR4", "RBMX", "TRDMT1", "NOP2", "HNRNPA2B1", "DUS1L", # tRNA dihydrouridine synthase "TYW5", # tRNA wybutosine-synthesizing protein (Supp. folder: # "Histone demethylation" — incorrect; TYW5 is a # tRNA modification enzyme, not a histone demethylase) "CBLL1") # Hakai, E3 ubiquitin ligase component of the m6A # writer complex ), "ncRNA biogenesis & silencing" = list( prefix = c("AGO", "TNRC6", "DICER", "DROSHA", "DGCR8", "PIWI", "EXPORTIN", "TARBP", "PRKRA", "LIN28", "KHSRP", "ZCCHC", "DIS3", "INTS", "EXOSC", "PAN2", "PAN3", "DDX5", "DDX17", "DDX6", "EDC", "DCP", "XRN", "CNOT", "UPF", "SMG", "TUT", "MOV10", "FMR", "PUM", "STAU", "NCBP", "PHAX", "RBM7", "ZC3H12"), # ZC3H12A-D are endoribonucleases (Regnase-1-4) # that degrade target mRNAs in the cytoplasm — # mRNA turnover/silencing machinery exact = c("RDRP", "MTREX", "SKIV2L2", "MTR4", "ZFC3H1", "SRRT", "ARS2", "XPO5", "XPO1", "CRM1", "RAN", "EWSR1", "FUS", "TARDBP", "HNRNPK", "PARN", "PABPC1", "PAPD5", "DZIP3", "NONO", "SFPQ", "PSPC1", "YBX1", "HSPA8", "HSP90AA1", "HSP90AB1", "PQBP1", "MCM3AP", "PRPF8", # pre-mRNA splicing factor 8 — core spliceosome # component involved in snRNA processing "EIF3H", # eIF3 subunit H — translation initiation factor # involved in mRNA silencing/regulation "EIF3F") # eIF3 subunit F — translation initiation factor; # has MPN deubiquitinase fold but functions in # translation regulation and miRNA-mediated silencing ), "Histone modification and variants" = list( prefix = c("KMT", "KDM", "SETD", "EZH", "SUV", "EHMT", "MLL", "ASH1", "DOT1", "JMJD", "HDAC", "SIRT", "KAT", "PRMT", "CARM", "NSD", "WHSC", "SMYD", "CDYL", "RIOX", "MINA", "NO66", "PRDM", "SUMO", "SENP", "H1-", "H2A", "H2B", "H3", "H4", # histone variants "NAT"), # NAT8, NAT8B, NAT8L — N-acetyltransferases (camello # family); Supp. folder groups all NATs under # "Histone acetylation" exact = c("CREBBP", "EP300", "HAT1", "ELP3", "PHF8", "PHF2", "JARID2", "PKN1", "HSPBAP1", "WDR5", "ASH2L", "SETMAR", "OGT", "OGA", "MACROH2A1", # histone variant macro-H2A.1 "ATAT1", # alpha-tubulin acetyltransferase; Supp. folder: # "Histone acetylation" (grouped with KATs) "GTF3C1", "GTF3C2", # TFIIIC subunits; Supp. folder: # "Histone acetylation" (have HAT activity) "NAA10", "NAA50", "NAA60", # N-alpha-acetyltransferases; # Supp. folder: "Histone acetylation" "BRD4", # bromodomain-containing protein 4; reads # acetylated histones; Supp. folder: # "Histone acetylation" "ESCO1", "ESCO2", # N-acetyltransferases (cohesin establishment); # Supp. folder: "Histone acetylation" "DLAT", # dihydrolipoyllysine-residue acetyltransferase; # Supp. folder: "Histone acetylation" "ATF2", # cyclic AMP-dependent TF; histone acetyltransferase # activity; Supp. folder: "Histone acetylation" "BLOC1S1", # biogenesis of lysosome-related organelles; # GCN5L1; Supp. folder: "Histone acetylation" "ACAT1", "ACAT2", # acetyl-CoA acetyltransferases; Supp. folder: # "Histone acetylation" "UTY", # inactive histone demethylase (KDM6 family homolog) "HR", # lysine demethylase and nuclear receptor corepressor "MAPT", # microtubule-associated protein tau; has histone # H3 binding / chromatin roles "BAZ2B") # bromodomain adjacent to zinc finger; chromatin # remodeling reader (Supp. folder: not present, # but BAZ1B is in "Histone phosphorylation") ), "Ubiquitin signaling" = list( prefix = c("USP", "UBE", "UBR", "UBA", "UCHL", "UCH", "RNF", "TRIM", "HUWE", "HERC", "MARCH", "RBX", "CUL", "FBX", "SKP", "SOCS", "OTU", "UBQLN", "PSMD", "JOSD", "DDB", "MIB"), exact = c("VCPIP1", "BAP1", "ATXN3", "ATXN3L", # ataxin-3-like protein; Josephin domain DUB "COPS6", "COPS5", "PARK2", "ZRANB1", "BRCA1", "RAG1", "TOM1", "DTX3L", "BRCC3", # Lys-63-specific deubiquitinase; Supp. folder: # "Histone Ubiquitination" (Deubiquitination) "MYSM1", # histone H2A deubiquitinase; Supp. folder: # "Histone Ubiquitination" (Deubiquitination) "STAMBPL1", # AMSH-LP; deubiquitinase; Supp. folder: # "Histone Ubiquitination" (Deubiquitination) "MPND", # MPN domain-containing protein; DUB fold; # associated with ubiquitin signaling "ALG13") # UDP-N-acetylglucosamine transferase subunit; # has deubiquitinase-like fold; Supp. folder: # not present, but classified here per its # role in ubiquitin-related pathways ), "Chromatin signaling" = list( prefix = c("MAP3K", "MAP2K", "MAPK", "EYA", "PPP4", "PPP6"), exact = c("CHUK", "IKBKB", "IKBKG", "PPP1CB", "PPP1CA", "PPP1CC", "PPP1R7", "PPP1R10", "PPP1R8", "PPP2CA", "AURKB", "MSK1", "MSK2", "RPS6KA4", "RPS6KA5", "PARK7", "BUB1", "DLK1", "DAPK3", "PKM", "JAK2", "DBF4", "DYRK1A", "YWHAQ", "PRKCB", "PRKCD", "AURKA", "AURKC", "STK4", "MST1", "MDC1", "MCPH1", "PAK2", "HIRA", "ATM", "ATR", "CSNK2B", "PRKAA1", "PRKAA2", "PRKDC", "PIM1", "CDK8", "CDK1", "CHEK1", "PPM1D", "PPP6C", "PPP4C", "CDCA2", "APBB1", "BAZ1B", "JKAMP", "ARID1A", "ARID1B", "SLK", # STE20-like kinase; Supp. folder: # "Histone phosphorylation" (slkb) "WEE1", # Wee1-like kinase; Supp. folder: # "Histone phosphorylation" "HASPIN", # histone H3T3 kinase (GSG2); chromatin signaling "SGO1", "SGO2") # Shugoshin proteins; protect centromeric # cohesion by recruiting PP2A; chromatin- # associated signaling at kinetochores ) ) for (cat_name in names(cats)) { cat_def <- cats[[cat_name]] if (gene %in% cat_def$exact) return(cat_name) for (p in cat_def$prefix) { if (str_starts(gene, fixed(p))) return(cat_name) } } return("Unclassified") } # Apply classification and add category column machinery_categorized <- machinery %>% mutate(epimachinery_category = vapply(gene_name, classify_gene, character(1))) machinery_categorized write_csv(machinery_categorized, "../output/20-supplementary-files/Machinery_categorized.csv") ``` # Alignment summary tables Compiling per-sample alignment and processing statistics for all three species ## RNA-seq Set up the metadata ```{r} species_meta <- tibble( species = c("D-Apul", "E-Peve", "F-Ptuh"), species_full = c("Acropora pulchra", "Porites evermanni", "Pocillopora tuahiniensis"), genome = c("Apulcra-genome.fa", "Porites_evermanni_v1.fa", "Pocillopora_meandrina_HIv1.assembly.fasta"), rna_prefix = c("../../D-Apul", "../../E-Peve", "../../F-Ptuh") ) species_meta ``` ### 1.1 Raw reads from MultiQC FastQC on raw reads ```{r} raw_reads <- rbind( read_tsv("../../D-Apul/output/01-Apul-RNA-trimming-FastQC/raw-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences'), read_tsv("../../E-Peve/output/01-Peve-RNA-trimming-FastQC/raw-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences'), read_tsv("../../F-Ptuh/output/01-Ptuh-RNA-trimming-FastQC/raw-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences') ) raw_reads ``` ### 1.2 Trimmed reads from MultiQC FastQC on trimmed reads ```{r} trimmed_reads <- rbind( read_tsv("../../D-Apul/output/01-Apul-RNA-trimming-FastQC/trimmed-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences'), read_tsv("../../E-Peve/output/01-Peve-RNA-trimming-FastQC/trimmed-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences'), read_tsv("../../F-Ptuh/output/01-Ptuh-RNA-trimming-FastQC/trimmed-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences') ) trimmed_reads ``` ### 1.3 Retained after trimming (%) ```{r} read_counts <- raw_reads %>% left_join(trimmed_reads, by = c("Sample")) %>% mutate(retained_pct = round(trimmed_reads$'Total Sequences' / raw_reads$'Total Sequences' * 100, 2)) read_counts ``` ### 1.4 Alignment rate, unique %, multi % ```{r} parse_hisat2_log <- function(lines) { text <- paste(lines, collapse = "\n") total_pairs <- as.numeric(str_match(text, "^(\\d+) reads")[2]) total_reads <- total_pairs * 2 # Overall alignment rate overall_rate <- as.numeric(str_match(text, "(\\S+)% overall alignment rate")[2]) # Concordant pairs conc_0 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned concordantly 0 times")[2]) conc_1 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned concordantly exactly 1 time")[2]) conc_gt1 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned concordantly >1 times")[2]) # Discordant pairs (from concordant-0-times) disc_1 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned discordantly 1 time")[2]) disc_1 <- if (is.na(disc_1)) 0 else disc_1 # Individual mates (from pairs that aligned 0 times concordantly AND discordantly) mates_total <- as.numeric(str_match(text, "(\\d+) mates make up the pairs")[2]) mates_total <- if (is.na(mates_total)) 0 else mates_total mates_0 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned 0 times\n")[2]) mates_0 <- if (is.na(mates_0)) 0 else mates_0 mates_1 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned exactly 1 time\n")[2]) mates_1 <- if (is.na(mates_1)) 0 else mates_1 mates_gt1 <- as.numeric(str_match(text, "(\\d+) \\(\\S+%\\) aligned >1 times\n")[2]) mates_gt1 <- if (is.na(mates_gt1)) 0 else mates_gt1 # Compute READ-level unique/multi (combining all alignment levels) # Each concordant pair = 2 reads; each discordant pair = 2 reads unique_reads <- (conc_1 * 2) + (disc_1 * 2) + mates_1 multi_reads <- (conc_gt1 * 2) + mates_gt1 unmapped_reads <- mates_0 mapped_reads <- unique_reads + multi_reads mapped_read_pairs <- mapped_reads/2 # Sanity check: unique + multi + unmapped should ≈ total_reads # (may not be exact due to mates from discordant pairs that also aligned individually) tibble( total_pairs = total_pairs, total_reads = total_reads, overall_align_rate = overall_rate, unique_pct = round(unique_reads / total_reads * 100, 2), multi_pct = round(multi_reads / total_reads * 100, 2), mapped_read_pairs = mapped_read_pairs ) } # D-Apul: consolidated log (all 5 samples in one file) parse_hisat2_consolidated <- function(path, samples) { lines <- readLines(path) block_starts <- which(str_detect(lines, "^\\d+ reads; of these:")) results <- list() for (i in seq_along(block_starts)) { end <- if (i < length(block_starts)) block_starts[i + 1] - 1 else length(lines) block <- lines[block_starts[i]:end] results[[i]] <- parse_hisat2_log(block) %>% mutate(sample = samples[i]) } bind_rows(results) } # E-Peve / F-Ptuh: one stderr file per sample parse_hisat2_per_file <- function(dir_path) { files <- list.files(dir_path, pattern = "_hisat\\.stderr$", full.names = TRUE) results <- list() for (f in files) { sname <- basename(f) %>% str_remove("_hisat.stderr") %>% str_remove("^RNA-") results[[sname]] <- parse_hisat2_log(readLines(f)) %>% mutate(sample = sname) } bind_rows(results) } hisat2_stats <- bind_rows( parse_hisat2_consolidated( "../../D-Apul/output/07-Apul-Hisat/hisat.out", samples = c("ACR-140", "ACR-145", "ACR-150", "ACR-173", "ACR-178") ), parse_hisat2_per_file("../../E-Peve/output/06-Peve-Hisat"), parse_hisat2_per_file("../../F-Ptuh/output/06-Ptuh-Hisat") ) %>% mutate(species = case_when( str_detect(sample, "^ACR") ~ "D-Apul", str_detect(sample, "^POR") ~ "E-Peve", str_detect(sample, "^POC") ~ "F-Ptuh" )) %>% select(species, sample, everything()) hisat2_stats ``` ### 1.5 Mismatch % This requires collecting info from the full *.sorted.bam files, which are very large and live on Gannet (https://gannet.fish.washington.edu/seashell/bu-github/deep-dive-expression/). ```{r, engine='bash', eval=FALSE} # The files are backed up to Steven's backup directory: # /volume/web/seashell/bu-github/deep-dive-expression # The BAM paths are: # D-Apul: D-Apul/output/07-Apul-Hisat/RNA-ACR-{140,145,150,173,178}.sorted.bam # E-Peve: E-Peve/output/06-Peve-Hisat/RNA-POR-{71,73,76,79,82}.sorted.bam # F-Ptuh: F-Ptuh/output/06-Ptuh-Hisat/RNA-POC-{47,48,50,53,57}.sorted.bam mkdir -p ../output/20-supplementary-files/RNAseq_HISAT_samtools_stats # D-Apul BAMs wget -r -l1 -nd -nc -A "*.sorted.bam" \ -P ../output/20-supplementary-files/RNAseq_HISAT_samtools_stats \ https://gannet.fish.washington.edu/seashell/bu-github/deep-dive-expression/D-Apul/output/07-Apul-Hisat/ # E-Peve BAMs wget -r -l1 -nd -nc -A "*.sorted.bam" \ -P ../output/20-supplementary-files/RNAseq_HISAT_samtools_stats \ https://gannet.fish.washington.edu/seashell/bu-github/deep-dive-expression/E-Peve/output/06-Peve-Hisat/ # F-Ptuh BAMs wget -r -l1 -nd -nc -A "*.sorted.bam" \ -P ../output/20-supplementary-files/RNAseq_HISAT_samtools_stats \ https://gannet.fish.washington.edu/seashell/bu-github/deep-dive-expression/F-Ptuh/output/06-Ptuh-Hisat/ # Then run samtools on each BAM for bam in ../output/20-supplementary-files/RNAseq_HISAT_samtools_stats/*.sorted.bam; do base=$(basename "$bam" .sorted.bam) /home/shared/samtools-1.12/samtools stats "$bam" > "../output/20-supplementary-files/RNAseq_HISAT_samtools_stats/${base}.stats.txt" /home/shared/samtools-1.12/samtools flagstat "$bam" > "../output/20-supplementary-files/RNAseq_HISAT_samtools_stats/${base}.flagstat.txt" done ``` ```{r, engine='bash', eval=FALSE} # Remove BAMs once done (they're quite large) rm ../output/20-supplementary-files/RNAseq_HISAT_samtools_stats/*.sorted.bam ``` Parse samtools stats output for mismatch rate ```{r} samtools_dir <- "../output/20-supplementary-files/RNAseq_HISAT_samtools_stats" # Sample configuration samples_cfg <- list( D_Apul = c("ACR-140", "ACR-145", "ACR-150", "ACR-173", "ACR-178"), E_Peve = c("POR-71", "POR-73", "POR-76", "POR-79", "POR-82"), F_Ptuh = c("POC-47", "POC-48", "POC-50", "POC-53", "POC-57") ) species_map <- c(D_Apul = "D-Apul", E_Peve = "E-Peve", F_Ptuh = "F-Ptuh") # Parse a single samtools stats file parse_samtools_stats <- function(path, sample_name, species_code) { if (!file.exists(path)) { return(tibble(species = species_code, sample = sample_name, mapped_reads = NA, mismatch_rate_pct = NA)) } lines <- readLines(path) # Mapped reads: SN\treads mapped:\t28412345 mapped <- str_match(lines[grepl("^SN\\treads mapped:", lines)], "^SN\\treads mapped:\\t(\\d+)")[2] mapped <- if (!is.na(mapped)) as.numeric(mapped) else NA # Mismatch rate (reported as "error rate"): # SN\terror rate:\t7.151265e-03\t# mismatches / bases mapped (cigar) error_line <- lines[grepl("^SN\\terror rate:", lines)] error_rate <- str_match(error_line, "^SN\\terror rate:\\t([\\d.eE+-]+)")[2] mismatch_rate_pct <- if (!is.na(error_rate)) as.numeric(error_rate) * 100 else NA tibble(species = species_code, sample = sample_name, mapped_reads = mapped, mismatch_rate_pct = round(mismatch_rate_pct, 4)) } # Build results from single directory samtools_results <- bind_rows(lapply(names(samples_cfg), function(sp) { species_code <- species_map[sp] bind_rows(lapply(samples_cfg[[sp]], function(s) { parse_samtools_stats( file.path(samtools_dir, paste0("RNA-", s, ".stats.txt")), s, species_code ) })) })) # Rename for final output samtools_summary <- samtools_results %>% rename( Species = species, Sample = sample, `Mapped reads` = mapped_reads, `Mismatch rate (%)` = mismatch_rate_pct ) samtools_summary ``` ### 1.6 Number of lncRNAs identified ```{r} count_lncRNAs <- function(filepath) { df <- read_tsv(filepath, col_names = TRUE) sample_cols <- grep(".sorted.bam$", colnames(df), value = TRUE) sample_names <- str_remove(sample_cols, ".*pipeline\\.RNA\\.") %>% str_remove("\\.sorted\\.bam$") # Count non-zero entries per sample column nonzero_counts <- sapply(sample_cols, function(col) { sum(df[[col]] > 0, na.rm = TRUE) }) # Build a summary tibble summary <- tibble( sample = sample_names, n_unique_lncRNA = nonzero_counts ) print(summary) } lncRNA_counts <- rbind( count_lncRNAs("../output/01.6-lncRNA-pipeline/Apul-lncRNA-counts-filtered.txt"), count_lncRNAs("../output/01.6-lncRNA-pipeline/Peve-lncRNA-counts-filtered.txt"), count_lncRNAs("../output/01.6-lncRNA-pipeline/Ptuh-lncRNA-counts-filtered.txt") ) lncRNA_counts ``` ## sRNA-seq alignment statistics ### 2.1 Species metadata (genome) ```{r} # Genome files are referenced in the ShortStack Rmd scripts: srna_genomes <- tibble( species = c("D-Apul", "E-Peve", "F-Ptuh"), genome = c("Apulcra-genome.fa", "Porites_evermanni_v1.fa", "Pocillopora_meandrina_HIv1.assembly.fasta") ) srna_genomes ``` ### 2.1 Raw reads ```{r} s_raw_reads <- rbind( read_tsv("~/deep-dive/D-Apul/output/08-Apul-sRNAseq-trimming/raw-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences'), read_tsv("~/deep-dive/E-Peve/output/06-Peve-sRNAseq-trimming/raw-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences'), read_tsv("~/deep-dive/F-Pmea/output/08-Pmea-sRNAseq-trimming/raw-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "RNA-")) %>% filter(str_detect(Sample, "_R1_001$")) %>% select(Sample, 'Total Sequences') ) s_raw_reads$Sample <- str_remove(s_raw_reads$Sample, "_R1_001") s_raw_reads ``` ### 2.2 Trimmed reads ```{r} s_trimmed_reads <- rbind( read_tsv("~/deep-dive/D-Apul/output/08.2-Apul-sRNAseq-trimming-31bp-fastp-merged/trimmed-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "sRNA-")) %>% select(Sample, 'Total Sequences'), read_tsv("~/deep-dive/E-Peve/output/06.2-Peve-sRNAseq-trimming-31bp-fastp-merged/trimmed-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% select(Sample, 'Total Sequences') %>% mutate(Sample = paste0("sRNA-", Sample)), read_tsv("~/deep-dive/F-Pmea/output/08.2-Pmea-sRNAseq-trimming-31bp-fastp-merged/trimmed-fastqc/multiqc_data/multiqc_fastqc.txt", col_names = TRUE) %>% filter(str_detect(Sample, "sRNA-")) %>% select(Sample, 'Total Sequences') ) s_trimmed_reads$Sample <- str_remove(s_trimmed_reads$Sample, "-fastp-adapters-polyG-31bp-merged") s_trimmed_reads ``` ### 2.3 Retained after trimming (%) ```{r} s_read_counts <- s_raw_reads %>% left_join(s_trimmed_reads, by = c("Sample")) %>% mutate(retained_pct = round(s_trimmed_reads$'Total Sequences' / s_raw_reads$'Total Sequences' * 100, 2)) s_read_counts ``` ### 2.4 ShortStack alignment stats from alignment_details.tsv NOTE: the Github repo README.md states that the H tag means "H: Very highly multi-mapped read (>=50 hits)". However,the log from running ShortStack 4.1.0 (e.g., the file ~/deep-dive-expression/D-Apul/output/11-Apul-sRNA-ShortStack_4.1.0-pulchra_genome/shortstack.log) states that the tag represents "Very highly multi-mapped (>=20 hits)(H)". It looks like the >=50 threshold was used in previous versions of ShortStack (confirmed by looking at deep-dive ShortStack logs), but the updated 4.1.0 version switched to a >=20 threshold. The Github repo's README is out of date. ```{r} # alignment_details.tsv columns: # readfile mapping_type sequence_length sequence_count read_count # mapping_type (from https://github.com/MikeAxtell/ShortStack) # U: Uniquely mapped (not a multimapper). # P: Multimapper placed using the method set by option --mmap. # R: Multimapper placed at random. # H: Very highly multi-mapped read (>=20 hits). # N: Unmapped reads. read_alignment_details <- function(path, species_code) { df <- read_tsv(path, show_col_types = FALSE) # Extract sample name from readfile path df <- df %>% mutate( sample = str_extract(readfile, "(ACR|POR|POC)-\\d+"), species = species_code ) # Aggregate by sample × mapping_type, sum sequence_count agg <- df %>% group_by(species, sample, mapping_type) %>% summarise(seq_count = sum(read_count), .groups = "drop") # Pivot to wide format agg_wide <- agg %>% pivot_wider(names_from = mapping_type, values_from = seq_count, values_fill = 0) # totals and percentages # U=unique, P=guided, R=random, H=very-high, N=not mapped agg_wide %>% mutate( total_reads = U + P + R + H + N, mapped_reads = U + P + R + H, overall_align_pct = round((mapped_reads) / total_reads * 100, 2), unique_pct = round(U / total_reads * 100, 2), guided_pct = round(P / total_reads * 100, 2), random_pct = round(R / total_reads * 100, 2), very_high_pct = round(H / total_reads * 100, 2), ) %>% select( species, sample, total_reads, overall_align_pct, unique_pct, guided_pct, random_pct, very_high_pct, mapped_reads ) %>% rename( `Overall alignment (%)` = overall_align_pct, `Uniquely mapped reads (%)` = unique_pct, `Multi-mapped reads placed with guidance (%)` = guided_pct, `Multi-mapped reads placed randomly (%)` = random_pct, `Very high multi-mapping reads (>=50 hits; %)` = very_high_pct, `Mapped sRNA reads` = mapped_reads ) } srna_align <- bind_rows( read_alignment_details("../../D-Apul/output/11-Apul-sRNA-ShortStack_4.1.0-pulchra_genome/ShortStack_out/alignment_details.tsv", "D-Apul"), read_alignment_details("../../E-Peve/output/05-Peve-sRNA-ShortStack_4.1.0/ShortStack_out/alignment_details.tsv", "E-Peve"), read_alignment_details("../../F-Ptuh/output/05-Ptuh-sRNA-ShortStack_4.1.0/ShortStack_out/alignment_details.tsv", "F-Ptuh") ) srna_align ``` ### 2.5 Mapped miRNA reads from ShortStack Counts.txt ```{r} # Counts.txt is a TSV: Coords, Name, MIRNA, _condensed, ... # Filter rows where MIRNA == "Y", then sum each sample column. read_mirna_counts <- function(path, species_code) { df <- read_tsv(path, show_col_types = FALSE) sample_cols <- colnames(df)[-(1:3)] # skip Coords, Name, MIRNA # Extract clean sample names from column names # e.g. "sRNA-ACR-140-S1-TP2-..._condensed" -> "ACR-140" samples <- str_extract(sample_cols, "(ACR|POR|POC)-\\d+") mirna_only <- df %>% filter(MIRNA == "Y") # Sum each sample column across all miRNA rows mapped_mirna <- colSums(mirna_only[, sample_cols]) tibble( species = species_code, sample = samples, mapped_mirna_reads = mapped_mirna ) } srna_mirna_reads <- bind_rows( read_mirna_counts("../../D-Apul/output/11-Apul-sRNA-ShortStack_4.1.0-pulchra_genome/ShortStack_out/Counts.txt", "D-Apul"), read_mirna_counts("../../E-Peve/output/05-Peve-sRNA-ShortStack_4.1.0/ShortStack_out/Counts.txt", "E-Peve"), read_mirna_counts("../../F-Ptuh/output/05-Ptuh-sRNA-ShortStack_4.1.0/ShortStack_out/Counts.txt", "F-Ptuh") ) srna_mirna_reads ``` ### 2.6 Number of miRNAs identified ```{r} count_expressed_mirnas <- function(path) { df <- read_tsv(path, show_col_types = FALSE) sample_cols <- grep("merged_condensed", colnames(df), value = TRUE) as.data.frame(t(df %>% filter(MIRNA == "Y") %>% summarise(across(all_of(sample_cols), ~ sum(.x > 0, na.rm = TRUE))))) } miRNA_counts <- rbind( count_expressed_mirnas("../../D-Apul/output/11-Apul-sRNA-ShortStack_4.1.0-pulchra_genome/ShortStack_out/Counts.txt"), count_expressed_mirnas("../../E-Peve/output/05-Peve-sRNA-ShortStack_4.1.0/ShortStack_out/Counts.txt"), count_expressed_mirnas("../../F-Ptuh/output/05-Ptuh-sRNA-ShortStack_4.1.0/ShortStack_out/Counts.txt") ) colnames(miRNA_counts) <- c("miRNA_count") miRNA_counts$Sample <- rownames(miRNA_counts) miRNA_counts <- miRNA_counts %>% select(Sample, miRNA_count) %>% remove_rownames() miRNA_counts$Sample <- miRNA_counts$Sample %>% str_remove("-S1-TP2-fastp-adapters-polyG-31bp-merged_condensed") %>% str_remove("sRNA-") miRNA_counts ``` # Mature miRNA names and sequences for all species Already have per-species, just need to combine into single table ```{r} miRNA_details <- rbind( read.csv("../../D-Apul/output/11-Apul-sRNA-ShortStack_4.1.0-pulchra_genome/ShortStack_out/Apul_Results_mature_named_miRNAs.csv", header = TRUE) %>% dplyr::select(-X), # Including an NA filter here, because it looks like the Peve file has an extra row indicating miRNAs identified in # Ashey et al 2026 which were not also identified here using the new version of ShortStack (4.1.0) read.csv("../../E-Peve/output/05-Peve-sRNA-ShortStack_4.1.0/ShortStack_out/Peve_Results_mature_named_miRNAs.csv", header = TRUE) %>% dplyr::select(-X) %>% drop_na(Name), read.csv("../../F-Ptuh/output/05-Ptuh-sRNA-ShortStack_4.1.0/ShortStack_out/Ptuh_Results_mature_named_miRNAs.csv", header = TRUE) %>% dplyr::select(-X) ) miRNA_details write.csv(miRNA_details, "../output/20-supplementary-files/AllSpecies_Results_mature_named_miRNAs.csv", col.names = TRUE, row.names = FALSE) ``` # miRanda binding information + PCC coexpression (miRNA-mRNA and miRNA-lncRNA) ## miRNA-mRNA ```{r} consolidate_regions <- function(path_3UTR, path_CDS, path_5UTR){ # Load for each region data_3UTR <- read.csv(path_3UTR) %>% dplyr::select(-X, -X.1) data_CDS <- read.csv(path_CDS) %>% dplyr::select(-X) data_5UTR <- read.csv(path_5UTR) %>% dplyr::select(-X) # Format and combine into pooled targets data_3UTR$region <- "3UTR" colnames(data_CDS) <- c("miRNA", "mRNA_coord", "score", "energy", "query_start", "query_end", "subject_start", "subject_end", "total_bp_shared", "query_similar", "subject_similar", "mRNA", "PCC.cor", "p_value", "adjusted_p_value") data_CDS$query_start_end <- paste0(data_CDS$query_start, " ", data_CDS$query_end) data_CDS$subject_start_end <- paste0(data_CDS$subject_start, " ", data_CDS$subject_end) data_CDS$region <- "CDS" data_CDS <- data_CDS %>% dplyr::select(colnames(data_3UTR)) colnames(data_5UTR) <- c("miRNA", "mRNA_coord", "score", "energy", "query_start", "query_end", "subject_start", "subject_end", "total_bp_shared", "query_similar", "subject_similar", "mRNA", "PCC.cor", "p_value", "adjusted_p_value") data_5UTR$query_start_end <- paste0(data_5UTR$query_start, " ", data_5UTR$query_end) data_5UTR$subject_start_end <- paste0(data_5UTR$subject_start, " ", data_5UTR$subject_end) data_5UTR$region <- "5UTR" data_5UTR <- data_5UTR %>% dplyr::select(colnames(data_3UTR)) # Confirm all 3 region tables have same column names and orders ifelse(identical(colnames(data_3UTR), colnames(data_CDS)) == TRUE, ifelse(identical(colnames(data_3UTR), colnames(data_5UTR)) == TRUE, data <- rbind(data_3UTR, data_CDS, data_5UTR), data <- NULL) ) data } apul_miRanda_PCC_mRNA <- consolidate_regions( "../../D-Apul/output/09-Apul-mRNA-miRNA-interactions/miranda_PCC_miRNA_mRNA.csv", "../../D-Apul/output/09.01-Apul-mRNA-miRNA-interactions-CDS_5UTR/miRanda_PCC_miRNA_CDS.csv", "../../D-Apul/output/09.01-Apul-mRNA-miRNA-interactions-CDS_5UTR/miRanda_PCC_miRNA_5UTR.csv") peve_miRanda_PCC_mRNA <- consolidate_regions( "../../E-Peve/output/10-Peve-mRNA-miRNA-interactions/Peve-miranda_PCC_miRNA_mRNA.csv", "../../E-Peve/output/10.01-Peve-mRNA-miRNA-interactions-CDS_5UTR/miRanda_PCC_miRNA_CDS.csv", "../../E-Peve/output/10.01-Peve-mRNA-miRNA-interactions-CDS_5UTR/miRanda_PCC_miRNA_5UTR.csv" ) ptuh_miRanda_PCC_mRNA <- consolidate_regions( "../../F-Ptuh/output/11-Ptuh-mRNA-miRNA-interactions/three_prime_interaction/Ptuh-miranda_PCC_miRNA_mRNA.csv", "../../F-Ptuh/output/11.01-Ptuh-mRNA-miRNA-interactions-CDS_5UTR/miRanda_PCC_miRNA_CDS.csv", "../../F-Ptuh/output/11.01-Ptuh-mRNA-miRNA-interactions-CDS_5UTR/miRanda_PCC_miRNA_5UTR.csv" ) # Combine miRanda_PCC_mRNA <- rbind(apul_miRanda_PCC_mRNA, peve_miRanda_PCC_mRNA, ptuh_miRanda_PCC_mRNA) # Filter to only significant PCC (this consitutes our putative interaction set, which is both # predicted to bind via miRanda and significantly coexpressed) miRanda_PCC_mRNA_sig <- miRanda_PCC_mRNA %>% filter(p_value < 0.05) # Save write.csv(miRanda_PCC_mRNA_sig, "../output/20-supplementary-files/miRanda_PCC_mRNA_sig.csv", row.names = FALSE) write.csv(miRanda_PCC_mRNA, "../output/20-supplementary-files/miRanda_PCC_mRNA.csv", row.names = FALSE) ``` ## miRNA-lncRNA ```{r} apul_miRanda_PCC_lncRNA <- read.csv("../../D-Apul/output/28-Apul-miRNA-lncRNA-interactions/miranda_PCC_miRNA_lncRNA.csv") %>% dplyr::select(-X, -X.1) peve_miRanda_PCC_lncRNA <- read.csv("../../E-Peve/output/15-Peve-miRNA-lncRNA-PCC/miranda_PCC_miRNA_lncRNA.csv") %>% dplyr::select(-X, -X.1) ptuh_miRanda_PCC_lncRNA <- read.csv("../../F-Ptuh/output/15-Ptuh-miRNA-lncRNA-PCC/miranda_PCC_miRNA_lncRNA.csv") %>% dplyr::select(-X, -X.1) # Ensure all three tables have matching column names and order identical(colnames(apul_miRanda_PCC_lncRNA), colnames(peve_miRanda_PCC_lncRNA)) identical(colnames(apul_miRanda_PCC_lncRNA), colnames(ptuh_miRanda_PCC_lncRNA)) # Combine miRanda_PCC_lncRNA <- rbind(apul_miRanda_PCC_lncRNA, peve_miRanda_PCC_lncRNA, ptuh_miRanda_PCC_lncRNA) # Reorder columns to have same order as miRNA-mRNA table (above) miRanda_PCC_lncRNA <- miRanda_PCC_lncRNA %>% dplyr::select(miRNA, lncRNA, PCC.cor, p_value, adjusted_p_value, score, energy, query_start_end, subject_start_end, total_bp_shared, query_similar, subject_similar) # Filter to only significant PCC (this constitutes our putative interaction set, which is both # predicted to bind via miRanda and significantly coexpressed) miRanda_PCC_lncRNA_sig <- miRanda_PCC_lncRNA %>% filter(p_value < 0.05) # Save write.csv(miRanda_PCC_lncRNA_sig, "../output/20-supplementary-files/miRanda_PCC_lncRNA_sig.csv", row.names = FALSE) write.csv(miRanda_PCC_lncRNA, "../output/20-supplementary-files/miRanda_PCC_lncRNA.csv", row.names = FALSE) ```