Hi,
I analysed my RNA-seq data for S. pombe once with Tuxedo pipeline and once with DESeq! Surprisingly I see one of the genes which is important for my experiment present in Tuxedo gene_exp.diff
output but not in deseq.res
!
In gene_exp.diff
I get nmt1 III:1837498-1844115
but when I look for it on pombase it has other coordination starting 1838335 and ending at 1839525 and I don't see this gene in deseq.res output at all. Am I missing something with annotations or there some other thing that I am not considering?
Your help is much appreciated.
/Parham
I am doing this from the RNA-seq workflow. I have all my
accepted_hits.bam
andaccepted_hits.bam.bai
intophat_all
. Does it help? Please let me know if you need more info!Yup, that helps. I would guess that either
exByGn
doesn't containnmt1
or for some reason the counts for it are 0 and it's getting filtered when you subsetcnts
. Both of those are things you should be able to check.I could check
cnts
and it had zero value for allbam
files! So it must beexByGn
not containingnmt1
as you said, but I couldn't figure out how to check it. Would you give me a hint on that?If it's all 0 in cnts then it'd be easier to just redo the counting with featureCounts.
Alright I did it again with
featureCounts
. Thanks for suggesting it! But here it still gives all zero for all BAM mapped reads! Here is counts.txt output and counts.txt.summary. Since I never used it before I am writing the code I used for this purpose. If you could kindly have a look on that if I am doing something wrong.This is due to the ncRNA on the opposite strand. Since you have a stranded protocol, it's impossible to discriminate between reads mapping to the ncRNA and nmt1 (well, you can and do get reads assigned to that ncRNA, but that's because part of it doesn't overlap another feature). I'm guessing that cufflinks is either ignoring the ncRNA or its EM algorithm prioritizes against it (or perhaps it uses coverage in gauging which gene to assign counts to, though I could see that back firing).
When I was running Tuxedo pipeline I was giving
--library-type -frsecondstrand
for all steps! Does that explain it?It does. If these really are stranded libraries then tell summarizeOverlaps() or featureCounts that.
Great now it has reads and makes sense! Thanks a lot!