SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



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

Reply
 
Thread Tools
Old 09-04-2012, 03:24 PM   #21
bwubb
Member
 
Location: Philadelphia

Join Date: Jan 2012
Posts: 55
Default

Ok this is really frustrating. Still is complaining I do not have header in my interval file. Here is a sample of my file...

Code:
@HD	VN:1.0	SO:coordinate
@SQ	SN:1	LN:249250621
@SQ	SN:2	LN:243199373
@SQ	SN:3	LN:198022430
@SQ	SN:4	LN:191154276
@SQ	SN:5	LN:180915260
@SQ	SN:6	LN:171115067
@SQ	SN:7	LN:159138663
@SQ	SN:8	LN:146364022
...

@SQ	SN:GL000192.1	LN:547496
@RG	ID:4346_TTAGGC_ID	PL:illumina	PU:TTAGGC	LB:4346_TTAGGC_LB	SM:4346_TTAGGC_SM
@PG	ID:bwa	PN:bwa	VN:0.5.9-r16
1	45787138	45787258	+	BI426105800_19859
1	45787178	45787298	+	BI426105800_19860
1	45790352	45790472	+	BI426105800_19867
1	45790392	45790512	+	BI426105800_19868
1	45791243	45791363	+	BI426105800_19875
1	45791283	45791403	+	BI426105800_19876
...
Am I doing it wrong? Does the Baits file also need to be in interval format? Thank you.
bwubb is offline   Reply With Quote
Old 09-16-2012, 07:35 AM   #22
pengchy
Senior Member
 
Location: China

Join Date: Feb 2009
Posts: 116
Default

Quote:
Originally Posted by bwubb View Post
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?
I think the “@HD VN:1.0 SO:coordinate” information can be added manually to the top of the header.
But why not "samtools sort" added itself?
pengchy is offline   Reply With Quote
Old 01-16-2013, 10:05 AM   #23
chongm
Member
 
Location: Canada

Join Date: Sep 2012
Posts: 21
Default

Has anyone used hsmetrics for Illumina truseq exome targets? I'm having issues running this (even though I'm using the same modified bed file - with header from samtools - for both target and bait intervals).

Here is what I've done

Quote:
samtools view -h sample.bam > header.txt
awk -F $'\t' 'BEGIN { OFS = FS } {print $1,$2,$3,$6,$4;}' truseq.bed > intervals
cat header.txt intervals > truseq_intervals
Then I ran the hsmetrics
Quote:
java -Xmx4g -Djava.io.tmpdir=/tmp/ \
-jar CalculateHsMetrics.jar \
INPUT=sample.bam \
OUTPUT=sample.hsmetrics \
TARGET_INTERVALS=truseq_intervals\
BAIT_INTERVALS=truseq_intervals
I run into the following error:
Quote:
Exception in thread "main" net.sf.picard.PicardException: Invalid interval record contains 22 fields: MISEQ:23:000000000-A34AH:1:1111:11133:10035 99 chrM 339 60 151M = 469 283 ACACATCTCTGCCAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAATCTCATCA AAAAABBBDDDDDDEEGGGGGGIIIHIIIIIIIIHIIIIIIHIIIIIIIIIIIIIIIIIIHHIIIIIIIIIIIHHHHEEHIIHHHHHHHHHHHHHHHHFHHGGGGHGGGGGGGGGGGGGGEGGGGGGGGGGGGGEGGGGGGEGGGGGGGGG X0:i:1 X1:i:0 MD:Z:71A79 RG:Z:Exome_2011-002 XG:i:0 AM:i:37 NM:i:1 SM:i:37 XM:i:1XO:i:0 XT:A:U
at net.sf.picard.util.IntervalList.fromReader(IntervalList.java:209)
at net.sf.picard.util.IntervalList.fromFile(IntervalList.java:169)
at net.sf.picard.analysis.directed.CollectTargetedMetrics.doWork(CollectTargetedMetrics.java:99)
at net.sf.picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:177)
at net.sf.picard.analysis.directed.CalculateHsMetrics.main(CalculateHsMetrics.java:73)
Can anyone help with this? I'm new to picard.
chongm is offline   Reply With Quote
Old 01-16-2013, 10:20 AM   #24
bwubb
Member
 
Location: Philadelphia

Join Date: Jan 2012
Posts: 55
Default

Can you show a bit of your intervals file?

Some of the header and some of the actual intervals.

Not sure if that is the issue, but its where I usually start.
bwubb is offline   Reply With Quote
Old 01-16-2013, 12:02 PM   #25
Jon_Keats
Senior Member
 
Location: Phoenix, AZ

Join Date: Mar 2010
Posts: 279
Default

Big H not little

samtools view -H

Not

samtools view -h
Jon_Keats is offline   Reply With Quote
Old 01-16-2013, 02:17 PM   #26
chongm
Member
 
Location: Canada

Join Date: Sep 2012
Posts: 21
Default Thanks!

Thank you Jon_Keats! The -H seemed to have done the trick!

chongm
chongm is offline   Reply With Quote
Old 01-17-2013, 03:13 PM   #27
chongm
Member
 
Location: Canada

Join Date: Sep 2012
Posts: 21
Default % reads on target does not equal % bp on target

Quote:
Originally Posted by ECO View Post
Another vote for Picard's CalculateHsMetrics. It's in the public Galaxy (http://main.g2.bx.psu.edu/ under "NGS: Picard (beta)").
I just wanted to note that Picard's hsmetrics will calculate the percentage of basepairs that map to target regions which isn't the same as the percentage of reads that map to target regions. e.g. if a read sits right on a target region, (correct me if i'm wrong), then picard's hsmetrics (PCT_USABLE_BASES_ON_BAIT) will count 50 bp out of 100 as mapping to the target region, whereas if you are taking % reads mapping to target region, you will count that 1 read as mapping to the target region.
chongm is offline   Reply With Quote
Old 04-08-2013, 01:43 PM   #28
rogic
Junior Member
 
Location: Vancouver BC

Join Date: Feb 2012
Posts: 1
Default

Hello all,

one follow-up question regarding the baits file - I downloaded the file from Agilent earray website (SureSelect Human All Exon 50Mb), however it does not contain the strand info which seem to be needed by Picard's CalculateHsMetrics tool. Did the bait files that other people used had the strand info? If yes, where did they get them? If not, how did they circumvent Picard's requirement? Is the strand info actually used by the tool?

Thanks,
Sanja
rogic 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:10 PM.


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