What type of counts data to import for performing Isoform analysis in edgeR
2
0
Entering edit mode
@mohammedtoufiq91-17679
Last seen 4 hours ago
United States

Hi,

I am interested in performing isoform analysis on short read data (150bp) using the edgeR package, following the example from the edgeR User's Guide, section "4.6 Differential transcript expression of human lung adenocarcinoma cell lines."

I ran my pipeline using the nf-core/rnasplice pipeline and obtained counts and TPM values: nf-core rnasplice output:

From the example 4.6 in the edgeR User's Guide, I tried importing the "quant.sf" files but experienced difficulties. Based on the example, I imported scaled counts as suggested:

# Define the path to your TSV file
file_path <- "/projects/salmon/tximport/salmon.merged.transcript_counts_scaled.tsv"

# Import the TSV file
scaled.counts <- read.delim(file_path, header = TRUE, sep = "\t", row.names = 1)

# Create DGEList object
y <- DGEList(counts = scaled.counts, samples = Samples_metadata)
dim(y)

Which of the following files makes the most sense to import for "4.6 Differential transcript expression of human lung adenocarcinoma cell lines"?

Counts from nf-core/rnasplice:

  • salmon.merged.transcript_counts.tsv: Matrix of isoform-level raw counts across all samples.
  • salmon.merged.transcript_counts_scaled.tsv: Matrix of isoform-level scaled raw counts across all samples.
  • salmon.merged.transcript_counts_length_scaled.tsv: Matrix of isoform-level length-scaled raw counts across all samples.
  • salmon.merged.transcript_counts_dtu_scaled.tsv: Matrix of isoform-level dtu scaled raw counts across all samples.

TPMs from nf-core/rnasplice:

  • salmon.merged.transcript_tpm.tsv: Matrix of isoform-level TPM values across all samples.
  • salmon.merged.transcript_tpm_scaled.tsv: Matrix of isoform-level scaled TPM values across all samples.
  • salmon.merged.transcript_tpm_length_scaled.tsv: Matrix of isoform-level length-scaled TPM values across all samples.
  • salmon.merged.transcript_tpm_dtu_scaled.tsv: Matrix of isoform-level dtu scaled TPM values across all samples.

limma isoform edgeR salmon R • 483 views
ADD COMMENT
1
Entering edit mode
ATpoint ★ 4.6k
@atpoint-13662
Last seen 10 hours ago
Germany

You don't need any of these files. You need the catchSalmon function which reads quant.sf and the inferential (Bootstrap/Gibbs) replicates. Where in nf-core you can find that I cannot tell, but this is what the user guide tells, see also ?catchSalmon.

ADD COMMENT
0
Entering edit mode

ATpoint Thank you. I did wanted to try using the quant.sf files from the salmon run by referring edgeR user guide example 4.6 Differential transcript expression of human lung adenocarcinoma cell lines I tried the below steps, but I see there are NA values in the Overdispersion, and after running scaled.counts, the table is populated with NA values (see below screenshots)

> setwd("/projects/RNA/salmon/")

> quant <- dirname(list.files(".", "quant.sf", recursive = TRUE, full.names = TRUE))

> print(quant)

[1] "./DD21_progenitors" "./DD39_progenitors" "./DD48_progenitors" "./DD98_progenitors"

> catch <- catchSalmon(paths = quant)

Reading ./DD21_progenitors, 427366 transcripts, 0 none samples
Reading ./DD39_progenitors, 427366 transcripts, 0 none samples
Reading ./DD48_progenitors, 427366 transcripts, 0 none samples
Reading ./DD98_progenitors, 427366 transcripts, 0 none samples

> scaled.counts <- catch$counts/catch$annotation$Overdispersion

ADD REPLY
1
Entering edit mode

Probably nf-core did not run with bootstrap/gibbs replicates. Check the code. See salmon docs on how to turn that on.

ADD REPLY
0
Entering edit mode

The poster needs to run Salmon with --numGibbsSamples=50. I suspect that nf-core just won't do that, so Salmon needs to be run directly.

ADD REPLY
1
Entering edit mode
@gordon-smyth
Last seen 36 minutes ago
WEHI, Melbourne, Australia

I haven't used nf-core but, from browsing through their documentation, they didn't seem to have any pipeline for differential transcript expression. To run the edgeR pipeline, you should run Salmon directly with --numGibbsSamples=50. Otherwise, edgeR has no way to measure the assignment uncertainty of each transcript count. See

Baldoni PL, Chen L, Smyth GK (2024). Faster and more accurate assessment of differential transcript expression with Gibbs sampling and edgeR v4. NAR Genomics and Bioinformatics 6(4), lqae151. https://doi.org/10.1093/nargab/lqae151

Note that the edgeR concept of "scaled counts" is quite different from scaled counts in nf-core.

ADD COMMENT
0
Entering edit mode

Gordon Smyth

I am exploring workflows for differential transcript expression analysis and have been using the nfcore rnasplice pipeline, which outputs a variety of files, including downstream results:

https://nf-co.re/rnasplice/1.0.4/ https://nf-co.re/rnasplice/1.0.4/docs/output/

In particular, I have been working with the salmon.merged.transcript_counts_length_scaled.tsv file, which provides a matrix of isoform level length-scaled raw counts across all samples. I wanted to ask if you think the following workflow (below) I am trying is appropriate for differential transcript expression analysis.

In the meantime, I will explore and re-run Salmon with the --numGibbsSamples=50 option to account for transcript count uncertainty before proceeding as well.

# Create DGEList object
y <- DGEList(counts = scaled.counts,
             samples = Samples_metadata)

# Filtering
keep <- filterByExpr(y, design)
table(keep)
y <- y[keep, , keep.lib.sizes=FALSE]
y <- normLibSizes(y, method = "TMM")

## Obtaining Log CPM values
logCPM <- cpm(y, log = TRUE, prior.count = 1)
class(logCPM)
logCPM <- as.data.frame(logCPM)

##  Dispersion estimation
y <- estimateDisp(y, design, robust=TRUE)
y$common.dispersion
plotBCV(y)

fit <- glmQLFit(y, design, robust=TRUE)
plotQLDisp(fit)

## Differential expression
contr <- makeContrasts(Treatment - Control, levels=design)
qlf <- glmQLFTest(fit, contrast=contr)
is.de <- decideTests(qlf, adjust.method = "none", p.value=0.005)
summary(is.de)
tt <- topTags(qlf,n = Inf)

# Get results
results <- topTags(qlf, n = Inf)

Your insights would be greatly appreciated.

ADD REPLY
0
Entering edit mode

ATPoint and I have both already told you that you cannot use any of the output files from nf-core, and why that is so. I explained to you the correct way to do the analysis, and I have referred you to a reference that explains the issues.

It is inappropriate to input the salmon.merged.transcript_counts_length_scaled.tsv file to an edgeR analysis. It will give poor results. It will poor statistical power overall but also with potential false discoveries for transcripts that are short or which overlap other transcripts.

ADD REPLY

Login before adding your answer.

Traffic: 531 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