SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
breakdancer bam2cfg error mei404 Bioinformatics 3 02-08-2017 02:36 AM
Misunderstanding breakdancer output ralonso Bioinformatics 3 04-11-2012 08:02 AM
breakdancer empty cfg (no output from bam2cfg) - Problem Solved jtjli Bioinformatics 0 09-13-2011 07:07 PM
breakdancer output question csoong Bioinformatics 6 12-23-2010 03:46 AM
how to make BreakDancer's bam2cfg.pl work? cliff Bioinformatics 5 09-07-2010 09:15 AM

Reply
 
Thread Tools
Old 12-14-2011, 01:42 AM   #1
zhongj
Member
 
Location: Beijing

Join Date: Dec 2011
Posts: 13
Unhappy BreakDancer empty cfg (no output from bam2cfg)

Hi, everyone, I had spent days to do trouble-shooting on bam2cfg.pl,
breakdancer-1.1_2011_02_21.zip. I had tried some ways to figure out the problem, but without any success. Owing to only one lib I have, I changed the header code of bam2cfg.pl,

Edited code:

open(BAM,"samtools view -h $fbam |") || die "unable to open $fbam\n";
while(<BAM>){
chomp;
if(/^\@PG/){ #getting RG=>LIB mapping from the bam header
my ($id)="bwa";
my ($lib)="1.2k";
my ($platform)="Illumina";
my ($sample)="88LN";
my ($insertsize)="1200";
#if(defined $insertsize && $insertsize>0){
#$lib=$sample . '_'. $lib;
$libs{$lib}=1;
$RGlib{$id}=$lib;
$RGplatform{$id}=$platform;
#}
}

this is the sorted bam file I have:

@SQ SN:scaffold00001 LN:10500
@SQ SN:scaffold00002 LN:2281
@SQ SN:scaffold00003 LN:27085
@SQ SN:scaffold00004 LN:12161
@SQ SN:scaffold00005 LN:2206
.
.
.
@PG ID:bwa PN:bwa VN:0.5.9-r16
.
.
.
HWI-ST833:6:4:18265:54791#0 145 scaffold00001 207 0 50M scaffold00079 243097 0 TTTACTAAAACCGATTGGNCCCGGACAATATTTCGATGTGGGCCGGCCCT ggggggfggg\[]]][K]B[e`gggggfggWgggggegcggggggggggg XT:A:R NM:i:1 SM:i:0 AM:i:0 X0:i:2 X1:i:2 XM:i:1 XO:i:0 XG:i:0 MD:Z:18T31 XA:Z:scaffold00100,+12080,50M,1;scaffold00042,-1233,50M,2;scaffold00007,-133793,50M,2;
.
.
.

"perl bam2cfg.pl *bam > jun.cfg"

At last: samtaols and perl run successfully, but no output from bam2cfg .

Could anyone help me?

Thanks in advance!
Best
Jun

Last edited by zhongj; 12-14-2011 at 01:50 AM.
zhongj is offline   Reply With Quote
Old 12-14-2011, 10:41 AM   #2
ddgenome
Junior Member
 
Location: St. Louis, MO

Join Date: Apr 2009
Posts: 7
Default

Your code does not seem quite right. Are you sure the conditional tests on @PG (program)? The standard version of bam2cfg.pl tests on @RG (read group). Does your BAM header contain a @PG entry? If not, your changes are not getting executed because they are inside an if block that never evaluates to true. Can you just add an RG group to your BAM using samtools reheader? Otherwise, you should probably just set the $libs{$lib}, $RGlib{$id}, and $RGplatform{$id} before you open the BAM.
__________________
David Dooling
ddgenome is offline   Reply With Quote
Old 12-14-2011, 05:45 PM   #3
zhongj
Member
 
Location: Beijing

Join Date: Dec 2011
Posts: 13
Default

Hi, ddgenome, thanks for your suggestion, I will have a try.
zhongj is offline   Reply With Quote
Old 12-14-2011, 10:29 PM   #4
zhongj
Member
 
Location: Beijing

Join Date: Dec 2011
Posts: 13
Default

It still do not work. My Bam header do contain a @PG entry. Any one could help me? BreakDancer is too annoying...

Quote:
Originally Posted by ddgenome View Post
Your code does not seem quite right. Are you sure the conditional tests on @PG (program)? The standard version of bam2cfg.pl tests on @RG (read group). Does your BAM header contain a @PG entry? If not, your changes are not getting executed because they are inside an if block that never evaluates to true. Can you just add an RG group to your BAM using samtools reheader? Otherwise, you should probably just set the $libs{$lib}, $RGlib{$id}, and $RGplatform{$id} before you open the BAM.
zhongj is offline   Reply With Quote
Old 12-15-2011, 08:24 AM   #5
ddgenome
Junior Member
 
Location: St. Louis, MO

Join Date: Apr 2009
Posts: 7
Default

It really seems your BAM is not well-formed (at least for sequence analysis). Looking closer at the alignment record from the example BAM you provided, the read does not have an RG tag. So you are creating an RG record, but then no reads are associated with it.
__________________
David Dooling
ddgenome is offline   Reply With Quote
Old 12-16-2011, 12:43 AM   #6
zhongj
Member
 
Location: Beijing

Join Date: Dec 2011
Posts: 13
Default

Yes, you are right! But how can I fix my BAM file to adapt bam2cfg.pl? Manually add "@RG" ahead every alignment line? Like the following

@R 1.2k GHWI-ST833:6:4:18265:54791#0 145 scaffold00001 207 0 50M scaffold00079 243097 0 TTTACTAAAACCGATTGGNCCCGGACAATATTTCGATGTGGGCCGGCCCT ggggggfggg\[]]][K]B[e`gggggfggWgggggegcggggggggggg XT:A:R NM:i:1 SM:i:0 AM:i:0 X0:i:2 X1:i:2 XM:i:1 XO:i:0 XG:i:0 MD:Z:18T31 XA:Z:scaffold00100,+12080,50M,1;scaffold00042,-1233,50M,2;scaffold00007,-133793,50M,2;


or How to fix bam2cfg.pl to adapt my BAM file?

Looking forward for your reply!
Thanks.

Quote:
Originally Posted by ddgenome View Post
It really seems your BAM is not well-formed (at least for sequence analysis). Looking closer at the alignment record from the example BAM you provided, the read does not have an RG tag. So you are creating an RG record, but then no reads are associated with it.

Last edited by zhongj; 12-16-2011 at 12:52 AM.
zhongj is offline   Reply With Quote
Old 12-19-2011, 01:00 PM   #7
ddgenome
Junior Member
 
Location: St. Louis, MO

Join Date: Apr 2009
Posts: 7
Default

You need to add the @RG entry to the header and then add an RG field pointing to that @RG entry for each alignment record. See the SAM/BAM specification for more information.

http://samtools.sourceforge.net/SAM1.pdf
__________________
David Dooling
ddgenome is offline   Reply With Quote
Old 01-25-2012, 06:30 AM   #8
Robby
Member
 
Location: Germany

Join Date: Mar 2011
Posts: 68
Default

Hello,

did you solve the problem? If yes, could you share the solution, please? I have no outpufile as well and I don't know, what is wrong with my bam file.

My header looks like the following lines:
Quote:
@HD VN:1.0 SO:coordinate
@SQ SN:chr1 LN:249250621
....
@SQ SN:chr22 LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566
@SQ SN:chrM LN:16571
@RG ID:sample CN:xxx PL:ILLUMINA SM:sample LB:sample PI:50
@PG ID:bwa PN:bwa VN:0.5.9-r16

My reads look like
Quote:
readid1 1123 chrM 1 60 101M = 160 260 GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTGG @@@FDFADFADFD>FH@FHEIIIIIIIFGGGGIIC@6BGHGCEEHIIIGIIH(B7@CHGGCCEE<CE/909)2:4>8<>B9>@>3@@4>@BB>@?9<9@ RG:Z:sample XT:A:U NM:i:0 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:101
HWI-ST758_0058:1:1104:17925:175341#CTTGTA 1187
chrM 1 60 101M = 237 337 GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTGG @@??DF?D;=<DBG<?F9E:AFG>HF99E9CCF=@6)?0D>D9??BDGDFFB<F8=CGCFGA@HIH3=:/5;BB<A4@8&055933>@>+39&059@@055 RG:Z:sample XT:A:U NM:i:0 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:101
HWI-ST758_0058:1:1105:15315:81776#CTTGTA 1123
The output is of bam2cfg_2.pl is:
Quote:
Processing bam:test.bam
Closing BAM file
Send TERM signal for 5164
samtools pid process 5164 is still there...
invoking kill -9 on 5164 ...
Closing samtools process : 5164
Could anyone help me?
Best regards
Robby
Robby is offline   Reply With Quote
Old 02-06-2012, 12:56 AM   #9
zhongj
Member
 
Location: Beijing

Join Date: Dec 2011
Posts: 13
Default

sorry,I am still working on it
zhongj is offline   Reply With Quote
Old 02-07-2012, 05:48 AM   #10
Robby
Member
 
Location: Germany

Join Date: Mar 2011
Posts: 68
Default

Dear zhongj,
thanks for your answer. My problem is already solved.

I used the following command
Quote:
perl bam2cfg.pl -g -h sample.bam > sample.cfg
I received the message mentioned before, but that is not an error message. For the first computations it is enough to use a specified amount of reads. The output of that script should look like:

Quote:
readgroup:xxx platform:ILLUMINA map:/path/to/bam/sample.bam readlen:xxx lib:xxx num:xxx lower:xxx upper:xxx mean:xxx std:xxx SWnormality:-xxx flag:0(xx.xx%)1(x.xx%)18(xx.xx%)2(x.xx%)32(x.xx%)4(x.xx%)8(x.xx%)30001 exe:samtools view
After that I proceeded with breakdancer-max1.2 and it worked. At least I received an output and try to analyze and understand that at present.

Best regards
Robby
Robby is offline   Reply With Quote
Old 02-22-2012, 12:10 PM   #11
slink
Junior Member
 
Location: Miami

Join Date: Jun 2010
Posts: 6
Default

Hi all,

I'm having a similar problem and I think it might be due to the way I'm originally aligning the file. I'm aligning with bowtie in the SAM format (-S), but the header is @PG and there are no @RG tags for the reads.

How are you aligning your files to get the Read Group Id?

Thanks,
Sara
slink is offline   Reply With Quote
Old 02-22-2012, 01:43 PM   #12
ddgenome
Junior Member
 
Location: St. Louis, MO

Join Date: Apr 2009
Posts: 7
Default

Google is your friend.

http://bowtie-bio.sourceforge.net/manual.shtml#sam
__________________
David Dooling
ddgenome is offline   Reply With Quote
Old 02-22-2012, 11:19 PM   #13
Robby
Member
 
Location: Germany

Join Date: Mar 2011
Posts: 68
Default

Hi slink,

if your mappings are already finished you can add the read group with Picard as well. For breakdancer ist is also enough, if the RG is mentioned in the header. So you can simply modify the header manually.

Best regards
Robby
Robby is offline   Reply With Quote
Old 10-18-2017, 10:12 AM   #14
etwatson
Member
 
Location: Los Angeles, CA

Join Date: Jun 2012
Posts: 18
Default

I'm going to bump this thread to mention to mention that BAMs merged with samtools merge will not produce output unless the -c flag is provided to combine @RG headers with colliding IDs - if appropriate.

Without the -c flag, samtools merge will create additional @RG IDs in the read mapping that may not be in the header, especially if a header is provided with the -h flag.
etwatson is offline   Reply With Quote
Reply

Tags
bam2cfg, breakdancer, structure variation, svs

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 04:54 AM.


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