Processing Data_Tam_RNAseq_2025_LB_vs_Mac_ATCC19606

gene_x 0 like s 229 view s

Tags: pipeline, RNA-seq

  1. Targets

    1. Could you please assist me with processing RNA-seq data? The reference genome is CP059040. I aim to analyze the data using PCA, a Venn diagram, and KEGG and GO annotation enrichment analysis.
    2. The samples are labeled as follows (where 'x' indicates the replicate number):
    4.     LB-AB-x
    5.     LB-IJ-x
    6.     LB-W1-x
    7.     LB-WT19606-x
    8.     LB-Y1-x
    9.     Mac-AB-x
    10.     Mac-IJ-x
    11.     Mac-W1-x
    12.     Mac-WT19606-x
    13.     Mac-Y1-x

  2. Download the raw data

    1. ./lnd login -u X101SC25015922-Z02-J002 -p m*********5
    2. ./lnd list
    3. ./lnd cp -d oss://  ./
    4. ./lnd cp oss://CP2024102300053 .  #Error
    5. ./lnd list oss://CP2024102300053
    6. ./lnd cp -d oss://CP2024102300053/H101SC25015922/RSMR00204 .
    7. #CP2024102300053/H101SC25015922/RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002

  3. Prepare raw data

    1. mkdir raw_data; cd raw_data
    3. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-AB-1/LB-AB-1_1.fq.gz LB-AB-r1_R1.fq.gz
    4. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-AB-1/LB-AB-1_2.fq.gz LB-AB-r1_R2.fq.gz
    5. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-AB-2/LB-AB-2_1.fq.gz LB-AB-r2_R1.fq.gz
    6. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-AB-2/LB-AB-2_2.fq.gz LB-AB-r2_R2.fq.gz
    7. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-AB-3/LB-AB-3_1.fq.gz LB-AB-r3_R1.fq.gz
    8. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-AB-3/LB-AB-3_2.fq.gz LB-AB-r3_R2.fq.gz
    10. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-IJ-1/LB-IJ-1_1.fq.gz LB-IJ-r1_R1.fq.gz
    11. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-IJ-1/LB-IJ-1_2.fq.gz LB-IJ-r1_R2.fq.gz
    12. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-IJ-2/LB-IJ-2_1.fq.gz LB-IJ-r2_R1.fq.gz
    13. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-IJ-2/LB-IJ-2_2.fq.gz LB-IJ-r2_R2.fq.gz
    14. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-IJ-4/LB-IJ-4_1.fq.gz LB-IJ-r4_R1.fq.gz
    15. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-IJ-4/LB-IJ-4_2.fq.gz LB-IJ-r4_R2.fq.gz
    17. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-W1-1/LB-W1-1_1.fq.gz LB-W1-r1_R1.fq.gz
    18. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-W1-1/LB-W1-1_2.fq.gz LB-W1-r1_R2.fq.gz
    19. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-W1-2/LB-W1-2_1.fq.gz LB-W1-r2_R1.fq.gz
    20. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-W1-2/LB-W1-2_2.fq.gz LB-W1-r2_R2.fq.gz
    21. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-W1-3/LB-W1-3_1.fq.gz LB-W1-r3_R1.fq.gz
    22. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-W1-3/LB-W1-3_2.fq.gz LB-W1-r3_R2.fq.gz
    24. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-WT19606-2/LB-WT19606-2_1.fq.gz LB-WT19606-r2_R1.fq.gz
    25. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-WT19606-2/LB-WT19606-2_2.fq.gz LB-WT19606-r2_R2.fq.gz
    26. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-WT19606-3/LB-WT19606-3_1.fq.gz LB-WT19606-r3_R1.fq.gz
    27. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-WT19606-3/LB-WT19606-3_2.fq.gz LB-WT19606-r3_R2.fq.gz
    28. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-WT19606-4/LB-WT19606-4_1.fq.gz LB-WT19606-r4_R1.fq.gz
    29. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-WT19606-4/LB-WT19606-4_2.fq.gz LB-WT19606-r4_R2.fq.gz
    31. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-Y1-2/LB-Y1-2_1.fq.gz LB-Y1-r2_R1.fq.gz
    32. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-Y1-2/LB-Y1-2_2.fq.gz LB-Y1-r2_R2.fq.gz
    33. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-Y1-3/LB-Y1-3_1.fq.gz LB-Y1-r3_R1.fq.gz
    34. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-Y1-3/LB-Y1-3_2.fq.gz LB-Y1-r3_R2.fq.gz
    35. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-Y1-4/LB-Y1-4_1.fq.gz LB-Y1-r4_R1.fq.gz
    36. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/LB-Y1-4/LB-Y1-4_2.fq.gz LB-Y1-r4_R2.fq.gz
    38. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-AB-1/Mac-AB-1_1.fq.gz Mac-AB-r1_R1.fq.gz
    39. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-AB-1/Mac-AB-1_2.fq.gz Mac-AB-r1_R2.fq.gz
    40. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-AB-2/Mac-AB-2_1.fq.gz Mac-AB-r2_R1.fq.gz
    41. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-AB-2/Mac-AB-2_2.fq.gz Mac-AB-r2_R2.fq.gz
    42. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-AB-3/Mac-AB-3_1.fq.gz Mac-AB-r3_R1.fq.gz
    43. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-AB-3/Mac-AB-3_2.fq.gz Mac-AB-r3_R2.fq.gz
    45. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-IJ-1/Mac-IJ-1_1.fq.gz Mac-IJ-r1_R1.fq.gz
    46. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-IJ-1/Mac-IJ-1_2.fq.gz Mac-IJ-r1_R2.fq.gz
    47. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-IJ-2/Mac-IJ-2_1.fq.gz Mac-IJ-r2_R1.fq.gz
    48. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-IJ-2/Mac-IJ-2_2.fq.gz Mac-IJ-r2_R2.fq.gz
    49. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-IJ-4/Mac-IJ-4_1.fq.gz Mac-IJ-r4_R1.fq.gz
    50. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-IJ-4/Mac-IJ-4_2.fq.gz Mac-IJ-r4_R2.fq.gz
    52. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-W1-1/Mac-W1-1_1.fq.gz Mac-W1-r1_R1.fq.gz
    53. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-W1-1/Mac-W1-1_2.fq.gz Mac-W1-r1_R2.fq.gz
    54. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-W1-2/Mac-W1-2_1.fq.gz Mac-W1-r2_R1.fq.gz
    55. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-W1-2/Mac-W1-2_2.fq.gz Mac-W1-r2_R2.fq.gz
    56. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-W1-3/Mac-W1-3_1.fq.gz Mac-W1-r3_R1.fq.gz
    57. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-W1-3/Mac-W1-3_2.fq.gz Mac-W1-r3_R2.fq.gz
    59. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-WT19606-2/Mac-WT19606-2_1.fq.gz Mac-WT19606-r2_R1.fq.gz
    60. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-WT19606-2/Mac-WT19606-2_2.fq.gz Mac-WT19606-r2_R2.fq.gz
    61. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-WT19606-3/Mac-WT19606-3_1.fq.gz Mac-WT19606-r3_R1.fq.gz
    62. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-WT19606-3/Mac-WT19606-3_2.fq.gz Mac-WT19606-r3_R2.fq.gz
    63. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-WT19606-4/Mac-WT19606-4_1.fq.gz Mac-WT19606-r4_R1.fq.gz
    64. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-WT19606-4/Mac-WT19606-4_2.fq.gz Mac-WT19606-r4_R2.fq.gz
    66. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-Y1-2/Mac-Y1-2_1.fq.gz Mac-Y1-r2_R1.fq.gz
    67. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-Y1-2/Mac-Y1-2_2.fq.gz Mac-Y1-r2_R2.fq.gz
    68. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-Y1-3/Mac-Y1-3_1.fq.gz Mac-Y1-r3_R1.fq.gz
    69. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-Y1-3/Mac-Y1-3_2.fq.gz Mac-Y1-r3_R2.fq.gz
    70. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-Y1-4/Mac-Y1-4_1.fq.gz Mac-Y1-r4_R1.fq.gz
    71. ln -s ../RSMR00204/X101SC25015922-Z02/X101SC25015922-Z02-J002/01.RawData/Mac-Y1-4/Mac-Y1-4_2.fq.gz Mac-Y1-r4_R2.fq.gz

  4. Preparing the directory trimmed

    1. mkdir trimmed trimmed_unpaired;
    2. for sample_id in LB-AB-r1 LB-AB-r2 LB-AB-r3  LB-IJ-r1 LB-IJ-r2 LB-IJ-r4  LB-W1-r1 LB-W1-r2 LB-W1-r3  LB-WT19606-r2 LB-WT19606-r3 LB-WT19606-r4  LB-Y1-r2 LB-Y1-r3 LB-Y1-r4    Mac-AB-r1 Mac-AB-r2 Mac-AB-r3  Mac-IJ-r1 Mac-IJ-r2 Mac-IJ-r4  Mac-W1-r1 Mac-W1-r2 Mac-W1-r3  Mac-WT19606-r2 Mac-WT19606-r3 Mac-WT19606-r4  Mac-Y1-r2 Mac-Y1-r3 Mac-Y1-r4; do
    3.         java -jar /home/jhuang/Tools/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 100 raw_data/${sample_id}_R1.fq.gz raw_data/${sample_id}_R2.fq.gz trimmed/${sample_id}_R1.fq.gz trimmed_unpaired/${sample_id}_R1.fq.gz trimmed/${sample_id}_R2.fq.gz trimmed_unpaired/${sample_id}_R2.fq.gz ILLUMINACLIP:/home/jhuang/Tools/Trimmomatic-0.36/adapters/TruSeq3-PE-2.fa:2:30:10:8:TRUE LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36 AVGQUAL:20; done 2> trimmomatic_pe.log;
    4. done

  5. Preparing samplesheet.csv

    1. sample,fastq_1,fastq_2,strandedness
    2. LB-AB-r1,LB-AB-r1_R1.fq.gz,LB-AB-r1_R2.fq.gz,auto
    3. LB-AB-r2,LB-AB-r2_R1.fq.gz,LB-AB-r2_R2.fq.gz,auto
    4. LB-AB-r3,LB-AB-r3_R1.fq.gz,LB-AB-r3_R2.fq.gz,auto
    5. LB-IJ-r1,LB-IJ-r1_R1.fq.gz,LB-IJ-r1_R2.fq.gz,auto
    6. LB-IJ-r2,LB-IJ-r2_R1.fq.gz,LB-IJ-r2_R2.fq.gz,auto
    7. LB-IJ-r4,LB-IJ-r4_R1.fq.gz,LB-IJ-r4_R2.fq.gz,auto
    8. LB-W1-r1,LB-W1-r1_R1.fq.gz,LB-W1-r1_R2.fq.gz,auto
    9. LB-W1-r2,LB-W1-r2_R1.fq.gz,LB-W1-r2_R2.fq.gz,auto
    10. LB-W1-r3,LB-W1-r3_R1.fq.gz,LB-W1-r3_R2.fq.gz,auto
    11. LB-WT19606-r2,LB-WT19606-r2_R1.fq.gz,LB-WT19606-r2_R2.fq.gz,auto
    12. LB-WT19606-r3,LB-WT19606-r3_R1.fq.gz,LB-WT19606-r3_R2.fq.gz,auto
    13. LB-WT19606-r4,LB-WT19606-r4_R1.fq.gz,LB-WT19606-r4_R2.fq.gz,auto
    14. LB-Y1-r2,LB-Y1-r2_R1.fq.gz,LB-Y1-r2_R2.fq.gz,auto
    15. LB-Y1-r3,LB-Y1-r3_R1.fq.gz,LB-Y1-r3_R2.fq.gz,auto
    16. LB-Y1-r4,LB-Y1-r4_R1.fq.gz,LB-Y1-r4_R2.fq.gz,auto
    17. Mac-AB-r1,Mac-AB-r1_R1.fq.gz,Mac-AB-r1_R2.fq.gz,auto
    18. Mac-AB-r2,Mac-AB-r2_R1.fq.gz,Mac-AB-r2_R2.fq.gz,auto
    19. Mac-AB-r3,Mac-AB-r3_R1.fq.gz,Mac-AB-r3_R2.fq.gz,auto
    20. Mac-IJ-r1,Mac-IJ-r1_R1.fq.gz,Mac-IJ-r1_R2.fq.gz,auto
    21. Mac-IJ-r2,Mac-IJ-r2_R1.fq.gz,Mac-IJ-r2_R2.fq.gz,auto
    22. Mac-IJ-r4,Mac-IJ-r4_R1.fq.gz,Mac-IJ-r4_R2.fq.gz,auto
    23. Mac-W1-r1,Mac-W1-r1_R1.fq.gz,Mac-W1-r1_R2.fq.gz,auto
    24. Mac-W1-r2,Mac-W1-r2_R1.fq.gz,Mac-W1-r2_R2.fq.gz,auto
    25. Mac-W1-r3,Mac-W1-r3_R1.fq.gz,Mac-W1-r3_R2.fq.gz,auto
    26. Mac-WT19606-r2,Mac-WT19606-r2_R1.fq.gz,Mac-WT19606-r2_R2.fq.gz,auto
    27. Mac-WT19606-r3,Mac-WT19606-r3_R1.fq.gz,Mac-WT19606-r3_R2.fq.gz,auto
    28. Mac-WT19606-r4,Mac-WT19606-r4_R1.fq.gz,Mac-WT19606-r4_R2.fq.gz,auto
    29. Mac-Y1-r2,Mac-Y1-r2_R1.fq.gz,Mac-Y1-r2_R2.fq.gz,auto
    30. Mac-Y1-r3,Mac-Y1-r3_R1.fq.gz,Mac-Y1-r3_R2.fq.gz,auto
    31. Mac-Y1-r4,Mac-Y1-r4_R1.fq.gz,Mac-Y1-r4_R2.fq.gz,auto
    33. #mv trimmed/* .

  6. nextflow run

    1. #Example1: http://xgenes.com/article/article-content/157/prepare-virus-gtf-for-nextflow-run/
    2. #docker pull nfcore/rnaseq
    3. ln -s /home/jhuang/Tools/nf-core-rnaseq-3.12.0/ rnaseq
    5. # ---- SUCCESSFUL with directly downloaded gff3 and fasta from NCBI using docker after replacing 'CDS' with 'exon' ----
    6. (host_env) /usr/local/bin/nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results    --fasta "/home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040.fasta" --gff "/home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_m.gff"        -profile docker -resume  --max_cpus 55 --max_memory 512.GB --max_time 2400.h    --save_align_intermeds --save_unaligned --save_reference    --aligner 'star_salmon'    --gtf_group_features 'gene_id'  --gtf_extra_attributes 'gene_name' --featurecounts_group_type 'gene_biotype' --featurecounts_feature_type 'transcript'

  7. Import data and pca-plot

    1.  #mamba activate r_env
    3. #install.packages("ggfun")
    4. # Import the required libraries
    5. library("AnnotationDbi")
    6. library("clusterProfiler")
    7. library("ReactomePA")
    8. library(gplots)
    9. library(tximport)
    10. library(DESeq2)
    11. #library("org.Hs.eg.db")
    12. library(dplyr)
    13. library(tidyverse)
    14. #install.packages("devtools")
    15. #devtools::install_version("gtable", version = "0.3.0")
    16. library(gplots)
    17. library("RColorBrewer")
    18. #install.packages("ggrepel")
    19. library("ggrepel")
    20. # install.packages("openxlsx")
    21. library(openxlsx)
    22. library(EnhancedVolcano)
    23. library(DESeq2)
    24. library(edgeR)
    26. setwd("~/DATA/Data_Tam_RNAseq_2025_LB_vs_Mac_ATCC19606/results/star_salmon")
    27. # Define paths to your Salmon output quantification files
    29. files <- c("LB-AB_r1" = "./LB-AB-r1/quant.sf",
    30.         "LB-AB_r2" = "./LB-AB-r2/quant.sf",
    31.         "LB-AB_r3" = "./LB-AB-r3/quant.sf",
    32.         "LB-IJ_r1" = "./LB-IJ-r1/quant.sf",
    33.         "LB-IJ_r2" = "./LB-IJ-r2/quant.sf",
    34.         "LB-IJ_r4" = "./LB-IJ-r4/quant.sf",
    35.         "LB-W1_r1" = "./LB-W1-r1/quant.sf",
    36.         "LB-W1_r2" = "./LB-W1-r2/quant.sf",
    37.         "LB-W1_r3" = "./LB-W1-r3/quant.sf",
    38.         "LB-WT19606_r2" = "./LB-WT19606-r2/quant.sf",
    39.         "LB-WT19606_r3" = "./LB-WT19606-r3/quant.sf",
    40.         "LB-WT19606_r4" = "./LB-WT19606-r4/quant.sf",
    41.         "LB-Y1_r2" = "./LB-Y1-r2/quant.sf",
    42.         "LB-Y1_r3" = "./LB-Y1-r3/quant.sf",
    43.         "LB-Y1_r4" = "./LB-Y1-r4/quant.sf",
    44.         "Mac-AB_r1" = "./Mac-AB-r1/quant.sf",
    45.         "Mac-AB_r2" = "./Mac-AB-r2/quant.sf",
    46.         "Mac-AB_r3" = "./Mac-AB-r3/quant.sf",
    47.         "Mac-IJ_r1" = "./Mac-IJ-r1/quant.sf",
    48.         "Mac-IJ_r2" = "./Mac-IJ-r2/quant.sf",
    49.         "Mac-IJ_r4" = "./Mac-IJ-r4/quant.sf",
    50.         "Mac-W1_r1" = "./Mac-W1-r1/quant.sf",
    51.         "Mac-W1_r2" = "./Mac-W1-r2/quant.sf",
    52.         "Mac-W1_r3" = "./Mac-W1-r3/quant.sf",
    53.         "Mac-WT19606_r2" = "./Mac-WT19606-r2/quant.sf",
    54.         "Mac-WT19606_r3" = "./Mac-WT19606-r3/quant.sf",
    55.         "Mac-WT19606_r4" = "./Mac-WT19606-r4/quant.sf",
    56.         "Mac-Y1_r2" = "./Mac-Y1-r2/quant.sf",
    57.         "Mac-Y1_r3" = "./Mac-Y1-r3/quant.sf",
    58.         "Mac-Y1_r4" = "./Mac-Y1-r4/quant.sf")
    60. # Import the transcript abundance data with tximport
    61. txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE)
    62. # Define the replicates and condition of the samples
    63. #replicate <- factor(c("r1", "r2", "r3", "r1", "r2", "r3", "r1", "r2", "r3"))
    64. #adeA and adeB encode a membrane fusion protein that is part of the AdeABC efflux pump, which contributes to multidrug resistance.
    65. #System: Part of the AdeIJK efflux pump, which includes: adeI — membrane fusion protein, adeJ — RND transporter, adeK — outer membrane factor
    66. condition <- factor(c("LB-AB","LB-AB","LB-AB", "LB-IJ","LB-IJ","LB-IJ", "LB-W1","LB-W1","LB-W1","LB-WT19606","LB-WT19606","LB-WT19606","LB-Y1","LB-Y1","LB-Y1","Mac-AB","Mac-AB","Mac-AB","Mac-IJ","Mac-IJ","Mac-IJ","Mac-W1","Mac-W1","Mac-W1","Mac-WT19606","Mac-WT19606","Mac-WT19606","Mac-Y1","Mac-Y1","Mac-Y1"))
    67. # Define the colData for DESeq2
    68. colData <- data.frame(condition=condition, row.names=names(files))
    70. # ------------------------
    71. # 1️⃣ Setup and input files
    72. # ------------------------
    74. # Read in transcript-to-gene mapping
    75. tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE)
    76. colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name")
    78. # Prepare tx2gene for gene-level summarization (remove gene_name if needed)
    79. tx2gene_geneonly <- tx2gene[, c("transcript_id", "gene_id")]
    81. # -------------------------------
    82. # 2️⃣ Transcript-level counts
    83. # -------------------------------
    84. # Create DESeqDataSet directly from tximport (transcript-level)
    85. dds_tx <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
    86. write.csv(counts(dds_tx), file="transcript_counts.csv")
    88. # --------------------------------
    89. # 3️⃣ Gene-level summarization
    90. # --------------------------------
    91. # Re-import Salmon data summarized at gene level
    92. txi_gene <- tximport(files, type="salmon", tx2gene=tx2gene_geneonly, txOut=FALSE)
    94. # Create DESeqDataSet for gene-level counts
    95. #dds <- DESeqDataSetFromTximport(txi_gene, colData=colData, design=~condition+replicate)
    96. dds <- DESeqDataSetFromTximport(txi_gene, colData=colData, design=~condition)
    98. # --------------------------------
    99. # 4️⃣ Raw counts table (with gene names)
    100. # --------------------------------
    101. # Extract raw gene-level counts
    102. counts_data <- as.data.frame(counts(dds, normalized=FALSE))
    103. counts_data$gene_id <- rownames(counts_data)
    105. # Add gene names
    106. tx2gene_unique <- unique(tx2gene[, c("gene_id", "gene_name")])
    107. counts_data <- merge(counts_data, tx2gene_unique, by="gene_id", all.x=TRUE)
    109. # Reorder columns: gene_id, gene_name, then counts
    110. count_cols <- setdiff(colnames(counts_data), c("gene_id", "gene_name"))
    111. counts_data <- counts_data[, c("gene_id", "gene_name", count_cols)]
    113. # --------------------------------
    114. # 5️⃣ Calculate CPM
    115. # --------------------------------
    116. library(edgeR)
    117. library(openxlsx)
    119. # Prepare count matrix for CPM calculation
    120. count_matrix <- as.matrix(counts_data[, !(colnames(counts_data) %in% c("gene_id", "gene_name"))])
    122. # Calculate CPM
    123. #cpm_matrix <- cpm(count_matrix, normalized.lib.sizes=FALSE)
    124. total_counts <- colSums(count_matrix)
    125. cpm_matrix <- t(t(count_matrix) / total_counts) * 1e6
    126. cpm_matrix <- as.data.frame(cpm_matrix)
    128. # Add gene_id and gene_name back to CPM table
    129. cpm_counts <- cbind(counts_data[, c("gene_id", "gene_name")], cpm_matrix)
    131. # --------------------------------
    132. # 6️⃣ Save outputs
    133. # --------------------------------
    134. write.csv(counts_data, "gene_raw_counts.csv", row.names=FALSE)
    135. write.xlsx(counts_data, "gene_raw_counts.xlsx", row.names=FALSE)
    136. write.xlsx(cpm_counts, "gene_cpm_counts.xlsx", row.names=FALSE)
    138. # -- (Optional) Save the rlog-transformed counts --
    139. dim(counts(dds))
    140. head(counts(dds), 10)
    141. rld <- rlogTransformation(dds)
    142. rlog_counts <- assay(rld)
    143. write.xlsx(as.data.frame(rlog_counts), "gene_rlog_transformed_counts.xlsx")
    145. # ---- (Optional for NACHREIHEN) split the factos media and strain from condition (for comparison Mac vs LB) ----
    146. # AdeIJK vs. AdeABC Efflux Pumps
    147. #     * AdeIJK is the "housekeeping" pump — always active, broadly expressed, contributing to background resistance.
    148. #     * AdeABC is the "emergency" pump — induced under stress or mutations, more potent in contributing to clinical multidrug resistance.
    149. #LB = Luria-Bertani broth (a standard rich growth medium)
    150. #Mac = MacConkey agar or broth (selective for Gram-negative bacteria)
    151. # - Growth medium   Media or Condition, GrowthMedium
    152. # - Bacterial strain/genotype   Strain or Isolate, Genotype, SampleType
    153. media <- factor(c("LB","LB","LB", "LB","LB","LB", "LB","LB","LB","LB","LB","LB","LB","LB","LB","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac","Mac"))
    154. strain <- factor(c("AB","AB","AB", "IJ","IJ","IJ", "W1","W1","W1","WT19606","WT19606","WT19606","Y1","Y1","Y1","AB","AB","AB","IJ","IJ","IJ","W1","W1","W1","WT19606","WT19606","WT19606","Y1","Y1","Y1"))
    155. # Define the colData for DESeq2
    156. colData <- data.frame(media=media, strain=strain, row.names=names(files))
    157. # -- transcript-level count data (x2) --
    158. # Create DESeqDataSet object
    159. dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~media+strain)
    160. #write.csv(counts(dds), file="transcript_counts_media_strain.csv")  #check correctness, it should be identical to transcript_counts.csv
    161. # -- gene-level count data (x2) --
    162. # Read in the tx2gene map from salmon_tx2gene.tsv
    163. tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE)
    164. # Set the column names
    165. colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name")
    166. # Remove the gene_name column if not needed
    167. tx2gene <- tx2gene[,1:2]
    168. # Import and summarize the Salmon data with tximport
    169. txi <- tximport(files, type = "salmon", tx2gene = tx2gene, txOut = FALSE)
    170. # Continue with the DESeq2 workflow as before...
    171. colData <- data.frame(media=media, strain=strain, row.names=names(files))
    172. dds <- DESeqDataSetFromTximport(txi, colData=colData, design=~media+strain)
    173. #dds <- dds[rowSums(counts(dds) > 3) > 2, ]    #3796->????
    174. #write.csv(counts(dds, normalized=FALSE), file="gene_counts_media_strain.csv")  #check correctness, it should be identical to gene_counts.csv
    175. # ---- (Optional for NACHREIHEN) END ----
    177. # -- pca --
    178. png("pca2.png", 1200, 800)
    179. plotPCA(rld, intgroup=c("condition"))
    180. dev.off()
    181. # -- heatmap --
    182. png("heatmap2.png", 1200, 800)
    183. distsRL <- dist(t(assay(rld)))
    184. mat <- as.matrix(distsRL)
    185. hc <- hclust(distsRL)
    186. hmcol <- colorRampPalette(brewer.pal(9,"GnBu"))(100)
    187. heatmap.2(mat, Rowv=as.dendrogram(hc),symm=TRUE, trace="none",col = rev(hmcol), margin=c(13, 13))
    188. dev.off()
    190. # -- pca_media_strain --
    191. png("pca_media.png", 1200, 800)
    192. plotPCA(rld, intgroup=c("media"))
    193. dev.off()
    194. png("pca_strain.png", 1200, 800)
    195. plotPCA(rld, intgroup=c("strain"))
    196. dev.off()

  8. (Optional; ERROR-->need to be debugged!) ) estimate size factors and dispersion values.

    1. #Size Factors: These are used to normalize the read counts across different samples. The size factor for a sample accounts for differences in sequencing depth (i.e., the total number of reads) and other technical biases between samples. After normalization with size factors, the counts should be comparable across samples. Size factors are usually calculated in a way that they reflect the median or mean ratio of gene expression levels between samples, assuming that most genes are not differentially expressed.
    2. #Dispersion: This refers to the variability or spread of gene expression measurements. In RNA-seq data analysis, each gene has its own dispersion value, which reflects how much the counts for that gene vary between different samples, more than what would be expected just due to the Poisson variation inherent in counting. Dispersion is important for accurately modeling the data and for detecting differentially expressed genes.
    3. #So in summary, size factors are specific to samples (used to make counts comparable across samples), and dispersion values are specific to genes (reflecting variability in gene expression).
    5. sizeFactors(dds)
    6. #NULL
    7. # Estimate size factors
    8. dds <- estimateSizeFactors(dds)
    9. # Estimate dispersions
    10. dds <- estimateDispersions(dds)
    11. #> sizeFactors(dds)
    13. #control_r1 control_r2  HSV.d2_r1  HSV.d2_r2  HSV.d4_r1  HSV.d4_r2  HSV.d6_r1
    14. #2.3282468  2.0251928  1.8036883  1.3767551  0.9341929  1.0911693  0.5454526
    15. #HSV.d6_r2  HSV.d8_r1  HSV.d8_r2
    16. #0.4604461  0.5799834  0.6803681
    18. # (DEBUG) If avgTxLength is Necessary
    19. #To simplify the computation and ensure sizeFactors are calculated:
    20. assays(dds)$avgTxLength <- NULL
    21. dds <- estimateSizeFactors(dds)
    22. sizeFactors(dds)
    23. #If you want to retain avgTxLength but suspect it is causing issues, you can explicitly instruct DESeq2 to compute size factors without correcting for library size with average transcript lengths:
    24. dds <- estimateSizeFactors(dds, controlGenes = NULL, use = FALSE)
    25. sizeFactors(dds)
    27. # If alone with virus data, the following BUG occured:
    28. #Still NULL --> BUG --> using manual calculation method for sizeFactor calculation!
    29.                     HeLa_TO_r1                      HeLa_TO_r2
    30.                     0.9978755                       1.1092227
    31. data.frame(genes = rownames(dds), dispersions = dispersions(dds))
    33. #Given the raw counts, the control_r1 and control_r2 samples seem to have a much lower sequencing depth (total read count) than the other samples. Therefore, when normalization methods are applied, the normalization factors for these control samples will be relatively high, boosting the normalized counts.
    34. 1/0.9978755=1.002129023
    35. 1/1.1092227=
    36. #bamCoverage --bam ../markDuplicates/${sample}Aligned.sortedByCoord.out.bam -o ${sample}_norm.bw --binSize 10 --scaleFactor  --effectiveGenomeSize 2864785220
    37. bamCoverage --bam ../markDuplicates/HeLa_TO_r1Aligned.sortedByCoord.out.markDups.bam -o HeLa_TO_r1.bw --binSize 10 --scaleFactor 1.002129023     --effectiveGenomeSize 2864785220
    38. bamCoverage --bam ../markDuplicates/HeLa_TO_r2Aligned.sortedByCoord.out.markDups.bam -o HeLa_TO_r2.bw --binSize 10 --scaleFactor  0.901532217        --effectiveGenomeSize 2864785220
    40. raw_counts <- counts(dds)
    41. normalized_counts <- counts(dds, normalized=TRUE)
    42. #write.table(raw_counts, file="raw_counts.txt", sep="\t", quote=F, col.names=NA)
    43. #write.table(normalized_counts, file="normalized_counts.txt", sep="\t", quote=F, col.names=NA)
    44. #convert bam to bigwig using deepTools by feeding inverse of DESeq’s size Factor
    45. estimSf <- function (cds){
    46.     # Get the count matrix
    47.     cts <- counts(cds)
    48.     # Compute the geometric mean
    49.     geomMean <- function(x) prod(x)^(1/length(x))
    50.     # Compute the geometric mean over the line
    51.     gm.mean  <-  apply(cts, 1, geomMean)
    52.     # Zero values are set to NA (avoid subsequentcdsdivision by 0)
    53.     gm.mean[gm.mean == 0] <- NA
    54.     # Divide each line by its corresponding geometric mean
    55.     # sweep(x, MARGIN, STATS, FUN = "-", check.margin = TRUE, ...)
    56.     # MARGIN: 1 or 2 (line or columns)
    57.     # STATS: a vector of length nrow(x) or ncol(x), depending on MARGIN
    58.     # FUN: the function to be applied
    59.     cts <- sweep(cts, 1, gm.mean, FUN="/")
    60.     # Compute the median over the columns
    61.     med <- apply(cts, 2, median, na.rm=TRUE)
    62.     # Return the scaling factor
    63.     return(med)
    64. }
    65. #https://dputhier.github.io/ASG/practicals/rnaseq_diff_Snf2/rnaseq_diff_Snf2.html
    66. #http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#data-transformations-and-visualization
    67. #https://hbctraining.github.io/DGE_workshop/lessons/02_DGE_count_normalization.html
    68. #https://hbctraining.github.io/DGE_workshop/lessons/04_DGE_DESeq2_analysis.html
    69. #https://genviz.org/module-04-expression/0004/02/01/DifferentialExpression/
    70. #DESeq2’s median of ratios [1]
    71. #EdgeR’s trimmed mean of M values (TMM) [2]
    72. #http://www.nathalievialaneix.eu/doc/html/TP1_normalization.html  #very good website!
    73. test_normcount <- sweep(raw_counts, 2, sizeFactors(dds), "/")
    74. sum(test_normcount != normalized_counts)

  9. Select the differentially expressed genes

    1. #https://galaxyproject.eu/posts/2020/08/22/three-steps-to-galaxify-your-tool/
    2. #https://www.biostars.org/p/282295/
    3. #https://www.biostars.org/p/335751/
    4. #> dds$condition
    5. #LB-AB       LB-IJ       LB-W1       LB-WT19606  LB-Y1       Mac-AB     Mac-IJ      Mac-W1      Mac-WT19606 Mac-Y1
    6. #CONSOLE: mkdir star_salmon/degenes
    8. setwd("degenes")
    9. #---- relevel to control ----
    10. dds$condition <- relevel(dds$condition, "LB-WT19606")
    11. dds = DESeq(dds, betaPrior=FALSE)
    12. resultsNames(dds)
    13. clist <- c("LB.AB_vs_LB.WT19606","LB.IJ_vs_LB.WT19606","LB.W1_vs_LB.WT19606","LB.Y1_vs_LB.WT19606")
    15. dds$condition <- relevel(dds$condition, "Mac-WT19606")
    16. dds = DESeq(dds, betaPrior=FALSE)
    17. resultsNames(dds)
    18. clist <- c("Mac.AB_vs_Mac.WT19606","Mac.IJ_vs_Mac.WT19606","Mac.W1_vs_Mac.WT19606","Mac.Y1_vs_Mac.WT19606")
    20. # - 如果你的实验是关注细菌在没有选择性压力下的生长、基因表达或一般行为，LB 是更好的对照。
    21. # - 如果你希望研究细菌在选择性压力下的行为（例如，针对革兰氏阴性细菌、测试抗生素耐药性或区分乳糖发酵菌），那么 MacConkey 更适合作为对照。
    22. dds$media <- relevel(dds$media, "LB")
    23. dds = DESeq(dds, betaPrior=FALSE)
    24. resultsNames(dds)
    25. clist <- c("Mac_vs_LB")
    27. dds$media <- relevel(dds$media, "Mac")
    28. dds = DESeq(dds, betaPrior=FALSE)
    29. resultsNames(dds)
    30. clist <- c("LB_vs_Mac")
    32. for (i in clist) {
    33.   #contrast = paste("condition", i, sep="_")
    34.   contrast = paste("media", i, sep="_")
    35.   res = results(dds, name=contrast)
    36.   res <- res[!is.na(res$log2FoldChange),]
    37.   res_df <- as.data.frame(res)
    39.   write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-"))
    40.   up <- subset(res_df, padj<=0.05 & log2FoldChange>=2)
    41.   down <- subset(res_df, padj<=0.05 & log2FoldChange<=-2)
    42.   write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-"))
    43.   write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-"))
    44. }
    46. # -- Under host-env --
    47. grep -P "\tgene\t" CP059040.gff > CP059040_gene.gff
    48. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.AB_vs_LB.WT19606-all.txt LB.AB_vs_LB.WT19606-all.csv
    49. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.AB_vs_LB.WT19606-up.txt LB.AB_vs_LB.WT19606-up.csv
    50. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.AB_vs_LB.WT19606-down.txt LB.AB_vs_LB.WT19606-down.csv
    51. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.IJ_vs_LB.WT19606-all.txt LB.IJ_vs_LB.WT19606-all.csv
    52. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.IJ_vs_LB.WT19606-up.txt LB.IJ_vs_LB.WT19606-up.csv
    53. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.IJ_vs_LB.WT19606-down.txt LB.IJ_vs_LB.WT19606-down.csv
    54. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.W1_vs_LB.WT19606-all.txt LB.W1_vs_LB.WT19606-all.csv
    55. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.W1_vs_LB.WT19606-up.txt LB.W1_vs_LB.WT19606-up.csv
    56. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.W1_vs_LB.WT19606-down.txt LB.W1_vs_LB.WT19606-down.csv
    57. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.Y1_vs_LB.WT19606-all.txt LB.Y1_vs_LB.WT19606-all.csv
    58. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.Y1_vs_LB.WT19606-up.txt LB.Y1_vs_LB.WT19606-up.csv
    59. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB.Y1_vs_LB.WT19606-down.txt LB.Y1_vs_LB.WT19606-down.csv
    60. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.AB_vs_Mac.WT19606-all.txt Mac.AB_vs_Mac.WT19606-all.csv
    61. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.AB_vs_Mac.WT19606-up.txt Mac.AB_vs_Mac.WT19606-up.csv
    62. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.AB_vs_Mac.WT19606-down.txt Mac.AB_vs_Mac.WT19606-down.csv
    63. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.IJ_vs_Mac.WT19606-all.txt Mac.IJ_vs_Mac.WT19606-all.csv
    64. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.IJ_vs_Mac.WT19606-up.txt Mac.IJ_vs_Mac.WT19606-up.csv
    65. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.IJ_vs_Mac.WT19606-down.txt Mac.IJ_vs_Mac.WT19606-down.csv
    66. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.W1_vs_Mac.WT19606-all.txt Mac.W1_vs_Mac.WT19606-all.csv
    67. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.W1_vs_Mac.WT19606-up.txt Mac.W1_vs_Mac.WT19606-up.csv
    68. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.W1_vs_Mac.WT19606-down.txt Mac.W1_vs_Mac.WT19606-down.csv
    69. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.Y1_vs_Mac.WT19606-all.txt Mac.Y1_vs_Mac.WT19606-all.csv
    70. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.Y1_vs_Mac.WT19606-up.txt Mac.Y1_vs_Mac.WT19606-up.csv
    71. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac.Y1_vs_Mac.WT19606-down.txt Mac.Y1_vs_Mac.WT19606-down.csv
    72. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac_vs_LB-all.txt Mac_vs_LB-all.csv
    73. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac_vs_LB-up.txt Mac_vs_LB-up.csv
    74. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff Mac_vs_LB-down.txt Mac_vs_LB-down.csv
    75. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB_vs_Mac-all.txt LB_vs_Mac-all.csv
    76. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB_vs_Mac-up.txt LB_vs_Mac-up.csv
    77. python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine/CP059040_gene.gff LB_vs_Mac-down.txt LB_vs_Mac-down.csv
    79. # ---- Mac_vs_LB ----
    80. res <- read.csv("Mac_vs_LB-all.csv")
    81. # Replace empty GeneName with modified GeneID
    82. res$GeneName <- ifelse(
    83.   res$GeneName == "" | is.na(res$GeneName),
    84.   gsub("gene-", "", res$GeneID),
    85.   res$GeneName
    86. )
    87. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    88. #print(duplicated_genes)
    89. # [1] "bfr"  "lipA" "ahpF" "pcaF" "alr"  "pcaD" "cydB" "lpdA" "pgaC" "ppk1"
    90. #[11] "pcaF" "tuf"  "galE" "murI" "yccS" "rrf"  "rrf"  "arsB" "ptsP" "umuD"
    91. #[21] "map"  "pgaB" "rrf"  "rrf"  "rrf"  "pgaD" "uraH" "benE"
    92. #res[res$GeneName == "bfr", ]
    94. #1st_strategy First occurrence is kept and Subsequent duplicates are removed
    95. #res <- res[!duplicated(res$GeneName), ]
    96. #2nd_strategy keep the row with the smallest padj value for each GeneName
    97. res <- res %>%
    98.   group_by(GeneName) %>%
    99.   slice_min(padj, with_ties = FALSE) %>%
    100.   ungroup()
    101. res <- as.data.frame(res)
    102. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    103. res <- res[order(res$padj, -res$log2FoldChange), ]
    105. # Assuming res is your dataframe and already processed
    106. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    107. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    108. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    109. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    110. # Create a new workbook
    111. wb <- createWorkbook()
    112. # Add the complete dataset as the first sheet
    113. addWorksheet(wb, "Complete_Data")
    114. writeData(wb, "Complete_Data", res)
    115. # Add the up-regulated genes as the second sheet
    116. addWorksheet(wb, "Up_Regulated")
    117. writeData(wb, "Up_Regulated", up_regulated)
    118. # Add the down-regulated genes as the third sheet
    119. addWorksheet(wb, "Down_Regulated")
    120. writeData(wb, "Down_Regulated", down_regulated)
    121. # Save the workbook to a file
    122. saveWorkbook(wb, "Gene_Expression_Mac_vs_LB.xlsx", overwrite = TRUE)
    124. # Set the 'GeneName' column as row.names
    125. rownames(res) <- res$GeneName
    126. # Drop the 'GeneName' column since it's now the row names
    127. res$GeneName <- NULL
    128. head(res)
    130. ## Ensure the data frame matches the expected format
    131. ## For example, it should have columns: log2FoldChange, padj, etc.
    132. #res <- as.data.frame(res)
    133. ## Remove rows with NA in log2FoldChange (if needed)
    134. #res <- res[!is.na(res$log2FoldChange),]
    136. # Replace padj = 0 with a small value
    137. res$padj[res$padj == 0] <- 1e-150
    139. #library(EnhancedVolcano)
    140. # Assuming res is already sorted and processed
    141. png("Mac_vs_LB.png", width=1200, height=2000)
    142. #max.overlaps = 10
    143. EnhancedVolcano(res,
    144.                 lab = rownames(res),
    145.                 x = 'log2FoldChange',
    146.                 y = 'padj',
    147.                 pCutoff = 1e-2,
    148.                 FCcutoff = 2,
    149.                 title = '',
    150.                 subtitleLabSize = 18,
    151.                 pointSize = 3.0,
    152.                 labSize = 5.0,
    153.                 colAlpha = 1,
    154.                 legendIconSize = 4.0,
    155.                 drawConnectors = TRUE,
    156.                 widthConnectors = 0.5,
    157.                 colConnectors = 'black',
    158.                 subtitle = expression("Mac versus LB"))
    159. dev.off()
    161. # ---- LB.AB_vs_LB.WT19606 ----
    162. res <- read.csv("LB.AB_vs_LB.WT19606-all.csv")
    163. # Replace empty GeneName with modified GeneID
    164. res$GeneName <- ifelse(
    165.   res$GeneName == "" | is.na(res$GeneName),
    166.   gsub("gene-", "", res$GeneID),
    167.   res$GeneName
    168. )
    169. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    171. res <- res %>%
    172.   group_by(GeneName) %>%
    173.   slice_min(padj, with_ties = FALSE) %>%
    174.   ungroup()
    175. res <- as.data.frame(res)
    176. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    177. res <- res[order(res$padj, -res$log2FoldChange), ]
    179. # Assuming res is your dataframe and already processed
    180. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    181. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    182. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    183. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    184. # Create a new workbook
    185. wb <- createWorkbook()
    186. # Add the complete dataset as the first sheet
    187. addWorksheet(wb, "Complete_Data")
    188. writeData(wb, "Complete_Data", res)
    189. # Add the up-regulated genes as the second sheet
    190. addWorksheet(wb, "Up_Regulated")
    191. writeData(wb, "Up_Regulated", up_regulated)
    192. # Add the down-regulated genes as the third sheet
    193. addWorksheet(wb, "Down_Regulated")
    194. writeData(wb, "Down_Regulated", down_regulated)
    195. # Save the workbook to a file
    196. saveWorkbook(wb, "Gene_Expression_LB.AB_vs_LB.WT19606.xlsx", overwrite = TRUE)
    198. # Set the 'GeneName' column as row.names
    199. rownames(res) <- res$GeneName
    200. # Drop the 'GeneName' column since it's now the row names
    201. res$GeneName <- NULL
    202. head(res)
    204. ## Ensure the data frame matches the expected format
    205. ## For example, it should have columns: log2FoldChange, padj, etc.
    206. #res <- as.data.frame(res)
    207. ## Remove rows with NA in log2FoldChange (if needed)
    208. #res <- res[!is.na(res$log2FoldChange),]
    210. # Replace padj = 0 with a small value
    211. res$padj[res$padj == 0] <- 1e-12
    213. #library(EnhancedVolcano)
    214. # Assuming res is already sorted and processed
    215. png("LB.AB_vs_LB.WT19606.png", width=1200, height=1200)
    216. #max.overlaps = 10
    217. EnhancedVolcano(res,
    218.                 lab = rownames(res),
    219.                 x = 'log2FoldChange',
    220.                 y = 'padj',
    221.                 pCutoff = 1e-2,
    222.                 FCcutoff = 2,
    223.                 title = '',
    224.                 subtitleLabSize = 18,
    225.                 pointSize = 3.0,
    226.                 labSize = 5.0,
    227.                 colAlpha = 1,
    228.                 legendIconSize = 4.0,
    229.                 drawConnectors = TRUE,
    230.                 widthConnectors = 0.5,
    231.                 colConnectors = 'black',
    232.                 subtitle = expression("LB.AB versus LB.WT19606"))
    233. dev.off()
    235. # ---- LB.IJ_vs_LB.WT19606 ----
    236. res <- read.csv("LB.IJ_vs_LB.WT19606-all.csv")
    237. # Replace empty GeneName with modified GeneID
    238. res$GeneName <- ifelse(
    239.   res$GeneName == "" | is.na(res$GeneName),
    240.   gsub("gene-", "", res$GeneID),
    241.   res$GeneName
    242. )
    243. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    245. res <- res %>%
    246.   group_by(GeneName) %>%
    247.   slice_min(padj, with_ties = FALSE) %>%
    248.   ungroup()
    249. res <- as.data.frame(res)
    250. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    251. res <- res[order(res$padj, -res$log2FoldChange), ]
    253. # Assuming res is your dataframe and already processed
    254. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    255. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    256. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    257. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    258. # Create a new workbook
    259. wb <- createWorkbook()
    260. # Add the complete dataset as the first sheet
    261. addWorksheet(wb, "Complete_Data")
    262. writeData(wb, "Complete_Data", res)
    263. # Add the up-regulated genes as the second sheet
    264. addWorksheet(wb, "Up_Regulated")
    265. writeData(wb, "Up_Regulated", up_regulated)
    266. # Add the down-regulated genes as the third sheet
    267. addWorksheet(wb, "Down_Regulated")
    268. writeData(wb, "Down_Regulated", down_regulated)
    269. # Save the workbook to a file
    270. saveWorkbook(wb, "Gene_Expression_LB.IJ_vs_LB.WT19606.xlsx", overwrite = TRUE)
    272. # Set the 'GeneName' column as row.names
    273. rownames(res) <- res$GeneName
    274. # Drop the 'GeneName' column since it's now the row names
    275. res$GeneName <- NULL
    276. head(res)
    278. ## Ensure the data frame matches the expected format
    279. ## For example, it should have columns: log2FoldChange, padj, etc.
    280. #res <- as.data.frame(res)
    281. ## Remove rows with NA in log2FoldChange (if needed)
    282. #res <- res[!is.na(res$log2FoldChange),]
    284. # Replace padj = 0 with a small value
    285. res$padj[res$padj == 0] <- 1e-12
    287. #library(EnhancedVolcano)
    288. # Assuming res is already sorted and processed
    289. png("LB.IJ_vs_LB.WT19606.png", width=1200, height=1200)
    290. #max.overlaps = 10
    291. EnhancedVolcano(res,
    292.                 lab = rownames(res),
    293.                 x = 'log2FoldChange',
    294.                 y = 'padj',
    295.                 pCutoff = 1e-2,
    296.                 FCcutoff = 2,
    297.                 title = '',
    298.                 subtitleLabSize = 18,
    299.                 pointSize = 3.0,
    300.                 labSize = 5.0,
    301.                 colAlpha = 1,
    302.                 legendIconSize = 4.0,
    303.                 drawConnectors = TRUE,
    304.                 widthConnectors = 0.5,
    305.                 colConnectors = 'black',
    306.                 subtitle = expression("LB.IJ versus LB.WT19606"))
    307. dev.off()
    309. # ---- LB.W1_vs_LB.WT19606 ----
    310. res <- read.csv("LB.W1_vs_LB.WT19606-all.csv")
    311. # Replace empty GeneName with modified GeneID
    312. res$GeneName <- ifelse(
    313.   res$GeneName == "" | is.na(res$GeneName),
    314.   gsub("gene-", "", res$GeneID),
    315.   res$GeneName
    316. )
    317. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    319. res <- res %>%
    320.   group_by(GeneName) %>%
    321.   slice_min(padj, with_ties = FALSE) %>%
    322.   ungroup()
    323. res <- as.data.frame(res)
    324. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    325. res <- res[order(res$padj, -res$log2FoldChange), ]
    327. # Assuming res is your dataframe and already processed
    328. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    329. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    330. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    331. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    332. # Create a new workbook
    333. wb <- createWorkbook()
    334. # Add the complete dataset as the first sheet
    335. addWorksheet(wb, "Complete_Data")
    336. writeData(wb, "Complete_Data", res)
    337. # Add the up-regulated genes as the second sheet
    338. addWorksheet(wb, "Up_Regulated")
    339. writeData(wb, "Up_Regulated", up_regulated)
    340. # Add the down-regulated genes as the third sheet
    341. addWorksheet(wb, "Down_Regulated")
    342. writeData(wb, "Down_Regulated", down_regulated)
    343. # Save the workbook to a file
    344. saveWorkbook(wb, "Gene_Expression_LB.W1_vs_LB.WT19606.xlsx", overwrite = TRUE)
    346. # Set the 'GeneName' column as row.names
    347. rownames(res) <- res$GeneName
    348. # Drop the 'GeneName' column since it's now the row names
    349. res$GeneName <- NULL
    350. head(res)
    352. ## Ensure the data frame matches the expected format
    353. ## For example, it should have columns: log2FoldChange, padj, etc.
    354. #res <- as.data.frame(res)
    355. ## Remove rows with NA in log2FoldChange (if needed)
    356. #res <- res[!is.na(res$log2FoldChange),]
    358. # Replace padj = 0 with a small value
    359. res$padj[res$padj == 0] <- 1e-12
    361. #library(EnhancedVolcano)
    362. # Assuming res is already sorted and processed
    363. png("LB.W1_vs_LB.WT19606.png", width=1200, height=1200)
    364. #max.overlaps = 10
    365. EnhancedVolcano(res,
    366.                 lab = rownames(res),
    367.                 x = 'log2FoldChange',
    368.                 y = 'padj',
    369.                 pCutoff = 1e-2,
    370.                 FCcutoff = 2,
    371.                 title = '',
    372.                 subtitleLabSize = 18,
    373.                 pointSize = 3.0,
    374.                 labSize = 5.0,
    375.                 colAlpha = 1,
    376.                 legendIconSize = 4.0,
    377.                 drawConnectors = TRUE,
    378.                 widthConnectors = 0.5,
    379.                 colConnectors = 'black',
    380.                 subtitle = expression("LB.W1 versus LB.WT19606"))
    381. dev.off()
    383. # ---- LB.Y1_vs_LB.WT19606 ----
    384. res <- read.csv("LB.Y1_vs_LB.WT19606-all.csv")
    385. # Replace empty GeneName with modified GeneID
    386. res$GeneName <- ifelse(
    387.   res$GeneName == "" | is.na(res$GeneName),
    388.   gsub("gene-", "", res$GeneID),
    389.   res$GeneName
    390. )
    391. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    393. res <- res %>%
    394.   group_by(GeneName) %>%
    395.   slice_min(padj, with_ties = FALSE) %>%
    396.   ungroup()
    397. res <- as.data.frame(res)
    398. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    399. res <- res[order(res$padj, -res$log2FoldChange), ]
    401. # Assuming res is your dataframe and already processed
    402. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    403. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    404. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    405. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    406. # Create a new workbook
    407. wb <- createWorkbook()
    408. # Add the complete dataset as the first sheet
    409. addWorksheet(wb, "Complete_Data")
    410. writeData(wb, "Complete_Data", res)
    411. # Add the up-regulated genes as the second sheet
    412. addWorksheet(wb, "Up_Regulated")
    413. writeData(wb, "Up_Regulated", up_regulated)
    414. # Add the down-regulated genes as the third sheet
    415. addWorksheet(wb, "Down_Regulated")
    416. writeData(wb, "Down_Regulated", down_regulated)
    417. # Save the workbook to a file
    418. saveWorkbook(wb, "Gene_Expression_LB.Y1_vs_LB.WT19606.xlsx", overwrite = TRUE)
    420. # Set the 'GeneName' column as row.names
    421. rownames(res) <- res$GeneName
    422. # Drop the 'GeneName' column since it's now the row names
    423. res$GeneName <- NULL
    424. head(res)
    426. ## Ensure the data frame matches the expected format
    427. ## For example, it should have columns: log2FoldChange, padj, etc.
    428. #res <- as.data.frame(res)
    429. ## Remove rows with NA in log2FoldChange (if needed)
    430. #res <- res[!is.na(res$log2FoldChange),]
    432. # Replace padj = 0 with a small value
    433. res$padj[res$padj == 0] <- 1e-12
    435. #library(EnhancedVolcano)
    436. # Assuming res is already sorted and processed
    437. png("LB.Y1_vs_LB.WT19606.png", width=1200, height=1200)
    438. #max.overlaps = 10
    439. EnhancedVolcano(res,
    440.                 lab = rownames(res),
    441.                 x = 'log2FoldChange',
    442.                 y = 'padj',
    443.                 pCutoff = 1e-2,
    444.                 FCcutoff = 2,
    445.                 title = '',
    446.                 subtitleLabSize = 18,
    447.                 pointSize = 3.0,
    448.                 labSize = 5.0,
    449.                 colAlpha = 1,
    450.                 legendIconSize = 4.0,
    451.                 drawConnectors = TRUE,
    452.                 widthConnectors = 0.5,
    453.                 colConnectors = 'black',
    454.                 subtitle = expression("LB.Y1 versus LB.WT19606"))
    455. dev.off()
    457. # ---- Mac.AB_vs_Mac.WT19606 ----
    458. res <- read.csv("Mac.AB_vs_Mac.WT19606-all.csv")
    459. # Replace empty GeneName with modified GeneID
    460. res$GeneName <- ifelse(
    461.   res$GeneName == "" | is.na(res$GeneName),
    462.   gsub("gene-", "", res$GeneID),
    463.   res$GeneName
    464. )
    465. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    467. res <- res %>%
    468.   group_by(GeneName) %>%
    469.   slice_min(padj, with_ties = FALSE) %>%
    470.   ungroup()
    471. res <- as.data.frame(res)
    472. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    473. res <- res[order(res$padj, -res$log2FoldChange), ]
    475. # Assuming res is your dataframe and already processed
    476. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    477. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    478. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    479. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    480. # Create a new workbook
    481. wb <- createWorkbook()
    482. # Add the complete dataset as the first sheet
    483. addWorksheet(wb, "Complete_Data")
    484. writeData(wb, "Complete_Data", res)
    485. # Add the up-regulated genes as the second sheet
    486. addWorksheet(wb, "Up_Regulated")
    487. writeData(wb, "Up_Regulated", up_regulated)
    488. # Add the down-regulated genes as the third sheet
    489. addWorksheet(wb, "Down_Regulated")
    490. writeData(wb, "Down_Regulated", down_regulated)
    491. # Save the workbook to a file
    492. saveWorkbook(wb, "Gene_Expression_Mac.AB_vs_Mac.WT19606.xlsx", overwrite = TRUE)
    494. # Set the 'GeneName' column as row.names
    495. rownames(res) <- res$GeneName
    496. # Drop the 'GeneName' column since it's now the row names
    497. res$GeneName <- NULL
    498. head(res)
    500. ## Ensure the data frame matches the expected format
    501. ## For example, it should have columns: log2FoldChange, padj, etc.
    502. #res <- as.data.frame(res)
    503. ## Remove rows with NA in log2FoldChange (if needed)
    504. #res <- res[!is.na(res$log2FoldChange),]
    506. # Replace padj = 0 with a small value
    507. res$padj[res$padj == 0] <- 1e-12
    509. #library(EnhancedVolcano)
    510. # Assuming res is already sorted and processed
    511. png("Mac.AB_vs_Mac.WT19606.png", width=1200, height=1200)
    512. #max.overlaps = 10
    513. EnhancedVolcano(res,
    514.                 lab = rownames(res),
    515.                 x = 'log2FoldChange',
    516.                 y = 'padj',
    517.                 pCutoff = 1e-2,
    518.                 FCcutoff = 2,
    519.                 title = '',
    520.                 subtitleLabSize = 18,
    521.                 pointSize = 3.0,
    522.                 labSize = 5.0,
    523.                 colAlpha = 1,
    524.                 legendIconSize = 4.0,
    525.                 drawConnectors = TRUE,
    526.                 widthConnectors = 0.5,
    527.                 colConnectors = 'black',
    528.                 subtitle = expression("Mac.AB versus Mac.WT19606"))
    529. dev.off()
    531. # ---- Mac.IJ_vs_Mac.WT19606 ----
    532. res <- read.csv("Mac.IJ_vs_Mac.WT19606-all.csv")
    533. # Replace empty GeneName with modified GeneID
    534. res$GeneName <- ifelse(
    535.   res$GeneName == "" | is.na(res$GeneName),
    536.   gsub("gene-", "", res$GeneID),
    537.   res$GeneName
    538. )
    539. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    541. res <- res %>%
    542.   group_by(GeneName) %>%
    543.   slice_min(padj, with_ties = FALSE) %>%
    544.   ungroup()
    545. res <- as.data.frame(res)
    546. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    547. res <- res[order(res$padj, -res$log2FoldChange), ]
    549. # Assuming res is your dataframe and already processed
    550. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    551. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    552. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    553. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    554. # Create a new workbook
    555. wb <- createWorkbook()
    556. # Add the complete dataset as the first sheet
    557. addWorksheet(wb, "Complete_Data")
    558. writeData(wb, "Complete_Data", res)
    559. # Add the up-regulated genes as the second sheet
    560. addWorksheet(wb, "Up_Regulated")
    561. writeData(wb, "Up_Regulated", up_regulated)
    562. # Add the down-regulated genes as the third sheet
    563. addWorksheet(wb, "Down_Regulated")
    564. writeData(wb, "Down_Regulated", down_regulated)
    565. # Save the workbook to a file
    566. saveWorkbook(wb, "Gene_Expression_Mac.IJ_vs_Mac.WT19606.xlsx", overwrite = TRUE)
    568. # Set the 'GeneName' column as row.names
    569. rownames(res) <- res$GeneName
    570. # Drop the 'GeneName' column since it's now the row names
    571. res$GeneName <- NULL
    572. head(res)
    574. ## Ensure the data frame matches the expected format
    575. ## For example, it should have columns: log2FoldChange, padj, etc.
    576. #res <- as.data.frame(res)
    577. ## Remove rows with NA in log2FoldChange (if needed)
    578. #res <- res[!is.na(res$log2FoldChange),]
    580. # Replace padj = 0 with a small value
    581. res$padj[res$padj == 0] <- 1e-12
    583. #library(EnhancedVolcano)
    584. # Assuming res is already sorted and processed
    585. png("Mac.IJ_vs_Mac.WT19606.png", width=1200, height=1200)
    586. #max.overlaps = 10
    587. EnhancedVolcano(res,
    588.                 lab = rownames(res),
    589.                 x = 'log2FoldChange',
    590.                 y = 'padj',
    591.                 pCutoff = 1e-2,
    592.                 FCcutoff = 2,
    593.                 title = '',
    594.                 subtitleLabSize = 18,
    595.                 pointSize = 3.0,
    596.                 labSize = 5.0,
    597.                 colAlpha = 1,
    598.                 legendIconSize = 4.0,
    599.                 drawConnectors = TRUE,
    600.                 widthConnectors = 0.5,
    601.                 colConnectors = 'black',
    602.                 subtitle = expression("Mac.IJ versus Mac.WT19606"))
    603. dev.off()
    605. # ---- Mac.W1_vs_Mac.WT19606 ----
    606. res <- read.csv("Mac.W1_vs_Mac.WT19606-all.csv")
    607. # Replace empty GeneName with modified GeneID
    608. res$GeneName <- ifelse(
    609.   res$GeneName == "" | is.na(res$GeneName),
    610.   gsub("gene-", "", res$GeneID),
    611.   res$GeneName
    612. )
    613. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    615. res <- res %>%
    616.   group_by(GeneName) %>%
    617.   slice_min(padj, with_ties = FALSE) %>%
    618.   ungroup()
    619. res <- as.data.frame(res)
    620. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    621. res <- res[order(res$padj, -res$log2FoldChange), ]
    623. # Assuming res is your dataframe and already processed
    624. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    625. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    626. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    627. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    628. # Create a new workbook
    629. wb <- createWorkbook()
    630. # Add the complete dataset as the first sheet
    631. addWorksheet(wb, "Complete_Data")
    632. writeData(wb, "Complete_Data", res)
    633. # Add the up-regulated genes as the second sheet
    634. addWorksheet(wb, "Up_Regulated")
    635. writeData(wb, "Up_Regulated", up_regulated)
    636. # Add the down-regulated genes as the third sheet
    637. addWorksheet(wb, "Down_Regulated")
    638. writeData(wb, "Down_Regulated", down_regulated)
    639. # Save the workbook to a file
    640. saveWorkbook(wb, "Gene_Expression_Mac.W1_vs_Mac.WT19606.xlsx", overwrite = TRUE)
    642. # Set the 'GeneName' column as row.names
    643. rownames(res) <- res$GeneName
    644. # Drop the 'GeneName' column since it's now the row names
    645. res$GeneName <- NULL
    646. head(res)
    648. ## Ensure the data frame matches the expected format
    649. ## For example, it should have columns: log2FoldChange, padj, etc.
    650. #res <- as.data.frame(res)
    651. ## Remove rows with NA in log2FoldChange (if needed)
    652. #res <- res[!is.na(res$log2FoldChange),]
    654. # Replace padj = 0 with a small value
    655. res$padj[res$padj == 0] <- 1e-12
    657. #library(EnhancedVolcano)
    658. # Assuming res is already sorted and processed
    659. png("Mac.W1_vs_Mac.WT19606.png", width=1200, height=1200)
    660. #max.overlaps = 10
    661. EnhancedVolcano(res,
    662.                 lab = rownames(res),
    663.                 x = 'log2FoldChange',
    664.                 y = 'padj',
    665.                 pCutoff = 1e-2,
    666.                 FCcutoff = 2,
    667.                 title = '',
    668.                 subtitleLabSize = 18,
    669.                 pointSize = 3.0,
    670.                 labSize = 5.0,
    671.                 colAlpha = 1,
    672.                 legendIconSize = 4.0,
    673.                 drawConnectors = TRUE,
    674.                 widthConnectors = 0.5,
    675.                 colConnectors = 'black',
    676.                 subtitle = expression("Mac.W1 versus Mac.WT19606"))
    677. dev.off()
    679. # ---- Mac.Y1_vs_Mac.WT19606 ----
    680. res <- read.csv("Mac.Y1_vs_Mac.WT19606-all.csv")
    681. # Replace empty GeneName with modified GeneID
    682. res$GeneName <- ifelse(
    683.   res$GeneName == "" | is.na(res$GeneName),
    684.   gsub("gene-", "", res$GeneID),
    685.   res$GeneName
    686. )
    687. duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    689. res <- res %>%
    690.   group_by(GeneName) %>%
    691.   slice_min(padj, with_ties = FALSE) %>%
    692.   ungroup()
    693. res <- as.data.frame(res)
    694. # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    695. res <- res[order(res$padj, -res$log2FoldChange), ]
    697. # Assuming res is your dataframe and already processed
    698. # Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
    699. up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
    700. # Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
    701. down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
    702. # Create a new workbook
    703. wb <- createWorkbook()
    704. # Add the complete dataset as the first sheet
    705. addWorksheet(wb, "Complete_Data")
    706. writeData(wb, "Complete_Data", res)
    707. # Add the up-regulated genes as the second sheet
    708. addWorksheet(wb, "Up_Regulated")
    709. writeData(wb, "Up_Regulated", up_regulated)
    710. # Add the down-regulated genes as the third sheet
    711. addWorksheet(wb, "Down_Regulated")
    712. writeData(wb, "Down_Regulated", down_regulated)
    713. # Save the workbook to a file
    714. saveWorkbook(wb, "Gene_Expression_Mac.Y1_vs_Mac.WT19606.xlsx", overwrite = TRUE)
    716. # Set the 'GeneName' column as row.names
    717. rownames(res) <- res$GeneName
    718. # Drop the 'GeneName' column since it's now the row names
    719. res$GeneName <- NULL
    720. head(res)
    722. ## Ensure the data frame matches the expected format
    723. ## For example, it should have columns: log2FoldChange, padj, etc.
    724. #res <- as.data.frame(res)
    725. ## Remove rows with NA in log2FoldChange (if needed)
    726. #res <- res[!is.na(res$log2FoldChange),]
    728. # Replace padj = 0 with a small value
    729. res$padj[res$padj == 0] <- 1e-12
    731. #library(EnhancedVolcano)
    732. # Assuming res is already sorted and processed
    733. png("Mac.Y1_vs_Mac.WT19606.png", width=1200, height=1200)
    734. #max.overlaps = 10
    735. EnhancedVolcano(res,
    736.                 lab = rownames(res),
    737.                 x = 'log2FoldChange',
    738.                 y = 'padj',
    739.                 pCutoff = 1e-2,
    740.                 FCcutoff = 2,
    741.                 title = '',
    742.                 subtitleLabSize = 18,
    743.                 pointSize = 3.0,
    744.                 labSize = 5.0,
    745.                 colAlpha = 1,
    746.                 legendIconSize = 4.0,
    747.                 drawConnectors = TRUE,
    748.                 widthConnectors = 0.5,
    749.                 colConnectors = 'black',
    750.                 subtitle = expression("Mac.Y1 versus Mac.WT19606"))
    751. dev.off()
    753. #TODO: annotate the Gene_Expression_xxx_vs_yyy.xlsx

  10. Clustering the genes and draw heatmap

    1. #http://xgenes.com/article/article-content/150/draw-venn-diagrams-using-matplotlib/
    2. #http://xgenes.com/article/article-content/276/go-terms-for-s-epidermidis/
    4. # save the Up-regulated and Down-regulated genes into -up.id and -down.id
    5. for i in Mac_vs_LB LB.AB_vs_LB.WT19606 LB.IJ_vs_LB.WT19606 LB.W1_vs_LB.WT19606 LB.Y1_vs_LB.WT19606 Mac.AB_vs_Mac.WT19606 Mac.IJ_vs_Mac.WT19606 Mac.W1_vs_Mac.WT19606 Mac.Y1_vs_Mac.WT19606; do
    6.   echo "cut -d',' -f1-1 ${i}-up.txt > ${i}-up.id";
    7.   echo "cut -d',' -f1-1 ${i}-down.txt > ${i}-down.id";
    8. done
    9. #5 LB.AB_vs_LB.WT19606-down.id
    10. #20 LB.AB_vs_LB.WT19606-up.id
    11. #64 LB.IJ_vs_LB.WT19606-down.id
    12. #69 LB.IJ_vs_LB.WT19606-up.id
    13. #23 LB.W1_vs_LB.WT19606-down.id
    14. #97 LB.W1_vs_LB.WT19606-up.id
    15. #9 LB.Y1_vs_LB.WT19606-down.id
    16. #20 LB.Y1_vs_LB.WT19606-up.id
    17. #20 Mac.AB_vs_Mac.WT19606-down.id
    18. #29 Mac.AB_vs_Mac.WT19606-up.id
    19. #65 Mac.IJ_vs_Mac.WT19606-down.id
    20. #197 Mac.IJ_vs_Mac.WT19606-up.id
    21. #359 Mac_vs_LB-down.id
    22. #308 Mac_vs_LB-up.id
    23. #290 Mac.W1_vs_Mac.WT19606-down.id
    24. #343 Mac.W1_vs_Mac.WT19606-up.id
    25. #75 Mac.Y1_vs_Mac.WT19606-down.id
    26. #0 Mac.Y1_vs_Mac.WT19606.png-down.id
    27. #0 Mac.Y1_vs_Mac.WT19606.png-up.id
    28. #68 Mac.Y1_vs_Mac.WT19606-up.id
    29. #2061 total
    31. cat *.id | sort -u > ids
    32. #Delete "GeneName"
    33. #add Gene_Id in the first line, delete the ""  #Note that using GeneID as index, rather than GeneName, since .txt contains only GeneID.
    34. GOI <- read.csv("ids")$Gene_Id    #1329
    35. RNASeq.NoCellLine <- assay(rld)
    36. #install.packages("gplots")
    37. library("gplots")
    38. #clustering methods: "ward.D", "ward.D2", "single", "complete", "average" (= UPGMA), "mcquitty" (= WPGMA), "median" (= WPGMC) or "centroid" (= UPGMC).  pearson or spearman
    39. datamat = RNASeq.NoCellLine[GOI, ]
    40. #datamat = RNASeq.NoCellLine
    41. write.csv(as.data.frame(datamat), file ="DEGs_heatmap_expression_data.txt")
    43. constant_rows <- apply(datamat, 1, function(row) var(row) == 0)
    44. if(any(constant_rows)) {
    45.   cat("Removing", sum(constant_rows), "constant rows.\n")
    46.   datamat <- datamat[!constant_rows, ]
    47. }
    48. hr <- hclust(as.dist(1-cor(t(datamat), method="pearson")), method="complete")
    49. hc <- hclust(as.dist(1-cor(datamat, method="spearman")), method="complete")
    50. mycl = cutree(hr, h=max(hr$height)/1.15)
    51. mycol = c("YELLOW", "BLUE", "ORANGE", "MAGENTA", "CYAN", "RED", "GREEN", "MAROON", "LIGHTBLUE", "PINK", "MAGENTA", "LIGHTCYAN", "LIGHTRED", "LIGHTGREEN");
    52. mycol = mycol[as.vector(mycl)]
    53. #png("DEGs_heatmap.png", width=900, height=800)
    54. #cex.lab=10, labRow="",
    55. png("DEGs_heatmap.png", width=1200, height=1000)
    56. heatmap.2(as.matrix(datamat),Rowv=as.dendrogram(hr),Colv = NA, dendrogram = 'row',labRow="",
    57.             scale='row',trace='none',col=bluered(75), cexCol=1.8,
    58.             RowSideColors = mycol, margins=c(10,2), cexRow=1.5, srtCol=30, lhei = c(1, 8), lwid=c(2, 8))  #rownames(datamat)
    59. #heatmap.2(datamat, Rowv=as.dendrogram(hr), col=bluered(75), scale="row", RowSideColors=mycol, trace="none", margin=c(5,5), sepwidth=c(0,0), dendrogram = 'row', Colv = 'false', density.info='none', labRow="", srtCol=30, lhei=c(0.1,2))
    60. dev.off()
    61. #### cluster members #####
    62. write.csv(names(subset(mycl, mycl == '1')),file='cluster1_YELLOW.txt')
    63. write.csv(names(subset(mycl, mycl == '2')),file='cluster2_DARKBLUE.txt')
    64. write.csv(names(subset(mycl, mycl == '3')),file='cluster3_DARKORANGE.txt')
    65. write.csv(names(subset(mycl, mycl == '4')),file='cluster4_DARKMAGENTA.txt')
    66. write.csv(names(subset(mycl, mycl == '5')),file='cluster5_DARKCYAN.txt')
    67. #~/Tools/csv2xls-0.4/csv_to_xls.py cluster*.txt -d',' -o DEGs_heatmap_cluster_members.xls
    68. #~/Tools/csv2xls-0.4/csv_to_xls.py DEGs_heatmap_expression_data.txt -d',' -o DEGs_heatmap_expression_data.xls;
    70. #### (NOT_WORKING) cluster members (adding annotations, note that it does not work for the bacteria, since it is not model-speices and we cannot use mart=ensembl) #####
    71. subset_1<-names(subset(mycl, mycl == '1'))
    72. data <- as.data.frame(datamat[rownames(datamat) %in% subset_1, ])  #2575
    73. subset_2<-names(subset(mycl, mycl == '2'))
    74. data <- as.data.frame(datamat[rownames(datamat) %in% subset_2, ])  #1855
    75. subset_3<-names(subset(mycl, mycl == '3'))
    76. data <- as.data.frame(datamat[rownames(datamat) %in% subset_3, ])  #217
    77. subset_4<-names(subset(mycl, mycl == '4'))
    78. data <- as.data.frame(datamat[rownames(datamat) %in% subset_4, ])  #
    79. subset_5<-names(subset(mycl, mycl == '5'))
    80. data <- as.data.frame(datamat[rownames(datamat) %in% subset_5, ])  #
    81. # Initialize an empty data frame for the annotated data
    82. annotated_data <- data.frame()
    83. # Determine total number of genes
    84. total_genes <- length(rownames(data))
    85. # Loop through each gene to annotate
    86. for (i in 1:total_genes) {
    87.     gene <- rownames(data)[i]
    88.     result <- getBM(attributes = c('ensembl_gene_id', 'external_gene_name', 'gene_biotype', 'entrezgene_id', 'chromosome_name', 'start_position', 'end_position', 'strand', 'description'),
    89.                     filters = 'ensembl_gene_id',
    90.                     values = gene,
    91.                     mart = ensembl)
    92.     # If multiple rows are returned, take the first one
    93.     if (nrow(result) > 1) {
    94.         result <- result[1, ]
    95.     }
    96.     # Check if the result is empty
    97.     if (nrow(result) == 0) {
    98.         result <- data.frame(ensembl_gene_id = gene,
    99.                             external_gene_name = NA,
    100.                             gene_biotype = NA,
    101.                             entrezgene_id = NA,
    102.                             chromosome_name = NA,
    103.                             start_position = NA,
    104.                             end_position = NA,
    105.                             strand = NA,
    106.                             description = NA)
    107.     }
    108.     # Transpose expression values
    109.     expression_values <- t(data.frame(t(data[gene, ])))
    110.     colnames(expression_values) <- colnames(data)
    111.     # Combine gene information and expression data
    112.     combined_result <- cbind(result, expression_values)
    113.     # Append to the final dataframe
    114.     annotated_data <- rbind(annotated_data, combined_result)
    115.     # Print progress every 100 genes
    116.     if (i %% 100 == 0) {
    117.         cat(sprintf("Processed gene %d out of %d\n", i, total_genes))
    118.     }
    119. }
    120. # Save the annotated data to a new CSV file
    121. write.csv(annotated_data, "cluster1_YELLOW.csv", row.names=FALSE)
    122. write.csv(annotated_data, "cluster2_DARKBLUE.csv", row.names=FALSE)
    123. write.csv(annotated_data, "cluster3_DARKORANGE.csv", row.names=FALSE)
    124. write.csv(annotated_data, "cluster4_DARKMAGENTA.csv", row.names=FALSE)
    125. write.csv(annotated_data, "cluster5_DARKCYAN.csv", row.names=FALSE)
    126. #~/Tools/csv2xls-0.4/csv_to_xls.py cluster*.csv -d',' -o DEGs_heatmap_clusters.xls

KEGG and GO annotations in non-model organisms

https://www.biobam.com/functional-analysis/

Blast2GO_workflow

  1. Assign KEGG and GO Terms (see diagram above)

    Since your organism is non-model, standard R databases (org.Hs.eg.db, etc.) won’t work. You’ll need to manually retrieve KEGG and GO annotations.

    Option 1 (KEGG Terms): EggNog based on orthology and phylogenies

    1. EggNOG-mapper assigns both KEGG Orthology (KO) IDs and GO terms.
    3. Install EggNOG-mapper:
    5.     mamba create -n eggnog_env python=3.8 eggnog-mapper -c conda-forge -c bioconda  #eggnog-mapper_2.1.12
    6.     mamba activate eggnog_env
    8. Run annotation:
    10.     #diamond makedb --in eggnog6.prots.faa -d eggnog_proteins.dmnd
    11.     mkdir /home/jhuang/mambaforge/envs/eggnog_env/lib/python3.8/site-packages/data/
    12.     download_eggnog_data.py --dbname eggnog.db -y --data_dir /home/jhuang/mambaforge/envs/eggnog_env/lib/python3.8/site-packages/data/
    13.     #NOT_WORKING: emapper.py -i CP059040_gene.fasta -o eggnog_dmnd_out --cpu 60 -m diamond[hmmer,mmseqs] --dmnd_db /home/jhuang/REFs/eggnog_data/data/eggnog_proteins.dmnd
    14.     python ~/Scripts/update_fasta_header.py CP059040_protein_.fasta CP059040_protein.fasta
    15.     emapper.py -i CP059040_protein.fasta -o eggnog_out --cpu 60 --resume
    16.     #----> result annotations.tsv: Contains KEGG, GO, and other functional annotations.
    17.     #---->  470.IX87_14445:
    18.         * 470 likely refers to the organism or strain (e.g., Acinetobacter baumannii ATCC 19606 or another related strain).
    19.         * IX87_14445 would refer to a specific gene or protein within that genome.
    21. Extract KEGG KO IDs from annotations.emapper.annotations.

    Option 2 (GO Terms from 'Blast2GO 5 Basic', saved in blast2go_annot.annot): Using Blast/Diamond + Blast2GO_GUI based on sequence alignment + GO mapping

    • 'Load protein sequences' (Tags: NONE, generated columns: Nr, SeqName) -->
    • Buttons 'blast' (Tags: BLASTED, generated columns: Description, Length, #Hits, e-Value, sim mean),
    • Button 'mapping' (Tags: MAPPED, generated columns: #GO, GO IDs, GO Names), "Mapping finished - Please proceed now to annotation."
    • Button 'annot' (Tags: ANNOTATED, generated columns: Enzyme Codes, Enzyme Names), "Annotation finished." * Used parameter 'Annotation CutOff': The Blast2GO Annotation Rule seeks to find the most specific GO annotations with a certain level of reliability. An annotation score is calculated for each candidate GO which is composed by the sequence similarity of the Blast Hit, the evidence code of the source GO and the position of the particular GO in the Gene Ontology hierarchy. This annotation score cutoff select the most specific GO term for a given GO branch which lies above this value. * Used parameter 'GO Weight' is a value which is added to Annotation Score of a more general/abstract Gene Ontology term for each of its more specific, original source GO terms. In this case, more general GO terms which summarise many original source terms (those ones directly associated to the Blast Hits) will have a higher Annotation Score.

    or blast2go_cli_v1.5.1 (NOT_USED)

    1.     #https://help.biobam.com/space/BCD/2250407989/Installation
    2.     #see ~/Scripts/blast2go_pipeline.sh

    Option 3 (GO Terms from 'Blast2GO 5 Basic', saved in blast2go_annot.annot2): Interpro based protein families / domains --> Button interpro * Button 'interpro' (Tags: INTERPRO, generated columns: InterPro IDs, InterPro GO IDs, InterPro GO Names) --> "InterProScan Finished - You can now merge the obtained GO Annotations."

    MERGE the results of InterPro GO IDs (Option 3) to GO IDs (Option 2) and generate final GO IDs * Button 'interpro'/'Merge InterProScan GOs to Annotation' --> "Merge (add and validate) all GO terms retrieved via InterProScan to the already existing GO annotation." --> "Finished merging GO terms from InterPro with annotations. Maybe you want to run ANNEX (Annotation Augmentation)." #* Button 'annot'/'ANNEX' --> "ANNEX finished. Maybe you want to do the next step: Enzyme Code Mapping."

    1. #-- before merging (blast2go_annot.annot) --
    2. #H0N29_18790     GO:0004842      ankyrin repeat domain-containing protein
    3. #H0N29_18790     GO:0085020
    4. #-- after merging (blast2go_annot.annot2) -->
    5. #H0N29_18790     GO:0031436      ankyrin repeat domain-containing protein
    6. #H0N29_18790     GO:0070531
    7. #H0N29_18790     GO:0004842
    8. #H0N29_18790     GO:0005515
    9. #H0N29_18790     GO:0085020

    Option 4 (NOT_USED): RFAM for non-colding RNA

    Option 5 (NOT_USED): PSORTb for subcellular localizations

    Option 6 (NOT_USED): KAAS (KEGG Automatic Annotation Server)

    • Go to KAAS
    • Upload your FASTA file.
    • Select an appropriate gene set.
    • Download the KO assignments.

  2. Find the Closest KEGG Organism Code (NOT_USED)

    Since your species isn't directly in KEGG, use a closely related organism.

    • Check available KEGG organisms:
      1. library(clusterProfiler)
      2. library(KEGGREST)
      4. kegg_organisms <- keggList("organism")
      6. Pick the closest relative (e.g., zebrafish "dre" for fish, Arabidopsis "ath" for plants).
      8. # Search for Acinetobacter in the list
      9. grep("Acinetobacter", kegg_organisms, ignore.case = TRUE, value = TRUE)
      10. # Gammaproteobacteria
      11. #Extract KO IDs from the eggnog results for  "Acinetobacter baumannii strain ATCC 19606"

  3. Find the Closest KEGG Organism for a Non-Model Species

    If your organism is not in KEGG, search for the closest relative:

    1.     grep("fish", kegg_organisms, ignore.case = TRUE, value = TRUE)  # Example search

    For KEGG pathway enrichment in non-model species, use "ko" instead of a species code (the code has been intergrated in the point 4):

    1.     kegg_enrich <- enrichKEGG(gene = gene_list, organism = "ko")  # "ko" = KEGG Orthology

  4. Perform KEGG and GO Enrichment in R (under dir ~/DATA/ata_Tam_RNAseq_2025_LB_vs_Mac_ATCC19606/results/star_salmon/degenes)

    1.     #BiocManager::install("GO.db")
    2.     #BiocManager::install("AnnotationDbi")
    4.     # Load required libraries
    5.     library(openxlsx)  # For Excel file handling
    6.     library(dplyr)     # For data manipulation
    7.     library(tidyr)
    8.     library(stringr)
    9.     library(clusterProfiler)  # For KEGG and GO enrichment analysis
    10.     #library(org.Hs.eg.db)  # Replace with appropriate organism database
    11.     library(GO.db)
    12.     library(AnnotationDbi)
    14.     setwd("~/DATA/Data_Tam_RNAseq_2025_LB_vs_Mac_ATCC19606/results/star_salmon/degenes")
    15.     # PREPARING go_terms and ec_terms: annot_* file: cut -f1-2 -d$'\t' blast2go_annot.annot2 > blast2go_annot.annot2_
    16.     # Step 1: Load the blast2go annotation file with a check for missing columns
    17.     annot_df <- read.table("/home/jhuang/b2gWorkspace_Tam_RNAseq_2024/blast2go_annot.annot2_",
    18.                         header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE)
    20.     # If the structure is inconsistent, we can make sure there are exactly 3 columns:
    21.     colnames(annot_df) <- c("GeneID", "Term")
    22.     # Step 2: Filter and aggregate GO and EC terms as before
    23.     go_terms <- annot_df %>%
    24.     filter(grepl("^GO:", Term)) %>%
    25.     group_by(GeneID) %>%
    26.     summarize(GOs = paste(Term, collapse = ","), .groups = "drop")
    27.     ec_terms <- annot_df %>%
    28.     filter(grepl("^EC:", Term)) %>%
    29.     group_by(GeneID) %>%
    30.     summarize(EC = paste(Term, collapse = ","), .groups = "drop")
    32.     # Load the results
    33.     #res <- read.csv("Mac_vs_LB-all.csv")     #up307, down358
    34.     #res <- read.csv("LB.AB_vs_LB.WT19606-all.csv")     #up307, down358
    35.     #res <- read.csv("LB.IJ_vs_LB.WT19606-all.csv")     #up307, down358
    36.     #res <- read.csv("LB.W1_vs_LB.WT19606-all.csv")     #up307, down358
    37.     #res <- read.csv("LB.Y1_vs_LB.WT19606-all.csv")     #up307, down358
    38.     #res <- read.csv("Mac.AB_vs_Mac.WT19606-all.csv")     #up307, down358
    39.     #res <- read.csv("Mac.IJ_vs_Mac.WT19606-all.csv")     #up307, down358
    40.     #res <- read.csv("Mac.W1_vs_Mac.WT19606-all.csv")     #up307, down358
    41.     res <- read.csv("Mac.Y1_vs_Mac.WT19606-all.csv")     #up307, down358
    43.     # Replace empty GeneName with modified GeneID
    44.     res$GeneName <- ifelse(
    45.         res$GeneName == "" | is.na(res$GeneName),
    46.         gsub("gene-", "", res$GeneID),
    47.         res$GeneName
    48.     )
    50.     # Remove duplicated genes by selecting the gene with the smallest padj
    51.     duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
    53.     res <- res %>%
    54.     group_by(GeneName) %>%
    55.     slice_min(padj, with_ties = FALSE) %>%
    56.     ungroup()
    58.     res <- as.data.frame(res)
    59.     # Sort res first by padj (ascending) and then by log2FoldChange (descending)
    60.     res <- res[order(res$padj, -res$log2FoldChange), ]
    61.     # Read eggnog annotations
    62.     eggnog_data <- read.delim("~/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/eggnog_out.emapper.annotations.txt", header = TRUE, sep = "\t")
    63.     # Remove the "gene-" prefix from GeneID in res to match eggnog 'query' format
    64.     res$GeneID <- gsub("gene-", "", res$GeneID)
    65.     # Merge eggnog data with res based on GeneID
    66.     res <- res %>% left_join(eggnog_data, by = c("GeneID" = "query"))
    68.     # Merge with the res dataframe
    69.     # Perform the left joins and rename columns
    70.     res_updated <- res %>%
    71.     left_join(go_terms, by = "GeneID") %>%
    72.     left_join(ec_terms, by = "GeneID") %>% dplyr::select(-EC.x, -GOs.x) %>% dplyr::rename(EC = EC.y, GOs = GOs.y)
    74.     # Filter up-regulated genes
    75.     up_regulated <- res_updated[res_updated$log2FoldChange > 2 & res_updated$padj < 0.01, ]
    76.     # Filter down-regulated genes
    77.     down_regulated <- res_updated[res_updated$log2FoldChange < -2 & res_updated$padj < 0.01, ]
    79.     # Create a new workbook
    80.     wb <- createWorkbook()
    81.     # Add the complete dataset as the first sheet (with annotations)
    82.     addWorksheet(wb, "Complete_Data")
    83.     writeData(wb, "Complete_Data", res_updated)
    84.     # Add the up-regulated genes as the second sheet (with annotations)
    85.     addWorksheet(wb, "Up_Regulated")
    86.     writeData(wb, "Up_Regulated", up_regulated)
    87.     # Add the down-regulated genes as the third sheet (with annotations)
    88.     addWorksheet(wb, "Down_Regulated")
    89.     writeData(wb, "Down_Regulated", down_regulated)
    90.     # Save the workbook to a file
    91.     saveWorkbook(wb, "Gene_Expression_with_Annotations_Urine_vs_MHB.xlsx", overwrite = TRUE)
    93.     # Set GeneName as row names after the join
    94.     rownames(res_updated) <- res_updated$GeneName
    95.     res_updated <- res_updated %>% dplyr::select(-GeneName)
    96.     ## Set the 'GeneName' column as row.names
    97.     #rownames(res_updated) <- res_updated$GeneName
    98.     ## Drop the 'GeneName' column since it's now the row names
    99.     #res_updated$GeneName <- NULL
    100.     # -- BREAK_1 --
    102.     # ---- Perform KEGG enrichment analysis (up_regulated) ----
    103.     gene_list_kegg_up <- up_regulated$KEGG_ko
    104.     gene_list_kegg_up <- gsub("ko:", "", gene_list_kegg_up)
    105.     kegg_enrichment_up <- enrichKEGG(gene = gene_list_kegg_up, organism = 'ko')
    106.     # -- convert the GeneID (Kxxxxxx) to the true GeneID --
    107.     # Step 0: Create KEGG to GeneID mapping
    108.     kegg_to_geneid_up <- up_regulated %>%
    109.     dplyr::select(KEGG_ko, GeneID) %>%
    110.     filter(!is.na(KEGG_ko)) %>%  # Remove missing KEGG KO entries
    111.     mutate(KEGG_ko = str_remove(KEGG_ko, "ko:"))  # Remove 'ko:' prefix if present
    112.     # Step 1: Clean KEGG_ko values (separate multiple KEGG IDs)
    113.     kegg_to_geneid_clean <- kegg_to_geneid_up %>%
    114.     mutate(KEGG_ko = str_remove_all(KEGG_ko, "ko:")) %>%  # Remove 'ko:' prefixes
    115.     separate_rows(KEGG_ko, sep = ",") %>%  # Ensure each KEGG ID is on its own row
    116.     filter(KEGG_ko != "-") %>%  # Remove invalid KEGG IDs ("-")
    117.     distinct()  # Remove any duplicate mappings
    118.     # Step 2.1: Expand geneID column in kegg_enrichment_up
    119.     expanded_kegg <- kegg_enrichment_up %>%
    120.     as.data.frame() %>%
    121.     separate_rows(geneID, sep = "/") %>%  # Split multiple KEGG IDs (Kxxxxx)
    122.     left_join(kegg_to_geneid_clean, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>%  # Explicitly handle many-to-many
    123.     distinct() %>%  # Remove duplicate matches
    124.     group_by(ID) %>%
    125.     summarise(across(everything(), ~ paste(unique(na.omit(.)), collapse = "/")), .groups = "drop")  # Re-collapse results
    126.     #dplyr::glimpse(expanded_kegg)
    127.     # Step 3.1: Replace geneID column in the original dataframe
    128.     kegg_enrichment_up_df <- as.data.frame(kegg_enrichment_up)
    129.     # Remove old geneID column and merge new one
    130.     kegg_enrichment_up_df <- kegg_enrichment_up_df %>%
    131.     dplyr::select(-geneID) %>%  # Remove old geneID column
    132.     left_join(expanded_kegg %>% dplyr::select(ID, GeneID), by = "ID") %>%  # Merge new GeneID column
    133.     dplyr::rename(geneID = GeneID)  # Rename column back to geneID
    135.     # ---- Perform KEGG enrichment analysis (down_regulated) ----
    136.     # Step 1: Extract KEGG KO terms from down-regulated genes
    137.     gene_list_kegg_down <- down_regulated$KEGG_ko
    138.     gene_list_kegg_down <- gsub("ko:", "", gene_list_kegg_down)
    139.     # Step 2: Perform KEGG enrichment analysis
    140.     kegg_enrichment_down <- enrichKEGG(gene = gene_list_kegg_down, organism = 'ko')
    141.     # --- Convert KEGG gene IDs (Kxxxxxx) to actual GeneIDs ---
    142.     # Step 3: Create KEGG to GeneID mapping from down_regulated dataset
    143.     kegg_to_geneid_down <- down_regulated %>%
    144.     dplyr::select(KEGG_ko, GeneID) %>%
    145.     filter(!is.na(KEGG_ko)) %>%  # Remove missing KEGG KO entries
    146.     mutate(KEGG_ko = str_remove(KEGG_ko, "ko:"))  # Remove 'ko:' prefix if present
    147.     # -- BREAK_2 --
    149.     # Step 4: Clean KEGG_ko values (handle multiple KEGG IDs)
    150.     kegg_to_geneid_down_clean <- kegg_to_geneid_down %>%
    151.     mutate(KEGG_ko = str_remove_all(KEGG_ko, "ko:")) %>%  # Remove 'ko:' prefixes
    152.     separate_rows(KEGG_ko, sep = ",") %>%  # Ensure each KEGG ID is on its own row
    153.     filter(KEGG_ko != "-") %>%  # Remove invalid KEGG IDs ("-")
    154.     distinct()  # Remove duplicate mappings
    155.     # Step 5: Expand geneID column in kegg_enrichment_down
    156.     expanded_kegg_down <- kegg_enrichment_down %>%
    157.     as.data.frame() %>%
    158.     separate_rows(geneID, sep = "/") %>%  # Split multiple KEGG IDs (Kxxxxx)
    159.     left_join(kegg_to_geneid_down_clean, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>%  # Handle many-to-many mappings
    160.     distinct() %>%  # Remove duplicate matches
    161.     group_by(ID) %>%
    162.     summarise(across(everything(), ~ paste(unique(na.omit(.)), collapse = "/")), .groups = "drop")  # Re-collapse results
    163.     # Step 6: Replace geneID column in the original kegg_enrichment_down dataframe
    164.     kegg_enrichment_down_df <- as.data.frame(kegg_enrichment_down) %>%
    165.     dplyr::select(-geneID) %>%  # Remove old geneID column
    166.     left_join(expanded_kegg_down %>% dplyr::select(ID, GeneID), by = "ID") %>%  # Merge new GeneID column
    167.     dplyr::rename(geneID = GeneID)  # Rename column back to geneID
    168.     # View the updated dataframe
    169.     head(kegg_enrichment_down_df)
    171.     # Create a new workbook
    172.     wb <- createWorkbook()
    173.     # Save enrichment results to the workbook
    174.     addWorksheet(wb, "KEGG_Enrichment_Up")
    175.     writeData(wb, "KEGG_Enrichment_Up", as.data.frame(kegg_enrichment_up_df))
    176.     # Save enrichment results to the workbook
    177.     addWorksheet(wb, "KEGG_Enrichment_Down")
    178.     writeData(wb, "KEGG_Enrichment_Down", as.data.frame(kegg_enrichment_down_df))
    180.     # Define gene list (up-regulated genes)
    181.     gene_list_go_up <- up_regulated$GeneID  # Extract the 149 up-regulated genes
    182.     gene_list_go_down <- down_regulated$GeneID  # Extract the 65 down-regulated genes
    184.     # Define background gene set (all genes in res)
    185.     background_genes <- res_updated$GeneID  # Extract the 3646 background genes
    187.     # Prepare GO annotation data from res
    188.     go_annotation <- res_updated[, c("GOs","GeneID")]  # Extract relevant columns
    189.     go_annotation <- go_annotation %>%
    190.     tidyr::separate_rows(GOs, sep = ",")  # Split multiple GO terms into separate rows
    191.     # -- BREAK_3 --
    193.     go_enrichment_up <- enricher(
    194.         gene = gene_list_go_up,                # Up-regulated genes
    195.         TERM2GENE = go_annotation,       # Custom GO annotation
    196.         pvalueCutoff = 0.05,             # Significance threshold
    197.         pAdjustMethod = "BH",
    198.         universe = background_genes      # Define the background gene set
    199.     )
    200.     go_enrichment_up <- as.data.frame(go_enrichment_up)
    202.     go_enrichment_down <- enricher(
    203.         gene = gene_list_go_down,                # Up-regulated genes
    204.         TERM2GENE = go_annotation,       # Custom GO annotation
    205.         pvalueCutoff = 0.05,             # Significance threshold
    206.         pAdjustMethod = "BH",
    207.         universe = background_genes      # Define the background gene set
    208.     )
    209.     go_enrichment_down <- as.data.frame(go_enrichment_down)
    211.     ## Remove the 'p.adjust' column since no adjusted methods have been applied --> In this version we have used pvalue filtering (see above)!
    212.     #go_enrichment_up <- go_enrichment_up[, !names(go_enrichment_up) %in% "p.adjust"]
    213.     # Update the Description column with the term descriptions
    214.     go_enrichment_up$Description <- sapply(go_enrichment_up$ID, function(go_id) {
    215.     # Using select to get the term description
    216.     term <- tryCatch({
    217.         AnnotationDbi::select(GO.db, keys = go_id, columns = "TERM", keytype = "GOID")
    218.     }, error = function(e) {
    219.         message(paste("Error for GO term:", go_id))  # Print which GO ID caused the error
    220.         return(data.frame(TERM = NA))  # In case of error, return NA
    221.     })
    223.     if (nrow(term) > 0) {
    224.         return(term$TERM)
    225.     } else {
    226.         return(NA)  # If no description found, return NA
    227.     }
    228.     })
    229.     ## Print the updated data frame
    230.     #print(go_enrichment_up)
    232.     ## Remove the 'p.adjust' column since no adjusted methods have been applied --> In this version we have used pvalue filtering (see above)!
    233.     #go_enrichment_down <- go_enrichment_down[, !names(go_enrichment_down) %in% "p.adjust"]
    234.     # Update the Description column with the term descriptions
    235.     go_enrichment_down$Description <- sapply(go_enrichment_down$ID, function(go_id) {
    236.     # Using select to get the term description
    237.     term <- tryCatch({
    238.         AnnotationDbi::select(GO.db, keys = go_id, columns = "TERM", keytype = "GOID")
    239.     }, error = function(e) {
    240.         message(paste("Error for GO term:", go_id))  # Print which GO ID caused the error
    241.         return(data.frame(TERM = NA))  # In case of error, return NA
    242.     })
    244.     if (nrow(term) > 0) {
    245.         return(term$TERM)
    246.     } else {
    247.         return(NA)  # If no description found, return NA
    248.     }
    249.     })
    251.     addWorksheet(wb, "GO_Enrichment_Up")
    252.     writeData(wb, "GO_Enrichment_Up", as.data.frame(go_enrichment_up))
    254.     addWorksheet(wb, "GO_Enrichment_Down")
    255.     writeData(wb, "GO_Enrichment_Down", as.data.frame(go_enrichment_down))
    257.     # Save the workbook with enrichment results
    258.     saveWorkbook(wb, "KEGG_and_GO_Enrichments_Urine_vs_MHB.xlsx", overwrite = TRUE)
    260.     #Error for GO term: GO:0006807: replace "GO:0006807 obsolete nitrogen compound metabolic process"
    261.     #obsolete nitrogen compound metabolic process #https://www.ebi.ac.uk/QuickGO/term/GO:0006807
    262.     #TODO: marked the color as yellow if the p.adjusted <= 0.05 in GO_enrichment!
    264.     #mv KEGG_and_GO_Enrichments_Urine_vs_MHB.xlsx KEGG_and_GO_Enrichments_Mac_vs_LB.xlsx
    265.     #Mac_vs_LB
    266.     #LB.AB_vs_LB.WT19606
    267.     #LB.IJ_vs_LB.WT19606
    268.     #LB.W1_vs_LB.WT19606
    269.     #LB.Y1_vs_LB.WT19606
    270.     #Mac.AB_vs_Mac.WT19606
    271.     #Mac.IJ_vs_Mac.WT19606
    272.     #Mac.W1_vs_Mac.WT19606
    273.     #Mac.Y1_vs_Mac.WT19606

  5. (DEBUG) Draw the Venn diagram to compare the total DEGs across AUM, Urine, and MHB, irrespective of up- or down-regulation.

    1.         library(openxlsx)
    3.         # Function to read and clean gene ID files
    4.         read_gene_ids <- function(file_path) {
    5.         # Read the gene IDs from the file
    6.         gene_ids <- readLines(file_path)
    8.         # Remove any quotes and trim whitespaces
    9.         gene_ids <- gsub('"', '', gene_ids)  # Remove quotes
    10.         gene_ids <- trimws(gene_ids)  # Trim whitespaces
    12.         # Remove empty entries or NAs
    13.         gene_ids <- gene_ids[gene_ids != "" & !is.na(gene_ids)]
    15.         return(gene_ids)
    16.         }
    18.         # Example list of LB files with both -up.id and -down.id for each condition
    19.         lb_files_up <- c("LB.AB_vs_LB.WT19606-up.id", "LB.IJ_vs_LB.WT19606-up.id",
    20.                         "LB.W1_vs_LB.WT19606-up.id", "LB.Y1_vs_LB.WT19606-up.id")
    21.         lb_files_down <- c("LB.AB_vs_LB.WT19606-down.id", "LB.IJ_vs_LB.WT19606-down.id",
    22.                         "LB.W1_vs_LB.WT19606-down.id", "LB.Y1_vs_LB.WT19606-down.id")
    24.         # Combine both up and down files for each condition
    25.         lb_files <- c(lb_files_up, lb_files_down)
    27.         # Read gene IDs for each file in LB group
    28.         #lb_degs <- setNames(lapply(lb_files, read_gene_ids), gsub("-(up|down).id", "", lb_files))
    29.         lb_degs <- setNames(lapply(lb_files, read_gene_ids), make.unique(gsub("-(up|down).id", "", lb_files)))
    31.         lb_degs_ <- list()
    32.         combined_set <- c(lb_degs[["LB.AB_vs_LB.WT19606"]], lb_degs[["LB.AB_vs_LB.WT19606.1"]])
    33.         #unique_combined_set <- unique(combined_set)
    34.         lb_degs_$AB <- combined_set
    35.         combined_set <- c(lb_degs[["LB.IJ_vs_LB.WT19606"]], lb_degs[["LB.IJ_vs_LB.WT19606.1"]])
    36.         lb_degs_$IJ <- combined_set
    37.         combined_set <- c(lb_degs[["LB.W1_vs_LB.WT19606"]], lb_degs[["LB.W1_vs_LB.WT19606.1"]])
    38.         lb_degs_$W1 <- combined_set
    39.         combined_set <- c(lb_degs[["LB.Y1_vs_LB.WT19606"]], lb_degs[["LB.Y1_vs_LB.WT19606.1"]])
    40.         lb_degs_$Y1 <- combined_set
    42.         # Example list of Mac files with both -up.id and -down.id for each condition
    43.         mac_files_up <- c("Mac.AB_vs_Mac.WT19606-up.id", "Mac.IJ_vs_Mac.WT19606-up.id",
    44.                         "Mac.W1_vs_Mac.WT19606-up.id", "Mac.Y1_vs_Mac.WT19606-up.id")
    45.         mac_files_down <- c("Mac.AB_vs_Mac.WT19606-down.id", "Mac.IJ_vs_Mac.WT19606-down.id",
    46.                         "Mac.W1_vs_Mac.WT19606-down.id", "Mac.Y1_vs_Mac.WT19606-down.id")
    48.         # Combine both up and down files for each condition in Mac group
    49.         mac_files <- c(mac_files_up, mac_files_down)
    51.         # Read gene IDs for each file in Mac group
    52.         mac_degs <- setNames(lapply(mac_files, read_gene_ids), make.unique(gsub("-(up|down).id", "", mac_files)))
    54.         mac_degs_ <- list()
    55.         combined_set <- c(mac_degs[["Mac.AB_vs_Mac.WT19606"]], mac_degs[["Mac.AB_vs_Mac.WT19606.1"]])
    56.         mac_degs_$AB <- combined_set
    57.         combined_set <- c(mac_degs[["Mac.IJ_vs_Mac.WT19606"]], mac_degs[["Mac.IJ_vs_Mac.WT19606.1"]])
    58.         mac_degs_$IJ <- combined_set
    59.         combined_set <- c(mac_degs[["Mac.W1_vs_Mac.WT19606"]], mac_degs[["Mac.W1_vs_Mac.WT19606.1"]])
    60.         mac_degs_$W1 <- combined_set
    61.         combined_set <- c(mac_degs[["Mac.Y1_vs_Mac.WT19606"]], mac_degs[["Mac.Y1_vs_Mac.WT19606.1"]])
    62.         mac_degs_$Y1 <- combined_set
    64.         # Function to clean sheet names to ensure no sheet name exceeds 31 characters
    65.         truncate_sheet_name <- function(names_list) {
    66.         sapply(names_list, function(name) {
    67.         if (nchar(name) > 31) {
    68.         return(substr(name, 1, 31))  # Truncate sheet name to 31 characters
    69.         }
    70.         return(name)
    71.         })
    72.         }
    74.         # Assuming lb_degs_ is already a list of gene sets (LB.AB, LB.IJ, etc.)
    76.         # Define intersections between different conditions for LB
    77.         inter_lb_ab_ij <- intersect(lb_degs_$AB, lb_degs_$IJ)
    78.         inter_lb_ab_w1 <- intersect(lb_degs_$AB, lb_degs_$W1)
    79.         inter_lb_ab_y1 <- intersect(lb_degs_$AB, lb_degs_$Y1)
    80.         inter_lb_ij_w1 <- intersect(lb_degs_$IJ, lb_degs_$W1)
    81.         inter_lb_ij_y1 <- intersect(lb_degs_$IJ, lb_degs_$Y1)
    82.         inter_lb_w1_y1 <- intersect(lb_degs_$W1, lb_degs_$Y1)
    84.         # Define intersections between three conditions for LB
    85.         inter_lb_ab_ij_w1 <- Reduce(intersect, list(lb_degs_$AB, lb_degs_$IJ, lb_degs_$W1))
    86.         inter_lb_ab_ij_y1 <- Reduce(intersect, list(lb_degs_$AB, lb_degs_$IJ, lb_degs_$Y1))
    87.         inter_lb_ab_w1_y1 <- Reduce(intersect, list(lb_degs_$AB, lb_degs_$W1, lb_degs_$Y1))
    88.         inter_lb_ij_w1_y1 <- Reduce(intersect, list(lb_degs_$IJ, lb_degs_$W1, lb_degs_$Y1))
    90.         # Define intersection between all four conditions for LB
    91.         inter_lb_ab_ij_w1_y1 <- Reduce(intersect, list(lb_degs_$AB, lb_degs_$IJ, lb_degs_$W1, lb_degs_$Y1))
    93.         # Now remove the intersected genes from each original set for LB
    94.         venn_list_lb <- list()
    96.         # For LB.AB, remove genes that are also in other conditions
    97.         venn_list_lb[["LB.AB_only"]] <- setdiff(lb_degs_$AB, union(inter_lb_ab_ij, union(inter_lb_ab_w1, inter_lb_ab_y1)))
    99.         # For LB.IJ, remove genes that are also in other conditions
    100.         venn_list_lb[["LB.IJ_only"]] <- setdiff(lb_degs_$IJ, union(inter_lb_ab_ij, union(inter_lb_ij_w1, inter_lb_ij_y1)))
    102.         # For LB.W1, remove genes that are also in other conditions
    103.         venn_list_lb[["LB.W1_only"]] <- setdiff(lb_degs_$W1, union(inter_lb_ab_w1, union(inter_lb_ij_w1, inter_lb_ab_w1_y1)))
    105.         # For LB.Y1, remove genes that are also in other conditions
    106.         venn_list_lb[["LB.Y1_only"]] <- setdiff(lb_degs_$Y1, union(inter_lb_ab_y1, union(inter_lb_ij_y1, inter_lb_ab_w1_y1)))
    108.         # Add the intersections for LB (same as before)
    109.         venn_list_lb[["LB.AB_AND_LB.IJ"]] <- inter_lb_ab_ij
    110.         venn_list_lb[["LB.AB_AND_LB.W1"]] <- inter_lb_ab_w1
    111.         venn_list_lb[["LB.AB_AND_LB.Y1"]] <- inter_lb_ab_y1
    112.         venn_list_lb[["LB.IJ_AND_LB.W1"]] <- inter_lb_ij_w1
    113.         venn_list_lb[["LB.IJ_AND_LB.Y1"]] <- inter_lb_ij_y1
    114.         venn_list_lb[["LB.W1_AND_LB.Y1"]] <- inter_lb_w1_y1
    116.         # Define intersections between three conditions for LB
    117.         venn_list_lb[["LB.AB_AND_LB.IJ_AND_LB.W1"]] <- inter_lb_ab_ij_w1
    118.         venn_list_lb[["LB.AB_AND_LB.IJ_AND_LB.Y1"]] <- inter_lb_ab_ij_y1
    119.         venn_list_lb[["LB.AB_AND_LB.W1_AND_LB.Y1"]] <- inter_lb_ab_w1_y1
    120.         venn_list_lb[["LB.IJ_AND_LB.W1_AND_LB.Y1"]] <- inter_lb_ij_w1_y1
    122.         # Define intersection between all four conditions for LB
    123.         venn_list_lb[["LB.AB_AND_LB.IJ_AND_LB.W1_AND_LB.Y1"]] <- inter_lb_ab_ij_w1_y1
    125.         # Assuming mac_degs_ is already a list of gene sets (Mac.AB, Mac.IJ, etc.)
    127.         # Define intersections between different conditions
    128.         inter_mac_ab_ij <- intersect(mac_degs_$AB, mac_degs_$IJ)
    129.         inter_mac_ab_w1 <- intersect(mac_degs_$AB, mac_degs_$W1)
    130.         inter_mac_ab_y1 <- intersect(mac_degs_$AB, mac_degs_$Y1)
    131.         inter_mac_ij_w1 <- intersect(mac_degs_$IJ, mac_degs_$W1)
    132.         inter_mac_ij_y1 <- intersect(mac_degs_$IJ, mac_degs_$Y1)
    133.         inter_mac_w1_y1 <- intersect(mac_degs_$W1, mac_degs_$Y1)
    135.         # Define intersections between three conditions
    136.         inter_mac_ab_ij_w1 <- Reduce(intersect, list(mac_degs_$AB, mac_degs_$IJ, mac_degs_$W1))
    137.         inter_mac_ab_ij_y1 <- Reduce(intersect, list(mac_degs_$AB, mac_degs_$IJ, mac_degs_$Y1))
    138.         inter_mac_ab_w1_y1 <- Reduce(intersect, list(mac_degs_$AB, mac_degs_$W1, mac_degs_$Y1))
    139.         inter_mac_ij_w1_y1 <- Reduce(intersect, list(mac_degs_$IJ, mac_degs_$W1, mac_degs_$Y1))
    141.         # Define intersection between all four conditions
    142.         inter_mac_ab_ij_w1_y1 <- Reduce(intersect, list(mac_degs_$AB, mac_degs_$IJ, mac_degs_$W1, mac_degs_$Y1))
    144.         # Now remove the intersected genes from each original set
    145.         venn_list_mac <- list()
    147.         # For Mac.AB, remove genes that are also in other conditions
    148.         venn_list_mac[["Mac.AB_only"]] <- setdiff(mac_degs_$AB, union(inter_mac_ab_ij, union(inter_mac_ab_w1, inter_mac_ab_y1)))
    150.         # For Mac.IJ, remove genes that are also in other conditions
    151.         venn_list_mac[["Mac.IJ_only"]] <- setdiff(mac_degs_$IJ, union(inter_mac_ab_ij, union(inter_mac_ij_w1, inter_mac_ij_y1)))
    153.         # For Mac.W1, remove genes that are also in other conditions
    154.         venn_list_mac[["Mac.W1_only"]] <- setdiff(mac_degs_$W1, union(inter_mac_ab_w1, union(inter_mac_ij_w1, inter_mac_ab_w1_y1)))
    156.         # For Mac.Y1, remove genes that are also in other conditions
    157.         venn_list_mac[["Mac.Y1_only"]] <- setdiff(mac_degs_$Y1, union(inter_mac_ab_y1, union(inter_mac_ij_y1, inter_mac_ab_w1_y1)))
    159.         # Add the intersections (same as before)
    160.         venn_list_mac[["Mac.AB_AND_Mac.IJ"]] <- inter_mac_ab_ij
    161.         venn_list_mac[["Mac.AB_AND_Mac.W1"]] <- inter_mac_ab_w1
    162.         venn_list_mac[["Mac.AB_AND_Mac.Y1"]] <- inter_mac_ab_y1
    163.         venn_list_mac[["Mac.IJ_AND_Mac.W1"]] <- inter_mac_ij_w1
    164.         venn_list_mac[["Mac.IJ_AND_Mac.Y1"]] <- inter_mac_ij_y1
    165.         venn_list_mac[["Mac.W1_AND_Mac.Y1"]] <- inter_mac_w1_y1
    167.         # Define intersections between three conditions
    168.         venn_list_mac[["Mac.AB_AND_Mac.IJ_AND_Mac.W1"]] <- inter_mac_ab_ij_w1
    169.         venn_list_mac[["Mac.AB_AND_Mac.IJ_AND_Mac.Y1"]] <- inter_mac_ab_ij_y1
    170.         venn_list_mac[["Mac.AB_AND_Mac.W1_AND_Mac.Y1"]] <- inter_mac_ab_w1_y1
    171.         venn_list_mac[["Mac.IJ_AND_Mac.W1_AND_Mac.Y1"]] <- inter_mac_ij_w1_y1
    173.         # Define intersection between all four conditions
    174.         venn_list_mac[["Mac.AB_AND_Mac.IJ_AND_Mac.W1_AND_Mac.Y1"]] <- inter_mac_ab_ij_w1_y1
    176.         # Save the gene IDs to Excel for further inspection (optional)
    177.         write.xlsx(lb_degs, file = "LB_DEGs.xlsx")
    178.         write.xlsx(mac_degs, file = "Mac_DEGs.xlsx")
    180.         # Clean sheet names and write the Venn intersection sets for LB and Mac groups into Excel files
    181.         write.xlsx(venn_list_lb, file = "Venn_LB_Genes_Intersect.xlsx", sheetName = truncate_sheet_name(names(venn_list_lb)), rowNames = FALSE)
    182.         write.xlsx(venn_list_mac, file = "Venn_Mac_Genes_Intersect.xlsx", sheetName = truncate_sheet_name(names(venn_list_mac)), rowNames = FALSE)
    184.         # Venn Diagram for LB group
    185.         venn1 <- ggvenn(lb_degs_,
    186.                         fill_color = c("skyblue", "tomato", "gold", "orchid"),
    187.                         stroke_size = 0.4,
    188.                         set_name_size = 5)
    189.         ggsave("Venn_LB_Genes.png", plot = venn1, width = 7, height = 7, dpi = 300)
    191.         # Venn Diagram for Mac group
    192.         venn2 <- ggvenn(mac_degs_,
    193.                         fill_color = c("lightgreen", "slateblue", "plum", "orange"),
    194.                         stroke_size = 0.4,
    195.                         set_name_size = 5)
    196.         ggsave("Venn_Mac_Genes.png", plot = venn2, width = 7, height = 7, dpi = 300)
    198.         cat("✅ All Venn intersection sets exported to Excel successfully.\n")

like unlike

点赞本文的读者

还没有人对此文章表态

本文有评论

没有评论

看文章，发评论，不要沉默

最受欢迎文章

  1. Why Do Significant Gene Lists Change After Adding Additional Conditions in Differential Gene Expression Analysis?
  2. Motif Discovery in Biological Sequences: A Comparison of MEME and HOMER
  3. Calling peaks using findPeaks of HOMER
  4. Updating Human Gene Identifiers using Ensembl BioMart: A Step-by-Step Guide
  5. pheatmap vs heatmap.2
  6. Setup conda environments
  7. Should the inputs for GSVA be normalized or raw?
  8. File format for single channel analysis of Agilent microarray data with Limma?
  9. PiCRUST2 Pipeline for Functional Prediction and Pathway Analysis in Metagenomics
  10. Kraken2 Installation and Usage Guide

最新文章

最多评论文章

推荐相似文章

How to correlate RNA-seq Data with Mass Spectrometry Proteomics Data?

RNA-seq Tam on Acinetobacter baumannii strain ATCC 19606 CP059040.1 (Data_Tam_RNAseq_2024)

Deciphering S. aureus with metatranscriptomics (Marc's project)

RNA-seq 2024 Ute from raw counts

© 2023 XGenes.com Impressum