How to isolate and examine discordant pairs?

[mcrawford]

Hello,
I am examining yeast that have undergone meiotic recombination, using paired-end sequencing on the Illumina MiSeq and aligned with bowtie2.
I want to look in my sequencing data for ectopic recombination, e.g. if one read pair mate aligns to one chromosome and the other mate to another chromosome. I think I can do this by isolating 'discordant' pairs, e.g. in the command line output you get something like:

650 pairs aligned concordantly 0 times; of these:
34 (5.23%) aligned discordantly 1 time

I want to be able to isolate and take a closer look at the discordant pairs, but I can't figure out how they are marked in the SAM file - I have been reading about the CIGAR string and FLAG field, but I don't think any of the values match what I want.

Many thanks for any help with this

[Brian Bushnell]

The FLAG field has bit 0x2 set if reads are properly paired. So, you want to filter for reads that have bit 0x2 unset. For example, with Reformat:

reformat.sh in=mapped.sam out=filtered.sam requiredbits=1 filterbits=14

That will require reads to have a mate (1) and forbid them from being properly paired, unmapped, or having an unmapped mate (2 & 4 & 8 = 14).

You can do the same thing with samtools -f and -F flags.

[mcrawford]

The FLAG field has bit 0x2 set if reads are properly paired. So, you want to filter for reads that have bit 0x2 unset. For example, with Reformat:

reformat.sh in=mapped.sam out=filtered.sam requiredbits=1 filterbits=14

That will require reads to have a mate (1) and forbid them from being properly paired, unmapped, or having an unmapped mate (2 & 4 & 8 = 14).

You can do the same thing with samtools -f and -F flags.

Thank you very much for your reply!

I am now trying to run bbmap reformat on Windows via the command line by entering:

java -ea -Xmx200m C:\bbmap\current\ jgi.ReformatReads in=reads.sam out=discordant.sam requiredbits=1 filterbits=14

However, I get the error: "Could not find or load main class C:\bbmap\current\". From googling the error it seems there is probably something wrong with the filepath but I tried a few different things and it still doesn't work. I'm not very experienced with java!

[GenoMax]

Is your BBMap directory under c:\? If you put it someplace else then replace the c:\ with that path (e.g. c:\download\bbmap\current).

[mcrawford]

Is your BBMap directory under c:\? If you put it someplace else then replace the c:\ with that path (e.g. c:\download\bbmap\current).

Yes, it is under C:\ (I deliberately put it there to make the path nice and short )

[GenoMax]

There should be a space between C:\bbmap\current\ and jgi.ReformatReads. Is it there? Hard to tell from post. Hopefully you did not move any files under bbmap folder.

[mcrawford]

There should be a space between C:\bbmap\current\ and jgi.ReformatReads. Is it there? Hard to tell from post. Hopefully you did not move any files under bbmap folder.

Yes, there is a space. I didn't move any files either, but I might try re-downloading bbmap just in case.

[westerman]

I am not that good at Java either but I do think you need a 'cp' before your path. When I run 'reformat.sh' from a Linux command line it says that it is running the program like so (single line broken down into multiple lines for clarity)

java
-ea
-Xmx200m
-cp /group/bioinfo/apps/apps/BBMap-34.94/current/
jgi.ReformatReads
PhiX_NoIndexfastq.gz
out=rick.tmp

[mcrawford]

I am not that good at Java either but I do think you need a 'cp' before your path. When I run 'reformat.sh' from a Linux command line it says that it is running the program like so (single line broken down into multiple lines for clarity)

java
-ea
-Xmx200m
-cp /group/bioinfo/apps/apps/BBMap-34.94/current/
jgi.ReformatReads
PhiX_NoIndexfastq.gz
out=rick.tmp

Thanks so much! Adding '-cp' works. However... it now can't find my .sam file. Which folder is it looking in? I tried putting the file in the 'bbmap', 'current' and 'jgi' folders. Or do I need the full path in front of each filename as well?

[westerman]

The full path couldn't hurt. Quoting for spaces, etc. if you have any.

[mcrawford]

OK, I added the full path to each filename and it works now! Thanks for all the help in this thread