Samtools Error Message, Please Help !

Question

[madkitty]

I'm having troubles with this simples query since one week already .. it says that no @SQ lines are in the header. I checked both input and output file, both have the header in full from @HQ to @SQ lines.

samtools view -h file1.bam | awk -F'\t' 'BEGIN { dict[65] dict[147]} $2 in dict' | samtools view -h -bS - > ~/header.bam

[samopen] no @SQ lines in the header.
[sam_read1] missing header? Abort!

I tried to add substr($0,1,1) == "@" to grab the header .. same issue :

samtools view -h file1.bam | awk -F'\t' 'substr($0,1,1) == "@" ||  BEGIN { dict[65] dict[147]} $2 in dict' | samtools view -h -bS - > ~/header.bam

awk: substr($0,1,1) == "@" || BEGIN { dict[65] dict[147]} $2 in dict
awk:                          ^ syntax error
[samopen] no @SQ lines in the header.
[sam_read1] missing header? Abort!

Here is the header :

@HD     VN:1.0 GO:none SO:coordinate
@SQ     SN:chr1 LN:197195432
@SQ     SN:chr2 LN:181748087
@SQ     SN:chr3 LN:159599783
@SQ     SN:chr4 LN:155630120
@SQ     SN:chr5 LN:152537259
@SQ     SN:chr6 LN:149517037
@SQ     SN:chr7 LN:152524553
@SQ     SN:chr8 LN:131738871
@SQ     SN:chr9 LN:124076172
@SQ     SN:chr10        LN:129993255
@SQ     SN:chr11        LN:121843856
@SQ     SN:chr12        LN:121257530
@SQ     SN:chr13        LN:120284312
@SQ     SN:chr14        LN:125194864
@SQ     SN:chr15        LN:103494974
@SQ     SN:chr16        LN:98319150
@SQ     SN:chr17        LN:95272651
@SQ     SN:chr18        LN:90772031
@SQ     SN:chr19        LN:61342430
@SQ     SN:chrX LN:166650296
@SQ     SN:chrY LN:15902555
@SQ     SN:chrM LN:16299
@PG     ID:bwa  PN:bwa  VN:0.5.9-r16

I'm clueless .. if you have another way to grab multiple values in column 2 please let me know (here only value 65 and 177 are represented, I actually need to grab about 10 different ones.

[madkitty]

Please note that this question is the main problem I'm encountering right now. I posted an sub-related questions yesterday that has been partially solved by adding the first missing line of the header @HQ. I think the main problem comes from the Syntax of my query, especially the BEGIN .. thinggy part.

[Leonor Palmeira]

Does replacing the quotes around the "\t" change anything : samtools view -h file1.bam | awk -F"\t" 'BEGIN { dict[65] dict[147]} $2 in dict' | samtools view -h -bS - > ~/header.bam ?

[madkitty]

no .. that doesn't do anything unfortunately..

[Vikas Bansal]

Do you want to output only a header from your file or you want to use reheader on some other file? Could you please mention which files do you have and what do you want to do?

[madkitty]

I have whole genome sequencing files from one sample. I want to do a Quality check on the reads and select only well paired reads. (then I will remove reads with lower quality mapping etc..). So I decided to start with a first step and select only reads that have a FLAG (column 2) of certain values : 65, 177 and so on.. I don't know any picard tool that can do that so that's why I'm trying with AWK

Answer 1

[Vikas Bansal]

I think , in your command, after using awk there is no header in pipe output. The output of " awk -F'\t' 'BEGIN { dict[65] dict[147]} $2 in dict' " will be just like a tab separated file. So while you are using "samtools view -h -bS -" , you will get this error.

Just use - (I did not test it)

samtools view -h file1.bam | awk -F'\t' 'BEGIN { dict[65] dict[147]} $2 in dict' > out.sam

and then you can add header first.

If you want to have reads with specific flags (as you mention in comments), then you can just use "-f" in "view" of samtools.

[madkitty]

Thanks for your answer, in your example the output file is a SAM file .. does it make a difference if I output to SAM instead of BAM ?

[Vikas Bansal]

The output in my example is sam file but without header. So if you want to have header again then you have to use reheader. BAM is just binary format of SAM and you can always convert SAM to BAM or vice versa using samtools. Please have a look at samtools.

[madkitty]

Okay now I understand Thanks a bunch :) But the reheader doesn't work for me so .. I will output directly to bam with -bS

[H.Hasani]

When it comes to header and correcting or filling the missing values, Picard suite could become very helpful Picard Standard Options