weird dispersion plot

[scseekers]

Hi @hartleys

I have RNASeq data for 4 tissue types A,B,C,D. I wish to perform a pairwise comparison of these samples (A-B, B-C, C-D, B-D, A-D) and therefore used the following code:

library(JunctionSeq)
gtf.file <- "./withNovel.forJunctionSeq.gff.gz"
samplesheet <- read.csv("samplesheet.csv", header = T, stringsAsFactors=FALSE)
jscs <- runJunctionSeqAnalyses(sample.files = countFiles,
                               sample.names = samplesheet$sample.names,
                               condition=factor(samplesheet$condition),
                               flat.gff.file = gtf.file,
                               nCores = 12,
                               analysis.type = "junctionsAndExons"
)

Also, in results, I get LFC values of merely A-B, A-C, A-D in results. Is there a way to set constrasts like in DESeq2.
However, I get a weird Dispersion plot. What could be the possible reason?

When I run two samples at a time like (B-D) or say (A-B) the exons and junctions follows the trend line:

[scseekers]

Well!
I also observed one weird thing. On comparing all samples together (A, B, C, D), I get a comparison (AvsB, AvsC, AvsD), and a single "genewiseResults.txt.gz" containing 213 Differentially Used Features (exons+junctions) at FDR of 0.01. However, when I run the Junctionseq to achieve, all vs all comparison using a FOR loop as shown below, I don't get any significant results at FDR of 0.05 or 0.01. Why such a discrepancy?
My samples_dxd looks like (.csv file)
Arep1,A
Arep2,A
Arep3,A
Brep1,B
Brep2,B
Brep3,B
Crep1,C
Crep2,C
Crep3,C
Drep1,D
Drep2,D
Drep3,D

and my jseq_group looks like (.csv file)
pair,group1,group2
A_B,A,B
A_C,A,B
C_B,C,B
C_D,C,D
B_D,B,D

for (i in 1:nrow(jseq_group)) {
samples_dxd <- subset(samplesheet, grepl(jseq_group$group1[i],samplesheet$condition) | grepl(jseq_group$group2[i], samplesheet$condition))
countFiles <- paste0("F:/",samples_dxd$sample.names,"/QC.spliceJunctionAndExonCounts.withNovel.forJunctionSeq.txt.gz")
jscs <- runJunctionSeqAnalyses(sample.files = countFiles,
                               sample.names = samples_dxd$sample.names,
                               condition=factor(samples_dxd$condition),
                               flat.gff.file = gtf.file,
                               nCores = 1,
                               analysis.type = "junctionsAndExons")
save(jscs, file = paste0(jseq_group$Pair[i],".Rdata"))
rm(jscs)
}
     
for (i in 1:nrow(jseq_group)) {
load(paste0(jseq_group$Pair[i],".Rdata"))
     writeCompleteResults(jscs,
                     outfile.prefix=jseq_group$Pair[i],
                     save.jscs = TRUE, FDR.threshold = 0.05)
     rm(jscs)
}