DSEQ2 analysis
1
0
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

DESEQ2 logfoldchange • 337 views
ADD COMMENT
1
Entering edit mode
7 months ago
ATpoint 86k

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.

ADD COMMENT

Login before adding your answer.

Traffic: 2033 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6