Question: Samtools Error Message, Please Help !
0
gravatar for madkitty
6.1 years ago by
madkitty570
Canada
madkitty570 wrote:

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.

bam samtools sam • 4.3k views
ADD COMMENTlink modified 2.4 years ago by Biostar ♦♦ 20 • written 6.1 years ago by madkitty570

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.

ADD REPLYlink written 6.1 years ago by madkitty570

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 ?

ADD REPLYlink written 6.1 years ago by Leonor Palmeira3.6k

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

ADD REPLYlink written 6.1 years ago by madkitty570

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?

ADD REPLYlink written 6.1 years ago by Vikas Bansal2.3k

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

ADD REPLYlink modified 6.1 years ago • written 6.1 years ago by madkitty570
2
gravatar for Vikas Bansal
6.1 years ago by
Vikas Bansal2.3k
Berlin, Germany
Vikas Bansal2.3k wrote:

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.

ADD COMMENTlink modified 6.1 years ago • written 6.1 years ago by Vikas Bansal2.3k

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 ?

ADD REPLYlink written 6.1 years ago by madkitty570

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.

ADD REPLYlink written 6.1 years ago by Vikas Bansal2.3k

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

ADD REPLYlink modified 6.1 years ago • written 6.1 years ago by madkitty570

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

ADD REPLYlink written 2.4 years ago by H.Hasani590
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1155 users visited in the last hour