Is it really CONSIDERED a batch effect when I extracted information from FASTQ file, reading first lines of it, and got the run_number and flowcell_ID from those lines? Or I am unintentionally reading too much by extracting such information for my RNA-sequence dataset which actually is not a batch effect?
•
written
by
Umair •
1 answer
Hey,
It is not 'over-interpreting' - the information that you have extracted can indeed be used to identify potential batches. In RNA-seq, the sequencing run and flowcell are well known sources of technical variation / batch effects, and, for this reason, are sometimes explicitly included in the statistical model. The flowcell ID, in particular, can be important.
To check if these are actually driving a batch effect in your data, I would advise to generate a PCA bi-plot (or heatmap) of your normalised counts and colour the samples by flowcell / run. If you see a clear separation, then, yes, there is a batch effect. In that case, you can use this information as a covariate in your model, e.g., in DESeq2's design formula.
Kevin
I would not consider the flow cell meaningfully and clearly impacting things based on your non-batch corrected plots. I'd leave it out of your model design personally. Note that checking additional PCs can also be helpful (PC3/4), but it can quickly become a ghost hunt. True batch effects are typically pretty obvious.
Looking at your PCA bi-plot, which was generated by my Bioc package, PCAtools, I would remove the sample at the top-right, then re-do everything, and then re-assess.
The thing with low sample sizes like n=3 is that it is hard to tell what really is an outlier. Removing points reduces power even more (which at n=3 is already low). Overall, there is no strong group separation, suggesting modest numbers of DEGs. What you can do is to use the limma framework for DE testing and either use arrayWeights() together with limma-trend, or either voomWithQualityWeights() or voomLmFit() (the latter with sample.weights = TRUE, see user guide for details) to automatically down-weight samples that show outlier behaviour for many genes according to the design. That is a relatively simple workaround to avoid manually removing samples while dampening the effect of "data-driven-defined" outliers. While this is not necessarily "better" than removing samples, it at least is automated and reproducible, while PCA interpretation is not, unless you put code to identify outliers by some cutoffs such as SD or MAD. See number of DEGs after using these weighting strategies, check MA- and volcano plots to get an idea whether there is any evidence for the treatment effect. Display DEGs in a heatmap, see whether some samples obviously do not match the group pattern. That all will give an idea of the complete dataset, which is more intuitive than relying on PCA alone.
That is definitely true.
I'm new to bioinformatics and this is why asking your guidance on every step. ATpoint and Kevin Blighe thank you for the reply. I dropped idea of running EBSeq. And I ran Sample weight diagnostics in R with following sequence (Please, see table)
- Pipeline 1 – limma-trend + arrayWeights on logCPM This is the limma-trend workflow on logCPM (no voom), plus empirical sample quality weights via arrayWeights. Mean–variance trend is handled via eBayes(trend=TRUE).
- Pipeline 2 – voomWithQualityWeights (automatic) This is the “official” voom + quality-weights combo in one call. It applies: voom mean–variance modeling --> observation-level precision weights. arrayWeights-style sample weights --> per-sample quality. Combined weights used in subsequent lmFit.
- Pipeline 3 – voom + separate arrayWeights + manual combo Here we: Run voom to get logCPM + observation weights. Run arrayWeights on the voom-transformed expression matrix. Multiply sample weights into the observation weights, mimicking voomWithQualityWeights but more explicit.
And I am attaching "VOOM with Quality Weights" and "Simple total DEG counts per contrast × method" charts too for your advice. Please, guide me after considering all these tables and graph, the next sequence of steps to conduct.
I'm the limma author. I have to say, you are making a lot of unnecessary work for yourself by running limma in all these different ways. The different limma calling sequences should in practice give broadly similar results, and the fact that you are getting such different results is a sure sign that you have implemented them in ways that are not really equivalent. The differences you show are too large to arise from differences in the functions themselves. In particular, voomWithQualityWeights is simply calling voom and arrayWeights iteratively, so you must get completely identical results from voom and arrayWeights if you do it right.
The newest member of the voom family is edgeR::voomLmFit, and it really makes all the earlier pipelines obsolete, although the other functions are kept for backward compatibility. I would recommend that you just use voomLmFit with sample.weights=TRUE.
In terms of batch effects, I would be asking myself, what is different about the samples on the right hand side of your PCA plots, well away from the majority of the other samples? Is there any sample annotation that might explain why they look so different? Do they have unusual normalization factors?
I would also be using robust=TRUE when calling eBayes.
Your PCA plots don't suggest any clear differences between the treatments, so getting few or even no DE genes might be a correct result.
Nice to meet you Gordon Smyth I did not know that one day I will be encountering someone like you. Kindly, can you let me know whether the steps I am performing are fine and I attach number of DEGs and header results for HS vs C. I ran following script in R-studio:
############################################################
## 0. Setup: packages and file paths
############################################################
suppressPackageStartupMessages({
library(edgeR) # DGEList, filterByExpr, calcNormFactors, voomLmFit
library(limma) # makeContrasts, eBayes, topTable
library(dplyr) # data wrangling
library(tibble) # rownames <-> columns
library(readr) # read_tsv
})
## ---- Edit these paths for your project ----
counts_file <- "/Users/umair/Desktop/First_Project_paper_material/EBSeq/data/gene_count_description.txt"
meta_file <- "/Users/umair/Desktop/First_Project_paper_material/EBSeq/data/metadata.txt"
stopifnot(file.exists(counts_file), file.exists(meta_file))
############################################################
## 1. Load counts and metadata, and align samples
############################################################
## 1.1 Load count matrix
## Assumptions:
## - Tab-delimited
## - First column = gene_id
## - Columns C1..HT3 = sample columns
raw_counts_df <- readr::read_tsv(counts_file, show_col_types = FALSE)
gene_col <- names(raw_counts_df)[1] # "gene_id"
all_cols <- names(raw_counts_df)
## Identify sample columns: all columns that look like sample names
## For my dataset, these are: C1..C3, S1..S3, H1..H3, T1..T3, HS1..HS3, HT1..HT3
## We can detect them by excluding annotation columns.
annot_cols <- c("gene_name", "gene_chr", "gene_start", "gene_end",
"gene_strand", "gene_length", "gene_biotype",
"gene_description", "tf_family")
sample_cols <- setdiff(all_cols, c(gene_col, annot_cols))
## Build a clean counts dataframe: gene_id as rownames, only sample columns kept
counts_df <- raw_counts_df %>%
dplyr::select(all_of(c(gene_col, sample_cols))) %>%
dplyr::rename(gene_id = !!gene_col) %>%
tibble::column_to_rownames("gene_id")
## 1.2 Load metadata
## For my current metadata:
## - "treatments" column = sample labels (C1..HT3)
## - "group" column = biological treatment (C, H, S, T, HS, HT)
meta_df <- readr::read_tsv(meta_file, show_col_types = FALSE)
## Standardize column names so the rest of the code is generic
meta_df <- meta_df %>%
dplyr::rename(
sample_raw = sample, # keep original sample ID if needed
sample = treatments, # this is C1..HT3, matches counts columns
treatment = group # C, H, S, T, HS, HT
)
meta_df$sample <- as.character(meta_df$sample)
## 1.3 Align counts and metadata by sample name
common_samples <- intersect(colnames(counts_df), meta_df$sample)
if (length(common_samples) < 2) {
stop("Too few overlapping samples between counts and metadata.
Check sample names in both files.")
}
## Sort samples to a consistent order
common_samples <- sort(common_samples)
counts_df <- counts_df[, common_samples, drop = FALSE]
meta_df <- meta_df %>%
dplyr::filter(sample %in% common_samples) %>%
dplyr::arrange(match(sample, common_samples))
stopifnot(identical(colnames(counts_df), meta_df$sample))
## 1.4 Final counts matrix and treatment factor
counts <- as.matrix(counts_df)
if (!is.numeric(counts)) {
stop("Counts matrix is not numeric. Check that count columns contain only numbers.")
}
## Explicit treatment factor order
meta_df$treatment <- factor(
meta_df$treatment,
levels = c("C", "H", "S", "T", "HS", "HT")
)
############################################################
## 2. edgeR container: filtering + TMM normalization
############################################################
## 2.1 Construct DGEList
dge <- DGEList(counts = counts, group = meta_df$treatment)
## 2.2 Filter lowly expressed genes using design-aware rule
design_for_filter <- model.matrix(~ meta_df$treatment)
keep <- filterByExpr(dge, design = design_for_filter)
dge <- dge[keep, , keep.lib.sizes = FALSE]
## 2.3 TMM normalization
dge <- calcNormFactors(dge, method = "TMM")
############################################################
## 3. Design matrix and contrasts
############################################################
## 3.1 Design: no intercept, one column per treatment
design <- model.matrix(~ 0 + treatment, data = meta_df)
colnames(design) <- sub("treatment", "", colnames(design))
colnames(design)
# Expect: "C" "H" "S" "T" "HS" "HT"
## 3.2 Contrasts: all treatments vs control C
contrast_matrix <- makeContrasts(
HvsC = H - C,
SvsC = S - C,
TvsC = T - C,
HSvsC = HS - C,
HTvsC = HT - C,
levels = design
)
############################################################
## 4. Canonical voomLmFit with sample.weights = TRUE
############################################################
## voomLmFit does:
## - voom mean–variance modelling
## - limma linear modelling
## - when sample.weights = TRUE: estimates sample quality weights
## (i.e. automatically down-weights outlier samples)
fit_voom <- voomLmFit(
dge,
design = design,
sample.weights = TRUE, # <-- automatic sample quality weights
plot = FALSE
)
## Apply contrasts and empirical Bayes moderation
fit_voom <- contrasts.fit(fit_voom, contrast_matrix)
fit_voom <- eBayes(fit_voom, robust = TRUE)
############################################################
## 5. Helpers to extract DEGs with thresholds
############################################################
## 5.1 Annotate DE status for one contrast
annotate_de <- function(tt,
lfc_cutoff = 1,
fdr_cutoff = 0.05,
lfc_col = "logFC",
fdr_col = "adj.P.Val") {
tt$gene_id <- rownames(tt)
tt <- dplyr::relocate(tt, gene_id)
tt <- tt %>%
dplyr::mutate(
DE = dplyr::case_when(
.data[[fdr_col]] <= fdr_cutoff & .data[[lfc_col]] >= lfc_cutoff ~ "up",
.data[[fdr_col]] <= fdr_cutoff & .data[[lfc_col]] <= -lfc_cutoff ~ "down",
TRUE ~ "ns"
)
)
tt
}
## 5.2 Extract full topTable + DEG list for all contrasts
get_top_and_deg <- function(fit,
lfc_cutoff = 1,
fdr_cutoff = 0.05) {
coef_names <- colnames(fit$coefficients)
top_list <- vector("list", length(coef_names))
names(top_list) <- coef_names
deg_list <- vector("list", length(coef_names))
names(deg_list) <- coef_names
for (cn in coef_names) {
tt <- limma::topTable(fit, coef = cn, number = Inf, sort.by = "P")
tt <- annotate_de(tt, lfc_cutoff = lfc_cutoff, fdr_cutoff = fdr_cutoff)
top_list[[cn]] <- tt
deg_list[[cn]] <- tt %>% dplyr::filter(DE != "ns")
}
list(top_tables = top_list, deg_lists = deg_list)
}
############################################################
## 6. Run DEG extraction and summarize counts
############################################################
res_voomLmFit <- get_top_and_deg(
fit_voom,
lfc_cutoff = 1, # |log2FC| >= 1
fdr_cutoff = 0.05 # FDR <= 0.05
)
## 6.1 DEG counts per contrast
deg_counts_voom <- tibble(
contrast = names(res_voomLmFit$deg_lists),
n_DEG = sapply(res_voomLmFit$deg_lists, nrow)
)
deg_counts_voom
## 6.2 Example: inspect results for HS vs C
head(res_voomLmFit$top_tables$HSvsC)
Looks ok.
There is no need to create different design matrices for filtering vs DE.
I strongly discourage the use of fold-change cutoffs for DE tables. They interfer with FDR control and give priority to low-expressed genes of potentially less interest. They do not make the results more reliable or confident, somewhat the opposite.
I understand that fold-changes cutoffs are popular in the literature, but they have no scientific or statistical basis. Why in the world would it help to add an additional DE filter when you have so few DE genes to start with?
Thank you Gordon Smyth for your reply and opening my eyes as I am novice to this field and now attaching new table after I removed the FC cutoffs. The number of genes detected have increased. 
for your reference Gordon Smyth , one of the objective of my study is to find difference of impact between single and their co-treatments i.e C is control seedling, H is a fungal pathogen inoculated on a seedling, and HS is combination of a fungal pathogen + a bacterial species on a seedling.
When I see threads like this with never-ending comments and answers, I realise that the onus is on you to take reflection on the information already received from experts like Gordon and I, and to then formulate your own solution.
Thank you and kind regards,
Kevin
Kevin Blighe Thank you for your comment. My questions may appear basic, but they come from genuine curiosity. You are, of course, free to skip them. I could ask an AI tool, but I value learning from the real-world experience of people in this community. That’s why I joined Biostars and why I ask questions here.
You are the best professional that I have ever seen - never forget that. Keep up the hard work.
The sequencing company appears to have used raw p-values instead of p-values adjusted for multiple testing, which makes the significance results complete nonsense.
In general, limma gives more rigorous error rate control than DESeq2, so will give slightly fewer DE genes than DESeq2 for most vanilla analyes where both packages fit exactly the the same model. However, limma becomes more powerful when sample weights are used, and will often give more DE genes than DESeq2 when sample weights are appropriate.
For your data, the HS and HT treatments are the most variable as well as well as the most DE, so a vanilla analysis that treats all the treatments as equally variable may be suspect. The limma sample weight approach takes care of that, but I would be tempted to run voomLmFit with var.group=treatment, which will tell limma to estimate weights for each treatment group rather than for each individual sample. With only three samples in each group, I think that may be more reliable.
The huge difference in DE results that you see here is, however, not due to any difference between limma vs DESeq2 but rather to a naive and inappropriate use of DESeq2. If I had to guess, I suspect the sequencing company might have found no DE results with a regular analysis so decided to report a raw p-value cutoff instead. Most journal reviewers would jump on that however.
I'm in the business of returning reliable analyses. It is an intended feature of limma that it does always give significant DE results.
limma is truly the best
Thank you Professor Smyth, and good morning to you there in Oz
Gordon Smyth thank you for your detailed reply and I will acknolwedge you in my paper. I will stop asking you questions as Kevin Blighe seems angry on me for my curiosity.
Not angry, dear Sir - just encouraging you to be the best Professional that you could ever be.
I and Professor Smyth are here for you.
See my suggestion about var.group, which I added to my previous comment, not noticing that you have already replied.
Log in to answer this question.
-
missing values for batch in RNA-seq data
written by donan •
-
Extracting CDS using python
written by Davi.marcon •
-
Data processing order: normalisation or batch effect correction first?
written by Sebastian Hesse
-
How can I extract the sequencing reads containing a specific linker/tag?
written by naeem40thju •
-
HOMER error -bash: ./ctcf.motif: Permission denied
written by gerardo_za_94 •
-
Extra adjustment variables during batch effect correction
written by n,n
-
Qs on merging microarray datasets and different ways on building a classifier on them.
written by arronar
-
Best way to merge RNA-seq data from two sequencing runs of the same samples
written by unsupervised_learner •
-
RNA-seq unbalanced batch effect correction
written by endikavarela •
-
News: New simulation dataset for 'A benchmark for RNA-seq quantification pipelines'
written by Michael Love
