[bob-loblaw]

Also I don't want to send you down a confusing path. I don't mind providing you with some help to get that pipeline working. One more thing to consider - do you expect insertions/deletions to be important? If so then RSEM may not be what you want since it uses bowtie1 for alignments. eXpress is a similar solution and with eXpress you can use alignments from bowtie1, bowtie2, bwa (with some tweaking) and really any aligner that can output all possible alignments for a given read. These tools attempt to disambiguate the alignments to a set of gene/protein/transcript sequences giving you "unique" mappings for even reads that can align equally well to several references. I've done a bit of benchmarking and honestly I haven't seen great results from eXpress but RSEM does pretty well. Both work VERY well if you are able to sum counts of sequences together for sequences that share exons or share sequence (as in multi-copy genes or alternatively spliced genes). They work OK in terms of per-sequence level counts - certainly better than what the aligners can do on their own - but certainly not perfect. Just keep in mind that you're per-sequence expressions will likely contain some false positives (maybe a lot...) and will also likely be missing a few true positives. In the end you're knowledge of which sequences in your database share sequence or share exons will help you immensely in getting stable and reliable read counts.

Well alternatively spliced genes won't be a problem, it's all bacteria I'm mapping to. RSEM I have been playing around with, but I have a pretty large sample size, and realigning all of them would be very time consuming, obviously I'll do it if necessary but I'd prefer not to have to. I haven't tried eXpress yet, but I will.

When I was mapping with bowtie2 I left it's reporting mode in default (i.e. report only the best alignment) but eXpress wants to be able to select the best alignment itself. Do you think this will be a big issue?

[sdriscoll]

Also I don't want to send you down a confusing path. I don't mind providing you with some help to get that pipeline working. One more thing to consider - do you expect insertions/deletions to be important? If so then RSEM may not be what you want since it uses bowtie1 for alignments. eXpress is a similar solution and with eXpress you can use alignments from bowtie1, bowtie2, bwa (with some tweaking) and really any aligner that can output all possible alignments for a given read. These tools attempt to disambiguate the alignments to a set of gene/protein/transcript sequences giving you "unique" mappings for even reads that can align equally well to several references. I've done a bit of benchmarking and honestly I haven't seen great results from eXpress but RSEM does pretty well. Both work VERY well if you are able to sum counts of sequences together for sequences that share exons or share sequence (as in multi-copy genes or alternatively spliced genes). They work OK in terms of per-sequence level counts - certainly better than what the aligners can do on their own - but certainly not perfect. Just keep in mind that you're per-sequence expressions will likely contain some false positives (maybe a lot...) and will also likely be missing a few true positives. In the end you're knowledge of which sequences in your database share sequence or share exons will help you immensely in getting stable and reliable read counts.

[bob-loblaw]

If you really want to align to this database you're using I suggest trying RSEM.

I'll look into it. Thank you.

[sdriscoll]

If you really want to align to this database you're using I suggest trying RSEM.

[sdriscoll]

Exactly as Simon said. If you were mapping to a genome reference then idxstats would return read counts per chromosome. It's absolutely more complicated to map to a transcriptome reference. A couple tools for that are eXpress and RSEM but neither of those will help you get counts at the gene level without you providing some knowledge of which references are from the same gene.

Probably the most straightforward approach is to align your reads to a genome reference (full chromosome sequences) with Tophat or STAR, if you have the RAM for it, then to count hits to genes with something like htseq-count which can find overlaps of genomic coordinates with gene features annotated in a GTF file.

[bob-loblaw]

This is what sdriscoll meant, when he said that you are mapping to a transcriptome. You have provided your aligner with a FASTA file which did not contain one sequence for each chromosome but one sequence for each transcript. (Otherwise, how would samtools know where the transcripts are, as you haven't supplied a GFF file.) This is known as "mapping against the transcriptome" and it is "bad" if you don't know exactly what you are doing, for various reasons that you'll find in old threads here.

Ah okay. I misunderstood, I thought sdriscoll was asking if I was mapping against a reference assembled from transcriptomics data, as opposed to the DNA sequences of predicted proteins from sequenced genomes (Which is what I'm using). Sorry for being a "student" and for for making a "mistake". I'll correct that post

[Simon Anders]

This is what sdriscoll meant, when he said that you are mapping to a transcriptome. You have provided your aligner with a FASTA file which did not contain one sequence for each chromosome but one sequence for each transcript. (Otherwise, how would samtools know where the transcripts are, as you haven't supplied a GFF file.) This is known as "mapping against the transcriptome" and it is "bad" if you don't know exactly what you are doing, for various reasons that you'll find in old threads here.

[bob-loblaw]

Actually, I don't get it. How do you get per-gene count with idxstat? I thought it only tells you the number of reads mapped to each reference sequence, i.e., to each chromosome. How do you get individual genes?

Maybe that depends on the reference that you've mapped to?

For me it throws out a table which has every gene ID from the genome I'm mapping to, along with it's length, and how many reads are mapped to that reference. For some of them it just gives me and organism and coordinates for the genome, but I can get more information about that from a GFF file I have of the annotations.

[Simon Anders]

Actually, I don't get it. How do you get per-gene count with idxstat? I thought it only tells you the number of reads mapped to each reference sequence, i.e., to each chromosome. How do you get individual genes?

[bob-loblaw]

This is all sounds as if some mistake happens during counting. Your approach of using samtools idxstats is rathe runorthodox, and I wonder if it is correct. It might be much safer to use some well-tested tool to obtain a count table instead of using some home-brewn solution.

From what I can see idxstats works perfectly, it needs to a bit of tweeking in R so that DESeq can read it, but nothing major.

I am open to suggestions though, can you give some examples of those well tested tools? Thanks

[Simon Anders]

This is all sounds as if some mistake happens during counting. Your approach of using samtools idxstats is rathe runorthodox, and I wonder if it is correct. It might be much safer to use some well-tested tool to obtain a count table instead of using some home-brewn solution.

[bob-loblaw]

So, to rephrase, you're mapping to a transcriptome reference and not a genome, correct? When you do this its very important that the rows you merge be all features that share exonic sequence - as in all alternative isoforms of a gene and even multi copy genes that reside in separate loci. Otherwise you're going to run into some confusion for sure. Bowtie is not designed nor capable of making alignment decisions between features with shared sequence beyond total random selection. It does make the same decision each time you run it but that's by design...they do a random number seeding trick to ensure this happens. So even if you're using an Ensemble annotation for mouse or human simply merging counts based on any of the provided ids isn't enough to remove the ambiguity problem.

I'm mapping against the DNA sequences of predicted proteins from sequenced genomes. The only rows that I merge are the ones which have identical IDs I hadn't even considered merging isoforms...

We originally did it to cut down on computational time more than anything (although it actually doesn't make that big of a difference in terms of the size of the count table)

[sdriscoll]

So, to rephrase, you're mapping to a transcriptome reference and not a genome, correct? When you do this its very important that the rows you merge be all features that share exonic sequence - as in all alternative isoforms of a gene and even multi copy genes that reside in separate loci. Otherwise you're going to run into some confusion for sure. Bowtie is not designed nor capable of making alignment decisions between features with shared sequence beyond total random selection. It does make the same decision each time you run it but that's by design...they do a random number seeding trick to ensure this happens. So even if you're using an Ensemble annotation for mouse or human simply merging counts based on any of the provided ids isn't enough to remove the ambiguity problem.

[bob-loblaw]

Dear bob loblaw

thank you. The behaviour you report is not reasonable, and somewhere in your workflow or tool chain there must be a mistake. Can you report the sequence of steps (script) that you perform, to reproduce your observation?

Best wishes
Wolfgang

Well my workflow is very simple. I'm using bowtie2 to map my reads to a reference database, and then I'm using samtools idxstats function to create a count table. Then I merge rows with duplicate IDs, then I put it into R and run DESeq on it. It's a very simple workflow.

I've also noticed that for the random subset of reads and the whole file the exact same proportion of reads is mapped (i.e. if 67.4% is mapped for the whole file, then in the subset file 67.4% is also mapped) which at least indicates that this difference in the number of differentially expressed genes between the whole file and subset isn't down to some non random effect in the subsetting process.

Any ideas?
Thanks

[sdriscoll]

i good rule of thumb when doing random subsetting of data is to do it more than once so you can observe that the results are stable and not jumping around due to the random sampling. So...sure go ahead and subset down to 9,000,000 reads but run a few iterations of the entire pipeline. if the results are stable you're good. if not then subsetting may not be appropriate, for whatever reason.