NexteraXT MiSeq Bowtie2 strange insert size distribution

[M4TTN]

Dear all, we are trying to reproducibly determine the size distribution of our inserts obtained from NexteraXT-processed MiSeq (2x300bp) multiplexed samples.

We have used Bowtie2 for alignment and SeqMonk for visualising the apparent insert lengths.

For unknown reasons we observe aberrant peaks at around 300bp, and another at around 550 bp, the latter of which disappears if we alter the insert length cutoff (-X option in bowtie) to 1000 (as was suggested elsewhere). The default being 500bp.

However, the peak at approx 300 bp persists. Any idea what this is?

In the following image, the left image is mapped with default Bowtie parameters. The right image with -X 1000.

https://www.dropbox.com/s/brwnp2onpfwkafl/SAM500vs1000.png?dl=0

We suspect it to be some sort of aberration during mapping due to the 300bp read length (possibly from the heavily overlapping reads?), but we can't really make sense of it.

[mikep]

How much QC have you done on the fragmentation prior to sequencing? It could be the kit, and I say that because a flaky nextera kit on Miseq has turned into a living nightmare for us on a genotyping project.

[M4TTN]

Personally I don't think this has anything explicitly to do with the fragmentation. However good/bad it is, I cannot see how this could create a peak in insert size of almost exactly 300bp.

Nevertheless, to give the details: we are using NexteraXT which only uses the inbuilt library normalisation steps. Starting DNA amount is 1-1.5ng, cut at some point with AMPure and finally with Illumina's NexteraXT library normalisation beads (possibly also just Ampure?) prior to pooling for MiSeq (2x300bp). Multiplexing 16-24 samples.

We see the same peaks in all of the samples we've analysed so far.

Overall genome coverage is OK (averages about 30x, but with some regions of very low coverage, and some very high).

I really think this has something to do with the way Bowtie analyses/reports fragments where the two reads fully overlap (which will be the case for any fragment shorter than 300bp).

Does anyone know how Bowtie deals with such a situation?
i.e. where the start position of Read1 = end position of Read2?

What about where the insert is <300bp long, but each read is actually >300 bp long?

i.e. for inserts that are 290-299 bp long, yet the read read length remains as 300bp in each in the FastQ file due to an inability for Illumina software to remove the 10 bp of putative adapter. Would this enable Bowtie to record the insert lengths as 300 (because the read lengths in the FastQ are 300bp long), even thought the actual insert was slightly shorter than 300bp.

Indeed, any information about how Illumina recognises adapter sequences at the ends and trims them when making the FastQ? Presumably it does trim, but requires a certain length of adapter sequence before determining that it can be trimmed off?

[Brian Bushnell]

You can get a reference-free insert size distribution histogram with BBMerge, like this:

bbmerge.sh in1=r1.fq in2=r2.fq ihist=ihist.txt loose

That will be based on read overlapping rather than mapping coordinates. It's very fast, and you can cap it at the first 1m pairs with the flag "reads=1000000" to get it to finish in a few seconds, which is plenty for a good histogram.

You can generate insert size histograms based on mapping with BBMap, again with the ihist flag. I have no idea how Bowtie deals with inserts shorter than read length, but both BBMerge and BBMap process them correctly; though with the default parameters and no "local" flag BBMap won't map them when the insert size is under roughly 75% of the read length. This can be fixed by doing adapter trimming prior to mapping, which I highly recommend for any aligner. Since adapter trimming (for fragment libraries) only removes bases on the right end of reads, it strictly increases the accuracy of insert size calculations, when done correctly.

Incidentally, running BBDuk with the "tbo" (trim by overlap) flag will allow the removal of adapters down to 1bp, because that removes adapters based on apparent insert size as calculated by overlap, rather than matching adapter sequences.

[M4TTN]

Thanks Brian. Please forgive my Unix ignorance, but how would I go about installing/using BBMerge?

[Brian Bushnell]

You just need to unzip and untar it, and then it will work as long as you have Java installed. You can extract it in one command with tar:

tar -zxvf BBMap_34.41.tar.gz

Then just run bbmerge.sh, which will be in the bbmap directory. For example:

/path/to/bbmap/bbmerge.sh in1=r1.fq in2=r2.fq ihist=ihist.txt loose

...where "/path/to/" is the location you extracted BBMap, and your working directory is wherever the reads are located (though you can provide absolute or relative paths to the read files, too).

[M4TTN]

But how does one run the bbmerge.sh file?

If I double-click it opens in Xcode.

Edit: OK. Got it!

[M4TTN]

OK. I am in the bbmap directory, and I can see the bbmerge.sh file.

But when I enter something similar to your code (but with my fastq files replacing r1.fq and r2.fq, I just get:

-bash: bbmerge.sh: command not found

[GenoMax]

You may need to add execute permissions to the shell scripts

$ chmod u+x bbmerge.sh

$ ./bbmerge.sh in1=r1.fq in2=r2.fq ihist=ihist.txt loose

[M4TTN]

ls -l the bbmap directory, gives:


-rwxr-xr-x 1 mn98 staff 5286 23 Jan 18:15 bbmerge.sh

And when trying to run the code, I still get:

-bash: bbmerge.sh: command not found

[M4TTN]

OK. Adding ./ ahead of bbmerge.sh etc

Now gives me this output:


Exception in thread "main" java.lang.NoClassDefFoundError: 2/jni/
Caused by: java.lang.ClassNotFoundException: 2.jni.
at java.net.URLClassLoader$1.run(URLClassLoader.java:202)
at java.security.AccessController.doPrivileged(Native Method)
at java.net.URLClassLoader.findClass(URLClassLoader.java:190)
at java.lang.ClassLoader.loadClass(ClassLoader.java:306)
at sun.misc.Launcher$AppClassLoader.loadClass(Launcher.java:301)
at java.lang.ClassLoader.loadClass(ClassLoader.java:247)

[Unfortunately I have next to no experience with Terminal/Unix...]

[GenoMax]

Just to make sure. You have not moved any of BBMap files/directories after you uncompressed the tar file? Are you running the command from within the BBMap directory?

[M4TTN]

I did not move the files after download.

In terminal I typed:

cd

Then dragged in the bbmap folder and pressed return.

I can check the content of this directory with ls -l

(That is about the limit of my Unix...)

To attempt to execute the command, I type it in as written by brian, but starting directly at "bbmerge.sh" and dragged fast.q files in where the r1.fq and r2.fq were.

I added the reads=1000000 before the hist command.

To get the second (longer) error message I started command with ./bbmerge.sh

[M4TTN]

OK..something happening now...all 8 cores active...

OUTPUT:


BBMerge version 5.5
Finished reading
Total time: 20.117 seconds.

Pairs: 1000000
Joined: 789709 78.971%
Ambiguous: 66486 6.649%
No Solution: 143805 14.381%
Too Short: 0 0.000%
Avg Insert: 279.9
Standard Deviation: 111.5
Mode: 249

Insert range: 36 - 590
90th percentile: 437
75th percentile: 348
50th percentile: 267
25th percentile: 202
10th percentile: 145

[GenoMax]

Success!

Looks like 79% of your reads could be merged.