% On-Off target

[tonio100680]

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

[Dameon]

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.

[laura]

You may find picards CalculateHsMetrics useful

http://picard.sourceforge.net/comman...ulateHsMetrics

[swbarnes2]

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:

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

[gwilymh]

Originally Posted by laura

You may find picards CalculateHsMetrics useful

http://picard.sourceforge.net/comman...ulateHsMetrics

What is the difference between the BAIT_INTERVALS and TARGET_INTERVALS files?

[laura]

Originally Posted by gwilymh

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

[Jon_Keats]

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:

       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?

[gwilymh]

Originally Posted by Jon_Keats

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:

       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]

Originally Posted by laura

You may find picards CalculateHsMetrics useful

http://picard.sourceforge.net/comman...ulateHsMetrics

Can Picard be used on a Windows system?

[gwilymh]

Originally Posted by swbarnes2

...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)

[swbarnes2]

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

[ECO]

Another vote for Picard's CalculateHsMetrics. It's in the public Galaxy (http://main.g2.bx.psu.edu/ under "NGS: Picard (beta)").

[Heisman]

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.

[gwilymh]

Originally Posted by Heisman

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)?

[Heisman]

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.

[gwilymh]

Originally Posted by Heisman

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

[bwubb]

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.

[ECO]

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):

#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

[bwubb]

Oh wonderful. Thank you so much.

[bwubb]

Originally Posted by ECO

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):

#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?