Entering edit mode
7 months ago
adi.gershon1
▴
10
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