compression level in samtools and gzip

samtools fastq can convert bam to fastq format, e.g. samtools fastq input.bam -o output.fastq

The output file will be automatically compressed if the file names have a .gz or .bgzf extension, e.g. 

samtools fastq input.bam -o output1.fastq.gz

Alternatively, you can also pipe the stdout to compressor explicitly, e.g.

samtools fastq input.bam | gzip > output2.fastq.gz

Interestingly, I noticed that output2.fastq.gz is significantly smaller than output1.fastq.gz, even though the uncompressed file content is the same. 

Actually, this is because of the different default compression ratios used in samtools and gzip. 

In samtools fastq, its default compression level is 1 (out of [0..9]) while gzip's default compression level is 6 (out of [1..9]). 

Explaination of "samtools flagstat" output

==> CRR119891.sam.flagstat <==

192013910 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 secondary

38830640 + 0 supplementary

0 + 0 duplicates

191681231 + 0 mapped (99.83% : N/A)

153183270 + 0 paired in sequencing

76591635 + 0 read1

76591635 + 0 read2

118230182 + 0 properly paired (77.18% : N/A)

152555312 + 0 with itself and mate mapped

295279 + 0 singletons (0.19% : N/A)

10160332 + 0 with mate mapped to a different chr

5532624 + 0 with mate mapped to a different chr (mapQ>=5)

Here are the notes from our careful check:

1. What's secondary, supplementary, and singleton? I was using "bwa mem" with default parameters. In that case, it only reports one alignment if the read can be mapped to multiple locations with equal quality. See discussion here: https://www.biostars.org/p/304614/#304620. That's why you don't see the secondary alignment in above example. Supplementary alignments are those alignments that are part of a chimeric alignment. For example, in a circRNA case, one read is broke into two pieces, mapped to different locations in back-splicing order. Then one of the two fragment will be supplementary (the other is not). Singleton is a mapped read whose mate is unmapped. So its flag has 0x1 and 0x8 but 0x4 is unset. Note that a singleton's mate (which is unmapped) can have a flag of 0110.  

2. What's the total number of reads? The number in the output is NOT for reads, but rather for alignment in the SAM files. The samtools flagstat only check the FLAG, not the read ID. So, in the above example, the total number of reads (R1+R2) should be 192013910 - 38830640. (If there is secondary or duplicate, it should also be excluded). It's also indicated by "paired in sequencing", 153183270 (R1 +R2). This can also be calculated by checking the read ID column e.g. cat CRR119891.sam | cut -f1 | sort -u | wc -l, then times 2 for paired-end reads.

3. How many reads are unmapped? The total alignment is sum of "mapped" + "singletons" + unmapped. So, the number of unmapped reads should be calcualted as 192013910-191681231-295279 = 37400 in above case. The unmapped reads can also be counted by checking the RNAME column in the SAM file (e.g. cat CRR119891.sam | awk '$3=="*"' | wc -l). 

4. Mappability? The above mapped % is based on "mapped" / "total" alignment, not based on reads. If for reads, it should be calculated based on above two numbers. 

samtools flagstat won't output the number of reads mapped

I have been using samtools flagstat to get statistics summary of my BAM/SAM file, like many tutorial suggested (e.g. Dave Tang's note).

But then I noticed that flagstat reports the number of mapping (or alignments/hits, whatever you like to call), but not the number of reads mapped. This can be very different if your alignment contains multiple mappers (and this is especially true if you use default setting or without setting "-g 1" in Tophat). To get the number of reads mapped (e.g. in order to get mappability), you can use the following command:

samtools view -cF 0x100 accepted_hits.bam

Note: Even with multiple mappers, only one hit has flag of 0x100 in the output SAM file. I've already explained this in one of my previous posts.

samtools tview only show alignments with non-* QUAL field

If your input file is fasta (or fastq without quality score), then your alignment from Tophat or Bowtie2 will generate SAM like:


$ samtools view alignment.bam  | fgrep -wi T0730373
t0730373 256 scaffold00006 1146736 0 27M * 0 0 ACTGGGAACACCGCGTGATGTTGGCTA * AS:i:-12 XS:i:-12 XN:i:0 XM:i:2 XO:i:0XG:i:0 NM:i:2 MD:Z:0C25T0 YT:Z:UU
t0730373 272 scaffold00006 1351262 0 27M * 0 0 TAGCCAACATCACGCGGTGTTCCCAGT * AS:i:-12 XS:i:-12 XN:i:0 XM:i:2 XO:i:0XG:i:0 NM:i:2 MD:Z:0A25G0 YT:Z:UU
t0730373 0 scaffold00006 1460535 0 27M * 0 0 ACTGGGAACACCGCGTGATGTTGGCTA IIIIIIIIIIIIIIIIIIIIIIIIIII AS:i:-12 XS:i:-12 XN:i:0 XM:i:2 XO:i:0 XG:i:0 NM:i:2 MD:Z:0C25T0 YT:Z:UU
t0730373 256 scaffold00144 669608 0 27M * 0 0 ACTGGGAACACCGCGTGATGTTGGCTA * AS:i:-12 XS:i:-12 XN:i:0 XM:i:2 XO:i:0XG:i:0 NM:i:2 MD:Z:0C25T0 YT:Z:UU
t0730373 256 scaffold00550 85141 0 27M * 0 0 ACTGGGAACACCGCGTGATGTTGGCTA * AS:i:-12 XS:i:-12 XN:i:0 XM:i:2 XO:i:0XG:i:0 NM:i:2 MD:Z:0C25T0 YT:Z:UU

You see that only the first alignment has the QUAL field, the other secondary alignments have a "*" in the field. And the samtools tview cannot display those alignment.  That's the tip!

sorting BAM makes the file smaller

Yes, it's not just me to notice this change. See here: http://seqanswers.com/forums/archive/index.php/t-13652.html

As Heng said, "BAM is compressed. Sorting helps to give a better compression ratio because similar sequences are grouped together."... So it's not because of the removal of the unmapped reads (which are put at the end).

The tips is - always sort the output BAM after converting a SAM, e.g.

samtools view -Sbu in.sam | samtools sort - in.sorted

mv in.sorted.bam in.bam

Sorted BAM is smaller and better for searching. 

---------------------------

Other tips are:

1. SAM->BAM does not require a sorted header, nor a header. 

If there is header, samtools view -Sb in.sam > in.bam

if there is no header, samtools view -Sbt genome.fa.fai in.sam > in.bam

But the in.bam will not follow the order in the sequence dictionary (genome.fa.fai) unless you sort it by samtools sort. 

2. If input BAM, cufflinks require a proper header for the BAM file, esp. the line of 

@HD VN:1.0 SO:coordinate

Without the line, even if your BAM (or SAM) is sorted, but cufflinks cannot tell it by the file, only if you provide the info through the @HD line. So, I guess 

@HD VN:1.0 SO:unsorted 

won't work.