SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Target enrichment performance pfrommolt Bioinformatics 65 03-17-2016 07:46 AM
Target enrichment cub103 Illumina/Solexa 19 04-11-2011 05:37 AM
Target Capture Technology Expo Illumina/Solexa 0 02-25-2011 10:08 AM
Target re-sequencing Muraya Introductions 0 11-09-2010 01:20 AM
Target Enrichmnet In-Solution mestro2 Sample Prep / Library Generation 10 07-28-2010 08:27 AM

Reply
 
Thread Tools
Old 12-20-2011, 06:48 AM   #1
tonio100680
Member
 
Location: FRANCE / Caen

Join Date: Apr 2010
Posts: 25
Default % On-Off target

hello,
I'm looking for a tool (or command line) to determine the % On-Off target + or - 50 bp of exon from my capture file but not annotated (!!!). Capture SureSelect agilent home.

Current pipeline:
GAIIx Illumina
CASAVA1.8
IGV
CNV-seq
SAMtools
BEDtools
GALAXY
NextGENe

Please HELP
tonio100680 is offline   Reply With Quote
Old 12-20-2011, 06:59 AM   #2
Dameon
Member
 
Location: St. Louis, MO - USA

Join Date: Dec 2011
Posts: 14
Default

Create a bed file of your Agilent SureSelect targets and use BEDtools to merge adjacent targets and then slopBed to add 50 bps to either side of your merged targets. Then use Bedtools BedtoBam to convert your bam file to a bed file and then use intersectBed to create an intersection of your bam.bed and the target.bed. This will create a bed file illustrating the target regions covered by your bam file which you can then parse for percent on and off target. I think there may be examples of this workflow in the BEDtools manual available online.
Dameon is offline   Reply With Quote
Old 12-20-2011, 11:02 PM   #3
laura
Senior Member
 
Location: Cambridge UK

Join Date: Sep 2008
Posts: 151
Default

You may find picards CalculateHsMetrics useful

http://picard.sourceforge.net/comman...ulateHsMetrics
laura is offline   Reply With Quote
Old 12-21-2011, 10:09 AM   #4
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

BEDTools intersectBed will work with bam files, and output .bam files. I use it that way all the time. Your command line look something like:

Quote:
intersectBed -abam yourbam.bam -b paddedExometarget.bed > intersect.bam
I just used Excel to change the coordinates in the target .bed file, to pad them. I'm not sure that's necessry, since intersectBed will get reads that are hanging off of your target regions.

So then use samtools flagstat to count the number of mapped reads of the original .bam, and then of the intersect.bam
swbarnes2 is offline   Reply With Quote
Old 01-20-2012, 12:30 PM   #5
gwilymh
Member
 
Location: Milwaukee

Join Date: Dec 2011
Posts: 72
Default

Quote:
Originally Posted by laura View Post
You may find picards CalculateHsMetrics useful

http://picard.sourceforge.net/comman...ulateHsMetrics
What is the difference between the BAIT_INTERVALS and TARGET_INTERVALS files?
gwilymh is offline   Reply With Quote
Old 01-20-2012, 12:36 PM   #6
laura
Senior Member
 
Location: Cambridge UK

Join Date: Sep 2008
Posts: 151
Default

Quote:
Originally Posted by gwilymh View Post
What is the difference between the BAIT_INTERVALS and TARGET_INTERVALS files?
To be honest we use the same file for both

I suspect that when you have files from your pull down supplier there may be some subtle differences but I don't think it matters to hugely
laura is offline   Reply With Quote
Old 01-20-2012, 07:15 PM   #7
Jon_Keats
Senior Member
 
Location: Phoenix, AZ

Join Date: Mar 2010
Posts: 279
Default

Picard HsMetrics is designed by the Broad, which helped develop the Agilent in solution capture method. In the case of Agilent they provide positions for both the baits and the target regions.

Think:
Code:
       1---------Target------------1
------IIIIIIIIIIII  exon  IIIIIIIIIIIIII--------
     xbaitx         ybaity          zbaitz

Unfortunately, Illumina only provides the target regions not the actual bait locations. So it is harder to decide if some of the non-exonic reads are uncaptured flow-through or captured regions that are not in the "official" targets. Clearly, in my opinion Illumina has captured entire 5kb regions that include 3 exons totaling 1kb resulting in 4kb of high coverage intronic region. Not sure if the designer was just lazy, has ulterior motives, some internal data that this makes target recovery the best?

Last edited by Jon_Keats; 01-20-2012 at 07:18 PM.
Jon_Keats is offline   Reply With Quote
Old 01-23-2012, 09:18 AM   #8
gwilymh
Member
 
Location: Milwaukee

Join Date: Dec 2011
Posts: 72
Default

Quote:
Originally Posted by Jon_Keats View Post
Picard HsMetrics is designed by the Broad, which helped develop the Agilent in solution capture method. In the case of Agilent they provide positions for both the baits and the target regions.

Think:
Code:
       1---------Target------------1
------IIIIIIIIIIII  exon  IIIIIIIIIIIIII--------
     xbaitx         ybaity          zbaitz

Unfortunately, Illumina only provides the target regions not the actual bait locations. So it is harder to decide if some of the non-exonic reads are uncaptured flow-through or captured regions that are not in the "official" targets. Clearly, in my opinion Illumina has captured entire 5kb regions that include 3 exons totaling 1kb resulting in 4kb of high coverage intronic region. Not sure if the designer was just lazy, has ulterior motives, some internal data that this makes target recovery the best?
A colleague of mine directed me to the UCSC Genome Browser Tables page (http://genome.ucsc.edu/cgi-bin/hgTables?org=human), which can be used to search gene names against genome and chromosomes to identify where on a given genome assembly a specific gene is located.
gwilymh is offline   Reply With Quote
Old 01-23-2012, 09:37 AM   #9
gwilymh
Member
 
Location: Milwaukee

Join Date: Dec 2011
Posts: 72
Default

Quote:
Originally Posted by laura View Post
You may find picards CalculateHsMetrics useful

http://picard.sourceforge.net/comman...ulateHsMetrics
Can Picard be used on a Windows system?
gwilymh is offline   Reply With Quote
Old 02-23-2012, 11:36 AM   #10
gwilymh
Member
 
Location: Milwaukee

Join Date: Dec 2011
Posts: 72
Default

Quote:
Originally Posted by swbarnes2 View Post
...So then use samtools flagstat to count the number of mapped reads of the original .bam, and then of the intersect.bam
Does anyone know the specific commands to compile and execute flagstat in samtools? The samtools literature is frustratingly vague!

(I am using samtools in Cygwin in a Windows 7 system)
gwilymh is offline   Reply With Quote
Old 02-23-2012, 12:04 PM   #11
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Quote:
samtools flagstat my_data.bam
You have to get samtools installed. I think you just do a make command. On my install, I think I had a slight problem with the curses file, someone admin at my work suggested I tweak that line in the make file slightly, and it works fine now.l
swbarnes2 is offline   Reply With Quote
Old 02-23-2012, 12:25 PM   #12
ECO
--Site Admin--
 
Location: SF Bay Area, CA, USA

Join Date: Oct 2007
Posts: 1,358
Default

Another vote for Picard's CalculateHsMetrics. It's in the public Galaxy (http://main.g2.bx.psu.edu/ under "NGS: Picard (beta)").
ECO is offline   Reply With Quote
Old 02-23-2012, 01:55 PM   #13
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Picard by default does +/- 250 bp. I recommend using the same file for baits and targets: you can have baits that extend past targets and hence get more coverage for baits than for target if you had all of your target sequence covered by baits. If you had say only 80% of your targets covered by baits, then you already know this, and it just complicates things to try to consider it again. So, again, I recommend using the actual bait intervals if possible.
Heisman is offline   Reply With Quote
Old 07-08-2012, 01:22 PM   #14
gwilymh
Member
 
Location: Milwaukee

Join Date: Dec 2011
Posts: 72
Default

Quote:
Originally Posted by Heisman View Post
Picard by default does +/- 250 bp. I recommend using the same file for baits and targets: you can have baits that extend past targets and hence get more coverage for baits than for target if you had all of your target sequence covered by baits. If you had say only 80% of your targets covered by baits, then you already know this, and it just complicates things to try to consider it again. So, again, I recommend using the actual bait intervals if possible.
I want to verify that picard does indeed use a +/- 250 bp around each target/bait, but have so far not been able to find this written down anywhere. Where did you come by this information?

Also, can the interval be modified (to, say, +/- 300bp)?
gwilymh is offline   Reply With Quote
Old 07-08-2012, 01:38 PM   #15
Heisman
Senior Member
 
Location: St. Louis

Join Date: Dec 2010
Posts: 535
Default

Somewhere there is a web page that has the code for each of the programs... I stumbled on it before and am not really a computer guy so don't know where to find it. But it had 250 set for that metric.
Heisman is offline   Reply With Quote
Old 07-08-2012, 01:46 PM   #16
gwilymh
Member
 
Location: Milwaukee

Join Date: Dec 2011
Posts: 72
Default

Quote:
Originally Posted by Heisman View Post
Somewhere there is a web page that has the code for each of the programs... I stumbled on it before and am not really a computer guy so don't know where to find it. But it had 250 set for that metric.
Thank you for your explanation
gwilymh is offline   Reply With Quote
Old 07-18-2012, 08:47 AM   #17
bwubb
Member
 
Location: Philadelphia

Join Date: Jan 2012
Posts: 58
Default

Does anyone know what format the INTERVALS files need to be? Im using simple bed files but I keep getting this error

Exception in thread "main" java.lang.IllegalStateException: Interval list file must contain header.

I tried adding a Chr Start End header but it doesnt like this either. The simplicity of this is confusing me I guess.

Thanks.
bwubb is offline   Reply With Quote
Old 07-18-2012, 10:14 AM   #18
ECO
--Site Admin--
 
Location: SF Bay Area, CA, USA

Join Date: Oct 2007
Posts: 1,358
Default

That intervals file is annoying to make... here's how I do it (basically you need to add a SAM header and rearrange some columns):

Code:
#Input files to CalculateHsMetrics need SAM header on an interval file ("picard interval file")
#example here: ftp://ftp.broadinstitute.org/pub/gsa/exampleFiles/thousand_genomes_alpha_redesign.targets.interval_list

#put header from bam file at the top of the BI file above "baits.txt"
samtools view -H aligned_reads.bam > header.txt

#interval file needs to look like this:
#1    1104841    1104940    +    target_1
#1    1105283    1105599    +    target_2
#1    1105712    1105860    +    target_3

#rearrange columns of baits bed file, and add SAM header 
awk '{print $1,$2,$3,$6,$4;}' SureSelect_baits.bed > bi.txt
cat header.txt bi.txt > baits.txt
ECO is offline   Reply With Quote
Old 07-18-2012, 10:41 AM   #19
bwubb
Member
 
Location: Philadelphia

Join Date: Jan 2012
Posts: 58
Default

Oh wonderful. Thank you so much.
bwubb is offline   Reply With Quote
Old 08-24-2012, 12:14 PM   #20
bwubb
Member
 
Location: Philadelphia

Join Date: Jan 2012
Posts: 58
Default

Quote:
Originally Posted by ECO View Post
That intervals file is annoying to make... here's how I do it (basically you need to add a SAM header and rearrange some columns):

Code:
#Input files to CalculateHsMetrics need SAM header on an interval file ("picard interval file")
#example here: ftp://ftp.broadinstitute.org/pub/gsa/exampleFiles/thousand_genomes_alpha_redesign.targets.interval_list

#put header from bam file at the top of the BI file above "baits.txt"
samtools view -H aligned_reads.bam > header.txt

#interval file needs to look like this:
#1    1104841    1104940    +    target_1
#1    1105283    1105599    +    target_2
#1    1105712    1105860    +    target_3

#rearrange columns of baits bed file, and add SAM header 
awk '{print $1,$2,$3,$6,$4;}' SureSelect_baits.bed > bi.txt
cat header.txt bi.txt > baits.txt
So this method more or less worked. Picard seems to demand the columns be tab separated so my awk was more like:

awk -F $'\t' 'BEGIN { OFS = FS } {print $1,$2,$3,$6,$4;}' SureSelect_baits.bed > bi.txt

Figured Id share.

I had a question about the header still though, and I expect this is something I just dont understand about the conversion to bam process or something with picard.

The CalculateHSMetrics still yells at me that interval file needs a header.
It seems my aligned_reads.bam files are lacking "@HD VN:1.0 SO:coordinate" at the very top. Is this abnormal?

If I use a bam file that has gone through Picard Read group assignment it does have the @HD etc., but it also will have a @RG line as well.

So do these interval files need to be made for each sample after RG assignment?
bwubb is offline   Reply With Quote
Reply

Tags
% on-off target

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 09:41 AM.


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