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.
1 answer
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.
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.
When it comes to header and correcting or filling the missing values, Picard suite could become very helpful Picard Standard Options
Log in to answer this question.
-
Having trouble with samtools cat
written by matthew.dambrosio.1 •
-
Reformat sam header issue
written by pg_canada •
-
GNU parallel on multiple files to re-header bam files
written by pierre.justeau •
-
Convert an altered sam file into a bam file
written by Jautis
-
StringTie "missing header error" [for both HISAT2 and STAR alignment files]
written by EagleEye
-
converting sam to bam format samtools
written by Hans
-
Sam To Bam Conversion
written by Sajan Franco •
-
Breakdancer Options, Unconsistency Please Help!
written by madkitty
-
Samtools Error Message : No @Sq Lines In The Header
written by madkitty
-
Bam Files And Reads Quality Check
written by 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.
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 ?
no .. that doesn't do anything unfortunately..
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?
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