MISO: Best-guessing insert size and stddev versus running single-ended for PE data

[fpbarthel]

I am conducting a miso analysis for a single gene on a large number of Paired End RNAseq bams that is spread out on several external HDs (twenty 4tb drives). Computing insert sizes for these is extremely tedious and slow. However, running MISO itself (in single-end mode or paired end mode with best guess insert sizes) for a single gene runs extremely fast. Because of the speed, calculating exact insert sizes is not an option. For several samples, the mean insert size is around 150-250 and the stddev is more or less 50. Therefore, in my opinion I have two options:

1. Run miso in single-end mode
2. Run miso in paired-end mode and best guessing insert size at 250 and sd at 50

I wonder if option 2) may be superior because the data is paired end. Can anyone comment on their thoughts? Would be extremely helpful.

[Brian Bushnell]

I would expect most analyses to be better with paired data, though I have never personally used MISO. You can calculate insert-size distributions quite rapidly with BBMerge. How long are your reads, and what kind of organism are the from?

[fpbarthel]

Hi Brain,

Does BBMerge insert-size distribution calculation work efficiently with transcriptome mapped reads, where reads may span exons that are kilobases apart?

At this point I only have miniBAMs for a single gene, generated from the large BAMs stored on external drive. These minibams may have as little as 20 reads. All is human data (hg19 aligned). Read lengths are 48 or 75, depending on the sample.

Thanks!

[Brian Bushnell]

Oh, that could be a bit of a problem... BBMerge does not care about the presence of introns, but it does require the reads to be overlapping or near-overlapping (so, those reads are probably too short for a 250bp average insert size). Although as long as you map to the transcriptome, an insert size calculated from mapping will be be unaffected by introns. There will still be a bit of uncertainty due to differential splicing, but I think you'll still get a pretty accurate value.

For example, if you convert a large bam to fastq, and the reads are in their original order or name-sorted, you can run BBMap like this:

bbmap.sh ref=transcriptome.fasta in=reads.fq reads=1m ihist=ihist.txt interleaved

That will map the first 1 million pairs to the transcriptome and calculate the insert size distribution, which should only take a minute or so per bam file.

If BBTools and samtools are installed, you can do the conversion like this:

reformat.sh in=x.bam out=x.fq reads=2m

...which will just convert the first 2 million reads (1 million pairs) of the bam file, so you don't have to convert the whole thing. Again, though, the bam file must be unsorted (original order) or name-sorted.