This is a test version of Biostars. For the public version, visit https://www.biostars.org.
no header in SAM file

I use the minimap2 to do alignment for my pacbio long-reads DNA data against reference genome. Here is the command I use below:

minimap2 -t 30 -2 -I 100g -ax map-pb genome.fa sample_name.fq > aln.sample_name.sam

to align sample_name.fa to reference genome.fa. And use samtools read the sam file with

samtools view aln.sample_name.sam | head

Here is the content of the first line:

m64291e_211224_041609/0/42782_42995 16 chr6 116495378 18 4S18M1D4M1D8M1I3M2I26M1D54M5I27M2I14M2I3M2I6M1I8M1I3M3I10M6S * 0 0 GTGCGCCTTTTTTATGGTTGGGTTATAGAGGGAGCCTATGAGATTTCTTGGTATTTTAGACTATTGAGTTGCCTTCCACCCTGCAGATCCCTCGCACTCCCCCGTACTGGATGTAACCCCGTGGTAGGACATGATTTTGCCTCCACCAAACTTCTCACTGTTTTCCTGGTGTATTAATCTCAGCCTCCATGGCGTCCGGTTCCAGTATGGTTT !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! NM:i:32 ms:i:226 AS:i:216 nn:i:0 tp:A:P cm:i:5 s1:i:52 s2:i:0 de:f:0.1122 rl:i:0

It seems that there is no header at all.

the raw data sample_name.fq is convert from raw bam file by samtools bam2fq, I wonder maybe some mistakes were taken in the fq in this step?

Has anyone had a similar situation to mine?

minimap samtools pacbio
• 4,668 views • link
updated by GenoMax GenoMax • written by Maxine • Maxine

1 answer

samtools view ... does not output the header by default. Try samtools view -h ....

• 0 views • link
written by Billy Rowell • Billy Rowell

Wow, It works! Thank you! I'm such an idiot.

Here's one more thing I noticed, the headers dominate 749 lines! Is that normal? Here are line 743 to 749 below:

@SQ     SN:original_scaffold_1837_pilon LN:9511
@SQ     SN:original_scaffold_906_pilon  LN:9172
@SQ     SN:original_scaffold_2215_pilon LN:8750
@SQ     SN:original_scaffold_2026_pilon LN:7171
@SQ     SN:original_scaffold_1113_pilon LN:4052
@PG     ID:minimap2     PN:minimap2     VN:2.24-r1122   CL:minimap2 -t 30 -2 -I 5g -ax map-pb /home/maxine91/projects/def-jfu/data/bufo_genome/genome.fa /home/maxine91/projects/def-jfu/data/Z03_fq.dir/cy201704.fq
@PG     ID:samtools     PN:samtools     PP:minimap2     VN:1.13 CL:samtools view -h aln.cy201704.sam
• 0 views • link
updated by GenoMax GenoMax • written by Maxine • Maxine

No problem!

It's hard to know exactly what to expect to see in your header without knowing the entire workflow, but it looks like you aligned against reference with many contigs. Each of those contigs will have an @SQ line in the header. I'm guessing that's the majority of these 749 header lines.

• 0 views • link
written by Billy Rowell • Billy Rowell

yes! the reference genome contains 749 contigs. Thank you so much!

• 0 views • link
written by Maxine • Maxine

Log in to answer this question.

More posts like this