Swish (fishpond) scaled counts result
1
0
Entering edit mode
675743 ▴ 10
Last seen 5 weeks ago

Dear Fishpond team.

I've been taking a look at the Swish (fishpond) package for my DTE/DGE analysis and it seem to fit my needs perfectly, however I have a small question about the scaled counts output.

As I understand, it is possible to plot the results of each transcript/gene with plotInfReps via idx but is it possible to also extract a matrix containing the scaled counts for all transcripts/genes? I would like to make my own plots of the scaled counts result.

Best regards, Fredrik

Swish fishpond • 164 views
0
Entering edit mode
@mikelove
Last seen 9 hours ago
United States

Thanks for your interest in using Swish.

So here is a typical run of swish:

library(SummarizedExperiment)
library(fishpond)
y <- makeSimSwishData()
y <- scaleInfReps(y)
## Progress: 20 on 20
y <- labelKeep(y)
set.seed(1)
y <- swish(y, x="condition")
## Generating test statistics over permutations
## Progress: 100 on 100


Then we can look at a particular feature (say a gene):

plotInfReps(y, which.min(mcols(y)\$qvalue), "condition")


This plots the distribution of the Gibbs samples (or bootstraps) for each sample.

The values that are plotted are stored in a series of matrices: the matrices named "infRep..." in assayNames(y).

If you want the values for a particular feature you can use the code:

infRepMatrix <- unlist( assays(y[idx,])[ grep("infRep", assayNames(y)) ] )


This matrix is nreps x nsamples. You could then compute the column median (the thick lines in the boxplot) or the column mean of these.

If you want to compute this for all features at once you can use the abind function to make an array and then use apply:

infReps <- assays(y)[ grep("infRep", assayNames(y)) ]
infArray <- abind::abind( as.list(infReps), along=3 )


infArray is nfeature x nsamples x nreps. You can compute summaries like so:

infMed <- apply(infArray, 1:2, median)
infMean <- apply(infArray, 1:2, mean)

1
Entering edit mode

Hello Michael.

Thank you for taking your time and answering the question in such a quick manner. The code you've supplied looks great, I will run the analysis in Swish and give it a go!

Thanks!

Regards Fredrik

0
Entering edit mode

Hi again Michael.

I tried the package and it looks great! I have one last question or request, if one would like to use the inbuilt plots (for portraying the inferential replicates) is it possible to change the name of the x-axis from sample number to sample name?

Thanks!

Best Regards, Fredrik

0
Entering edit mode

Sure, we went with numbers just because they are short, so work regardless of whether there are few or many replicates.

You can customize the x-axis by specifying: xaxis=FALSE, xlab="" in the call to plotInfReps().

Then you can add any labels you like, e.g.:

axis(1, 1:10, letters[1:10])

0
Entering edit mode

Hi Michael.

Thank you for the quick reply.

Best Regards, Fredrik