This is a test version of Biostars. For the public version, visit https://www.biostars.org.
DSEQ2 analysis

Hi I did gene expression analysis with dseq2

I'm trying to understand if positive logfoldchange mean mot expression in the mutant or in the control I'm a bit confused this is my code

#!/usr/bin/env Rscript

# Set library path and load required libraries
.libPaths("/sci/labs/maayan.salton/adi8897/R_libs")
library(DESeq2)

# Define the full paths to your count files for the first analysis (Israeli_WBP4_mutant)
countFiles <- c(
  "SRR9600556_sorted_counts.txt",
  "SRR9600557_sorted_counts.txt",
  "SRR9600558_sorted_counts.txt",
  "SRR9600559_sorted_counts.txt",
  "output_1024-3.counts"  # Israeli_WBP4_mutant
)

# Sample names extracted from file names
sampleNames <- gsub("_sorted_counts.txt|_counts", "", countFiles)

# Define conditions for each sample (first four are controls, last one is the mutant)
conditions <- c("control", "control", "control", "control", "mutant")

# Create a sample information dataframe
sampleInfo <- data.frame(row.names = sampleNames, condition = factor(conditions, levels = c("control", "mutant")))

# Initialize an empty list to store count data from each file
countDataList <- list()

# Loop through count files to read and store data
for (i in seq_along(countFiles)) {
  countFilePath <- file.path("/sci/labs/maayan.salton/adi8897/Alternative_Splicing_project/WBP4/gene_expression", countFiles[i])
  counts <- read.table(countFilePath, header = FALSE, col.names = c("gene", sampleNames[i]))
  countDataList[[i]] <- counts
}

# Merge all count data into a single data frame by gene
mergedCounts <- Reduce(function(x, y) merge(x, y, by = "gene", all = TRUE), countDataList)

# Convert merged data frame to a matrix for DESeq2, with genes as rows and samples as columns
countMatrix <- as.matrix(mergedCounts[-1])
rownames(countMatrix) <- mergedCounts$gene

# Make sure column names of the count matrix match the row names of the sample information dataframe
colnames(countMatrix) <- sampleNames  # Ensuring the names align exactly

# Create DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = countMatrix,
                              colData = sampleInfo,
                              design = ~ condition)

# Run the DESeq analysis
dds <- DESeq(dds)

# Getting results
results <- results(dds)

# Ordering the results by p-value
orderedResults <- results[order(results$pvalue), ]

# Write the results to a CSV file
write.csv(as.data.frame(orderedResults), "DESeq2_results_Israeli_WBP4_mutant_vs_control.csv")

# Optional: Generate an MA plot
plotMA(results, main="MA Plot for Israeli WBP4 Mutant vs Control", ylim=c(-2,2))

thanks for your help

deseq2 logfoldchange

1 answer

The tool is called DESeq2. Anyway, if your factor is factor(conditions, levels = c("control", "mutant")) then the first level is the reference which in default results() means that positive logFC indicate overexpression in mutant.

Log in to answer this question.