SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
High performance computing & Tophat & Cufflinks paolo.kunder Bioinformatics 1 01-24-2012 09:02 AM
questions about samtools mpileup & bcftools chenjy Bioinformatics 0 07-26-2011 05:21 AM
Bowtie & Samtools Questions with SOLiD data earisme Bioinformatics 1 09-16-2010 02:34 PM

Reply
 
Thread Tools
Old 08-18-2011, 03:35 AM   #1
tonio100680
Member
 
Location: FRANCE / Caen

Join Date: Apr 2010
Posts: 25
Default samtools & flagstat & awk

Hi,
Here is my command :
$ samtools merge L002_LBCO1.bam L002_LBCO1_chr*.bam

$ samtools view L002_LBCO1.bam | awk '$3== "chr1.fa" && $4>= 45787123 && $4<= 45787316 || $3== "chr1.fa" && $4>= 45790335 && $4<= 45790528' > essai_chr1.bam

$ samtools flagstat L002_LBCO1.bam
3933498 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
3868863 + 0 mapped (98.36%:nan%)
3933498 + 0 paired in sequencing
1966749 + 0 read1
1966749 + 0 read2
3787076 + 0 properly paired (96.28%:nan%)
3804228 + 0 with itself and mate mapped
64635 + 0 singletons (1.64%:nan%)
7262 + 0 with mate mapped to a different chr
7121 + 0 with mate mapped to a different chr (mapQ>=5)

$ samtools flagstat essai_chr1.bam
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[bam_flagstat_core] Truncated file? Continue anyway.
0 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
0 + 0 mapped (nan%:nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (nan%:nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (nan%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

I would compare files L002_LBCO1.bam (merge result) with essai_chr1. bam (awk result). Apparently he does not recognize the file essai_chr1.bam as bam file.

Do you have a solution?

Best regards
tonio100680 is offline   Reply With Quote
Old 08-18-2011, 05:08 AM   #2
chenyao
Member
 
Location: Beijing

Join Date: Jul 2011
Posts: 74
Default

Quote:
Originally Posted by tonio100680 View Post
Hi,
Here is my command :
$ samtools merge L002_LBCO1.bam L002_LBCO1_chr*.bam

$ samtools view L002_LBCO1.bam | awk '$3== "chr1.fa" && $4>= 45787123 && $4<= 45787316 || $3== "chr1.fa" && $4>= 45790335 && $4<= 45790528' > essai_chr1.bam

$ samtools flagstat L002_LBCO1.bam
3933498 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
3868863 + 0 mapped (98.36%:nan%)
3933498 + 0 paired in sequencing
1966749 + 0 read1
1966749 + 0 read2
3787076 + 0 properly paired (96.28%:nan%)
3804228 + 0 with itself and mate mapped
64635 + 0 singletons (1.64%:nan%)
7262 + 0 with mate mapped to a different chr
7121 + 0 with mate mapped to a different chr (mapQ>=5)

$ samtools flagstat essai_chr1.bam
[bam_header_read] EOF marker is absent. The input is probably truncated.
[bam_header_read] invalid BAM binary header (this is not a BAM file).
[bam_flagstat_core] Truncated file? Continue anyway.
0 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
0 + 0 mapped (nan%:nan%)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (nan%:nan%)
0 + 0 with itself and mate mapped
0 + 0 singletons (nan%:nan%)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

I would compare files L002_LBCO1.bam (merge result) with essai_chr1. bam (awk result). Apparently he does not recognize the file essai_chr1.bam as bam file.

Do you have a solution?

Best regards
truly, the awk command can not generate 'bam' file, you can try making 'san' file then convert to bam file by samtools, but I am not sure it can work.
chenyao is offline   Reply With Quote
Old 08-18-2011, 07:43 AM   #3
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,178
Default

Tonio,

As chenyao already correctly stated, the output of your samtools view | awk pipe is plain text. You could save that to essai_chr1.sam and then convert that sam file to bam.

There is an easier way though. You don't need to use awk to filter the alignments you want; samtools can do this for you, and output a bam file directly. You can use region specifiers in your samtools view command. These define regions using chromosome names and coordinates, you can give multiple regions on one command line. They are written in the format
Code:
<chromosome_name>:<start_position>-<end_position>
and should be placed after the input sam/bam file name. Your command pipe above could be rewritten:

Code:
samtools view -bh -o essai_chr1.bam  L002_LBCO1.bam chr1:45787123-45787316 chr1:45790335-45790528
-b tells sammtools to output in bam format
-h will include the header lines in the output file
-o gives the name of the output file
kmcarr is offline   Reply With Quote
Reply

Tags
awk, flagstat, merged bam, samtools

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:56 PM.


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