SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Convert merged BAM back to per lane BAM or FASTQ file danielsbrewer Bioinformatics 6 10-03-2013 07:29 AM
Editing a Kegg Ortholgy file for metagenomic analysis in Qiime Giorgio C Bioinformatics 0 11-19-2012 11:49 AM
Question about editing a tabbed text file.. shyam_la Bioinformatics 24 06-14-2012 09:06 PM
Editing the content of an ab1 file thedamian Sanger/Dye Terminator 3 04-12-2012 07:14 AM

Reply
 
Thread Tools
Old 01-04-2013, 07:34 AM   #1
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default Editing BAM file

The below is the BAM file produced by running tophat.

samtools view -h accepted_hits.bam
@HD VN:1.0 SO:coordinate
@SQ SN:Mt:3452:Mt3.5.1Chr1 LN:34297684
@SQ SN:Mt:3453:Mt3.5.1Chr2 LN:34376035
@SQ SN:Mt:3454:Mt3.5Chr3 LN:43359096
@SQ SN:Mt:3455:Mt3.5.1Chr4 LN:47685368
...
...
...
@SQ SN:Mt:3460:AC130806.58 LN:141498
@SQ SN:Mt:3461:AC229707.11 LN:3289
@SQ SN:Mt:3462:AC135504.61 LN:132492
...
...
...
@PG ID:TopHat VN:2.0.4 CL:/usr/local/share/tophat-2.0.4.Linux_x86_64/tophat -p 6 -o IR_t1_Mt -G Mt.gff Mt 121108_I7001139_FCD1FBHACXX_L4_GSL_0040.06_1.fq,121108_I7001139_FCD1FBHACXX_L5_GSL_0040.06_1.fq,121108_I7001139_FCD
1FBHACXX_L6_GSL_0040.06_1.fq,121108_I7001139_FCD1FBHACXX_L7_GSL_0040.06_1.fq,121108_I7001139_FCD1FBHACXX_L8_GSL_0040.06_1.fq
FCD1FBHACXX:5:1114:18853:38567#GCCAATA 0 Mt:3452:Mt3.5.1Chr1 150172 50 100M * 0 0 CCCGGGCGCGTATTGCCGAAAAGAAGGTTTAGCTTCAGTTGTGCTACTCTAGACCTGCGGAAGACAGGTGCTGCATGCAAAGGCCTAAGTCCCGCGATGT bbbeeeeegggggi
hhfigiiehdfffaegdfhiihhiiifcgegggeeeeeeddbddaaaYaa`bca__bbcbcbbcccccbaabb_bcdcbaX[aaa[ AS:i:0 X
...
...
...

However, because of too many colons existing between SN and LN, it seems that cufflinks does not work.
So, I think I will have to remove "Mt:3452:" part.

My question is that "Could I use sed command to edit BAM file?"
syintel87 is offline   Reply With Quote
Old 01-04-2013, 08:58 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Not directly. Something like

Code:
samtools view data.bam | sed whatever | samtools view -bSh > corrected.bam
Would work.
swbarnes2 is offline   Reply With Quote
Old 01-04-2013, 09:04 AM   #3
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by swbarnes2 View Post
Not directly. Something like

Code:
samtools view data.bam | sed whatever | samtools view -bSh > corrected.bam
Would work.
Dear swbarnes2,

Oh! Really helpful!! Thank you so much!!!!
syintel87 is offline   Reply With Quote
Old 01-04-2013, 02:31 PM   #4
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by swbarnes2 View Post
Not directly. Something like

Code:
samtools view data.bam | sed whatever | samtools view -bSh > corrected.bam
Would work.

When I put the command below,
samtools view acc_hits.bam | sed 's/Mt\:[0-9][0-9][0-9][0-9]\://' | samtools view -bSh > corrected.bam

the result is the follwing:
Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]

Options: -b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-? longer help


I don't think it works.
Could you see something wrong?
syintel87 is offline   Reply With Quote
Old 01-04-2013, 03:58 PM   #5
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

Quote:
Originally Posted by syintel87 View Post
However, because of too many colons existing between SN and LN, it seems that cufflinks does not work.
So, I think I will have to remove "Mt:3452:" part.
If you are right and editing the reference names like that fixes things, you've found a bug in cufflinks - If so please report it so that cufflinks can be fixed.
maubp is offline   Reply With Quote
Old 01-04-2013, 04:34 PM   #6
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

Quote:
Originally Posted by maubp View Post
If you are right and editing the reference names like that fixes things, you've found a bug in cufflinks - If so please report it so that cufflinks can be fixed.
Yes, while running cufflinks, the error message appeared:
Error: sort order of reads in BAMs must be the same

According to the tips given in another posting, I am trying editing BAM file.
(http://seqanswers.com/forums/showthread.php?t=9304)

Though I sent an email to cufflinks, I haven't received a reply yet.
syintel87 is offline   Reply With Quote
Old 01-07-2013, 03:38 AM   #7
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,060
Default

Quote:
Originally Posted by syintel87 View Post
When I put the command below,
samtools view acc_hits.bam | sed 's/Mt\:[0-9][0-9][0-9][0-9]\://' | samtools view -bSh > corrected.bam


I don't think it works.
Could you see something wrong?
Have you tried doing the three steps independently instead of a single command line?
GenoMax is offline   Reply With Quote
Old 01-07-2013, 05:01 AM   #8
TobyH
Junior Member
 
Location: York

Join Date: Aug 2012
Posts: 5
Default

Quote:
Originally Posted by syintel87 View Post
When I put the command below,
samtools view acc_hits.bam | sed 's/Mt\:[0-9][0-9][0-9][0-9]\://' | samtools view -bSh > corrected.bam

the result is the follwing:
Usage: samtools view [options] <in.bam>|<in.sam> [region1 [...]]

Options: -b output BAM
-h print header for the SAM output
-H print header only (no alignments)
-S input is SAM
-u uncompressed BAM output (force -b)
-x output FLAG in HEX (samtools-C specific)
-X output FLAG in string (samtools-C specific)
-c print only the count of matching records
-t FILE list of reference names and lengths (force -S) [null]
-T FILE reference sequence file (force -S) [null]
-o FILE output file name [stdout]
-R FILE list of read groups to be outputted [null]
-f INT required flag, 0 for unset [0]
-F INT filtering flag, 0 for unset [0]
-q INT minimum mapping quality [0]
-l STR only output reads in library STR [null]
-r STR only output reads in read group STR [null]
-? longer help


I don't think it works.
Could you see something wrong?
If you're still having a problem with this, I believe that you need to specify for samtools view to take input from STDIN at the third step. This is done by a single '-' IIRC. So like this...

samtools view acc_hits.bam | sed 's/Mt\:[0-9][0-9][0-9][0-9]\://' | samtools view -bSh - > corrected.bam
TobyH is offline   Reply With Quote
Old 01-07-2013, 06:38 AM   #9
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

It's likely the problem is due to trying to create a BAM file from input with no header (try a "samtools view -h acc_hits.bam" at the beginning of the command). You can't create a BAM file without a header.
dpryan is offline   Reply With Quote
Old 01-07-2013, 07:30 AM   #10
syintel87
Member
 
Location: Universe

Join Date: Dec 2012
Posts: 81
Default

It seems that errors do not occur by running the command below.

samtools view -h acc_hits.bam | sed 's/Mt\:[0-9][0-9][0-9][0-9]\://' | samtools view -bSh - > corrected.bam

However, seeing "corrected.bam", there was no change.
"sed" command does not seem to be applied to editing bam file.
syintel87 is offline   Reply With Quote
Old 01-07-2013, 09:35 AM   #11
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Yes, I accidentally left out the dash in the second samtools view command. If nothing is changing, there's something wrong with the sed command. In general, piping in and out of samtools is the way to edit .bam files, it's better than making gigantic intermediate sam files that you are going to recompress anyway.
swbarnes2 is offline   Reply With Quote
Old 01-08-2013, 12:08 AM   #12
syfo
Just a member
 
Location: Southern EU

Join Date: Nov 2012
Posts: 103
Default

Quote:
Originally Posted by swbarnes2 View Post
Yes, I accidentally left out the dash in the second samtools view command. If nothing is changing, there's something wrong with the sed command. In general, piping in and out of samtools is the way to edit .bam files, it's better than making gigantic intermediate sam files that you are going to recompress anyway.
Right.

Quote:
Originally Posted by syintel87 View Post
It seems that errors do not occur by running the command below.

samtools view -h acc_hits.bam | sed 's/Mt\:[0-9][0-9][0-9][0-9]\://' | samtools view -bSh - > corrected.bam

However, seeing "corrected.bam", there was no change.
"sed" command does not seem to be applied to editing bam file.
I am surprised that nothing changes. This "sed" command seems correct to me and it works on the text from your first post. Note that this "sed" command after the first pipe does not exactly edit the bam file but the text flow from the "samtools view" command.
syfo 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:38 AM.


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