| batch | condition | media | treatment | induced | exp_5 | |
|---|---|---|---|---|---|---|
| SK064 | 1 | LB-IPTG | LB | DMSO | yes | control |
| SK065 | 2 | LB-IPTG | LB | DMSO | yes | control |
| SK066 | 3 | LB-IPTG | LB | DMSO | yes | control |
| SK076 | 1 | LB-Eb-IPTG | LB | Eb | yes | treatment |
| SK077 | 2 | LB-Eb-IPTG | LB | Eb | yes | treatment |
| SK078 | 3 | LB-Eb-IPTG | LB | Eb | yes | treatment |
## batch condition media treatment induced exp_5
## SK064 1 LB-IPTG LB DMSO yes control
## SK065 2 LB-IPTG LB DMSO yes control
## SK066 3 LB-IPTG LB DMSO yes control
## SK076 1 LB-Eb-IPTG LB Eb yes treatment
## SK077 2 LB-Eb-IPTG LB Eb yes treatment
## SK078 3 LB-Eb-IPTG LB Eb yes treatment
## [1] "The control is: LB-IPTG"
## [1] "The treatment is: LB-Eb-IPTG"
# Load raw count data
rawcounts <- read.csv("../../raw-data/htseq_SK_unfiltered_reverse_strand_count.csv",
row.names = 1, header = TRUE)
# Append a prefix to the batch and run columns to turn to character class
metadata[["batch"]] <- paste0("b", metadata[["batch"]])
metadata[["run"]] <- paste0("r", metadata[["run"]])
# Subset the rawcounts dataframe to only include samples that are present in experiment
samples <- rownames(metadata)
rawcounts <- select(rawcounts, all_of(samples))
# Ensure that the order of samples in metadata matches the order of columns in rawcounts
all(rownames(metadata) == colnames(rawcounts)) # Should return TRUE if order is correct## [1] TRUE
# Define the design formula for DESeq2, including main condition and batch effect
design_formula <- ~ exp_5 + batch
# Create DESeq2 dataset object from the filtered count matrix and metadata
dds <- DESeqDataSetFromMatrix(
countData = rawcounts,
colData = metadata,
design = design_formula
)## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
# Sum up the counts for each sample
rawcounts_sums <- colSums(rawcounts)
# Turn it into a data frame
data <- data.frame(
Category = names(rawcounts_sums),
Count = as.numeric(rawcounts_sums)
)
# Make sure sample order is consistent
data <- data %>%
mutate(Category = factor(Category, levels = sort(unique(Category))))
# Plot bar chart of library sizes
p <- ggplot(data, aes(x = Category, y = Count, fill = Category)) +
geom_bar(stat = "identity") +
# coord_flip() # flip if too many samples to read labels easily
labs(y = "Library size (M)", x = NULL) + # remove x label, rename y axis
# scale_fill_brewer(palette = "Dark2") + # optional color palette
scale_y_continuous(
labels = scales::label_number(scale = 1e-6), # show counts in millions
limits = c(0, (max(data$Count) + 1e6)) # add a bit of space on top
) +
theme(
panel.grid = element_blank(),
panel.background = element_blank(),
axis.text = element_text(size = text.size, angle = 45, hjust = 1), # slant labels
axis.title = element_text(size = axis.text.size),
legend.text = element_text(size = text.size),
legend.title = element_text(size = axis.text.size),
panel.border = element_rect(fill = NA),
legend.position = "none", # no legend needed since each bar is a different sample
axis.ticks = element_line(linewidth = 0.2),
legend.key = element_blank()
)
# Show the plot
p# Save the plot as a PDF with today’s date and project name
ggsave(paste0(Sys.Date(), "_library_size_", exp, "_", treatment,".vs.", control, ".pdf"), device = "pdf")## Saving 7 x 5 in image
# Stack rawcounts into single vector
dat <- stack(rawcounts)
# Format expression distribution plot with the expression level and density
expression_dist <- ggplot(dat, aes(x = log2(values+1), fill=ind, label = ind )) +
geom_density(alpha = 0.5, position = "identity") +
ylab("Density") +
geom_textdensity(hjust = "ymax", vjust = 0,
text_only = TRUE, text_smoothing = 20) +
theme(panel.grid = element_blank(),
panel.background = element_blank(),
axis.text = element_text(size = text.size),
axis.title = element_text(size = axis.text.size),
legend.text = element_text(size = text.size),
legend.title = element_text(size = axis.text.size),
panel.border = element_rect(fill = NA),
legend.position = "none",
axis.ticks = element_line(linewidth = 0.2),
legend.key = element_blank()) +
coord_cartesian(xlim = c(0, 15))
expression_dist#save figure
cairo_pdf(paste0(Sys.Date(), "_expression_dist_", exp, "_", treatment, ".vs.", control, ".pdf"))
expression_dist
dev.off() ## quartz_off_screen
## 2
# Prepare data by getting library sizes and non-zero counts
x <- colSums(rawcounts)
y <- colSums(rawcounts != 0)
df <- data.frame(x,y)
df$ind <- rownames(df)
# Calculate axis limits
x_min <- min(df$x)
x_max <- max(df$x)
y_min <- min(df$y)
y_max <- max(df$y)
# Non-zero plot
nonzero <- ggplot(data = df, aes(x = x, y = y, color = ind, label = ind)) +
geom_point(size = 3) +
geom_text(hjust = -0.1, vjust = -0.1, color = "black", size = 5) +
scale_y_continuous(limits=c(y_min, y_max)) +
scale_x_continuous(
limits = c(x_min-1e3, x_max+9e5),
labels = scales::comma_format(scale = 1e-6, accuracy = 1) # Format labels without decimals
) +
xlab("Millions of reads") + ylab("Number of non-zero genes observed") +
labs(color = "Sample") +
theme(panel.grid = element_blank(),
panel.background = element_blank(),
axis.text = element_text(size = text.size),
axis.title = element_text(size = axis.text.size),
legend.text = element_text(size = text.size),
legend.title = element_text(size = axis.text.size),
panel.border = element_rect(fill = NA),
legend.position = "none",
axis.ticks = element_line(linewidth = 0.2),
legend.key = element_blank()) +
coord_cartesian()
# Show nonzero plot
nonzero# Save plot
cairo_pdf(paste0(Sys.Date(), "_non-zero_", exp, "_", treatment, ".vs.", control, ".pdf"))
nonzero
dev.off() ## quartz_off_screen
## 2
# Load GFF annotation
gff_annot <- load_gff_annotations("../../raw-data/paeruginosa_pa14.gff",
id_col = "gene_id", type = "gene")## Returning a df with 16 columns and 5979 rows.
## seqnames start end width strand source type score phase
## gene1650835 chromosome 483 2027 1545 + PseudoCAP gene NA 0
## gene1650837 chromosome 2056 3159 1104 + PseudoCAP gene NA 0
## gene1650839 chromosome 3169 4278 1110 + PseudoCAP gene NA 0
## gene1650841 chromosome 4275 6695 2421 + PseudoCAP gene NA 0
## gene1650843 chromosome 7018 7791 774 - PseudoCAP gene NA 0
## gene1650845 chromosome 7803 8339 537 - PseudoCAP gene NA 0
## gene_id Name Dbxref Alias name Parent locus
## gene1650835 gene1650835 GeneID:4384099 PA14_00010 dnaA <NA>
## gene1650837 gene1650837 GeneID:4384100 PA14_00020 dnaN <NA>
## gene1650839 gene1650839 GeneID:4384101 PA14_00030 recF <NA>
## gene1650841 gene1650841 GeneID:4384102 PA14_00050 gyrB <NA>
## gene1650843 gene1650843 GeneID:4385186 PA14_00060 PA14_00060 <NA>
## gene1650845 gene1650845 GeneID:4385187 PA14_00070 PA14_00070 <NA>
## SK064 SK065 SK066 SK076 SK077 SK078
## gene1650835 2234 2528 1719 1251 1173 1320
## gene1650837 1266 1754 1368 1120 1029 1066
## gene1650839 722 770 657 317 392 372
## gene1650841 3671 4307 3347 2698 2837 2620
## gene1650843 436 435 313 204 184 174
## gene1650845 257 299 177 141 94 170
## [1] 5979 16
## [1] 5979 6
## Mode FALSE TRUE
## logical 932 5047
## Mode TRUE
## logical 5979
# Create an expression object with counts, metadata, and gene info
expt <- create_expt(count_dataframe = rawcounts,
metadata = metadata,
gene_info = gff_annot)## Reading the sample metadata.
## Did not find the column: sampleid.
## The rownames do not appear numeric, using them.
## The sample definitions comprises: 6 rows(samples) and 8 columns(metadata fields).
## Matched 5979 annotations and counts.
## Bringing together the count matrix and gene information.
## Some annotations were lost in merging, setting them to 'undefined'.
## Saving the expressionset to 'expt.rda'.
## The final expressionset has 5979 features and 6 samples.
# Normalize counts log2 CPM
cpm_normalized_exp <- normalize_expt(expt,
transform = "log2",
norm = "raw",
convert = "cpm",
row_min = "10")## transform_counts: Found 270 values equal to 0, adding 1 to the matrix.
####### Calculate RPKM ########
# Grab gene lengths from GFF (should be in base pairs)
gene_lengths <- gff_annot$width
names(gene_lengths) <- rownames(gff_annot) # set gene names
# Make sure gene_lengths are in the same order as rawcounts rows
gene_lengths <- gene_lengths[rownames(rawcounts)]
# Make DGEList object for edgeR
dge <- DGEList(counts = rawcounts)
# Calculate RPKM (reads per kilobase per million)
rpkm <- rpkm(dge, gene.length = gene_lengths)
# Wrap RPKM in an experiment object too
expt_rpkm <- create_expt(count_dataframe = rpkm,
metadata = metadata,
gene_info = gff_annot)## Reading the sample metadata.
## Did not find the column: sampleid.
## The rownames do not appear numeric, using them.
## The sample definitions comprises: 6 rows(samples) and 8 columns(metadata fields).
## Matched 5979 annotations and counts.
## Bringing together the count matrix and gene information.
## Some annotations were lost in merging, setting them to 'undefined'.
## Saving the expressionset to 'expt.rda'.
## The final expressionset has 5979 features and 6 samples.
# log2 transform the RPKM (again, no normalization beyond log)
rpkm_normalized_exp <- normalize_expt(expt_rpkm,
transform = "log2",
norm = "raw")## transform_counts: Found 270 values equal to 0, adding 1 to the matrix.
# Get the actual expression matrices from the experiment objects
rpkm <- exprs(rpkm_normalized_exp)
cpm <- exprs(cpm_normalized_exp)
# Merge CPM and RPKM matrices with gene annotation info
cpm_annotated <- merge(gff_annot, cpm, by = 0)
rpkm_annotated <- merge(gff_annot, rpkm, by = 0)This generates a small heatmap of the alginate genes across samples to visualize expression.
# Loop through CPM and RPKM to plot both heatmaps
for (type in c("cpm", "rpkm")) {
# Define gene list
alg_operon_genes <- c(
"algA", "algF", "algJ", "algI", "algL", "algX", "algG", "algE",
"algK", "alg44", "alg8", "algD", "mucC", "mucB", "algU", "algW",
"algP", "algQ", "algR", "algZ", "algC", "algB"
)
# Get the right data per type of transformation
data_annot <- get(paste0(type, "_annotated"))
# Fix rownames using gene name (make unique)
rownames(data_annot) <- make.names(data_annot$name, unique = TRUE)
# Drop undesired metadata columns
cols_to_exclude <- c("Row.names","seqnames","start","end","width",
"strand","source","type","score","phase",
"gene_id","Name","Dbxref","Alias","name",
"Parent","locus")
keep <- !names(data_annot) %in% cols_to_exclude
expr_data <- data_annot[, keep]
# Extract matrix for alg genes
alg_genes <- expr_data[alg_operon_genes, , drop = FALSE]
mat <- alg_genes - rowMeans(alg_genes)
# Pull sample annotations
anno <- as.data.frame(metadata[, c("condition", "induced")])
colnames(anno) <- c("Condition", "Induced")
# Color settings for annotation
# ann_colors <- list(
# Condition = c("LB-DMSO" = "#D95F02", "LB-IPTG" = "#66A61E"),
# Induced = c("no" = "lightgrey", "yes" = "gray27")
# )
# Scale data by row (Z-score per gene)
scaled_mat <- t(scale(t(mat)))
# Build heatmap
ht <- Heatmap(
scaled_mat,
name = "Z-score",
top_annotation = HeatmapAnnotation(df = anno#,
#col = ann_colors
),
cluster_rows = TRUE,
cluster_columns = TRUE,
row_names_side = "left",
column_names_side = "bottom",
column_names_rot = 45,
row_names_gp = gpar(fontsize = 10),
column_names_gp = gpar(fontsize = 10),
show_row_dend = TRUE,
show_column_dend = TRUE
)
# Draw and save
draw(ht)
cairo_pdf(paste0(Sys.Date(), "_alg_heatmap_", type, "_", exp, "_", treatment, ".vs.", control, ".pdf"))
draw(ht)
dev.off()
}# Create DESeqDataSet
dds <- DESeqDataSetFromMatrix(
countData = rawcounts,
colData = metadata,
design = design_formula
)## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
## [1] "Number of genes before filtering:"
## [1] 5979
# Filter: keep genes with counts >=10 in at least 2 samples
print("Filtering genes with counts >=10 in at least 2 samples")## [1] "Filtering genes with counts >=10 in at least 2 samples"
keep <- rowSums(counts(dds) >= 10) >= 2
dds_filtered <- dds[keep, ]
# Check number of genes after filtering
print("Number of genes after filtering:")## [1] "Number of genes after filtering:"
## [1] 5723
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## [1] "Intercept" "exp_5_treatment_vs_control"
## [3] "batch_b2_vs_b1" "batch_b3_vs_b1"
res <- results(dds_final, contrast = c("exp_5", "treatment", "control"))
resdf <- as.data.frame(res)
resdf$gene_id <- rownames(resdf)
cpm <- as.data.frame(cpm)
rpkm <- as.data.frame(rpkm)
# append type of value to the dataframe
colnames(rpkm) <- paste0(colnames(rpkm), "_rpkm")
colnames(cpm) <- paste0(colnames(cpm), "_cpm")
rawcounts2 <- rawcounts
colnames(rawcounts2) <- paste0(colnames(rawcounts2), "_rawcounts")
# ensure column gene_id is there
rpkm$gene_id <- rownames(rpkm)
cpm$gene_id <- rownames(cpm)
rawcounts2$gene_id <- rownames(rawcounts2)
gff_annot <- gff_annot[rownames(resdf),]
dim(resdf)## [1] 5723 7
# merge with siggenes
resdf <- full_join(resdf, gff_annot, by = "gene_id", copy = FALSE)
head(resdf)## baseMean log2FoldChange lfcSE stat pvalue padj
## 1 1626.6686 -0.35245305 0.2052851 -1.71689571 0.085998229 0.21174657
## 2 1236.7503 -0.01285259 0.2148096 -0.05983246 0.952289075 0.97714378
## 3 512.7254 -0.57072684 0.2266342 -2.51827360 0.011793167 0.05086285
## 4 3163.4949 -0.04378508 0.2007550 -0.21810205 0.827349597 0.91257387
## 5 274.2181 -0.63849199 0.2390568 -2.67088021 0.007565264 0.03676861
## 6 179.5968 -0.43813691 0.2921709 -1.49959118 0.133720333 0.28575757
## gene_id seqnames start end width strand source type score phase
## 1 gene1650835 chromosome 483 2027 1545 + PseudoCAP gene NA 0
## 2 gene1650837 chromosome 2056 3159 1104 + PseudoCAP gene NA 0
## 3 gene1650839 chromosome 3169 4278 1110 + PseudoCAP gene NA 0
## 4 gene1650841 chromosome 4275 6695 2421 + PseudoCAP gene NA 0
## 5 gene1650843 chromosome 7018 7791 774 - PseudoCAP gene NA 0
## 6 gene1650845 chromosome 7803 8339 537 - PseudoCAP gene NA 0
## Name Dbxref Alias name Parent locus
## 1 GeneID:4384099 PA14_00010 dnaA <NA>
## 2 GeneID:4384100 PA14_00020 dnaN <NA>
## 3 GeneID:4384101 PA14_00030 recF <NA>
## 4 GeneID:4384102 PA14_00050 gyrB <NA>
## 5 GeneID:4385186 PA14_00060 PA14_00060 <NA>
## 6 GeneID:4385187 PA14_00070 PA14_00070 <NA>
## locusId accession GI scaffoldId start stop strand Alias
## 1 2194572 YP_788156.1 116053721 4582 483 2027 + PA14_00010
## 2 2194573 YP_788157.1 116053722 4582 2056 3159 + PA14_00020
## 3 2194574 YP_788158.1 116053723 4582 3169 4278 + PA14_00030
## 4 2194575 YP_788159.1 116053724 4582 4275 6695 + PA14_00050
## 5 2194576 YP_788160.1 116053725 4582 7791 7018 - PA14_00060
## 6 2194577 YP_788161.1 116053726 4582 8339 7803 - PA14_00070
## name desc COG
## 1 dnaA chromosomal replication initiator protein DnaA (NCBI) COG593
## 2 dnaN DNA polymerase III, beta chain (NCBI) COG592
## 3 recF DNA replication and repair protein RecF (NCBI) COG1195
## 4 gyrB DNA gyrase subunit B (NCBI) COG187
## 5 PA14_00060 putative acyltransferase (NCBI) COG204
## 6 PA14_00070 putative histidinol-phosphatase (NCBI) COG241
## COGFun
## 1 L
## 2 L
## 3 L
## 4 L
## 5 I
## 6 E
## COGDesc
## 1 ATPase involved in DNA replication initiation
## 2 DNA polymerase sliding clamp subunit (PCNA homolog)
## 3 Recombinational DNA repair ATPase (RecF pathway)
## 4 Type IIA topoisomerase (DNA gyrase/topo II, topoisomerase IV), B subunit
## 5 1-acyl-sn-glycerol-3-phosphate acyltransferase
## 6 Histidinol phosphatase and related phosphatases
## TIGRFam
## 1 TIGR00362 chromosomal replication initiator protein DnaA [dnaA]
## 2 TIGR00663 DNA polymerase III, beta subunit [dnaN]
## 3 TIGR00611 DNA replication and repair protein RecF [recF]
## 4 TIGR01059 DNA gyrase, B subunit [gyrB]
## 5
## 6 TIGR01656 histidinol-phosphate phosphatase domain,TIGR01662 HAD hydrolase, family IIIA
## TIGRRoles
## 1 DNA metabolism:DNA replication, recombination, and repair
## 2 DNA metabolism:DNA replication, recombination, and repair
## 3 DNA metabolism:DNA replication, recombination, and repair
## 4 DNA metabolism:DNA replication, recombination, and repair
## 5
## 6 Unknown function:Enzymes of unknown specificity
## GO
## 1 GO:0006270,GO:0006275,GO:0003688,GO:0017111,GO:0005524
## 2 GO:0006260,GO:0003677,GO:0003893,GO:0008408,GO:0016449,GO:0019984,GO:0003889,GO:0003894,GO:0015999,GO:0016450,GO:0003890,GO:0003895,GO:0016000,GO:0016451,GO:0003891,GO:0016448,GO:0016452
## 3 GO:0006281,GO:0005694,GO:0005524,GO:0017111,GO:0003697
## 4 GO:0006304,GO:0006265,GO:0005694,GO:0003918,GO:0005524
## 5 GO:0008152,GO:0003841
## 6 GO:0000105,GO:0004401
## EC ECDesc
## 1
## 2 2.7.7.7 DNA-directed DNA polymerase.
## 3
## 4 5.99.1.3 DNA topoisomerase (ATP-hydrolyzing).
## 5 2.3.1.51 1-acylglycerol-3-phosphate O-acyltransferase.
## 6 3.1.3.-
col_to_remove <- intersect(colnames(pa_go),colnames(resdf))
col_to_remove <- setdiff(col_to_remove, "Alias") # returns all except the second vector
pa_go <- pa_go %>% select(-all_of(col_to_remove))
resdf <- full_join(resdf, rpkm, by = "gene_id", copy = FALSE)
resdf <- full_join(resdf, cpm, by = "gene_id", copy = FALSE)
resdf <- full_join(resdf, rawcounts2, by = "gene_id", copy = FALSE)
#resdf <- full_join(resdf, pa_go, by = "Alias", copy = FALSE)
resdf <- merge(resdf, pa_go, by = "Alias", copy = FALSE)
head(resdf)## Alias baseMean log2FoldChange lfcSE stat pvalue
## 1 PA14_00010 1626.6686 -0.35245305 0.2052851 -1.71689571 0.085998229
## 2 PA14_00020 1236.7503 -0.01285259 0.2148096 -0.05983246 0.952289075
## 3 PA14_00030 512.7254 -0.57072684 0.2266342 -2.51827360 0.011793167
## 4 PA14_00050 3163.4949 -0.04378508 0.2007550 -0.21810205 0.827349597
## 5 PA14_00060 274.2181 -0.63849199 0.2390568 -2.67088021 0.007565264
## 6 PA14_00070 179.5968 -0.43813691 0.2921709 -1.49959118 0.133720333
## padj gene_id seqnames start end width strand source type
## 1 0.21174657 gene1650835 chromosome 483 2027 1545 + PseudoCAP gene
## 2 0.97714378 gene1650837 chromosome 2056 3159 1104 + PseudoCAP gene
## 3 0.05086285 gene1650839 chromosome 3169 4278 1110 + PseudoCAP gene
## 4 0.91257387 gene1650841 chromosome 4275 6695 2421 + PseudoCAP gene
## 5 0.03676861 gene1650843 chromosome 7018 7791 774 - PseudoCAP gene
## 6 0.28575757 gene1650845 chromosome 7803 8339 537 - PseudoCAP gene
## score phase Name Dbxref name Parent locus SK064_rpkm SK065_rpkm
## 1 NA 0 GeneID:4384099 dnaA <NA> 8.176179 8.165287
## 2 NA 0 GeneID:4384100 dnaN <NA> 7.842995 8.122956
## 3 NA 0 GeneID:4384101 recF <NA> 7.029755 6.934080
## 4 NA 0 GeneID:4384102 gyrB <NA> 8.244496 8.285574
## 5 NA 0 GeneID:4385186 PA14_00060 <NA> 6.823946 6.633153
## 6 NA 0 GeneID:4385187 PA14_00070 <NA> 6.591050 6.619830
## SK066_rpkm SK076_rpkm SK077_rpkm SK078_rpkm SK064_cpm SK065_cpm SK066_cpm
## 1 7.893669 7.376584 7.169888 7.503045 8.802026 8.791120 8.519135
## 2 8.048415 7.700114 7.463936 7.678665 7.985143 8.265208 8.190642
## 3 6.988440 5.888950 6.077132 6.165199 7.179220 7.083470 7.137873
## 4 8.205787 7.835021 7.792523 7.842403 9.517304 9.558461 9.478520
## 5 6.444121 5.775223 5.516477 5.598970 6.458065 6.267796 6.079357
## 6 6.152874 5.769860 5.086157 6.084175 5.706889 5.735417 5.273245
## SK076_cpm SK077_cpm SK078_cpm SK064_rawcounts SK065_rawcounts SK066_rawcounts
## 1 8.001125 7.793956 8.127844 2234 2528 1719
## 2 7.842201 7.605906 7.820742 1266 1754 1368
## 3 6.037095 6.225573 6.313765 722 770 657
## 4 9.106910 9.064301 9.114312 3671 4307 3347
## 5 5.413300 5.156056 5.238041 436 435 313
## 6 4.895473 4.225312 5.205387 257 299 177
## SK076_rawcounts SK077_rawcounts SK078_rawcounts locusId accession GI
## 1 1251 1173 1320 2194572 YP_788156.1 116053721
## 2 1120 1029 1066 2194573 YP_788157.1 116053722
## 3 317 392 372 2194574 YP_788158.1 116053723
## 4 2698 2837 2620 2194575 YP_788159.1 116053724
## 5 204 184 174 2194576 YP_788160.1 116053725
## 6 141 94 170 2194577 YP_788161.1 116053726
## scaffoldId stop desc COG
## 1 4582 2027 chromosomal replication initiator protein DnaA (NCBI) COG593
## 2 4582 3159 DNA polymerase III, beta chain (NCBI) COG592
## 3 4582 4278 DNA replication and repair protein RecF (NCBI) COG1195
## 4 4582 6695 DNA gyrase subunit B (NCBI) COG187
## 5 4582 7018 putative acyltransferase (NCBI) COG204
## 6 4582 7803 putative histidinol-phosphatase (NCBI) COG241
## COGFun
## 1 L
## 2 L
## 3 L
## 4 L
## 5 I
## 6 E
## COGDesc
## 1 ATPase involved in DNA replication initiation
## 2 DNA polymerase sliding clamp subunit (PCNA homolog)
## 3 Recombinational DNA repair ATPase (RecF pathway)
## 4 Type IIA topoisomerase (DNA gyrase/topo II, topoisomerase IV), B subunit
## 5 1-acyl-sn-glycerol-3-phosphate acyltransferase
## 6 Histidinol phosphatase and related phosphatases
## TIGRFam
## 1 TIGR00362 chromosomal replication initiator protein DnaA [dnaA]
## 2 TIGR00663 DNA polymerase III, beta subunit [dnaN]
## 3 TIGR00611 DNA replication and repair protein RecF [recF]
## 4 TIGR01059 DNA gyrase, B subunit [gyrB]
## 5
## 6 TIGR01656 histidinol-phosphate phosphatase domain,TIGR01662 HAD hydrolase, family IIIA
## TIGRRoles
## 1 DNA metabolism:DNA replication, recombination, and repair
## 2 DNA metabolism:DNA replication, recombination, and repair
## 3 DNA metabolism:DNA replication, recombination, and repair
## 4 DNA metabolism:DNA replication, recombination, and repair
## 5
## 6 Unknown function:Enzymes of unknown specificity
## GO
## 1 GO:0006270,GO:0006275,GO:0003688,GO:0017111,GO:0005524
## 2 GO:0006260,GO:0003677,GO:0003893,GO:0008408,GO:0016449,GO:0019984,GO:0003889,GO:0003894,GO:0015999,GO:0016450,GO:0003890,GO:0003895,GO:0016000,GO:0016451,GO:0003891,GO:0016448,GO:0016452
## 3 GO:0006281,GO:0005694,GO:0005524,GO:0017111,GO:0003697
## 4 GO:0006304,GO:0006265,GO:0005694,GO:0003918,GO:0005524
## 5 GO:0008152,GO:0003841
## 6 GO:0000105,GO:0004401
## EC ECDesc
## 1
## 2 2.7.7.7 DNA-directed DNA polymerase.
## 3
## 4 5.99.1.3 DNA topoisomerase (ATP-hydrolyzing).
## 5 2.3.1.51 1-acylglycerol-3-phosphate O-acyltransferase.
## 6 3.1.3.-
## [1] 5716 54
## Alias baseMean log2FoldChange lfcSE stat pvalue
## 1 PA14_00010 1626.6686 -0.35245305 0.2052851 -1.71689571 0.085998229
## 2 PA14_00020 1236.7503 -0.01285259 0.2148096 -0.05983246 0.952289075
## 3 PA14_00030 512.7254 -0.57072684 0.2266342 -2.51827360 0.011793167
## 4 PA14_00050 3163.4949 -0.04378508 0.2007550 -0.21810205 0.827349597
## 5 PA14_00060 274.2181 -0.63849199 0.2390568 -2.67088021 0.007565264
## 6 PA14_00070 179.5968 -0.43813691 0.2921709 -1.49959118 0.133720333
## padj gene_id seqnames start end width strand source type
## 1 0.21174657 gene1650835 chromosome 483 2027 1545 + PseudoCAP gene
## 2 0.97714378 gene1650837 chromosome 2056 3159 1104 + PseudoCAP gene
## 3 0.05086285 gene1650839 chromosome 3169 4278 1110 + PseudoCAP gene
## 4 0.91257387 gene1650841 chromosome 4275 6695 2421 + PseudoCAP gene
## 5 0.03676861 gene1650843 chromosome 7018 7791 774 - PseudoCAP gene
## 6 0.28575757 gene1650845 chromosome 7803 8339 537 - PseudoCAP gene
## score phase Name Dbxref name Parent locus SK064_rpkm SK065_rpkm
## 1 NA 0 GeneID:4384099 dnaA <NA> 8.176179 8.165287
## 2 NA 0 GeneID:4384100 dnaN <NA> 7.842995 8.122956
## 3 NA 0 GeneID:4384101 recF <NA> 7.029755 6.934080
## 4 NA 0 GeneID:4384102 gyrB <NA> 8.244496 8.285574
## 5 NA 0 GeneID:4385186 PA14_00060 <NA> 6.823946 6.633153
## 6 NA 0 GeneID:4385187 PA14_00070 <NA> 6.591050 6.619830
## SK066_rpkm SK076_rpkm SK077_rpkm SK078_rpkm SK064_cpm SK065_cpm SK066_cpm
## 1 7.893669 7.376584 7.169888 7.503045 8.802026 8.791120 8.519135
## 2 8.048415 7.700114 7.463936 7.678665 7.985143 8.265208 8.190642
## 3 6.988440 5.888950 6.077132 6.165199 7.179220 7.083470 7.137873
## 4 8.205787 7.835021 7.792523 7.842403 9.517304 9.558461 9.478520
## 5 6.444121 5.775223 5.516477 5.598970 6.458065 6.267796 6.079357
## 6 6.152874 5.769860 5.086157 6.084175 5.706889 5.735417 5.273245
## SK076_cpm SK077_cpm SK078_cpm SK064_rawcounts SK065_rawcounts SK066_rawcounts
## 1 8.001125 7.793956 8.127844 2234 2528 1719
## 2 7.842201 7.605906 7.820742 1266 1754 1368
## 3 6.037095 6.225573 6.313765 722 770 657
## 4 9.106910 9.064301 9.114312 3671 4307 3347
## 5 5.413300 5.156056 5.238041 436 435 313
## 6 4.895473 4.225312 5.205387 257 299 177
## SK076_rawcounts SK077_rawcounts SK078_rawcounts locusId accession GI
## 1 1251 1173 1320 2194572 YP_788156.1 116053721
## 2 1120 1029 1066 2194573 YP_788157.1 116053722
## 3 317 392 372 2194574 YP_788158.1 116053723
## 4 2698 2837 2620 2194575 YP_788159.1 116053724
## 5 204 184 174 2194576 YP_788160.1 116053725
## 6 141 94 170 2194577 YP_788161.1 116053726
## scaffoldId stop desc COG
## 1 4582 2027 chromosomal replication initiator protein DnaA (NCBI) COG593
## 2 4582 3159 DNA polymerase III, beta chain (NCBI) COG592
## 3 4582 4278 DNA replication and repair protein RecF (NCBI) COG1195
## 4 4582 6695 DNA gyrase subunit B (NCBI) COG187
## 5 4582 7018 putative acyltransferase (NCBI) COG204
## 6 4582 7803 putative histidinol-phosphatase (NCBI) COG241
## COGFun
## 1 L
## 2 L
## 3 L
## 4 L
## 5 I
## 6 E
## COGDesc
## 1 ATPase involved in DNA replication initiation
## 2 DNA polymerase sliding clamp subunit (PCNA homolog)
## 3 Recombinational DNA repair ATPase (RecF pathway)
## 4 Type IIA topoisomerase (DNA gyrase/topo II, topoisomerase IV), B subunit
## 5 1-acyl-sn-glycerol-3-phosphate acyltransferase
## 6 Histidinol phosphatase and related phosphatases
## TIGRFam
## 1 TIGR00362 chromosomal replication initiator protein DnaA [dnaA]
## 2 TIGR00663 DNA polymerase III, beta subunit [dnaN]
## 3 TIGR00611 DNA replication and repair protein RecF [recF]
## 4 TIGR01059 DNA gyrase, B subunit [gyrB]
## 5
## 6 TIGR01656 histidinol-phosphate phosphatase domain,TIGR01662 HAD hydrolase, family IIIA
## TIGRRoles
## 1 DNA metabolism:DNA replication, recombination, and repair
## 2 DNA metabolism:DNA replication, recombination, and repair
## 3 DNA metabolism:DNA replication, recombination, and repair
## 4 DNA metabolism:DNA replication, recombination, and repair
## 5
## 6 Unknown function:Enzymes of unknown specificity
## GO
## 1 GO:0006270,GO:0006275,GO:0003688,GO:0017111,GO:0005524
## 2 GO:0006260,GO:0003677,GO:0003893,GO:0008408,GO:0016449,GO:0019984,GO:0003889,GO:0003894,GO:0015999,GO:0016450,GO:0003890,GO:0003895,GO:0016000,GO:0016451,GO:0003891,GO:0016448,GO:0016452
## 3 GO:0006281,GO:0005694,GO:0005524,GO:0017111,GO:0003697
## 4 GO:0006304,GO:0006265,GO:0005694,GO:0003918,GO:0005524
## 5 GO:0008152,GO:0003841
## 6 GO:0000105,GO:0004401
## EC ECDesc
## 1
## 2 2.7.7.7 DNA-directed DNA polymerase.
## 3
## 4 5.99.1.3 DNA topoisomerase (ATP-hydrolyzing).
## 5 2.3.1.51 1-acylglycerol-3-phosphate O-acyltransferase.
## 6 3.1.3.-
## Mode FALSE
## logical 5716
## Mode FALSE TRUE
## logical 5272 444
sig_genes <- resdf[abs(resdf$log2FoldChange) >= logfc_cutoff & resdf$padj < pval_cutoff, ] %>% drop_na(log2FoldChange)
sig_genes$direction <- ifelse(sig_genes$log2FoldChange > 0, "up", "down")
summary(is.na(sig_genes$log2FoldChange))## Mode FALSE
## logical 29
## Mode FALSE
## logical 29
## [1] 29 55
## Alias baseMean log2FoldChange lfcSE stat pvalue
## 152 PA14_01860 60.46235 -2.236306 0.4947510 -4.520064 6.182105e-06
## 360 PA14_04650 14554.03916 2.029396 0.2229824 9.101148 8.938164e-20
## 570 PA14_07355 250.08453 3.270462 0.3580363 9.134443 6.574423e-20
## 1518 PA14_19560 102.09570 2.100120 0.3280252 6.402314 1.530397e-10
## 1520 PA14_19580 36.38396 2.024175 0.4675001 4.329786 1.492544e-05
## 1689 PA14_21680 3450.24358 2.057359 0.2441252 8.427474 3.532072e-17
## padj gene_id seqnames start end width strand source
## 152 1.059589e-04 gene1651139 chromosome 168754 169749 996 + PseudoCAP
## 360 7.864095e-17 gene1651577 chromosome 413653 414192 540 - PseudoCAP
## 570 6.941276e-17 gene1651999 chromosome 631670 632008 339 + PseudoCAP
## 1518 1.077196e-08 gene1653961 chromosome 1690298 1691446 1149 + PseudoCAP
## 1520 2.251183e-04 gene1653965 chromosome 1692252 1693076 825 + PseudoCAP
## 1689 1.695074e-14 gene1654319 chromosome 1878700 1879299 600 + PseudoCAP
## type score phase Name Dbxref name Parent locus SK064_rpkm
## 152 gene NA 0 GeneID:4383550 PA14_01860 <NA> 4.279016
## 360 gene NA 0 GeneID:4384279 pfpI <NA> 11.190724
## 570 gene NA 0 GeneID:4384180 PA14_07355 <NA> 5.007502
## 1518 gene NA 0 GeneID:4385230 ssuD <NA> 3.166867
## 1520 gene NA 0 GeneID:4385228 ssuB <NA> 1.963340
## 1689 gene NA 0 GeneID:4381273 PA14_21680 <NA> 8.728896
## SK065_rpkm SK066_rpkm SK076_rpkm SK077_rpkm SK078_rpkm SK064_cpm SK065_cpm
## 152 4.696355 4.589707 2.993446 0.8126278 2.655353 4.273532 4.690796
## 360 11.031055 11.484658 13.037210 12.7561007 12.858014 10.302281 10.142674
## 570 4.502150 5.099123 8.682287 7.5648775 6.845299 3.531763 3.060589
## 1518 2.695528 3.435322 4.970845 4.4145574 4.710203 3.346263 2.866733
## 1520 2.063104 2.846336 3.792644 3.5521295 3.838849 1.762222 1.857005
## 1689 8.835231 9.381954 10.814799 10.5290682 10.624484 7.994195 8.100370
## SK066_cpm SK076_cpm SK077_cpm SK078_cpm SK064_rawcounts SK065_rawcounts
## 152 4.584166 2.988391 0.8101404 2.650490 92 142
## 360 10.596118 12.148388 11.8673096 11.969211 6329 6460
## 570 3.618301 7.128476 6.0190153 5.308916 53 42
## 1518 3.618301 5.165245 4.6061369 4.903417 46 36
## 1520 2.610739 3.537025 3.3004517 3.582546 12 15
## 1689 8.646429 10.078367 9.7927534 9.888127 1274 1564
## SK066_rawcounts SK076_rawcounts SK077_rawcounts SK078_rawcounts locusId
## 152 108 34 4 25 2194724
## 360 7268 22249 19831 18984 2194942
## 570 53 681 339 183 2195152
## 1518 53 171 124 137 2196119
## 1520 24 52 47 52 2196121
## 1689 1878 5295 4704 4483 2196297
## accession GI scaffoldId stop
## 152 YP_788308.1 116053871 4582 169749
## 360 YP_788526.1 116054083 4582 413653
## 570 YP_788736.1 116054407 4582 632008
## 1518 YP_789703.1 116051464 4582 1691446
## 1520 YP_789705.1 116051462 4582 1693076
## 1689 YP_789881.1 116051288 4582 1879299
## desc COG
## 152 putative transmembrane sensor (NCBI) COG3712
## 360 protease PfpI (NCBI) COG693
## 570 hypothetical protein (NCBI) COG599
## 1518 putative sulfonate monooxygenase (NCBI) COG2141
## 1520 putative aliphatic sulfonate transport ATP-binding protein (NCBI) COG1116
## 1689 hypothetical protein (NCBI)
## COGFun
## 152 PT
## 360 R
## 570 S
## 1518 C
## 1520 P
## 1689
## COGDesc
## 152 Fe2+-dicitrate sensor, membrane component
## 360 Putative intracellular protease/amidase
## 570 Uncharacterized homolog of gamma-carboxymuconolactone decarboxylase subunit
## 1518 Coenzyme F420-dependent N5,N10-methylene tetrahydromethanopterin reductase and related flavin-dependent oxidoreductases
## 1520 ABC-type nitrate/sulfonate/bicarbonate transport system, ATPase component
## 1689
## TIGRFam
## 152
## 360 TIGR01382 intracellular protease, PfpI family
## 570 TIGR00778 alkylhydroperoxidase AhpD family core domain
## 1518 TIGR03565 alkanesulfonate monooxygenase, FMNH(2)-dependent [ssuD]
## 1520
## 1689
## TIGRRoles
## 152
## 360 Protein fate:Degradation of proteins, peptides, and glycopeptides
## 570 Unknown function:General
## 1518 Central intermediary metabolism:Sulfur metabolism
## 1520
## 1689
## GO EC ECDesc direction
## 152 down
## 360 GO:0016798 3.2.-.- up
## 570 up
## 1518 GO:0008726 1.14.14.5 Alkanesulfonate monooxygenase. up
## 1520 GO:0005524,GO:0016887 up
## 1689 up
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
# Variance stabilizing transformations
vsd <- vst(vst, blind = TRUE)
# PCA plot
pcaData <- plotPCA(vsd, intgroup = c("exp_5", "batch", "condition"), returnData = TRUE)## using ntop=500 top features by variance
percentVar <- round(100 * attr(pcaData, "percentVar"))
# Calculate the maximum absolute value among PC1 and PC2
max_abs_val <- max(c(abs(min(pcaData$PC1)), abs(min(pcaData$PC2)), abs(max(pcaData$PC1)), abs(max(pcaData$PC2))))
axis_size <- 10
theme_set(theme_minimal(base_family = "Arial"))
pca <- ggplot(pcaData, aes(PC1, PC2, color = condition, shape = condition)) +
geom_point(size = 3) +
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
scale_colour_brewer(palette = "Dark2") +
scale_shape_manual(values = c(1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16)) +
coord_fixed(ratio = 1, xlim = c(-max_abs_val, max_abs_val), ylim = c(-max_abs_val, max_abs_val)) +
labs(color = "condition", shape = "condition") +
guides(shape = guide_legend(title = "Condition"), # Adjust size of shape legend
color = guide_legend(title = "Condition")) +
theme(
panel.grid = element_blank(),
panel.background = element_blank(),
axis.text = element_text(size = axis_size),
axis.ticks = element_line(linewidth = 0.2),
axis.line = element_line(linewidth = 0, color = "black"), # Set the axis line linewidth
panel.border = element_rect(fill = NA, linewidth = 0.5), # Set the panel border linewidth
legend.key = element_rect(fill = "white", color = NA) # Set legend key background to white and remove border
)
print(pca)cairo_pdf(paste0(Sys.Date(),"_pca_",exp, "_", treatment,".vs.", control,".pdf", sep = ""))
pca
dev.off()## quartz_off_screen
## 2
sequence <- c("treatment")
for (i in sequence) {
set.seed(0)
res <- results(dds_final, contrast = c("exp_5", i, "control"))
gp2 <- wes_palettes$Zissou1
key <- read.csv("../../raw-data/pa14_gene_key.csv", header = TRUE)
rownames(key) <- key$gene_id
res <- merge(key, as.data.frame(res), by=0)
res <- as.data.frame(res)
res$Alias <- rownames(res)
# Define gene list
alg_genes <- c(
"algA", "algF", "algJ", "algI", "algL", "algX", "algG", "algE",
"algK", "alg44", "alg8", "algD", "mucC", "mucB", "algU", "algW",
"algP", "algQ", "algR", "algZ", "algC", "algB"
)
logfc_cutoff
keyvals <- ifelse(
res$log2FoldChange <= -logfc_cutoff & res$padj <= pval_cutoff, gp2[1],
ifelse(res$log2FoldChange >= logfc_cutoff & res$padj <= pval_cutoff, gp2[5],
'lightgrey'))
keyvals[is.na(keyvals)] <- 'lightgrey'
names(keyvals)[keyvals == gp2[5]] <- 'Up-regulated'
names(keyvals)[keyvals == 'lightgrey'] <- 'NS'
names(keyvals)[keyvals == gp2[1]] <- 'Down-regulated'
volcano <- EnhancedVolcano(res,
lab = res[["name"]],
selectLab = alg_genes,
x = 'log2FoldChange',
y = 'padj',
xlab = bquote(~Log[2]~ 'fold change'),
pCutoff = pval_cutoff,
FCcutoff = logfc_cutoff,
cutoffLineType = "dashed",
pointSize = 1,
labSize = 3.0,
axisLabSize = 14,
labCol = 'black',
labFace = 'bold',
#col = c('grey', gp2[2], gp2[3], gp2[1]),
colCustom = keyvals,
colAlpha = 5/5,
legendPosition = 'bottom',
legendLabSize = 12,
legendIconSize = 4.0,
drawConnectors = TRUE,
arrowheads = FALSE,
widthConnectors = 0.2,
gridlines.major = FALSE,
borderWidth = 0.3,
gridlines.minor = FALSE,
title = paste0(treatment," vs ", control, sep = ""),
subtitle = "",
#caption = paste0("total = ", nrow(toptable), " variables"),
caption = paste0("Volcano plot with log2 FC of ", logfc_cutoff," and P value of ",pval_cutoff, "."),
max.overlaps = 30,
#ylim = c(0,10),
colConnectors = 'black') +
theme(axis.ticks=element_line(linewidth =0.2)) # +
# scale_x_continuous(breaks = c(-5, -2.5, 0, 2.5, 5), limits = c(-6,6))
print(volcano)
cairo_pdf(paste0(Sys.Date(),"_volcano_", exp, "_", treatment,".vs.", control,".pdf", sep = ""))
print(volcano)
dev.off()
}sequence <- c("treatment")
for (i in sequence) {
set.seed(0)
res <- results(dds_final, contrast = c("exp_5", i, "control"))
gp2 <- wes_palettes$Zissou1
res$gene_id <- rownames(res)
res <- as.data.frame(res)
head(res)
res <- full_join(key, as.data.frame(res), by="gene_id")
rownames(res) <- res$Alias
any(is.na(res)) # Check if there are any NAs
res <- res %>% filter(!is.na(baseMean))
ma <- ggmaplot(res,
fdr = pval_cutoff,
fc = 2^(logfc_cutoff),
size = 1,
palette = c((gp2[5]), (gp2[1])),
select.top.method = "fc",
# ylim = c(-15,15),
legend = "bottom",
title = paste0(treatment," vs ", control, sep = ""),
subtitle = paste0("Volcano plot with log2 FC of ", logfc_cutoff," and P value of ",pval_cutoff, "."),
font.label = c( 12),
label.rectangle = FALSE,
genenames = res$name,
font.legend = c(12),
font.main = c("bold", 14),
alpha = 1.5,
ggtheme = ggplot2::theme(panel.grid =element_blank(),
axis.text = element_text(size = 12),
axis.title = element_text(size = 14),
panel.background = element_blank(),
panel.border = element_rect(fill = NA)))
print(ma)
cairo_pdf(paste0(Sys.Date(),"_ma_", exp, "_", treatment,".vs.", control,".pdf"))
print(ma)
dev.off()
}