SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Bam to Fastq on discordant mate pairs Gio79 Bioinformatics 3 11-09-2013 03:17 AM
FPKM and discordant pairs westeros Bioinformatics 0 09-26-2013 05:14 AM
How to examine paired reads ramiro2k Bioinformatics 3 10-25-2012 07:13 AM
Using DESeq to examine control variability z1dane Bioinformatics 0 08-10-2010 12:30 AM

Reply
 
Thread Tools
Old 06-11-2015, 05:00 AM   #1
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default How to isolate and examine discordant pairs?

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
mcrawford is offline   Reply With Quote
Old 06-11-2015, 09:09 AM   #2
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,707
Default

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.
Brian Bushnell is offline   Reply With Quote
Old 06-12-2015, 02:05 AM   #3
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default

Quote:
Originally Posted by Brian Bushnell View Post
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!
mcrawford is offline   Reply With Quote
Old 06-12-2015, 02:38 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,978
Default

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).
GenoMax is offline   Reply With Quote
Old 06-12-2015, 02:48 AM   #5
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default

Quote:
Originally Posted by GenoMax View Post
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 )
mcrawford is offline   Reply With Quote
Old 06-12-2015, 02:55 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,978
Default

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.
GenoMax is offline   Reply With Quote
Old 06-12-2015, 03:36 AM   #7
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default

Quote:
Originally Posted by GenoMax View Post
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.
mcrawford is offline   Reply With Quote
Old 06-12-2015, 04:47 AM   #8
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

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
westerman is offline   Reply With Quote
Old 06-12-2015, 04:55 AM   #9
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default

Quote:
Originally Posted by westerman View Post
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?
mcrawford is offline   Reply With Quote
Old 06-12-2015, 04:59 AM   #10
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

The full path couldn't hurt. Quoting for spaces, etc. if you have any.
westerman is offline   Reply With Quote
Old 06-12-2015, 05:05 AM   #11
mcrawford
Member
 
Location: Brighton

Join Date: Jun 2015
Posts: 10
Default

OK, I added the full path to each filename and it works now! Thanks for all the help in this thread
mcrawford is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 05:34 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO