gene_x 0 like s 13 view s
Tags: pipeline
Preparing raw data
They are wildtype strains grown in different medium.
Urine - human urine
AUM - artificial urine medium
MHB - Mueller-Hinton broth
Urine（人类尿液）：pH值、比重、温度、污染物、化学成分、微生物负荷。
AUM（人工尿液培养基）：pH值、营养成分、无菌性、渗透压、温度、污染物。
MHB（Mueller-Hinton培养基）：pH值、无菌性、营养成分、温度、渗透压、抗生素浓度。
mkdir raw_data; cd raw_data
ln -s ../X101SC24105589-Z01-J001/01.RawData/AUM-1/AUM-1_1.fq.gz AUM_r1_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/AUM-1/AUM-1_2.fq.gz AUM_r1_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/AUM-2/AUM-2_1.fq.gz AUM_r2_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/AUM-2/AUM-2_2.fq.gz AUM_r2_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/AUM-3/AUM-3_1.fq.gz AUM_r3_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/AUM-3/AUM-3_2.fq.gz AUM_r3_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/MHB-1/MHB-1_1.fq.gz MHB_r1_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/MHB-1/MHB-1_2.fq.gz MHB_r1_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/MHB-2/MHB-2_1.fq.gz MHB_r2_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/MHB-2/MHB-2_2.fq.gz MHB_r2_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/MHB-3/MHB-3_1.fq.gz MHB_r3_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/MHB-3/MHB-3_2.fq.gz MHB_r3_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/Urine-1/Urine-1_1.fq.gz Urine_r1_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/Urine-1/Urine-1_2.fq.gz Urine_r1_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/Urine-2/Urine-2_1.fq.gz Urine_r2_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/Urine-2/Urine-2_2.fq.gz Urine_r2_R2.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/Urine-3/Urine-3_1.fq.gz Urine_r3_R1.fq.gz
ln -s ../X101SC24105589-Z01-J001/01.RawData/Urine-3/Urine-3_2.fq.gz Urine_r3_R2.fq.gz
(Optional) using trinity to find the most closely reference
Trinity --seqType fq --max_memory 50G --left trimmed/wt_r1_R1.fastq.gz --right trimmed/wt_r1_R2.fastq.gz --CPU 12
#https://www.genome.jp/kegg/tables/br08606.html#prok
acb KGB Acinetobacter baumannii ATCC 17978 2007 GenBank
abm KGB Acinetobacter baumannii SDF 2008 GenBank
aby KGB Acinetobacter baumannii AYE 2008 GenBank
abc KGB Acinetobacter baumannii ACICU 2008 GenBank
abn KGB Acinetobacter baumannii AB0057 2008 GenBank
abb KGB Acinetobacter baumannii AB307-0294 2008 GenBank
abx KGB Acinetobacter baumannii 1656-2 2012 GenBank
abz KGB Acinetobacter baumannii MDR-ZJ06 2012 GenBank
abr KGB Acinetobacter baumannii MDR-TJ 2012 GenBank
abd KGB Acinetobacter baumannii TCDC-AB0715 2012 GenBank
abh KGB Acinetobacter baumannii TYTH-1 2012 GenBank
abad KGB Acinetobacter baumannii D1279779 2013 GenBank
abj KGB Acinetobacter baumannii BJAB07104 2013 GenBank
abab KGB Acinetobacter baumannii BJAB0715 2013 GenBank
abaj KGB Acinetobacter baumannii BJAB0868 2013 GenBank
abaz KGB Acinetobacter baumannii ZW85-1 2013 GenBank
abk KGB Acinetobacter baumannii AbH12O-A2 2014 GenBank
abau KGB Acinetobacter baumannii AB030 2014 GenBank
abaa KGB Acinetobacter baumannii AB031 2014 GenBank
abw KGB Acinetobacter baumannii AC29 2014 GenBank
abal KGB Acinetobacter baumannii LAC-4 2015 GenBank
#Note that the Acinetobacter baumannii strain ATCC 19606 chromosome, complete genome (GenBank: CP059040.1) was choosen as reference!
Downloading CP059040.fasta and CP059040.gff from GenBank
(Optional) Preparing CP059040.fasta, CP059040_gene.gff3 and CP059040.bed
#Reference genome: https://www.ncbi.nlm.nih.gov/nuccore/CP059040
cp /media/jhuang/Elements2/Data_Tam_RNASeq3/CP059040.fasta . # Elements (Anna C.arnes)
cp /media/jhuang/Elements2/Data_Tam_RNASeq3/CP059040_gene.gff3 .
cp /media/jhuang/Elements2/Data_Tam_RNASeq3/CP059040_gene.gtf .
cp /media/jhuang/Elements2/Data_Tam_RNASeq3/CP059040.bed .
rsync -a -P CP059040.fasta jhuang@hamm:~/DATA/Data_Tam_RNAseq_2024/
rsync -a -P CP059040_gene.gff3 jhuang@hamm:~/DATA/Data_Tam_RNAseq_2024/
rsync -a -P CP059040.bed jhuang@hamm:~/DATA/Data_Tam_RNAseq_2024/
(base) jhuang@WS-2290C:/media/jhuang/Elements2/Data_Tam_RNASeq3$ find . -name "CP059040*"
./CP059040.fasta
./CP059040.bed
./CP059040.gb
./CP059040.gff3
./CP059040.gff3_backup
./CP059040_full.gb
./CP059040_gene.gff3
./CP059040_gene.gtf
./CP059040_gene_old.gff3
./CP059040_rRNA.gff3
./CP059040_rRNA_v.gff3
# ---- REF: Acinetobacter baumannii ATCC 17978 (DEBUG, gene_name failed) ----
#gffread -E -F -T GCA_000015425.1_ASM1542v1_genomic.gff -o GCA_000015425.1_ASM1542v1_genomic.gtf_
#grep "CDS" GCA_000015425.1_ASM1542v1_genomic.gtf_ > GCA_000015425.1_ASM1542v1_genomic.gtf
#sed -i -e "s/\tCDS\t/\texon\t/g" GCA_000015425.1_ASM1542v1_genomic.gtf
#gffread -E -F --bed GCA_000015425.1_ASM1542v1_genomic.gtf -o GCA_000015425.1_ASM1542v1_genomic.bed
grep "locus_tag" GCA_000015425.1_ASM1542v1_genomic.gtf_ > GCA_000015425.1_ASM1542v1_genomic.gtf
sed -i -e "s/\ttranscript\t/\texon\t/g" GCA_000015425.1_ASM1542v1_genomic.gtf # or using fc_count_type=transcript
sed -i -e "s/\tgene_name\t/\tName\t/g" GCA_000015425.1_ASM1542v1_genomic.gtf
gffread -E -F --bed GCA_000015425.1_ASM1542v1_genomic.gtf -o GCA_000015425.1_ASM1542v1_genomic.bed
#grep "gene_name" GCA_000015425.1_ASM1542v1_genomic.gtf | wc -l #69=3887-3803
cp CP059040.gff3 CP059040_backup.gff3
sed -i -e "s/\tGenbank\tgene\t/\tGenbank_gene\t/g" CP059040.gff3
grep "Genbank_gene" CP059040.gff3 > CP059040_gene.gff3
sed -i -e "s/\tGenbank_gene\t/\tGenbank\tgene\t/g" CP059040_gene.gff3
#3796-3754=42--> they are pseudogene since grep "pseudogene" CP059040.gff3 | wc -l = 42
# --------------------------------------------------------------------------------------------------------------------------------------------------
# ---------- PREPARING gff3 file including gene_biotype=protein_coding+gene_biotype=tRNA = total(3754)) and gene_biotype=pseudogene(42) ------------
cp CP059040.gff3 CP059040_backup.gff3
sed -i -e "s/\tGenbank\tgene\t/\tGenbank_gene\t/g" CP059040.gff3
grep "Genbank_gene" CP059040.gff3 > CP059040_gene.gff3
sed -i -e "s/\tGenbank_gene\t/\tGenbank\tgene\t/g" CP059040_gene.gff3
grep "gene_biotype=pseudogene" CP059040.gff3_backup >> CP059040_gene.gff3 #-->3796
#The whole point of the GTF format was to standardise certain aspects that are left open in GFF. Hence, there are many different valid ways to encode the same information in a valid GFF format, and any parser or converter needs to be written specifically for the choices the author of the GFF file made. For example, a GTF file requires the gene ID attribute to be called "gene_id", while in GFF files, it may be "ID", "Gene", something different, or completely missing.
# from gff3 to gtf
sed -i -e "s/\tID=gene-/\tgene_id \"/g" CP059040_gene.gtf
sed -i -e "s/;/\"; /g" CP059040_gene.gtf
sed -i -e "s/=/=\"/g" CP059040_gene.gtf
#sed -i -e "s/\n/\"\n/g" CP059040_gene.gtf
#using editor instead!
#The following is GTF-format.
CP000521.1 Genbank exon 95 1492 . + . transcript_id "gene0"; gene_id "gene0"; Name "A1S_0001"; gbkey "Gene"; gene_biotype "protein_coding"; locus_tag "A1S_0001";
#NZ_MJHA01000001.1 RefSeq region 1 8663 . + . ID=id0;Dbxref=taxon:575584;Name=unnamed1;collected-by=IG Schaub;collection-date=1948;country=USA: Vancouver;culture-collection=ATCC:19606;gbkey=Src;genome=plasmid;isolation-source=urine;lat-lon=37.53 N 75.4 W;map=unlocalized;mol_type=genomic DNA;nat-host=Homo sapiens;plasmid-name=unnamed1;strain=ATCC 19606;type-material=type strain of Acinetobacter baumannii
#NZ_MJHA01000001.1 RefSeq gene 228 746 . - . ID=gene0;Name=BIT33_RS00005;gbkey=Gene;gene_biotype=protein_coding;locus_tag=BIT33_RS00005;old_locus_tag=BIT33_18795
#NZ_MJHA01000001.1 Protein Homology CDS 228 746 . - 0 ID=cds0;Parent=gene0;Dbxref=Genbank:WP_000839337.1;Name=WP_000839337.1;gbkey=CDS;inference=COORDINATES: similar to AA sequence:RefSeq:WP_000839337.1;product=hypothetical protein;protein_id=WP_000839337.1;transl_table=11
##gff-version 3
##sequence-region CP059040.1 1 3980852
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=470
gffread -E -F --bed CP059040.gff3 -o CP059040.bed #-->3796
##prepare the GTF-format (see above) --> ERROR! ----> using CP059040.gff3
##stringtie adeIJ.abx_r1.sorted.bam -o adeIJ.abx_r1.sorted_transcripts.gtf -v -G /media/jhuang/Elements/Data_Tam_RNASeq3/CP059040.gff3 -A adeIJ.abx_r1.sorted.gene_abund.txt -C adeIJ.abx_r1.sorted.bam.cov_refs.gtf -e -b adeIJ.abx_r1.sorted_ballgown
#[01/21 10:57:46] Loading reference annotation (guides)..
#GFF warning: merging adjacent/overlapping segments of gene-H0N29_00815 on CP059040.1 (179715-179786, 179788-180810)
#[01/21 10:57:46] 3796 reference transcripts loaded.
#Default stack size for threads: 8388608
#WARNING: no reference transcripts found for genomic sequence "gi|1906906720|gb|CP059040.1|"! (mismatched reference names?)
#WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!
#Please make sure the -G annotation file uses the same naming convention for the genome sequences.
#[01/21 10:58:30] All threads finished.
# ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
# The specified gene identifier attribute is 'Name'
# An example of attributes included in your GTF annotation is 'ID=exon-H0N29_00075-1;Parent=rna-H0N29_00075;gbkey=rRNA;locus_tag=H0N29_00075;product=16S ribosomal RNA'
# The program has to termin
# ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
# The specified gene identifier attribute is 'gene_biotype'
# An example of attributes included in your GTF annotation is 'ID=exon-H0N29_00075-1;Parent=rna-H0N29_00075;gbkey=rRNA;locus_tag=H0N29_00075;product=16S ribosomal RNA'
# The program has to terminate.
#grep "ID=cds-" CP059040.gff3 | wc -l
#grep "ID=exon-" CP059040.gff3 | wc -l
#grep "ID=gene-" CP059040.gff3 | wc -l #the same as H0N29_18980/5=3796
grep "gbkey=" CP059040.gff3 | wc -l 7695
grep "ID=id-" CP059040.gff3 | wc -l 5
grep "locus_tag=" CP059040.gff3 | wc -l 7689
#...
cds 3701 locus_tag=xxxx, no gene_biotype
exon 96 locus_tag=xxxx, no gene_biotype
gene 3796 locus_tag=xxxx, gene_biotype=xxxx,
id (riboswitch+direct_repeat,5) both no --> ignoring them!! # grep "ID=id-" CP059040.gff3
rna 96 locus_tag=xxxx, no gene_biotype
------------------
7694
cp CP059040.gff3_backup CP059040.gff3
grep "^##" CP059040.gff3 > CP059040_gene.gff3
grep "ID=gene" CP059040.gff3 >> CP059040_gene.gff3
#!!!!VERY_IMPORTANT!!!!: change type '\tCDS\t' to '\texon\t'!
sed -i -e "s/\tgene\t/\texon\t/g" CP059040_gene.gff3
Preparing the directory trimmed
mkdir trimmed trimmed_unpaired;
for sample_id in AUM_r1 AUM_r2 AUM_r3 Urine_r1 Urine_r2 Urine_r3 MHB_r1 MHB_r2 MHB_r3; do \
for sample_id in MHB_r1 MHB_r2 MHB_r3; do \
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;
done
Preparing samplesheet.csv
sample,fastq_1,fastq_2,strandedness
AUM_r1,AUM_r1_R1.fq.gz,AUM_r1_R2.fq.gz,auto
AUM_r2,AUM_r2_R1.fq.gz,AUM_r2_R2.fq.gz,auto
AUM_r3,AUM_r3_R1.fq.gz,AUM_r3_R2.fq.gz,auto
MHB_r1,MHB_r1_R1.fq.gz,MHB_r1_R2.fq.gz,auto
MHB_r2,MHB_r2_R1.fq.gz,MHB_r2_R2.fq.gz,auto
MHB_r3,MHB_r3_R1.fq.gz,MHB_r3_R2.fq.gz,auto
Urine_r1,Urine_r1_R1.fq.gz,Urine_r1_R2.fq.gz,auto
Urine_r2,Urine_r2_R1.fq.gz,Urine_r2_R2.fq.gz,auto
Urine_r3,Urine_r3_R1.fq.gz,Urine_r3_R2.fq.gz,auto
nextflow run
#Example1: http://xgenes.com/article/article-content/157/prepare-virus-gtf-for-nextflow-run/
docker pull nfcore/rnaseq
ln -s /home/jhuang/Tools/nf-core-rnaseq-3.12.0/ rnaseq
#Default: --gtf_group_features 'gene_id' --gtf_extra_attributes 'gene_name' --featurecounts_group_type 'gene_biotype' --featurecounts_feature_type 'exon'
#(host_env) !NOT_WORKING! jhuang@WS-2290C:~/DATA/Data_Tam_RNAseq_2024$ /usr/local/bin/nextflow run rnaseq/main.nf --input samplesheet.csv --outdir results --fasta "/home/jhuang/DATA/Data_Tam_RNAseq_2024/CP059040.fasta" --gff "/home/jhuang/DATA/Data_Tam_RNAseq_2024/CP059040.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'
# -- DEBUG_1 (CDS --> exon in CP059040.gff) --
#Checking the record (see below) in results/genome/CP059040.gtf
#In ./results/genome/CP059040.gtf e.g. "CP059040.1 Genbank transcript 1 1398 . + . transcript_id "gene-H0N29_00005"; gene_id "gene-H0N29_00005"; gene_name "dnaA"; Name "dnaA"; gbkey "Gene"; gene "dnaA"; gene_biotype "protein_coding"; locus_tag "H0N29_00005";"
#--featurecounts_feature_type 'transcript' returns only the tRNA results
#Since the tRNA records have "transcript and exon". In gene records, we have "transcript and CDS". replace the CDS with exon
grep -P "\texon\t" CP059040.gff | sort | wc -l #96
grep -P "cmsearch\texon\t" CP059040.gff | wc -l #=10 ignal recognition particle sRNA small typ, transfer-messenger RNA, 5S ribosomal RNA
grep -P "Genbank\texon\t" CP059040.gff | wc -l #=12 16S and 23S ribosomal RNA
grep -P "tRNAscan-SE\texon\t" CP059040.gff | wc -l #tRNA 74
wc -l star_salmon/AUM_r3/quant.genes.sf #--featurecounts_feature_type 'transcript' results in 96 records!
grep -P "\tCDS\t" CP059040.gff | wc -l #3701
sed 's/\tCDS\t/\texon\t/g' CP059040.gff > CP059040_m.gff
grep -P "\texon\t" CP059040_m.gff | sort | wc -l #3797
# -- DEBUG_2: combination of 'CP059040_m.gff' and 'exon' results in ERROR, using 'transcript' instead!
--gff "/home/jhuang/DATA/Data_Tam_RNAseq_2024/CP059040_m.gff" --featurecounts_feature_type 'transcript'
# ---- SUCCESSFUL with directly downloaded gff3 and fasta from NCBI using docker after replacing 'CDS' with 'exon' ----
(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'
# -- DEBUG_3: make sure the header of fasta is the same to the *_m.gff file
Import data and pca-plot
#mamba activate r_env
#install.packages("ggfun")
# Import the required libraries
library("AnnotationDbi")
library("clusterProfiler")
library("ReactomePA")
library(gplots)
library(tximport)
library(DESeq2)
#library("org.Hs.eg.db")
library(dplyr)
library(tidyverse)
#install.packages("devtools")
#devtools::install_version("gtable", version = "0.3.0")
library(gplots)
library("RColorBrewer")
#install.packages("ggrepel")
library("ggrepel")
# install.packages("openxlsx")
library(openxlsx)
library(EnhancedVolcano)
library(DESeq2)
library(edgeR)
setwd("~/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/results/star_salmon")
# Define paths to your Salmon output quantification files
files <- c("Urine_r1" = "./Urine_r1/quant.sf",
"Urine_r2" = "./Urine_r2/quant.sf",
"Urine_r3" = "./Urine_r3/quant.sf",
"MHB_r1" = "./MHB_r1/quant.sf",
"MHB_r2" = "./MHB_r2/quant.sf",
"MHB_r3" = "./MHB_r3/quant.sf")
# Import the transcript abundance data with tximport
txi <- tximport(files, type = "salmon", txIn = TRUE, txOut = TRUE)
# Define the replicates and condition of the samples
replicate <- factor(c("r1", "r2", "r3", "r1", "r2", "r3"))
condition <- factor(c("Urine","Urine","Urine", "MHB","MHB","MHB"))
# Define the colData for DESeq2
colData <- data.frame(condition=condition, replicate=replicate, row.names=names(files))
# ------------------------
# 1️⃣ Setup and input files
# ------------------------
# Read in transcript-to-gene mapping
tx2gene <- read.table("salmon_tx2gene.tsv", header=FALSE, stringsAsFactors=FALSE)
colnames(tx2gene) <- c("transcript_id", "gene_id", "gene_name")
# Prepare tx2gene for gene-level summarization (remove gene_name if needed)
tx2gene_geneonly <- tx2gene[, c("transcript_id", "gene_id")]
# -------------------------------
# 2️⃣ Transcript-level counts
# -------------------------------
# Create DESeqDataSet directly from tximport (transcript-level)
dds_tx <- DESeqDataSetFromTximport(txi, colData=colData, design=~condition)
write.csv(counts(dds_tx), file="transcript_counts.csv")
# --------------------------------
# 3️⃣ Gene-level summarization
# --------------------------------
# Re-import Salmon data summarized at gene level
txi_gene <- tximport(files, type="salmon", tx2gene=tx2gene_geneonly, txOut=FALSE)
# Create DESeqDataSet for gene-level counts
dds <- DESeqDataSetFromTximport(txi_gene, colData=colData, design=~condition+replicate)
# --------------------------------
# 4️⃣ Raw counts table (with gene names)
# --------------------------------
# Extract raw gene-level counts
counts_data <- as.data.frame(counts(dds, normalized=FALSE))
counts_data$gene_id <- rownames(counts_data)
# Add gene names
tx2gene_unique <- unique(tx2gene[, c("gene_id", "gene_name")])
counts_data <- merge(counts_data, tx2gene_unique, by="gene_id", all.x=TRUE)
# Reorder columns: gene_id, gene_name, then counts
count_cols <- setdiff(colnames(counts_data), c("gene_id", "gene_name"))
counts_data <- counts_data[, c("gene_id", "gene_name", count_cols)]
# --------------------------------
# 5️⃣ Calculate CPM
# --------------------------------
library(edgeR)
library(openxlsx)
# Prepare count matrix for CPM calculation
count_matrix <- as.matrix(counts_data[, !(colnames(counts_data) %in% c("gene_id", "gene_name"))])
# Calculate CPM
#cpm_matrix <- cpm(count_matrix, normalized.lib.sizes=FALSE)
total_counts <- colSums(count_matrix)
cpm_matrix <- t(t(count_matrix) / total_counts) * 1e6
cpm_matrix <- as.data.frame(cpm_matrix)
# Add gene_id and gene_name back to CPM table
cpm_counts <- cbind(counts_data[, c("gene_id", "gene_name")], cpm_matrix)
# --------------------------------
# 6️⃣ Save outputs
# --------------------------------
write.csv(counts_data, "gene_raw_counts.csv", row.names=FALSE)
write.xlsx(counts_data, "gene_raw_counts.xlsx", row.names=FALSE)
write.xlsx(cpm_counts, "gene_cpm_counts.xlsx", row.names=FALSE)
PCA dim(counts(dds)) head(counts(dds), 10) rld <- rlogTransformation(dds)
# draw simple pca and heatmap
#mat <- assay(rld)
#mm <- model.matrix(~condition, colData(rld))
#mat <- limma::removeBatchEffect(mat, batch=rld$batch, design=mm)
#assay(rld) <- mat
# -- pca --
png("pca.png", 1200, 800)
plotPCA(rld, intgroup=c("condition"))
dev.off()
# -- heatmap --
png("heatmap.png", 1200, 800)
distsRL <- dist(t(assay(rld)))
mat <- as.matrix(distsRL)
hc <- hclust(distsRL)
hmcol <- colorRampPalette(brewer.pal(9,"GnBu"))(100)
heatmap.2(mat, Rowv=as.dendrogram(hc),symm=TRUE, trace="none",col = rev(hmcol), margin=c(13, 13))
dev.off()
Select the differentially expressed genes
#https://galaxyproject.eu/posts/2020/08/22/three-steps-to-galaxify-your-tool/
#https://www.biostars.org/p/282295/
#https://www.biostars.org/p/335751/
#> dds$condition
#[1] Urine Urine Urine MHB MHB MHB
#Levels: MHB Urine
#CONSOLE: mkdir star_salmon/degenes
setwd("degenes")
#---- relevel to control ----
dds$condition <- relevel(dds$condition, "MHB")
dds = DESeq(dds, betaPrior=FALSE)
resultsNames(dds)
clist <- c("Urine_vs_MHB")
for (i in clist) {
contrast = paste("condition", i, sep="_")
res = results(dds, name=contrast)
res <- res[!is.na(res$log2FoldChange),]
res_df <- as.data.frame(res)
write.csv(as.data.frame(res_df[order(res_df$pvalue),]), file = paste(i, "all.txt", sep="-"))
up <- subset(res_df, padj<=0.05 & log2FoldChange>=1.35)
down <- subset(res_df, padj<=0.05 & log2FoldChange<=-1.35)
write.csv(as.data.frame(up[order(up$log2FoldChange,decreasing=TRUE),]), file = paste(i, "up.txt", sep="-"))
write.csv(as.data.frame(down[order(abs(down$log2FoldChange),decreasing=TRUE),]), file = paste(i, "down.txt", sep="-"))
}
# -- Under host-env --
grep -P "\tgene\t" CP059040.gff > CP059040_gene.gff
python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/CP059040_gene.gff Urine_vs_MHB-all.txt Urine_vs_MHB-all.csv
python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/CP059040_gene.gff Urine_vs_MHB-up.txt Urine_vs_MHB-up.csv
python3 ~/Scripts/replace_gene_names.py /home/jhuang/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/CP059040_gene.gff Urine_vs_MHB-down.txt Urine_vs_MHB-down.csv
res <- read.csv("Urine_vs_MHB-all.csv")
# Replace empty GeneName with modified GeneID
res$GeneName <- ifelse(
res$GeneName == "" | is.na(res$GeneName),
gsub("gene-", "", res$GeneID),
res$GeneName
)
duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
res <- res %>%
group_by(GeneName) %>%
slice_min(padj, with_ties = FALSE) %>%
ungroup()
res <- as.data.frame(res)
# Sort res first by padj (ascending) and then by log2FoldChange (descending)
res <- res[order(res$padj, -res$log2FoldChange), ]
# Assuming res is your dataframe and already processed
# Filter up-regulated genes: log2FoldChange > 2 and padj < 1e-2
up_regulated <- res[res$log2FoldChange > 2 & res$padj < 1e-2, ]
# Filter down-regulated genes: log2FoldChange < -2 and padj < 1e-2
down_regulated <- res[res$log2FoldChange < -2 & res$padj < 1e-2, ]
# Create a new workbook
wb <- createWorkbook()
# Add the complete dataset as the first sheet
addWorksheet(wb, "Complete_Data")
writeData(wb, "Complete_Data", res)
# Add the up-regulated genes as the second sheet
addWorksheet(wb, "Up_Regulated")
writeData(wb, "Up_Regulated", up_regulated)
# Add the down-regulated genes as the third sheet
addWorksheet(wb, "Down_Regulated")
writeData(wb, "Down_Regulated", down_regulated)
# Save the workbook to a file
saveWorkbook(wb, "Gene_Expression_Urine_vs_MHB.xlsx", overwrite = TRUE)
# Set the 'GeneName' column as row.names
rownames(res) <- res$GeneName
# Drop the 'GeneName' column since it's now the row names
res$GeneName <- NULL
head(res)
## Ensure the data frame matches the expected format
## For example, it should have columns: log2FoldChange, padj, etc.
#res <- as.data.frame(res)
## Remove rows with NA in log2FoldChange (if needed)
#res <- res[!is.na(res$log2FoldChange),]
# Replace padj = 0 with a small value
res$padj[res$padj == 0] <- 1e-305
#library(EnhancedVolcano)
# Assuming res is already sorted and processed
png("Urine_vs_MHB.png", width=1200, height=2000)
#max.overlaps = 10
EnhancedVolcano(res,
lab = rownames(res),
x = 'log2FoldChange',
y = 'padj',
pCutoff = 1e-2,
FCcutoff = 2,
title = '',
subtitleLabSize = 18,
pointSize = 3.0,
labSize = 5.0,
colAlpha = 1,
legendIconSize = 4.0,
drawConnectors = TRUE,
widthConnectors = 0.5,
colConnectors = 'black',
subtitle = expression("Urine versus MHB"))
dev.off()
KEGG and GO annotations in non-model organisms
https://www.biobam.com/functional-analysis/
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.
Preparing file 1 eggnog_out.emapper.annotations.txt for the R-code below: (KEGG Terms): EggNog based on orthology and phylogenies
EggNOG-mapper assigns both KEGG Orthology (KO) IDs and GO terms.
Install EggNOG-mapper:
mamba create -n eggnog_env python=3.8 eggnog-mapper -c conda-forge -c bioconda #eggnog-mapper_2.1.12
mamba activate eggnog_env
Run annotation:
#diamond makedb --in eggnog6.prots.faa -d eggnog_proteins.dmnd
mkdir /home/jhuang/mambaforge/envs/eggnog_env/lib/python3.8/site-packages/data/
download_eggnog_data.py --dbname eggnog.db -y --data_dir /home/jhuang/mambaforge/envs/eggnog_env/lib/python3.8/site-packages/data/
#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
python ~/Scripts/update_fasta_header.py CP059040_protein_.fasta CP059040_protein.fasta
emapper.py -i CP059040_protein.fasta -o eggnog_out --cpu 60 --resume
#----> result annotations.tsv: Contains KEGG, GO, and other functional annotations.
#----> 470.IX87_14445:
* 470 likely refers to the organism or strain (e.g., Acinetobacter baumannii ATCC 19606 or another related strain).
* IX87_14445 would refer to a specific gene or protein within that genome.
Extract KEGG KO IDs from annotations.emapper.annotations.
Preparing file 2 blast2go_annot.annot2_ for the R-code below:
Basic (GO Terms from 'Blast2GO 5 Basic', saved in blast2go_annot.annot): Using Blast/Diamond + Blast2GO_GUI based on sequence alignment + GO mapping
Advanced (GO Terms from 'Blast2GO 5 Basic'): Interpro based protein families / domains --> Button interpro
MERGE the results of InterPro GO IDs (advanced) to GO IDs (basic) and generate final GO IDs, saved in blast2go_annot.annot2
PREPARING go_terms and ec_terms: annot_* file:
cut -f1-2 -d$'\t' blast2go_annot.annot2 > blast2go_annot.annot2_
Perform KEGG and GO Enrichment in R
#BiocManager::install("GO.db")
#BiocManager::install("AnnotationDbi")
# Load required libraries
library(openxlsx) # For Excel file handling
library(dplyr) # For data manipulation
library(tidyr)
library(stringr)
library(clusterProfiler) # For KEGG and GO enrichment analysis
#library(org.Hs.eg.db) # Replace with appropriate organism database
library(GO.db)
library(AnnotationDbi)
setwd("~/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/results/star_salmon/degenes")
# Step 1: Load the blast2go annotation file with a check for missing columns
annot_df <- read.table("/home/jhuang/b2gWorkspace_Tam_RNAseq_2024/blast2go_annot.annot2_",
header = FALSE, sep = "\t", stringsAsFactors = FALSE, fill = TRUE)
# If the structure is inconsistent, we can make sure there are exactly 3 columns:
colnames(annot_df) <- c("GeneID", "Term")
# Step 2: Filter and aggregate GO and EC terms as before
go_terms <- annot_df %>%
filter(grepl("^GO:", Term)) %>%
group_by(GeneID) %>%
summarize(GOs = paste(Term, collapse = ","), .groups = "drop")
ec_terms <- annot_df %>%
filter(grepl("^EC:", Term)) %>%
group_by(GeneID) %>%
summarize(EC = paste(Term, collapse = ","), .groups = "drop")
# Load the results
res <- read.csv("Urine_vs_MHB-all.csv") #up259, down138
# Replace empty GeneName with modified GeneID
res$GeneName <- ifelse(
res$GeneName == "" | is.na(res$GeneName),
gsub("gene-", "", res$GeneID),
res$GeneName
)
# Remove duplicated genes by selecting the gene with the smallest padj
duplicated_genes <- res[duplicated(res$GeneName), "GeneName"]
res <- res %>%
group_by(GeneName) %>%
slice_min(padj, with_ties = FALSE) %>%
ungroup()
res <- as.data.frame(res)
# Sort res first by padj (ascending) and then by log2FoldChange (descending)
res <- res[order(res$padj, -res$log2FoldChange), ]
# Read eggnog annotations
eggnog_data <- read.delim("~/DATA/Data_Tam_RNAseq_2024_AUM_MHB_Urine_ATCC19606/eggnog_out.emapper.annotations.txt", header = TRUE, sep = "\t")
# Remove the "gene-" prefix from GeneID in res to match eggnog 'query' format
res$GeneID <- gsub("gene-", "", res$GeneID)
# Merge eggnog data with res based on GeneID
res <- res %>%
left_join(eggnog_data, by = c("GeneID" = "query"))
# DEBUG: NOT_NECESSARY, since res has already GeneName
##Convert row names to a new column 'GeneName' in res
#res_with_geneName <- res %>%
#mutate(GeneName = rownames(res)) %>%
#as.data.frame() # Ensure that it's a regular data frame without row names
## View the result
#head(res_with_geneName)
# Merge with the res dataframe
# Perform the left joins and rename columns
res_updated <- res %>%
left_join(go_terms, by = "GeneID") %>%
left_join(ec_terms, by = "GeneID") %>% dplyr::select(-EC.x, -GOs.x) %>% dplyr::rename(EC = EC.y, GOs = GOs.y)
# DEBUG: NOT_NECESSARY, since 'GeneName' is already the first column.
## Reorder columns to move 'GeneName' as the first column in res_updated
#res_updated <- res_updated %>%
#select(GeneName, everything())
## Count the number of rows in the KEGG_ko, GOs, EC columns that have non-missing values
#num_non_missing_KEGG_ko <- sum(res_updated$KEGG_ko != "-" & !is.na(res_updated$KEGG_ko))
#print(num_non_missing_KEGG_ko)
##[1] 2030
#num_non_missing_GOs <- sum(res_updated$GOs != "-" & !is.na(res_updated$GOs))
#print(num_non_missing_GOs)
##[1] 2865 --> 2875
#num_non_missing_EC <- sum(res_updated$EC != "-" & !is.na(res_updated$EC))
#print(num_non_missing_EC)
##[1] 1701
# Filter up-regulated genes
up_regulated <- res_updated[res_updated$log2FoldChange > 2 & res_updated$padj < 0.01, ]
# Filter down-regulated genes
down_regulated <- res_updated[res_updated$log2FoldChange < -2 & res_updated$padj < 0.01, ]
# Create a new workbook
wb <- createWorkbook()
# Add the complete dataset as the first sheet (with annotations)
addWorksheet(wb, "Complete_Data")
writeData(wb, "Complete_Data", res_updated)
# Add the up-regulated genes as the second sheet (with annotations)
addWorksheet(wb, "Up_Regulated")
writeData(wb, "Up_Regulated", up_regulated)
# Add the down-regulated genes as the third sheet (with annotations)
addWorksheet(wb, "Down_Regulated")
writeData(wb, "Down_Regulated", down_regulated)
# Save the workbook to a file
saveWorkbook(wb, "Gene_Expression_with_Annotations_Urine_vs_MHB.xlsx", overwrite = TRUE)
# Set GeneName as row names after the join
rownames(res_updated) <- res_updated$GeneName
res_updated <- res_updated %>% dplyr::select(-GeneName)
## Set the 'GeneName' column as row.names
#rownames(res_updated) <- res_updated$GeneName
## Drop the 'GeneName' column since it's now the row names
#res_updated$GeneName <- NULL
# ---- Perform KEGG enrichment analysis (up_regulated) ----
gene_list_kegg_up <- up_regulated$KEGG_ko
gene_list_kegg_up <- gsub("ko:", "", gene_list_kegg_up)
kegg_enrichment_up <- enrichKEGG(gene = gene_list_kegg_up, organism = 'ko')
# -- convert the GeneID (Kxxxxxx) to the true GeneID --
# Step 0: Create KEGG to GeneID mapping
kegg_to_geneid_up <- up_regulated %>%
dplyr::select(KEGG_ko, GeneID) %>%
filter(!is.na(KEGG_ko)) %>% # Remove missing KEGG KO entries
mutate(KEGG_ko = str_remove(KEGG_ko, "ko:")) # Remove 'ko:' prefix if present
# Step 1: Clean KEGG_ko values (separate multiple KEGG IDs)
kegg_to_geneid_clean <- kegg_to_geneid_up %>%
mutate(KEGG_ko = str_remove_all(KEGG_ko, "ko:")) %>% # Remove 'ko:' prefixes
separate_rows(KEGG_ko, sep = ",") %>% # Ensure each KEGG ID is on its own row
filter(KEGG_ko != "-") %>% # Remove invalid KEGG IDs ("-")
distinct() # Remove any duplicate mappings
# Step 2.1: Expand geneID column in kegg_enrichment_up
expanded_kegg <- kegg_enrichment_up %>%
as.data.frame() %>%
separate_rows(geneID, sep = "/") %>% # Split multiple KEGG IDs (Kxxxxx)
left_join(kegg_to_geneid_clean, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>% # Explicitly handle many-to-many
distinct() %>% # Remove duplicate matches
group_by(ID) %>%
summarise(across(everything(), ~ paste(unique(na.omit(.)), collapse = "/")), .groups = "drop") # Re-collapse results
#dplyr::glimpse(expanded_kegg)
# Step 3.1: Replace geneID column in the original dataframe
kegg_enrichment_up_df <- as.data.frame(kegg_enrichment_up)
# Remove old geneID column and merge new one
kegg_enrichment_up_df <- kegg_enrichment_up_df %>%
dplyr::select(-geneID) %>% # Remove old geneID column
left_join(expanded_kegg %>% dplyr::select(ID, GeneID), by = "ID") %>% # Merge new GeneID column
dplyr::rename(geneID = GeneID) # Rename column back to geneID
# ---- Perform KEGG enrichment analysis (down_regulated) ----
# Step 1: Extract KEGG KO terms from down-regulated genes
gene_list_kegg_down <- down_regulated$KEGG_ko
gene_list_kegg_down <- gsub("ko:", "", gene_list_kegg_down)
# Step 2: Perform KEGG enrichment analysis
kegg_enrichment_down <- enrichKEGG(gene = gene_list_kegg_down, organism = 'ko')
# --- Convert KEGG gene IDs (Kxxxxxx) to actual GeneIDs ---
# Step 3: Create KEGG to GeneID mapping from down_regulated dataset
kegg_to_geneid_down <- down_regulated %>%
dplyr::select(KEGG_ko, GeneID) %>%
filter(!is.na(KEGG_ko)) %>% # Remove missing KEGG KO entries
mutate(KEGG_ko = str_remove(KEGG_ko, "ko:")) # Remove 'ko:' prefix if present
# Step 4: Clean KEGG_ko values (handle multiple KEGG IDs)
kegg_to_geneid_down_clean <- kegg_to_geneid_down %>%
mutate(KEGG_ko = str_remove_all(KEGG_ko, "ko:")) %>% # Remove 'ko:' prefixes
separate_rows(KEGG_ko, sep = ",") %>% # Ensure each KEGG ID is on its own row
filter(KEGG_ko != "-") %>% # Remove invalid KEGG IDs ("-")
distinct() # Remove duplicate mappings
# Step 5: Expand geneID column in kegg_enrichment_down
expanded_kegg_down <- kegg_enrichment_down %>%
as.data.frame() %>%
separate_rows(geneID, sep = "/") %>% # Split multiple KEGG IDs (Kxxxxx)
left_join(kegg_to_geneid_down_clean, by = c("geneID" = "KEGG_ko"), relationship = "many-to-many") %>% # Handle many-to-many mappings
distinct() %>% # Remove duplicate matches
group_by(ID) %>%
summarise(across(everything(), ~ paste(unique(na.omit(.)), collapse = "/")), .groups = "drop") # Re-collapse results
# Step 6: Replace geneID column in the original kegg_enrichment_down dataframe
kegg_enrichment_down_df <- as.data.frame(kegg_enrichment_down) %>%
dplyr::select(-geneID) %>% # Remove old geneID column
left_join(expanded_kegg_down %>% dplyr::select(ID, GeneID), by = "ID") %>% # Merge new GeneID column
dplyr::rename(geneID = GeneID) # Rename column back to geneID
# View the updated dataframe
head(kegg_enrichment_down_df)
# Create a new workbook
wb <- createWorkbook()
# Save enrichment results to the workbook
addWorksheet(wb, "KEGG_Enrichment_Up")
writeData(wb, "KEGG_Enrichment_Up", as.data.frame(kegg_enrichment_up_df))
# Save enrichment results to the workbook
addWorksheet(wb, "KEGG_Enrichment_Down")
writeData(wb, "KEGG_Enrichment_Down", as.data.frame(kegg_enrichment_down_df))
saveWorkbook(wb, "KEGG_Enrichment.xlsx", overwrite = TRUE)
# ---- Perform GO enrichment analysis (TODO: extract the merged GO IDs from 'Blast2GO 5 Basic' and adapt the code below!)----
# Define gene list (up-regulated genes)
gene_list_go_up <- up_regulated$GeneID # Extract the 149 up-regulated genes
gene_list_go_down <- down_regulated$GeneID # Extract the 65 down-regulated genes
# Define background gene set (all genes in res)
background_genes <- res_updated$GeneID # Extract the 3646 background genes
# Prepare GO annotation data from res
go_annotation <- res_updated[, c("GOs","GeneID")] # Extract relevant columns
go_annotation <- go_annotation %>%
tidyr::separate_rows(GOs, sep = ",") # Split multiple GO terms into separate rows
# Perform GO enrichment analysis, where pAdjustMethod is one of "holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none"
go_enrichment_up <- enricher(
gene = gene_list_go_up, # Up-regulated genes
TERM2GENE = go_annotation, # Custom GO annotation
pvalueCutoff = 1.0, # Significance threshold
pAdjustMethod = "BH",
universe = background_genes # Define the background gene set
)
go_enrichment_up <- as.data.frame(go_enrichment_up)
go_enrichment_down <- enricher(
gene = gene_list_go_down, # Up-regulated genes
TERM2GENE = go_annotation, # Custom GO annotation
pvalueCutoff = 1.0, # Significance threshold
pAdjustMethod = "BH",
universe = background_genes # Define the background gene set
)
go_enrichment_down <- as.data.frame(go_enrichment_down)
## Remove the 'p.adjust' column since no adjusted methods have been applied!
#go_enrichment_up <- go_enrichment_up[, !names(go_enrichment_up) %in% "p.adjust"]
# Update the Description column with the term descriptions
go_enrichment_up$Description <- sapply(go_enrichment_up$ID, function(go_id) {
# Using select to get the term description
term <- tryCatch({
AnnotationDbi::select(GO.db, keys = go_id, columns = "TERM", keytype = "GOID")
}, error = function(e) {
message(paste("Error for GO term:", go_id)) # Print which GO ID caused the error
return(data.frame(TERM = NA)) # In case of error, return NA
})
if (nrow(term) > 0) {
return(term$TERM)
} else {
return(NA) # If no description found, return NA
}
})
## Print the updated data frame
#print(go_enrichment_up)
## Remove the 'p.adjust' column since no adjusted methods have been applied!
#go_enrichment_down <- go_enrichment_down[, !names(go_enrichment_down) %in% "p.adjust"]
# Update the Description column with the term descriptions
go_enrichment_down$Description <- sapply(go_enrichment_down$ID, function(go_id) {
# Using select to get the term description
term <- tryCatch({
AnnotationDbi::select(GO.db, keys = go_id, columns = "TERM", keytype = "GOID")
}, error = function(e) {
message(paste("Error for GO term:", go_id)) # Print which GO ID caused the error
return(data.frame(TERM = NA)) # In case of error, return NA
})
if (nrow(term) > 0) {
return(term$TERM)
} else {
return(NA) # If no description found, return NA
}
})
addWorksheet(wb, "GO_Enrichment_Up")
writeData(wb, "GO_Enrichment_Up", as.data.frame(go_enrichment_up))
addWorksheet(wb, "GO_Enrichment_Down")
writeData(wb, "GO_Enrichment_Down", as.data.frame(go_enrichment_down))
# Save the workbook with enrichment results
saveWorkbook(wb, "KEGG_and_GO_Enrichments_Urine_vs_MHB.xlsx", overwrite = TRUE)
#Error for GO term: GO:0006807: replace GO:0006807 obsolete nitrogen compound metabolic process
#TODO: marked the color as yellow if the p.adjusted <= 0.05 in GO_enrichment!
Finalizing the KEGG and GO Enrichment table
1. NOTE: geneIDs in KEGG_Enrichment have been already translated from ko to geneID in H0N29_*-format;
2. NEED_MANUAL_DELETION: p.adjust values have been calculated, we have to filter all records in GO_Enrichment-results by |p.adjust|<=0.05.
点赞本文的读者
还没有人对此文章表态
没有评论
Workflow using PICRUSt2 for Data_Karoline_16S_2025
Viral genome assembly and recombination analysis for Data_Sophie_HDV_Sequences
© 2023 XGenes.com Impressum