Hey Everyone, I have a concept question: batches. So, my input samplesheet is this:

sample	fastq_1	fastq_2	strandedness	condition
I1	/home/tud03125/01.RawData/01.RawData/I1/I1_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I1/I1_2.fq.gz	reverse	FASTED
I2	/home/tud03125/01.RawData/01.RawData/I2/I2_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I2/I2_2.fq.gz	reverse	FASTED
I3	/home/tud03125/01.RawData/01.RawData/I3/I3_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I3/I3_2.fq.gz	reverse	HCD
I4	/home/tud03125/01.RawData/01.RawData/I4/I4_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I4/I4_2.fq.gz	reverse	HCD
I5	/home/tud03125/01.RawData/01.RawData/I5/I5_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I5/I5_2.fq.gz	reverse	HCD_MK2206
I6	/home/tud03125/01.RawData/01.RawData/I6/I6_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I6/I6_2.fq.gz	reverse	HCD_MK2206
I7	/home/tud03125/01.RawData/01.RawData/I7/I7_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I7/I7_2.fq.gz	reverse	HCD_MK2206
I8	/home/tud03125/01.RawData/01.RawData/I8/I8_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I8/I8_2.fq.gz	reverse	HCD_MK2206
I9	/home/tud03125/01.RawData/01.RawData/I9/I9_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I9/I9_2.fq.gz	reverse	HCD_SPHINX31
I10	/home/tud03125/01.RawData/01.RawData/I10/I10_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I10/I10_2.fq.gz	reverse	HCD_SPHINX31
I11	/home/tud03125/01.RawData/01.RawData/I11/I11_1.fq.gz	/home/tud03125/01.RawData/01.RawData/I11/I11_2.fq.gz	reverse	HCD_SPHINX31
IA1	/home/tud03125/IA_dataset/FASTQ/fasted_rep1_R1.fq.gz	/home/tud03125/IA_dataset/FASTQ/fasted_rep1_R2.fq.gz	reverse	FASTED
IA2	/home/tud03125/IA_dataset/FASTQ/fasted_rep2_R1.fq.gz	/home/tud03125/IA_dataset/FASTQ/fasted_rep2_R2.fq.gz	reverse	FASTED
IA3	/home/tud03125/IA_dataset/FASTQ/fasted_rep3_R1.fq.gz	/home/tud03125/IA_dataset/FASTQ/fasted_rep3_R2.fq.gz	reverse	FASTED
IA4	/home/tud03125/IA_dataset/FASTQ/highcarb_rep1_R1.fq.gz	/home/tud03125/IA_dataset/FASTQ/highcarb_rep1_R2.fq.gz	reverse	HCD
IA5	/home/tud03125/IA_dataset/FASTQ/highcarb_rep2_R1.fq.gz	/home/tud03125/IA_dataset/FASTQ/highcarb_rep2_R2.fq.gz	reverse	HCD
IA6	/home/tud03125/IA_dataset/FASTQ/highcarb_rep3_R1.fq.gz	/home/tud03125/IA_dataset/FASTQ/highcarb_rep3_R2.fq.gz	reverse	HCD

There’s two different samples: one with

I

and the other with

IA

prefix characters. They’re two different datasets and thus two different batches. I know that

RNA-seq

pipeline does

DESeq2

, which addresses batches. But, I don’t think

RNAsplice

have

DESeq2

. Plus, my boss said that when she used

SUPPA

(not through this pipeline), they don’t consider batches. So, can

RNAsplice

address or do batch correction? Or, is it alright to run

RNAsplice

with datasets that potentially can have different batches?

I know that

RNA-seq

pipeline does

DESeq2

...

It technically does, but only for QC purposes -- you get a counts table (or equivalent) and use that to your statisitics outside the pipeline. #CE8SSJV3N receives no metadata and has no way to take batches into account, but for the counts table you don't need that. #C0493FQQZ3P, on the other hand, does quite a bit of statistics, and some tools are more sophisticated than others. I don't know enough about SUPPA, but I don't think rMATS is equipped to use batch information, but e.g. DEXSeq definitely can, provided you give it the proper model design and metadata. I don't think #C0493FQQZ3P has a mechanism of doing that -- its metadata is limited to the

condition

column in the sample sheet. If you wanted to do that, you would have to run DEXSeq (or any other tool with that capability) outside the pipeline with appropriate parameters, although you could use files and indices generated by the pipeline as input for that.

@Thomas Danhorn Got it. The other question: can different batches affect PSI values (whether in SUPPA or rMATS)? If they do, is there a way to incorporate batch corrections towards the calculated PSI values (or "IncLevel" for rMATS)? So when saying "I don't think rMATS is equipped to use batch information," is that like saying "rMATS does not consider batch information when calculating their PSI values?"

Yes, as far as I know rMATS does not consider batches. > can different batches affect PSI values? Presumably. The questions is, what do you call a "batch", and what are the differences between batches? You might use that to refer to samples being prepare on different days, possibly by a different person, but otherwise with the same kit, or maybe with different kits, or you might refer to samples cultivated at different times, which mean (slightly) different environmental conditions, possibly different batches of reagents, etc. How does e.g. a library prep kit affect the results of sequencing? I don't know -- maybe the buffer composition affects sequences with different base compositions (GC content) or RNA with different secondary structures differently. Maybe this leads to a bias in favor of certain genes or transcripts. Maybe something leads to more or less degradation of RNA, and this may affect short transcripts more than longer ones. When looking purely at gene expression, this may not matter that much, but for splicing it might. I don't know if these technical effects are string enough to meaningfully change the results. I'd probably be more worried about biological effects from (slightly) different conditions, because a small input difference can be amplified by regulatory cascades. How exactly this will affect the statistic also depends -- is the batch confounded with the experimental variable, e.g. were all controls grown and processed one day, and the treatments a month later? In that case you are screwed anyway, because there is just no way to tell if some change is from the treatment or the batch. But if you have e.g. three replicates of each, and a pair of control + treatment was done each day, then then the batch basically introduces "noise" -- the replicates become more different from each other. Since a difference is only significant if the between-group changes are clearly larger than the within-group variance (i.e. the "noise"), increased noise leads to less significance. Not great, but you'll have to live with it if there is no way to correct for it. > If they do, is there a way to incorporate batch corrections towards the calculated PSI values? I don't see a good way of doing that. With RNAseq, you get counts, and you can shrink or inflate them (e.g. with COMBATseq). For splicing this is much more complex, and especially for the JC files in rMATS, you have junction spanning reads, not just counts, and for JCEC you have to integrate those with non-junction reads that "support" a transcript. Given that the input are BAM files, there is no way to modify them in any way to "batch-correct" the resulting PSIs, and even if there is some way to incorporate covariates like batch, rMATS would have to handle that later, and given the complexities, I don't see a practical way of doing that (but someone smarter than me might). DEXSeq is a different beast, because in terms of methodology it is basically like DESeq(2) applied to exons -- you have counts, and you can add variables (including batch) to your model of how these counts come about.