Hi there,
I've got somewhat of a tricky situation in deciding how to go about setting up my differential expression analysis in DESeq2 for the experiment that was done and the questions being asked. In this experimental set up, a bacterial strain is either infected or not infected with a species specific phage.
I have three different genotypes, two conditions (infected/uninfected), and four time points. My main concern is whether I should/if it is appropriate to split up count matrices for the bacteria and the phage.
We want to look at which genes are differentially expressed between different genotypes when the bacteria has been infected by the phage, but of course if there is no phage then the number of reads mapping the each phage gene in the uninfected samples is zero across the board. Based on the DESeq2 vignette, I believe I can pre-filter based on a certain minimum value of reads, which would remove the phage genes anyway, and I think this solves my problem of zeros in this case. It seems to be that simply splitting the count matrix into bacterial/phage counts solves this as well, but I'm not sure how this would affect normalization/estimation of statistical parameters.
We also want to look at how the expression of viral genes changes between genotypes over time. In this case, it's only the viral counts that are of interest, and the "uninfected" condition is irrelevant. So this seems like another situation where simply splitting the counts into two separate phage and bacteria matrices also leads to a simplified scenario.
Is there anything inherently wrong with this approach? Or will I be totally messing up the estimation of parameters for the model to be accurate? Any input would be much appreciated as this is my first shot at RNAseq data analysis and it's a bit more complicated (as far I can tell) than I might have hoped for in a first attempt.
Thanks,
Sean
Can you tell us some stats which may help people trying to answer your question. Total numbers of reads for each sample (were these comparable), Overall % of reads that aligned (bacteria+virus) for each sample, out of the aligned reads what fraction of reads were aligned to bacteria + virus for each sample? Were these numbers consistent across samples? Do you have biological replicates for these samples for all categories?
Ah yes, of course. I have biological triplicate for all samples, two time points for each uninfected sample (t=0, and t=30), four time points for each infected sample (t=0, t=10, t=20, t=30), the and the three genotypes, giving a total of 54 samples.
The sequencing was done as paired-end 150 bp, alignments were carried out using bowtie2 in end to end mode (since splicing not a concern) after adapter removal with cutdadapt.
For the tables below c/t denotes uninfected or infected, the number after for timepoints, r1/r2/r3 for each biological replicate. It seems that the total reads aligning to samples is fairly consistent, some slight outliers here and there but no obvious pattern between genotypes, treatments, or time points.
I'm sorry if the formatting is horrendous, I just did markdown tables but not sure the best way to go about it otherwise on here.
TOTAL READS PER SAMPLE
The averages for genotype_1, genotype_2, and genotype_3 are 9346603, 10271899, 9187888, respectively. So a bit higher for genotype_2 but nothing that seems crazy?
For alignment percentage, all samples had >95% of read pairs aligning concordantly exactly 1 time, <5% aligning concordantly more than 1 time.
PERCENT READS ALIGNING TO BACTERA/VIRUS/PLASMID
I think things looks fairly consistent across samples. We see fewer reads aligning to the virus in genotype_2, but that is expected based on what is known about the genotype and this particular infection model. The plasmid in the bacteria also has a fairly consistent number of reads and we shouldn't expect it to vary much across samples. The percent of reads mapping is a bit higher in genotype_3, but that's almost certainly because there is a gene on the plasmid that is not present in the other two genotypes, so I would say that is expected as well.
Agreed. With a two genome experiment this can be a little tricky but you seem to have consistent data. Hopefully you will have some suggestions soon.