## -rw-r--r-- 1 root root 321M Mar 25 06:10 eg/ERR188273_chrX.sam
+## -rw-r--r-- 1 root root 321M Mar 25 07:18 eg/ERR188273_chrX.sam
Size of BAM file.
ls -lh eg/ERR188273_chrX.bam
-## -rw-r--r-- 1 root root 67M Mar 25 06:08 eg/ERR188273_chrX.bam
+## -rw-r--r-- 1 root root 67M Mar 25 07:16 eg/ERR188273_chrX.bam
We can use head to view a SAM file.
head eg/ERR188273_chrX.sam
@@ -312,9 +312,9 @@ Converting a BAM file to a CRAM fi
ls -lh eg/ERR188273_chrX.[sbcr]*am
-## -rw-r--r-- 1 root root 67M Mar 25 06:08 eg/ERR188273_chrX.bam
-## -rw-r--r-- 1 root root 40M Mar 25 06:10 eg/ERR188273_chrX.cram
-## -rw-r--r-- 1 root root 321M Mar 25 06:10 eg/ERR188273_chrX.sam
+## -rw-r--r-- 1 root root 67M Mar 25 07:16 eg/ERR188273_chrX.bam
+## -rw-r--r-- 1 root root 40M Mar 25 07:18 eg/ERR188273_chrX.cram
+## -rw-r--r-- 1 root root 321M Mar 25 07:18 eg/ERR188273_chrX.sam
You can use samtools view to view a CRAM file just as you would for a
BAM file.
@@ -343,8 +343,8 @@ Sorting a SAM/BAM file
ls -l eg/ERR188273_chrX.bam
ls -l eg/sorted.bam
-## -rw-r--r-- 1 root root 69983526 Mar 25 06:08 eg/ERR188273_chrX.bam
-## -rw-r--r-- 1 root root 69983599 Mar 25 06:10 eg/sorted.bam
+## -rw-r--r-- 1 root root 69983526 Mar 25 07:16 eg/ERR188273_chrX.bam
+## -rw-r--r-- 1 root root 69983599 Mar 25 07:18 eg/sorted.bam
You should use use additional threads (if they are available) to speed
up sorting; to use four threads, use -@ 4.
@@ -352,18 +352,18 @@ Sorting a SAM/BAM file
time samtools sort eg/ERR188273_chrX.sam -o eg/sorted.bam
##
-## real 0m8.762s
-## user 0m8.436s
-## sys 0m0.268s
+## real 0m8.755s
+## user 0m8.460s
+## sys 0m0.248s
Time taken using four threads.
time samtools sort -@ 4 eg/ERR188273_chrX.sam -o eg/sorted.bam
## [bam_sort_core] merging from 0 files and 4 in-memory blocks...
##
-## real 0m2.897s
-## user 0m10.521s
-## sys 0m0.483s
+## real 0m2.905s
+## user 0m10.518s
+## sys 0m0.505s
Many of the SAMtools subtools can use additional threads, so make use of
them if you have the resources!
@@ -736,26 +736,26 @@ Comparing BAM files
samtools mpileup -s -f test_ref.fa aln_bwa.bam aln_mm.bam | head -20
## [mpileup] 2 samples in 2 input files
-## 1 1694 C 1 ^]. D ] 1 ^]. D ]
-## 1 1695 G 1 . I ] 1 . I ]
-## 1 1696 T 1 . J ] 1 . J ]
-## 1 1697 T 1 . J ] 1 . J ]
-## 1 1698 A 1 . J ] 1 . J ]
-## 1 1699 A 1 . J ] 1 . J ]
-## 1 1700 T 1 . J ] 1 . J ]
-## 1 1701 C 1 . J ] 1 . J ]
-## 1 1702 A 1 . J ] 1 . J ]
-## 1 1703 A 1 . J ] 1 . J ]
-## 1 1704 T 1 . J ] 1 . J ]
-## 1 1705 C 1 . J ] 1 . J ]
-## 1 1706 C 1 . J ] 1 . J ]
-## 1 1707 C 1 . J ] 1 . J ]
-## 1 1708 C 1 . J ] 1 . J ]
-## 1 1709 C 1 . J ] 1 . J ]
-## 1 1710 A 1 . J ] 1 . J ]
-## 1 1711 A 1 . J ] 1 . J ]
-## 1 1712 C 1 . J ] 1 . J ]
-## 1 1713 A 1 . J ] 1 . J ]
+## 1 4020 C 1 ^]. E ] 1 ^]. E ]
+## 1 4021 T 1 . J ] 1 . J ]
+## 1 4022 C 1 . J ] 1 . J ]
+## 1 4023 T 1 . J ] 1 . J ]
+## 1 4024 G 1 . J ] 1 . J ]
+## 1 4025 T 1 . J ] 1 . J ]
+## 1 4026 T 1 . J ] 1 . J ]
+## 1 4027 A 1 . J ] 1 . J ]
+## 1 4028 T 1 . J ] 1 . J ]
+## 1 4029 A 1 . J ] 1 . J ]
+## 1 4030 G 1 . J ] 1 . J ]
+## 1 4031 C 1 . J ] 1 . J ]
+## 1 4032 G 1 . J ] 1 . J ]
+## 1 4033 G 1 . J ] 1 . J ]
+## 1 4034 G 1 . J ] 1 . J ]
+## 1 4035 A 1 . J ] 1 . J ]
+## 1 4036 T 1 . J ] 1 . J ]
+## 1 4037 T 1 . J ] 1 . J ]
+## 1 4038 A 1 . J ] 1 . J ]
+## 1 4039 C 1 . J ] 1 . J ]
Another approach is to use
deepTools and the
@@ -809,7 +809,8 @@
Coverage
of bases covered per chromosome/reference sequence.
samtools depth will return three columns: reference, position, and
coverage.
-samtools depth eg/ERR188273_chrX.bam | head
+samtools depth -@ 4 eg/ERR188273_chrX.bam > ERR188273_depth.tsv
+head ERR188273_depth.tsv
## chrX 21649 1
## chrX 21650 1
@@ -823,8 +824,10 @@ Coverage
## chrX 21658 1
The average depth can be calculated by summing the third column and
-dividing by the total number of bases (be sure to use -a).
-samtools depth -@ 2 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf "Bases covered:\t%.3f\nCoverage:\t%.3f\n", $., $cov}'
+dividing by the total number of bases (be sure to use -a with
+samtools depth as that will output all positions including zero
+depth).
+samtools depth -@ 4 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf "Bases covered:\t%.3f\nCoverage:\t%.3f\n", $., $cov}'
## Bases covered: 156040895.000
## Coverage: 0.532
@@ -843,7 +846,8 @@ Coverage
-samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam | head
+samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam > ERR188273_mpileup.tsv
+head ERR188273_mpileup.tsv
## [mpileup] 1 samples in 1 input files
## chrX 251271 g 1 ^]. @ ]
@@ -862,7 +866,8 @@ Coverage
some filtering by default. In the case of this example, read pairs that
are not both mapped will be ignored. To count these “orphan” reads, use
the --count-orphans argument.
-samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam | head
+samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam > ERR188273_mpileup_orphans.tsv
+head ERR188273_mpileup_orphans.tsv
## [mpileup] 1 samples in 1 input files
## chrX 21649 g 0 * * *
@@ -911,6 +916,30 @@ Coverage
average depth of each covered base = meandepth
percentage of bases covered = covbases
+The mosdepth tool can also
+calculate depth (and much faster than samtools depth) per base or
+within a given window. The output is given in a BED file, where the
+fourth column indicates the coverage.
+mosdepth ERR188275 eg/ERR188273_chrX.bam
+gunzip -c ERR188273.per-base.bed.gz | head
+
+## gzip: ERR188273.per-base.bed.gz: No such file or directory
+
+Coverage in using a 500 bp window.
+mosdepth -n --fast-mode --by 500 ERR188275_500 eg/ERR188273_chrX.bam
+gunzip -c ERR188275_500.regions.bed.gz | head
+
+## chrX 0 500 0.00
+## chrX 500 1000 0.00
+## chrX 1000 1500 0.00
+## chrX 1500 2000 0.00
+## chrX 2000 2500 0.00
+## chrX 2500 3000 0.00
+## chrX 3000 3500 0.00
+## chrX 3500 4000 0.00
+## chrX 4000 4500 0.00
+## chrX 4500 5000 0.00
+
Stargazers over time
@@ -956,5 +985,5 @@ Stargazers over time
diff --git a/search/search_index.json b/search/search_index.json
index 7cbc8de..5be7a8b 100644
--- a/search/search_index.json
+++ b/search/search_index.json
@@ -1 +1 @@
-{"config":{"indexing":"full","lang":["en"],"min_search_length":3,"prebuild_index":false,"separator":"[\\s\\-]+"},"docs":[{"location":"","text":"Learning the BAM format Introduction SAMtools provides various (sub)tools for manipulating alignments in the SAM/BAM format. The SAM (Sequence Alignment/Map) format (BAM is just the binary form of SAM) is currently the de facto standard for storing large nucleotide sequence alignments. If you are working with high-throughput sequencing data, at some point you will probably have to deal with SAM/BAM files, so familiarise yourself with them! For the latest information on SAMtools, please refer to the release notes . The examples in this README use the ERR188273_chrX.bam BAM file (stored in the eg folder) generated as per https://github.com/davetang/rnaseq using the HISAT2 + StringTie2 RNA-seq pipeline. This README is generated using the create_readme.sh script; if you want to generate this file yourself, please use this Docker image and the Makefile in this directory. For example: # clone this repo git clone https://github.com/davetang/learning_bam_file.git cd learning_bam_file docker pull davetang/r_build:4.1.2 docker run --rm -it -v $(pwd):/work davetang/r_build:4.1.2 /bin/bash # inside the Docker container make Installing SAMtools For installing SAMtools, I recommend using Conda and the Bioconda samtools package . I also recommend using Miniconda instead of Anaconda because Anaconda comes with a lot of tools/packages that you will probably not use. I wrote a short introduction to Conda if you want to find learn more. Once you have installed Miniconda, you can install SAMtools as follows: conda install -c bioconda samtools Otherwise you can download the source and compile it yourself; change dir to the location you want samtools to be installed. samtools will be installed in ${dir}/bin , so make sure this is in your $PATH . #!/usr/bin/env bash set -euo pipefail ver=1.15 tool=samtools url=https://github.com/samtools/${tool}/releases/download/${ver}/${tool}-${ver}.tar.bz2 dir=${HOME}/local wget ${url} tar xjf ${tool}-${ver}.tar.bz2 cd ${tool}-${ver} ./configure --prefix=${dir} make && make install cd .. rm -rf ${tool}-${ver} ${tool}-${ver}.tar.bz2 >&2 echo Done exit 0 Basic usage If you run samtools on the terminal without any parameters or with --help , all the available utilities are listed: samtools --help ## ## Program: samtools (Tools for alignments in the SAM format) ## Version: 1.19.2 (using htslib 1.19.1) ## ## Usage: samtools [options] ## ## Commands: ## -- Indexing ## dict create a sequence dictionary file ## faidx index/extract FASTA ## fqidx index/extract FASTQ ## index index alignment ## ## -- Editing ## calmd recalculate MD/NM tags and '=' bases ## fixmate fix mate information ## reheader replace BAM header ## targetcut cut fosmid regions (for fosmid pool only) ## addreplacerg adds or replaces RG tags ## markdup mark duplicates ## ampliconclip clip oligos from the end of reads ## ## -- File operations ## collate shuffle and group alignments by name ## cat concatenate BAMs ## consensus produce a consensus Pileup/FASTA/FASTQ ## merge merge sorted alignments ## mpileup multi-way pileup ## sort sort alignment file ## split splits a file by read group ## quickcheck quickly check if SAM/BAM/CRAM file appears intact ## fastq converts a BAM to a FASTQ ## fasta converts a BAM to a FASTA ## import Converts FASTA or FASTQ files to SAM/BAM/CRAM ## reference Generates a reference from aligned data ## reset Reverts aligner changes in reads ## ## -- Statistics ## bedcov read depth per BED region ## coverage alignment depth and percent coverage ## depth compute the depth ## flagstat simple stats ## idxstats BAM index stats ## cram-size list CRAM Content-ID and Data-Series sizes ## phase phase heterozygotes ## stats generate stats (former bamcheck) ## ampliconstats generate amplicon specific stats ## ## -- Viewing ## flags explain BAM flags ## head header viewer ## tview text alignment viewer ## view SAM<->BAM<->CRAM conversion ## depad convert padded BAM to unpadded BAM ## samples list the samples in a set of SAM/BAM/CRAM files ## ## -- Misc ## help [cmd] display this help message or help for [cmd] ## version detailed version information Viewing Use bioSyntax to prettify your output. samtools view aln.bam | sam-less Converting a SAM file to a BAM file A BAM file is just a SAM file but stored in binary format; you should always convert your SAM files into BAM format since they are smaller in size and are faster to manipulate. I don\u2019t have a SAM file in the example folder, so let\u2019s create one and check out the first ten lines. Note: remember to use -h to ensure the SAM file contains the sequence header information. Generally, I recommend storing only sorted BAM files as they use even less disk space and are faster to process. samtools view -h eg/ERR188273_chrX.bam > eg/ERR188273_chrX.sam Notice that the SAM file is much larger than the BAM file. Size of SAM file. ls -lh eg/ERR188273_chrX.sam ## -rw-r--r-- 1 root root 321M Mar 25 06:10 eg/ERR188273_chrX.sam Size of BAM file. ls -lh eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 67M Mar 25 06:08 eg/ERR188273_chrX.bam We can use head to view a SAM file. head eg/ERR188273_chrX.sam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -h eg/ERR188273_chrX.bam ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH aln.bam If the header information is available, we can convert a SAM file into BAM by using samtools view -b . In newer versions of SAMtools, the input format is auto-detected, so we no longer need the -S parameter. samtools view -b eg/ERR188273_chrX.sam > eg/my.bam Converting a BAM file to a CRAM file The CRAM format is even more compact. Use samtools view with the -T and -C arguments to convert a BAM file into CRAM. samtools view -T genome/chrX.fa -C -o eg/ERR188273_chrX.cram eg/ERR188273_chrX.bam ls -lh eg/ERR188273_chrX.[sbcr]*am ## -rw-r--r-- 1 root root 67M Mar 25 06:08 eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 40M Mar 25 06:10 eg/ERR188273_chrX.cram ## -rw-r--r-- 1 root root 321M Mar 25 06:10 eg/ERR188273_chrX.sam You can use samtools view to view a CRAM file just as you would for a BAM file. samtools view eg/ERR188273_chrX.cram | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:10 MD:Z:75 NM:i:0 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD< test_ref.fa script/random_paired_end.pl -f test_ref.fa -l 100 -n 10000 -m 300 -d 0.1 bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l100_n10000_d300_1984_1.fq.gz l100_n10000_d300_1984_2.fq.gz > aln.sam 2> /dev/null samtools flagstat will indicate that some reads (about 10%) mapped to different chromosomes. samtools flagstat aln.sam ## 20000 + 0 in total (QC-passed reads + QC-failed reads) ## 20000 + 0 primary ## 0 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 20000 + 0 mapped (100.00% : N/A) ## 20000 + 0 primary mapped (100.00% : N/A) ## 20000 + 0 paired in sequencing ## 10000 + 0 read1 ## 10000 + 0 read2 ## 18012 + 0 properly paired (90.06% : N/A) ## 20000 + 0 with itself and mate mapped ## 0 + 0 singletons (0.00% : N/A) ## 1988 + 0 with mate mapped to a different chr ## 1988 + 0 with mate mapped to a different chr (mapQ>=5) Flag of a proper pair. samtools flag $(samtools view -f 2 aln.sam | head -1 | cut -f2) ## 0x63 99 PAIRED,PROPER_PAIR,MREVERSE,READ1 Flag of a pair (that is not a proper pair). samtools flag $(samtools view -F 2 aln.sam | head -1 | cut -f2) ## 0x61 97 PAIRED,MREVERSE,READ1 Filtering unmapped reads Use -F 4 to filter out unmapped reads. samtools view -F 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.mapped.bam Use -f 4 to keep only unmapped reads. samtools view -f 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.unmapped.bam We can use the flags subcommand to confirm that a value of four represents an unmapped read. samtools flags 4 ## 0x4 4 UNMAP Extracting entries mapping to a specific loci Use samtools view and the ref:start-end syntax to extract reads mapping within a specific genomic loci; this requires a BAM index file. samtools view eg/ERR188273_chrX.bam chrX:20000-30000 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP Note that this takes into account the mapping of the entire read and not just the starting position. For example, if you specified chrX:20000-30000, a 75 bp long read that starts its mapping from position 19999 will also be returned. In addition, you can save the output as another BAM file if you want. samtools view -b eg/ERR188273_chrX.bam chrX:20000-30000 > eg/ERR188273_chrX_20000_30000.bam If you want reads mapped to a single reference (e.g. chromosome), just specify the ref and leave out the start and end values. samtools view eg/ERR188273_chrX.bam chrX | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YT:Z:UP NH:i:10 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD<DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH eg/first.bam Once again, you can use flags to verify this (it also accepts hexadecimal input). samtools flags 0x0040 ## 0x40 64 READ1 Stats For simple statistics use samtools flagstat . samtools flagstat eg/ERR188273_chrX.bam ## 1176360 + 0 in total (QC-passed reads + QC-failed reads) ## 1160084 + 0 primary ## 16276 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 1126961 + 0 mapped (95.80% : N/A) ## 1110685 + 0 primary mapped (95.74% : N/A) ## 1160084 + 0 paired in sequencing ## 580042 + 0 read1 ## 580042 + 0 read2 ## 1060858 + 0 properly paired (91.45% : N/A) ## 1065618 + 0 with itself and mate mapped ## 45067 + 0 singletons (3.88% : N/A) ## 0 + 0 with mate mapped to a different chr ## 0 + 0 with mate mapped to a different chr (mapQ>=5) For more stats, use samtools stats . samtools stats eg/ERR188273_chrX.bam | grep ^SN ## SN raw total sequences: 1160084 # excluding supplementary and secondary reads ## SN filtered sequences: 0 ## SN sequences: 1160084 ## SN is sorted: 1 ## SN 1st fragments: 580042 ## SN last fragments: 580042 ## SN reads mapped: 1110685 ## SN reads mapped and paired: 1065618 # paired-end technology bit set + both mates mapped ## SN reads unmapped: 49399 ## SN reads properly paired: 1060858 # proper-pair bit set ## SN reads paired: 1160084 # paired-end technology bit set ## SN reads duplicated: 0 # PCR or optical duplicate bit set ## SN reads MQ0: 905 # mapped and MQ=0 ## SN reads QC failed: 0 ## SN non-primary alignments: 16276 ## SN supplementary alignments: 0 ## SN total length: 87006300 # ignores clipping ## SN total first fragment length: 43503150 # ignores clipping ## SN total last fragment length: 43503150 # ignores clipping ## SN bases mapped: 83301375 # ignores clipping ## SN bases mapped (cigar): 83064942 # more accurate ## SN bases trimmed: 0 ## SN bases duplicated: 0 ## SN mismatches: 423271 # from NM fields ## SN error rate: 5.095663e-03 # mismatches / bases mapped (cigar) ## SN average length: 75 ## SN average first fragment length: 75 ## SN average last fragment length: 75 ## SN maximum length: 75 ## SN maximum first fragment length: 75 ## SN maximum last fragment length: 75 ## SN average quality: 36.0 ## SN insert size average: 182.7 ## SN insert size standard deviation: 176.0 ## SN inward oriented pairs: 530763 ## SN outward oriented pairs: 1042 ## SN pairs with other orientation: 1004 ## SN pairs on different chromosomes: 0 ## SN percentage of properly paired reads (%): 91.4 samtools calmd/fillmd The calmd or fillmd tool is useful for visualising mismatches and insertions in an alignment of a read to a reference genome. The -e argument changes identical bases between the read and reference into = . samtools view -b eg/ERR188273_chrX.bam | samtools fillmd -e - genome/chrX.fa > eg/ERR188273_chrX_fillmd.bam head eg/ERR188273_chrX_fillmd.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -b eg/ERR188273_chrX.bam ## @PG ID:samtools.1 PN:samtools PP:samtools VN:1.19.2 CL:samtools fillmd -e - genome/chrX.fa ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 =========================================================================== @@HHEGIIGAGIIIBGIIG@FECHDGHCHHHGHHFFFFFDEACC@ ## @ERR188273.14904746 ## GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT ## + ## @@HHEGIIGAGIIIBGIIG@FECH eg/ERR188273_chrX_rand.bam Count number of reads Use samtools idxstats to print stats on a BAM file; this requires an index file which is created by running samtools index . The output of idxstats is a file with four tab-delimited columns: Reference name Sequence length of reference Number of mapped reads Number of unmapped reads samtools idxstats eg/ERR188273_chrX.bam ## chrX 156040895 1126961 45067 ## * 0 0 4332 We can use this with awk to calculate: The number of mapped reads by summing the third column. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3} END {print s}' ## 1126961 The number of reads, which is the sum of mapped and unmapped reads. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3+$4} END {print s}' ## 1176360 Obtaining genomic sequence Use faidx to fetch genomic sequence; coordinates are 1-based. We need to first index the reference FASTA file that was used to map the reads. samtools faidx genome/chrX.fa Now we can obtain the sequence. samtools faidx genome/chrX.fa chrX:300000-300100 ## >chrX:300000-300100 ## ctgagatcgtgccactgcactccagcctgggcgacagagcgagactccatctcaaaaaaa ## aaaaaaaaaaaaaagaTggggtctctctatgttggccaggt Comparing BAM files The output from mpileup can be used to compare BAM files. The commands below generates alignments using bwa and minimap2 . len=100 n=10000 m=300 script/generate_random_seq.pl 30 1000000 1984 > test_ref.fa script/random_paired_end.pl -f test_ref.fa -l ${len} -n ${n} -m ${m} bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_bwa.bam minimap2 -ax sr test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_mm.bam The BAM files can be used with mpileup to compare the depths. samtools mpileup -s -f test_ref.fa aln_bwa.bam aln_mm.bam | head -20 ## [mpileup] 2 samples in 2 input files ## 1 1694 C 1 ^]. D ] 1 ^]. D ] ## 1 1695 G 1 . I ] 1 . I ] ## 1 1696 T 1 . J ] 1 . J ] ## 1 1697 T 1 . J ] 1 . J ] ## 1 1698 A 1 . J ] 1 . J ] ## 1 1699 A 1 . J ] 1 . J ] ## 1 1700 T 1 . J ] 1 . J ] ## 1 1701 C 1 . J ] 1 . J ] ## 1 1702 A 1 . J ] 1 . J ] ## 1 1703 A 1 . J ] 1 . J ] ## 1 1704 T 1 . J ] 1 . J ] ## 1 1705 C 1 . J ] 1 . J ] ## 1 1706 C 1 . J ] 1 . J ] ## 1 1707 C 1 . J ] 1 . J ] ## 1 1708 C 1 . J ] 1 . J ] ## 1 1709 C 1 . J ] 1 . J ] ## 1 1710 A 1 . J ] 1 . J ] ## 1 1711 A 1 . J ] 1 . J ] ## 1 1712 C 1 . J ] 1 . J ] ## 1 1713 A 1 . J ] 1 . J ] Another approach is to use deepTools and the bamCompare command. The bigWig output file shows the ratio of reads between b1 and b2 in 50 bp (default) windows. Converting reference names One of the most annoying bioinformatics problems is the use of different chromosome names, e.g. chr1 vs 1, in different references even when the sequences are identical. The GRCh38 reference downloaded from Ensembl has chromosome names without the chr : >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF Whereas the reference names from UCSC has the chr : >chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38 Luckily you can change the reference names using samtools reheader but just make sure your reference sequences are actually identical. samtools view eg/ERR188273_chrX.bam | head -2 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP View header samtools view -H eg/ERR188273_chrX.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -H eg/ERR188273_chrX.bam Substitute header with new name. samtools view -H eg/ERR188273_chrX.bam | sed 's/SN:chrX/SN:X/' > eg/my_header Save bam file with new ref and check it out. samtools reheader eg/my_header eg/ERR188273_chrX.bam > eg/ERR188273_X.bam samtools view eg/ERR188273_X.bam | head -2 ## ERR188273.4711308 73 X 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 X 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP Coverage Coverage can mean the: average depth of each covered base percentage of bases covered samtools depth and samtools mpileup can be used to indicate the depth of each covered base (and used to calculate the average depth. samtools coverage will provide both the average depth and percentage of bases covered per chromosome/reference sequence. samtools depth will return three columns: reference, position, and coverage. samtools depth eg/ERR188273_chrX.bam | head ## chrX 21649 1 ## chrX 21650 1 ## chrX 21651 1 ## chrX 21652 1 ## chrX 21653 1 ## chrX 21654 1 ## chrX 21655 1 ## chrX 21656 1 ## chrX 21657 1 ## chrX 21658 1 The average depth can be calculated by summing the third column and dividing by the total number of bases (be sure to use -a ). samtools depth -@ 2 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf \"Bases covered:\\t%.3f\\nCoverage:\\t%.3f\\n\", $., $cov}' ## Bases covered: 156040895.000 ## Coverage: 0.532 The samtools mpileup command also provides depth information (but not for reads that have a mapping quality of 0, by default) with some additional information: Sequence name 1-based coordinate Reference base (when used with -f ) Number of reads covering this position Read bases Base qualities Alignment mapping qualities (when used with -s ) samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam | head ## [mpileup] 1 samples in 1 input files ## chrX 251271 g 1 ^]. @ ] ## chrX 251272 a 1 . @ ] ## chrX 251273 a 1 . < ] ## chrX 251274 a 1 . D ] ## chrX 251275 a 1 . D ] ## chrX 251276 a 1 . D ] ## chrX 251277 t 1 . D ] ## chrX 251278 g 1 . D ] ## chrX 251279 g 1 . F ] ## chrX 251280 g 1 . B ] Note that the start of the samtools mpileup output differ from the start of the samtools depth output. This is because mpileup performs some filtering by default. In the case of this example, read pairs that are not both mapped will be ignored. To count these \u201corphan\u201d reads, use the --count-orphans argument. samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam | head ## [mpileup] 1 samples in 1 input files ## chrX 21649 g 0 * * * ## chrX 21650 a 1 . D ! ## chrX 21651 t 1 . F ! ## chrX 21652 c 1 . F ! ## chrX 21653 a 1 . H ! ## chrX 21654 c 1 . G ! ## chrX 21655 g 1 . H ! ## chrX 21656 a 1 . B ! ## chrX 21657 g 1 . H ! ## chrX 21658 g 1 . I ! In addition mpileup performs \u201cper-Base Alignment Quality\u201d (BAQ) by default and will adjust base quality scores. The default behaviour to to skip bases with baseQ/BAQ smaller than 13. If you are finding discrepancies between mpileup \u2019s coverage calculation with another coverage tool, you can either set --min-BQ to 0 or use --no-BAQ to disable BAQ. I have an old blog post on using mpileup . samtools coverage will provide the following coverage statistics: rname - Reference name / chromosome startpos - Start position endpos - End position (or sequence length) numreads - Number reads aligned to the region (after filtering) covbases - Number of covered bases with depth >= 1 coverage - Proportion of covered bases [0..1] meandepth - Mean depth of coverage meanbaseq - Mean base quality in covered region meanmapq - Mean mapping quality of selected reads samtools coverage eg/ERR188273_chrX.bam ## #rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq ## chrX 1 156040895 1110685 3402037 2.18022 0.532299 36.3 59.4 The example BAM file only contains reads for chrX hence the statistics are only returned for chrX . Returning to our coverage definition at the start of this section: average depth of each covered base = meandepth percentage of bases covered = covbases Stargazers over time","title":"Home"},{"location":"#learning-the-bam-format","text":"","title":"Learning the BAM format"},{"location":"#introduction","text":"SAMtools provides various (sub)tools for manipulating alignments in the SAM/BAM format. The SAM (Sequence Alignment/Map) format (BAM is just the binary form of SAM) is currently the de facto standard for storing large nucleotide sequence alignments. If you are working with high-throughput sequencing data, at some point you will probably have to deal with SAM/BAM files, so familiarise yourself with them! For the latest information on SAMtools, please refer to the release notes . The examples in this README use the ERR188273_chrX.bam BAM file (stored in the eg folder) generated as per https://github.com/davetang/rnaseq using the HISAT2 + StringTie2 RNA-seq pipeline. This README is generated using the create_readme.sh script; if you want to generate this file yourself, please use this Docker image and the Makefile in this directory. For example: # clone this repo git clone https://github.com/davetang/learning_bam_file.git cd learning_bam_file docker pull davetang/r_build:4.1.2 docker run --rm -it -v $(pwd):/work davetang/r_build:4.1.2 /bin/bash # inside the Docker container make","title":"Introduction"},{"location":"#installing-samtools","text":"For installing SAMtools, I recommend using Conda and the Bioconda samtools package . I also recommend using Miniconda instead of Anaconda because Anaconda comes with a lot of tools/packages that you will probably not use. I wrote a short introduction to Conda if you want to find learn more. Once you have installed Miniconda, you can install SAMtools as follows: conda install -c bioconda samtools Otherwise you can download the source and compile it yourself; change dir to the location you want samtools to be installed. samtools will be installed in ${dir}/bin , so make sure this is in your $PATH . #!/usr/bin/env bash set -euo pipefail ver=1.15 tool=samtools url=https://github.com/samtools/${tool}/releases/download/${ver}/${tool}-${ver}.tar.bz2 dir=${HOME}/local wget ${url} tar xjf ${tool}-${ver}.tar.bz2 cd ${tool}-${ver} ./configure --prefix=${dir} make && make install cd .. rm -rf ${tool}-${ver} ${tool}-${ver}.tar.bz2 >&2 echo Done exit 0","title":"Installing SAMtools"},{"location":"#basic-usage","text":"If you run samtools on the terminal without any parameters or with --help , all the available utilities are listed: samtools --help ## ## Program: samtools (Tools for alignments in the SAM format) ## Version: 1.19.2 (using htslib 1.19.1) ## ## Usage: samtools [options] ## ## Commands: ## -- Indexing ## dict create a sequence dictionary file ## faidx index/extract FASTA ## fqidx index/extract FASTQ ## index index alignment ## ## -- Editing ## calmd recalculate MD/NM tags and '=' bases ## fixmate fix mate information ## reheader replace BAM header ## targetcut cut fosmid regions (for fosmid pool only) ## addreplacerg adds or replaces RG tags ## markdup mark duplicates ## ampliconclip clip oligos from the end of reads ## ## -- File operations ## collate shuffle and group alignments by name ## cat concatenate BAMs ## consensus produce a consensus Pileup/FASTA/FASTQ ## merge merge sorted alignments ## mpileup multi-way pileup ## sort sort alignment file ## split splits a file by read group ## quickcheck quickly check if SAM/BAM/CRAM file appears intact ## fastq converts a BAM to a FASTQ ## fasta converts a BAM to a FASTA ## import Converts FASTA or FASTQ files to SAM/BAM/CRAM ## reference Generates a reference from aligned data ## reset Reverts aligner changes in reads ## ## -- Statistics ## bedcov read depth per BED region ## coverage alignment depth and percent coverage ## depth compute the depth ## flagstat simple stats ## idxstats BAM index stats ## cram-size list CRAM Content-ID and Data-Series sizes ## phase phase heterozygotes ## stats generate stats (former bamcheck) ## ampliconstats generate amplicon specific stats ## ## -- Viewing ## flags explain BAM flags ## head header viewer ## tview text alignment viewer ## view SAM<->BAM<->CRAM conversion ## depad convert padded BAM to unpadded BAM ## samples list the samples in a set of SAM/BAM/CRAM files ## ## -- Misc ## help [cmd] display this help message or help for [cmd] ## version detailed version information","title":"Basic usage"},{"location":"#viewing","text":"Use bioSyntax to prettify your output. samtools view aln.bam | sam-less","title":"Viewing"},{"location":"#converting-a-sam-file-to-a-bam-file","text":"A BAM file is just a SAM file but stored in binary format; you should always convert your SAM files into BAM format since they are smaller in size and are faster to manipulate. I don\u2019t have a SAM file in the example folder, so let\u2019s create one and check out the first ten lines. Note: remember to use -h to ensure the SAM file contains the sequence header information. Generally, I recommend storing only sorted BAM files as they use even less disk space and are faster to process. samtools view -h eg/ERR188273_chrX.bam > eg/ERR188273_chrX.sam Notice that the SAM file is much larger than the BAM file. Size of SAM file. ls -lh eg/ERR188273_chrX.sam ## -rw-r--r-- 1 root root 321M Mar 25 06:10 eg/ERR188273_chrX.sam Size of BAM file. ls -lh eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 67M Mar 25 06:08 eg/ERR188273_chrX.bam We can use head to view a SAM file. head eg/ERR188273_chrX.sam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -h eg/ERR188273_chrX.bam ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH aln.bam If the header information is available, we can convert a SAM file into BAM by using samtools view -b . In newer versions of SAMtools, the input format is auto-detected, so we no longer need the -S parameter. samtools view -b eg/ERR188273_chrX.sam > eg/my.bam","title":"Converting a SAM file to a BAM file"},{"location":"#converting-a-bam-file-to-a-cram-file","text":"The CRAM format is even more compact. Use samtools view with the -T and -C arguments to convert a BAM file into CRAM. samtools view -T genome/chrX.fa -C -o eg/ERR188273_chrX.cram eg/ERR188273_chrX.bam ls -lh eg/ERR188273_chrX.[sbcr]*am ## -rw-r--r-- 1 root root 67M Mar 25 06:08 eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 40M Mar 25 06:10 eg/ERR188273_chrX.cram ## -rw-r--r-- 1 root root 321M Mar 25 06:10 eg/ERR188273_chrX.sam You can use samtools view to view a CRAM file just as you would for a BAM file. samtools view eg/ERR188273_chrX.cram | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:10 MD:Z:75 NM:i:0 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD< test_ref.fa script/random_paired_end.pl -f test_ref.fa -l 100 -n 10000 -m 300 -d 0.1 bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l100_n10000_d300_1984_1.fq.gz l100_n10000_d300_1984_2.fq.gz > aln.sam 2> /dev/null samtools flagstat will indicate that some reads (about 10%) mapped to different chromosomes. samtools flagstat aln.sam ## 20000 + 0 in total (QC-passed reads + QC-failed reads) ## 20000 + 0 primary ## 0 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 20000 + 0 mapped (100.00% : N/A) ## 20000 + 0 primary mapped (100.00% : N/A) ## 20000 + 0 paired in sequencing ## 10000 + 0 read1 ## 10000 + 0 read2 ## 18012 + 0 properly paired (90.06% : N/A) ## 20000 + 0 with itself and mate mapped ## 0 + 0 singletons (0.00% : N/A) ## 1988 + 0 with mate mapped to a different chr ## 1988 + 0 with mate mapped to a different chr (mapQ>=5) Flag of a proper pair. samtools flag $(samtools view -f 2 aln.sam | head -1 | cut -f2) ## 0x63 99 PAIRED,PROPER_PAIR,MREVERSE,READ1 Flag of a pair (that is not a proper pair). samtools flag $(samtools view -F 2 aln.sam | head -1 | cut -f2) ## 0x61 97 PAIRED,MREVERSE,READ1","title":"Proper pair"},{"location":"#filtering-unmapped-reads","text":"Use -F 4 to filter out unmapped reads. samtools view -F 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.mapped.bam Use -f 4 to keep only unmapped reads. samtools view -f 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.unmapped.bam We can use the flags subcommand to confirm that a value of four represents an unmapped read. samtools flags 4 ## 0x4 4 UNMAP","title":"Filtering unmapped reads"},{"location":"#extracting-entries-mapping-to-a-specific-loci","text":"Use samtools view and the ref:start-end syntax to extract reads mapping within a specific genomic loci; this requires a BAM index file. samtools view eg/ERR188273_chrX.bam chrX:20000-30000 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP Note that this takes into account the mapping of the entire read and not just the starting position. For example, if you specified chrX:20000-30000, a 75 bp long read that starts its mapping from position 19999 will also be returned. In addition, you can save the output as another BAM file if you want. samtools view -b eg/ERR188273_chrX.bam chrX:20000-30000 > eg/ERR188273_chrX_20000_30000.bam If you want reads mapped to a single reference (e.g. chromosome), just specify the ref and leave out the start and end values. samtools view eg/ERR188273_chrX.bam chrX | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YT:Z:UP NH:i:10 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD<DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH eg/first.bam Once again, you can use flags to verify this (it also accepts hexadecimal input). samtools flags 0x0040 ## 0x40 64 READ1","title":"Extracting only the first read from paired end BAM files"},{"location":"#stats","text":"For simple statistics use samtools flagstat . samtools flagstat eg/ERR188273_chrX.bam ## 1176360 + 0 in total (QC-passed reads + QC-failed reads) ## 1160084 + 0 primary ## 16276 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 1126961 + 0 mapped (95.80% : N/A) ## 1110685 + 0 primary mapped (95.74% : N/A) ## 1160084 + 0 paired in sequencing ## 580042 + 0 read1 ## 580042 + 0 read2 ## 1060858 + 0 properly paired (91.45% : N/A) ## 1065618 + 0 with itself and mate mapped ## 45067 + 0 singletons (3.88% : N/A) ## 0 + 0 with mate mapped to a different chr ## 0 + 0 with mate mapped to a different chr (mapQ>=5) For more stats, use samtools stats . samtools stats eg/ERR188273_chrX.bam | grep ^SN ## SN raw total sequences: 1160084 # excluding supplementary and secondary reads ## SN filtered sequences: 0 ## SN sequences: 1160084 ## SN is sorted: 1 ## SN 1st fragments: 580042 ## SN last fragments: 580042 ## SN reads mapped: 1110685 ## SN reads mapped and paired: 1065618 # paired-end technology bit set + both mates mapped ## SN reads unmapped: 49399 ## SN reads properly paired: 1060858 # proper-pair bit set ## SN reads paired: 1160084 # paired-end technology bit set ## SN reads duplicated: 0 # PCR or optical duplicate bit set ## SN reads MQ0: 905 # mapped and MQ=0 ## SN reads QC failed: 0 ## SN non-primary alignments: 16276 ## SN supplementary alignments: 0 ## SN total length: 87006300 # ignores clipping ## SN total first fragment length: 43503150 # ignores clipping ## SN total last fragment length: 43503150 # ignores clipping ## SN bases mapped: 83301375 # ignores clipping ## SN bases mapped (cigar): 83064942 # more accurate ## SN bases trimmed: 0 ## SN bases duplicated: 0 ## SN mismatches: 423271 # from NM fields ## SN error rate: 5.095663e-03 # mismatches / bases mapped (cigar) ## SN average length: 75 ## SN average first fragment length: 75 ## SN average last fragment length: 75 ## SN maximum length: 75 ## SN maximum first fragment length: 75 ## SN maximum last fragment length: 75 ## SN average quality: 36.0 ## SN insert size average: 182.7 ## SN insert size standard deviation: 176.0 ## SN inward oriented pairs: 530763 ## SN outward oriented pairs: 1042 ## SN pairs with other orientation: 1004 ## SN pairs on different chromosomes: 0 ## SN percentage of properly paired reads (%): 91.4","title":"Stats"},{"location":"#samtools-calmdfillmd","text":"The calmd or fillmd tool is useful for visualising mismatches and insertions in an alignment of a read to a reference genome. The -e argument changes identical bases between the read and reference into = . samtools view -b eg/ERR188273_chrX.bam | samtools fillmd -e - genome/chrX.fa > eg/ERR188273_chrX_fillmd.bam head eg/ERR188273_chrX_fillmd.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -b eg/ERR188273_chrX.bam ## @PG ID:samtools.1 PN:samtools PP:samtools VN:1.19.2 CL:samtools fillmd -e - genome/chrX.fa ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 =========================================================================== @@HHEGIIGAGIIIBGIIG@FECHDGHCHHHGHHFFFFFDEACC@ ## @ERR188273.14904746 ## GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT ## + ## @@HHEGIIGAGIIIBGIIG@FECH eg/ERR188273_chrX_rand.bam","title":"Random subsampling of BAM file"},{"location":"#count-number-of-reads","text":"Use samtools idxstats to print stats on a BAM file; this requires an index file which is created by running samtools index . The output of idxstats is a file with four tab-delimited columns: Reference name Sequence length of reference Number of mapped reads Number of unmapped reads samtools idxstats eg/ERR188273_chrX.bam ## chrX 156040895 1126961 45067 ## * 0 0 4332 We can use this with awk to calculate: The number of mapped reads by summing the third column. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3} END {print s}' ## 1126961 The number of reads, which is the sum of mapped and unmapped reads. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3+$4} END {print s}' ## 1176360","title":"Count number of reads"},{"location":"#obtaining-genomic-sequence","text":"Use faidx to fetch genomic sequence; coordinates are 1-based. We need to first index the reference FASTA file that was used to map the reads. samtools faidx genome/chrX.fa Now we can obtain the sequence. samtools faidx genome/chrX.fa chrX:300000-300100 ## >chrX:300000-300100 ## ctgagatcgtgccactgcactccagcctgggcgacagagcgagactccatctcaaaaaaa ## aaaaaaaaaaaaaagaTggggtctctctatgttggccaggt","title":"Obtaining genomic sequence"},{"location":"#comparing-bam-files","text":"The output from mpileup can be used to compare BAM files. The commands below generates alignments using bwa and minimap2 . len=100 n=10000 m=300 script/generate_random_seq.pl 30 1000000 1984 > test_ref.fa script/random_paired_end.pl -f test_ref.fa -l ${len} -n ${n} -m ${m} bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_bwa.bam minimap2 -ax sr test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_mm.bam The BAM files can be used with mpileup to compare the depths. samtools mpileup -s -f test_ref.fa aln_bwa.bam aln_mm.bam | head -20 ## [mpileup] 2 samples in 2 input files ## 1 1694 C 1 ^]. D ] 1 ^]. D ] ## 1 1695 G 1 . I ] 1 . I ] ## 1 1696 T 1 . J ] 1 . J ] ## 1 1697 T 1 . J ] 1 . J ] ## 1 1698 A 1 . J ] 1 . J ] ## 1 1699 A 1 . J ] 1 . J ] ## 1 1700 T 1 . J ] 1 . J ] ## 1 1701 C 1 . J ] 1 . J ] ## 1 1702 A 1 . J ] 1 . J ] ## 1 1703 A 1 . J ] 1 . J ] ## 1 1704 T 1 . J ] 1 . J ] ## 1 1705 C 1 . J ] 1 . J ] ## 1 1706 C 1 . J ] 1 . J ] ## 1 1707 C 1 . J ] 1 . J ] ## 1 1708 C 1 . J ] 1 . J ] ## 1 1709 C 1 . J ] 1 . J ] ## 1 1710 A 1 . J ] 1 . J ] ## 1 1711 A 1 . J ] 1 . J ] ## 1 1712 C 1 . J ] 1 . J ] ## 1 1713 A 1 . J ] 1 . J ] Another approach is to use deepTools and the bamCompare command. The bigWig output file shows the ratio of reads between b1 and b2 in 50 bp (default) windows.","title":"Comparing BAM files"},{"location":"#converting-reference-names","text":"One of the most annoying bioinformatics problems is the use of different chromosome names, e.g. chr1 vs 1, in different references even when the sequences are identical. The GRCh38 reference downloaded from Ensembl has chromosome names without the chr : >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF Whereas the reference names from UCSC has the chr : >chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38 Luckily you can change the reference names using samtools reheader but just make sure your reference sequences are actually identical. samtools view eg/ERR188273_chrX.bam | head -2 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP View header samtools view -H eg/ERR188273_chrX.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -H eg/ERR188273_chrX.bam Substitute header with new name. samtools view -H eg/ERR188273_chrX.bam | sed 's/SN:chrX/SN:X/' > eg/my_header Save bam file with new ref and check it out. samtools reheader eg/my_header eg/ERR188273_chrX.bam > eg/ERR188273_X.bam samtools view eg/ERR188273_X.bam | head -2 ## ERR188273.4711308 73 X 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 X 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP","title":"Converting reference names"},{"location":"#coverage","text":"Coverage can mean the: average depth of each covered base percentage of bases covered samtools depth and samtools mpileup can be used to indicate the depth of each covered base (and used to calculate the average depth. samtools coverage will provide both the average depth and percentage of bases covered per chromosome/reference sequence. samtools depth will return three columns: reference, position, and coverage. samtools depth eg/ERR188273_chrX.bam | head ## chrX 21649 1 ## chrX 21650 1 ## chrX 21651 1 ## chrX 21652 1 ## chrX 21653 1 ## chrX 21654 1 ## chrX 21655 1 ## chrX 21656 1 ## chrX 21657 1 ## chrX 21658 1 The average depth can be calculated by summing the third column and dividing by the total number of bases (be sure to use -a ). samtools depth -@ 2 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf \"Bases covered:\\t%.3f\\nCoverage:\\t%.3f\\n\", $., $cov}' ## Bases covered: 156040895.000 ## Coverage: 0.532 The samtools mpileup command also provides depth information (but not for reads that have a mapping quality of 0, by default) with some additional information: Sequence name 1-based coordinate Reference base (when used with -f ) Number of reads covering this position Read bases Base qualities Alignment mapping qualities (when used with -s ) samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam | head ## [mpileup] 1 samples in 1 input files ## chrX 251271 g 1 ^]. @ ] ## chrX 251272 a 1 . @ ] ## chrX 251273 a 1 . < ] ## chrX 251274 a 1 . D ] ## chrX 251275 a 1 . D ] ## chrX 251276 a 1 . D ] ## chrX 251277 t 1 . D ] ## chrX 251278 g 1 . D ] ## chrX 251279 g 1 . F ] ## chrX 251280 g 1 . B ] Note that the start of the samtools mpileup output differ from the start of the samtools depth output. This is because mpileup performs some filtering by default. In the case of this example, read pairs that are not both mapped will be ignored. To count these \u201corphan\u201d reads, use the --count-orphans argument. samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam | head ## [mpileup] 1 samples in 1 input files ## chrX 21649 g 0 * * * ## chrX 21650 a 1 . D ! ## chrX 21651 t 1 . F ! ## chrX 21652 c 1 . F ! ## chrX 21653 a 1 . H ! ## chrX 21654 c 1 . G ! ## chrX 21655 g 1 . H ! ## chrX 21656 a 1 . B ! ## chrX 21657 g 1 . H ! ## chrX 21658 g 1 . I ! In addition mpileup performs \u201cper-Base Alignment Quality\u201d (BAQ) by default and will adjust base quality scores. The default behaviour to to skip bases with baseQ/BAQ smaller than 13. If you are finding discrepancies between mpileup \u2019s coverage calculation with another coverage tool, you can either set --min-BQ to 0 or use --no-BAQ to disable BAQ. I have an old blog post on using mpileup . samtools coverage will provide the following coverage statistics: rname - Reference name / chromosome startpos - Start position endpos - End position (or sequence length) numreads - Number reads aligned to the region (after filtering) covbases - Number of covered bases with depth >= 1 coverage - Proportion of covered bases [0..1] meandepth - Mean depth of coverage meanbaseq - Mean base quality in covered region meanmapq - Mean mapping quality of selected reads samtools coverage eg/ERR188273_chrX.bam ## #rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq ## chrX 1 156040895 1110685 3402037 2.18022 0.532299 36.3 59.4 The example BAM file only contains reads for chrX hence the statistics are only returned for chrX . Returning to our coverage definition at the start of this section: average depth of each covered base = meandepth percentage of bases covered = covbases","title":"Coverage"},{"location":"#stargazers-over-time","text":"","title":"Stargazers over time"}]}
\ No newline at end of file
+{"config":{"indexing":"full","lang":["en"],"min_search_length":3,"prebuild_index":false,"separator":"[\\s\\-]+"},"docs":[{"location":"","text":"Learning the BAM format Introduction SAMtools provides various (sub)tools for manipulating alignments in the SAM/BAM format. The SAM (Sequence Alignment/Map) format (BAM is just the binary form of SAM) is currently the de facto standard for storing large nucleotide sequence alignments. If you are working with high-throughput sequencing data, at some point you will probably have to deal with SAM/BAM files, so familiarise yourself with them! For the latest information on SAMtools, please refer to the release notes . The examples in this README use the ERR188273_chrX.bam BAM file (stored in the eg folder) generated as per https://github.com/davetang/rnaseq using the HISAT2 + StringTie2 RNA-seq pipeline. This README is generated using the create_readme.sh script; if you want to generate this file yourself, please use this Docker image and the Makefile in this directory. For example: # clone this repo git clone https://github.com/davetang/learning_bam_file.git cd learning_bam_file docker pull davetang/r_build:4.1.2 docker run --rm -it -v $(pwd):/work davetang/r_build:4.1.2 /bin/bash # inside the Docker container make Installing SAMtools For installing SAMtools, I recommend using Conda and the Bioconda samtools package . I also recommend using Miniconda instead of Anaconda because Anaconda comes with a lot of tools/packages that you will probably not use. I wrote a short introduction to Conda if you want to find learn more. Once you have installed Miniconda, you can install SAMtools as follows: conda install -c bioconda samtools Otherwise you can download the source and compile it yourself; change dir to the location you want samtools to be installed. samtools will be installed in ${dir}/bin , so make sure this is in your $PATH . #!/usr/bin/env bash set -euo pipefail ver=1.15 tool=samtools url=https://github.com/samtools/${tool}/releases/download/${ver}/${tool}-${ver}.tar.bz2 dir=${HOME}/local wget ${url} tar xjf ${tool}-${ver}.tar.bz2 cd ${tool}-${ver} ./configure --prefix=${dir} make && make install cd .. rm -rf ${tool}-${ver} ${tool}-${ver}.tar.bz2 >&2 echo Done exit 0 Basic usage If you run samtools on the terminal without any parameters or with --help , all the available utilities are listed: samtools --help ## ## Program: samtools (Tools for alignments in the SAM format) ## Version: 1.19.2 (using htslib 1.19.1) ## ## Usage: samtools [options] ## ## Commands: ## -- Indexing ## dict create a sequence dictionary file ## faidx index/extract FASTA ## fqidx index/extract FASTQ ## index index alignment ## ## -- Editing ## calmd recalculate MD/NM tags and '=' bases ## fixmate fix mate information ## reheader replace BAM header ## targetcut cut fosmid regions (for fosmid pool only) ## addreplacerg adds or replaces RG tags ## markdup mark duplicates ## ampliconclip clip oligos from the end of reads ## ## -- File operations ## collate shuffle and group alignments by name ## cat concatenate BAMs ## consensus produce a consensus Pileup/FASTA/FASTQ ## merge merge sorted alignments ## mpileup multi-way pileup ## sort sort alignment file ## split splits a file by read group ## quickcheck quickly check if SAM/BAM/CRAM file appears intact ## fastq converts a BAM to a FASTQ ## fasta converts a BAM to a FASTA ## import Converts FASTA or FASTQ files to SAM/BAM/CRAM ## reference Generates a reference from aligned data ## reset Reverts aligner changes in reads ## ## -- Statistics ## bedcov read depth per BED region ## coverage alignment depth and percent coverage ## depth compute the depth ## flagstat simple stats ## idxstats BAM index stats ## cram-size list CRAM Content-ID and Data-Series sizes ## phase phase heterozygotes ## stats generate stats (former bamcheck) ## ampliconstats generate amplicon specific stats ## ## -- Viewing ## flags explain BAM flags ## head header viewer ## tview text alignment viewer ## view SAM<->BAM<->CRAM conversion ## depad convert padded BAM to unpadded BAM ## samples list the samples in a set of SAM/BAM/CRAM files ## ## -- Misc ## help [cmd] display this help message or help for [cmd] ## version detailed version information Viewing Use bioSyntax to prettify your output. samtools view aln.bam | sam-less Converting a SAM file to a BAM file A BAM file is just a SAM file but stored in binary format; you should always convert your SAM files into BAM format since they are smaller in size and are faster to manipulate. I don\u2019t have a SAM file in the example folder, so let\u2019s create one and check out the first ten lines. Note: remember to use -h to ensure the SAM file contains the sequence header information. Generally, I recommend storing only sorted BAM files as they use even less disk space and are faster to process. samtools view -h eg/ERR188273_chrX.bam > eg/ERR188273_chrX.sam Notice that the SAM file is much larger than the BAM file. Size of SAM file. ls -lh eg/ERR188273_chrX.sam ## -rw-r--r-- 1 root root 321M Mar 25 07:18 eg/ERR188273_chrX.sam Size of BAM file. ls -lh eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 67M Mar 25 07:16 eg/ERR188273_chrX.bam We can use head to view a SAM file. head eg/ERR188273_chrX.sam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -h eg/ERR188273_chrX.bam ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH aln.bam If the header information is available, we can convert a SAM file into BAM by using samtools view -b . In newer versions of SAMtools, the input format is auto-detected, so we no longer need the -S parameter. samtools view -b eg/ERR188273_chrX.sam > eg/my.bam Converting a BAM file to a CRAM file The CRAM format is even more compact. Use samtools view with the -T and -C arguments to convert a BAM file into CRAM. samtools view -T genome/chrX.fa -C -o eg/ERR188273_chrX.cram eg/ERR188273_chrX.bam ls -lh eg/ERR188273_chrX.[sbcr]*am ## -rw-r--r-- 1 root root 67M Mar 25 07:16 eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 40M Mar 25 07:18 eg/ERR188273_chrX.cram ## -rw-r--r-- 1 root root 321M Mar 25 07:18 eg/ERR188273_chrX.sam You can use samtools view to view a CRAM file just as you would for a BAM file. samtools view eg/ERR188273_chrX.cram | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:10 MD:Z:75 NM:i:0 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD< test_ref.fa script/random_paired_end.pl -f test_ref.fa -l 100 -n 10000 -m 300 -d 0.1 bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l100_n10000_d300_1984_1.fq.gz l100_n10000_d300_1984_2.fq.gz > aln.sam 2> /dev/null samtools flagstat will indicate that some reads (about 10%) mapped to different chromosomes. samtools flagstat aln.sam ## 20000 + 0 in total (QC-passed reads + QC-failed reads) ## 20000 + 0 primary ## 0 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 20000 + 0 mapped (100.00% : N/A) ## 20000 + 0 primary mapped (100.00% : N/A) ## 20000 + 0 paired in sequencing ## 10000 + 0 read1 ## 10000 + 0 read2 ## 18012 + 0 properly paired (90.06% : N/A) ## 20000 + 0 with itself and mate mapped ## 0 + 0 singletons (0.00% : N/A) ## 1988 + 0 with mate mapped to a different chr ## 1988 + 0 with mate mapped to a different chr (mapQ>=5) Flag of a proper pair. samtools flag $(samtools view -f 2 aln.sam | head -1 | cut -f2) ## 0x63 99 PAIRED,PROPER_PAIR,MREVERSE,READ1 Flag of a pair (that is not a proper pair). samtools flag $(samtools view -F 2 aln.sam | head -1 | cut -f2) ## 0x61 97 PAIRED,MREVERSE,READ1 Filtering unmapped reads Use -F 4 to filter out unmapped reads. samtools view -F 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.mapped.bam Use -f 4 to keep only unmapped reads. samtools view -f 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.unmapped.bam We can use the flags subcommand to confirm that a value of four represents an unmapped read. samtools flags 4 ## 0x4 4 UNMAP Extracting entries mapping to a specific loci Use samtools view and the ref:start-end syntax to extract reads mapping within a specific genomic loci; this requires a BAM index file. samtools view eg/ERR188273_chrX.bam chrX:20000-30000 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP Note that this takes into account the mapping of the entire read and not just the starting position. For example, if you specified chrX:20000-30000, a 75 bp long read that starts its mapping from position 19999 will also be returned. In addition, you can save the output as another BAM file if you want. samtools view -b eg/ERR188273_chrX.bam chrX:20000-30000 > eg/ERR188273_chrX_20000_30000.bam If you want reads mapped to a single reference (e.g. chromosome), just specify the ref and leave out the start and end values. samtools view eg/ERR188273_chrX.bam chrX | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YT:Z:UP NH:i:10 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD<DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH eg/first.bam Once again, you can use flags to verify this (it also accepts hexadecimal input). samtools flags 0x0040 ## 0x40 64 READ1 Stats For simple statistics use samtools flagstat . samtools flagstat eg/ERR188273_chrX.bam ## 1176360 + 0 in total (QC-passed reads + QC-failed reads) ## 1160084 + 0 primary ## 16276 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 1126961 + 0 mapped (95.80% : N/A) ## 1110685 + 0 primary mapped (95.74% : N/A) ## 1160084 + 0 paired in sequencing ## 580042 + 0 read1 ## 580042 + 0 read2 ## 1060858 + 0 properly paired (91.45% : N/A) ## 1065618 + 0 with itself and mate mapped ## 45067 + 0 singletons (3.88% : N/A) ## 0 + 0 with mate mapped to a different chr ## 0 + 0 with mate mapped to a different chr (mapQ>=5) For more stats, use samtools stats . samtools stats eg/ERR188273_chrX.bam | grep ^SN ## SN raw total sequences: 1160084 # excluding supplementary and secondary reads ## SN filtered sequences: 0 ## SN sequences: 1160084 ## SN is sorted: 1 ## SN 1st fragments: 580042 ## SN last fragments: 580042 ## SN reads mapped: 1110685 ## SN reads mapped and paired: 1065618 # paired-end technology bit set + both mates mapped ## SN reads unmapped: 49399 ## SN reads properly paired: 1060858 # proper-pair bit set ## SN reads paired: 1160084 # paired-end technology bit set ## SN reads duplicated: 0 # PCR or optical duplicate bit set ## SN reads MQ0: 905 # mapped and MQ=0 ## SN reads QC failed: 0 ## SN non-primary alignments: 16276 ## SN supplementary alignments: 0 ## SN total length: 87006300 # ignores clipping ## SN total first fragment length: 43503150 # ignores clipping ## SN total last fragment length: 43503150 # ignores clipping ## SN bases mapped: 83301375 # ignores clipping ## SN bases mapped (cigar): 83064942 # more accurate ## SN bases trimmed: 0 ## SN bases duplicated: 0 ## SN mismatches: 423271 # from NM fields ## SN error rate: 5.095663e-03 # mismatches / bases mapped (cigar) ## SN average length: 75 ## SN average first fragment length: 75 ## SN average last fragment length: 75 ## SN maximum length: 75 ## SN maximum first fragment length: 75 ## SN maximum last fragment length: 75 ## SN average quality: 36.0 ## SN insert size average: 182.7 ## SN insert size standard deviation: 176.0 ## SN inward oriented pairs: 530763 ## SN outward oriented pairs: 1042 ## SN pairs with other orientation: 1004 ## SN pairs on different chromosomes: 0 ## SN percentage of properly paired reads (%): 91.4 samtools calmd/fillmd The calmd or fillmd tool is useful for visualising mismatches and insertions in an alignment of a read to a reference genome. The -e argument changes identical bases between the read and reference into = . samtools view -b eg/ERR188273_chrX.bam | samtools fillmd -e - genome/chrX.fa > eg/ERR188273_chrX_fillmd.bam head eg/ERR188273_chrX_fillmd.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -b eg/ERR188273_chrX.bam ## @PG ID:samtools.1 PN:samtools PP:samtools VN:1.19.2 CL:samtools fillmd -e - genome/chrX.fa ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 =========================================================================== @@HHEGIIGAGIIIBGIIG@FECHDGHCHHHGHHFFFFFDEACC@ ## @ERR188273.14904746 ## GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT ## + ## @@HHEGIIGAGIIIBGIIG@FECH eg/ERR188273_chrX_rand.bam Count number of reads Use samtools idxstats to print stats on a BAM file; this requires an index file which is created by running samtools index . The output of idxstats is a file with four tab-delimited columns: Reference name Sequence length of reference Number of mapped reads Number of unmapped reads samtools idxstats eg/ERR188273_chrX.bam ## chrX 156040895 1126961 45067 ## * 0 0 4332 We can use this with awk to calculate: The number of mapped reads by summing the third column. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3} END {print s}' ## 1126961 The number of reads, which is the sum of mapped and unmapped reads. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3+$4} END {print s}' ## 1176360 Obtaining genomic sequence Use faidx to fetch genomic sequence; coordinates are 1-based. We need to first index the reference FASTA file that was used to map the reads. samtools faidx genome/chrX.fa Now we can obtain the sequence. samtools faidx genome/chrX.fa chrX:300000-300100 ## >chrX:300000-300100 ## ctgagatcgtgccactgcactccagcctgggcgacagagcgagactccatctcaaaaaaa ## aaaaaaaaaaaaaagaTggggtctctctatgttggccaggt Comparing BAM files The output from mpileup can be used to compare BAM files. The commands below generates alignments using bwa and minimap2 . len=100 n=10000 m=300 script/generate_random_seq.pl 30 1000000 1984 > test_ref.fa script/random_paired_end.pl -f test_ref.fa -l ${len} -n ${n} -m ${m} bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_bwa.bam minimap2 -ax sr test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_mm.bam The BAM files can be used with mpileup to compare the depths. samtools mpileup -s -f test_ref.fa aln_bwa.bam aln_mm.bam | head -20 ## [mpileup] 2 samples in 2 input files ## 1 4020 C 1 ^]. E ] 1 ^]. E ] ## 1 4021 T 1 . J ] 1 . J ] ## 1 4022 C 1 . J ] 1 . J ] ## 1 4023 T 1 . J ] 1 . J ] ## 1 4024 G 1 . J ] 1 . J ] ## 1 4025 T 1 . J ] 1 . J ] ## 1 4026 T 1 . J ] 1 . J ] ## 1 4027 A 1 . J ] 1 . J ] ## 1 4028 T 1 . J ] 1 . J ] ## 1 4029 A 1 . J ] 1 . J ] ## 1 4030 G 1 . J ] 1 . J ] ## 1 4031 C 1 . J ] 1 . J ] ## 1 4032 G 1 . J ] 1 . J ] ## 1 4033 G 1 . J ] 1 . J ] ## 1 4034 G 1 . J ] 1 . J ] ## 1 4035 A 1 . J ] 1 . J ] ## 1 4036 T 1 . J ] 1 . J ] ## 1 4037 T 1 . J ] 1 . J ] ## 1 4038 A 1 . J ] 1 . J ] ## 1 4039 C 1 . J ] 1 . J ] Another approach is to use deepTools and the bamCompare command. The bigWig output file shows the ratio of reads between b1 and b2 in 50 bp (default) windows. Converting reference names One of the most annoying bioinformatics problems is the use of different chromosome names, e.g. chr1 vs 1, in different references even when the sequences are identical. The GRCh38 reference downloaded from Ensembl has chromosome names without the chr : >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF Whereas the reference names from UCSC has the chr : >chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38 Luckily you can change the reference names using samtools reheader but just make sure your reference sequences are actually identical. samtools view eg/ERR188273_chrX.bam | head -2 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP View header samtools view -H eg/ERR188273_chrX.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -H eg/ERR188273_chrX.bam Substitute header with new name. samtools view -H eg/ERR188273_chrX.bam | sed 's/SN:chrX/SN:X/' > eg/my_header Save bam file with new ref and check it out. samtools reheader eg/my_header eg/ERR188273_chrX.bam > eg/ERR188273_X.bam samtools view eg/ERR188273_X.bam | head -2 ## ERR188273.4711308 73 X 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 X 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP Coverage Coverage can mean the: average depth of each covered base percentage of bases covered samtools depth and samtools mpileup can be used to indicate the depth of each covered base (and used to calculate the average depth. samtools coverage will provide both the average depth and percentage of bases covered per chromosome/reference sequence. samtools depth will return three columns: reference, position, and coverage. samtools depth -@ 4 eg/ERR188273_chrX.bam > ERR188273_depth.tsv head ERR188273_depth.tsv ## chrX 21649 1 ## chrX 21650 1 ## chrX 21651 1 ## chrX 21652 1 ## chrX 21653 1 ## chrX 21654 1 ## chrX 21655 1 ## chrX 21656 1 ## chrX 21657 1 ## chrX 21658 1 The average depth can be calculated by summing the third column and dividing by the total number of bases (be sure to use -a with samtools depth as that will output all positions including zero depth). samtools depth -@ 4 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf \"Bases covered:\\t%.3f\\nCoverage:\\t%.3f\\n\", $., $cov}' ## Bases covered: 156040895.000 ## Coverage: 0.532 The samtools mpileup command also provides depth information (but not for reads that have a mapping quality of 0, by default) with some additional information: Sequence name 1-based coordinate Reference base (when used with -f ) Number of reads covering this position Read bases Base qualities Alignment mapping qualities (when used with -s ) samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam > ERR188273_mpileup.tsv head ERR188273_mpileup.tsv ## [mpileup] 1 samples in 1 input files ## chrX 251271 g 1 ^]. @ ] ## chrX 251272 a 1 . @ ] ## chrX 251273 a 1 . < ] ## chrX 251274 a 1 . D ] ## chrX 251275 a 1 . D ] ## chrX 251276 a 1 . D ] ## chrX 251277 t 1 . D ] ## chrX 251278 g 1 . D ] ## chrX 251279 g 1 . F ] ## chrX 251280 g 1 . B ] Note that the start of the samtools mpileup output differ from the start of the samtools depth output. This is because mpileup performs some filtering by default. In the case of this example, read pairs that are not both mapped will be ignored. To count these \u201corphan\u201d reads, use the --count-orphans argument. samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam > ERR188273_mpileup_orphans.tsv head ERR188273_mpileup_orphans.tsv ## [mpileup] 1 samples in 1 input files ## chrX 21649 g 0 * * * ## chrX 21650 a 1 . D ! ## chrX 21651 t 1 . F ! ## chrX 21652 c 1 . F ! ## chrX 21653 a 1 . H ! ## chrX 21654 c 1 . G ! ## chrX 21655 g 1 . H ! ## chrX 21656 a 1 . B ! ## chrX 21657 g 1 . H ! ## chrX 21658 g 1 . I ! In addition mpileup performs \u201cper-Base Alignment Quality\u201d (BAQ) by default and will adjust base quality scores. The default behaviour to to skip bases with baseQ/BAQ smaller than 13. If you are finding discrepancies between mpileup \u2019s coverage calculation with another coverage tool, you can either set --min-BQ to 0 or use --no-BAQ to disable BAQ. I have an old blog post on using mpileup . samtools coverage will provide the following coverage statistics: rname - Reference name / chromosome startpos - Start position endpos - End position (or sequence length) numreads - Number reads aligned to the region (after filtering) covbases - Number of covered bases with depth >= 1 coverage - Proportion of covered bases [0..1] meandepth - Mean depth of coverage meanbaseq - Mean base quality in covered region meanmapq - Mean mapping quality of selected reads samtools coverage eg/ERR188273_chrX.bam ## #rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq ## chrX 1 156040895 1110685 3402037 2.18022 0.532299 36.3 59.4 The example BAM file only contains reads for chrX hence the statistics are only returned for chrX . Returning to our coverage definition at the start of this section: average depth of each covered base = meandepth percentage of bases covered = covbases The mosdepth tool can also calculate depth (and much faster than samtools depth ) per base or within a given window. The output is given in a BED file, where the fourth column indicates the coverage. mosdepth ERR188275 eg/ERR188273_chrX.bam gunzip -c ERR188273.per-base.bed.gz | head ## gzip: ERR188273.per-base.bed.gz: No such file or directory Coverage in using a 500 bp window. mosdepth -n --fast-mode --by 500 ERR188275_500 eg/ERR188273_chrX.bam gunzip -c ERR188275_500.regions.bed.gz | head ## chrX 0 500 0.00 ## chrX 500 1000 0.00 ## chrX 1000 1500 0.00 ## chrX 1500 2000 0.00 ## chrX 2000 2500 0.00 ## chrX 2500 3000 0.00 ## chrX 3000 3500 0.00 ## chrX 3500 4000 0.00 ## chrX 4000 4500 0.00 ## chrX 4500 5000 0.00 Stargazers over time","title":"Home"},{"location":"#learning-the-bam-format","text":"","title":"Learning the BAM format"},{"location":"#introduction","text":"SAMtools provides various (sub)tools for manipulating alignments in the SAM/BAM format. The SAM (Sequence Alignment/Map) format (BAM is just the binary form of SAM) is currently the de facto standard for storing large nucleotide sequence alignments. If you are working with high-throughput sequencing data, at some point you will probably have to deal with SAM/BAM files, so familiarise yourself with them! For the latest information on SAMtools, please refer to the release notes . The examples in this README use the ERR188273_chrX.bam BAM file (stored in the eg folder) generated as per https://github.com/davetang/rnaseq using the HISAT2 + StringTie2 RNA-seq pipeline. This README is generated using the create_readme.sh script; if you want to generate this file yourself, please use this Docker image and the Makefile in this directory. For example: # clone this repo git clone https://github.com/davetang/learning_bam_file.git cd learning_bam_file docker pull davetang/r_build:4.1.2 docker run --rm -it -v $(pwd):/work davetang/r_build:4.1.2 /bin/bash # inside the Docker container make","title":"Introduction"},{"location":"#installing-samtools","text":"For installing SAMtools, I recommend using Conda and the Bioconda samtools package . I also recommend using Miniconda instead of Anaconda because Anaconda comes with a lot of tools/packages that you will probably not use. I wrote a short introduction to Conda if you want to find learn more. Once you have installed Miniconda, you can install SAMtools as follows: conda install -c bioconda samtools Otherwise you can download the source and compile it yourself; change dir to the location you want samtools to be installed. samtools will be installed in ${dir}/bin , so make sure this is in your $PATH . #!/usr/bin/env bash set -euo pipefail ver=1.15 tool=samtools url=https://github.com/samtools/${tool}/releases/download/${ver}/${tool}-${ver}.tar.bz2 dir=${HOME}/local wget ${url} tar xjf ${tool}-${ver}.tar.bz2 cd ${tool}-${ver} ./configure --prefix=${dir} make && make install cd .. rm -rf ${tool}-${ver} ${tool}-${ver}.tar.bz2 >&2 echo Done exit 0","title":"Installing SAMtools"},{"location":"#basic-usage","text":"If you run samtools on the terminal without any parameters or with --help , all the available utilities are listed: samtools --help ## ## Program: samtools (Tools for alignments in the SAM format) ## Version: 1.19.2 (using htslib 1.19.1) ## ## Usage: samtools [options] ## ## Commands: ## -- Indexing ## dict create a sequence dictionary file ## faidx index/extract FASTA ## fqidx index/extract FASTQ ## index index alignment ## ## -- Editing ## calmd recalculate MD/NM tags and '=' bases ## fixmate fix mate information ## reheader replace BAM header ## targetcut cut fosmid regions (for fosmid pool only) ## addreplacerg adds or replaces RG tags ## markdup mark duplicates ## ampliconclip clip oligos from the end of reads ## ## -- File operations ## collate shuffle and group alignments by name ## cat concatenate BAMs ## consensus produce a consensus Pileup/FASTA/FASTQ ## merge merge sorted alignments ## mpileup multi-way pileup ## sort sort alignment file ## split splits a file by read group ## quickcheck quickly check if SAM/BAM/CRAM file appears intact ## fastq converts a BAM to a FASTQ ## fasta converts a BAM to a FASTA ## import Converts FASTA or FASTQ files to SAM/BAM/CRAM ## reference Generates a reference from aligned data ## reset Reverts aligner changes in reads ## ## -- Statistics ## bedcov read depth per BED region ## coverage alignment depth and percent coverage ## depth compute the depth ## flagstat simple stats ## idxstats BAM index stats ## cram-size list CRAM Content-ID and Data-Series sizes ## phase phase heterozygotes ## stats generate stats (former bamcheck) ## ampliconstats generate amplicon specific stats ## ## -- Viewing ## flags explain BAM flags ## head header viewer ## tview text alignment viewer ## view SAM<->BAM<->CRAM conversion ## depad convert padded BAM to unpadded BAM ## samples list the samples in a set of SAM/BAM/CRAM files ## ## -- Misc ## help [cmd] display this help message or help for [cmd] ## version detailed version information","title":"Basic usage"},{"location":"#viewing","text":"Use bioSyntax to prettify your output. samtools view aln.bam | sam-less","title":"Viewing"},{"location":"#converting-a-sam-file-to-a-bam-file","text":"A BAM file is just a SAM file but stored in binary format; you should always convert your SAM files into BAM format since they are smaller in size and are faster to manipulate. I don\u2019t have a SAM file in the example folder, so let\u2019s create one and check out the first ten lines. Note: remember to use -h to ensure the SAM file contains the sequence header information. Generally, I recommend storing only sorted BAM files as they use even less disk space and are faster to process. samtools view -h eg/ERR188273_chrX.bam > eg/ERR188273_chrX.sam Notice that the SAM file is much larger than the BAM file. Size of SAM file. ls -lh eg/ERR188273_chrX.sam ## -rw-r--r-- 1 root root 321M Mar 25 07:18 eg/ERR188273_chrX.sam Size of BAM file. ls -lh eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 67M Mar 25 07:16 eg/ERR188273_chrX.bam We can use head to view a SAM file. head eg/ERR188273_chrX.sam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -h eg/ERR188273_chrX.bam ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH aln.bam If the header information is available, we can convert a SAM file into BAM by using samtools view -b . In newer versions of SAMtools, the input format is auto-detected, so we no longer need the -S parameter. samtools view -b eg/ERR188273_chrX.sam > eg/my.bam","title":"Converting a SAM file to a BAM file"},{"location":"#converting-a-bam-file-to-a-cram-file","text":"The CRAM format is even more compact. Use samtools view with the -T and -C arguments to convert a BAM file into CRAM. samtools view -T genome/chrX.fa -C -o eg/ERR188273_chrX.cram eg/ERR188273_chrX.bam ls -lh eg/ERR188273_chrX.[sbcr]*am ## -rw-r--r-- 1 root root 67M Mar 25 07:16 eg/ERR188273_chrX.bam ## -rw-r--r-- 1 root root 40M Mar 25 07:18 eg/ERR188273_chrX.cram ## -rw-r--r-- 1 root root 321M Mar 25 07:18 eg/ERR188273_chrX.sam You can use samtools view to view a CRAM file just as you would for a BAM file. samtools view eg/ERR188273_chrX.cram | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:2 MD:Z:70 NM:i:0 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 YT:Z:UP NH:i:10 MD:Z:75 NM:i:0 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD< test_ref.fa script/random_paired_end.pl -f test_ref.fa -l 100 -n 10000 -m 300 -d 0.1 bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l100_n10000_d300_1984_1.fq.gz l100_n10000_d300_1984_2.fq.gz > aln.sam 2> /dev/null samtools flagstat will indicate that some reads (about 10%) mapped to different chromosomes. samtools flagstat aln.sam ## 20000 + 0 in total (QC-passed reads + QC-failed reads) ## 20000 + 0 primary ## 0 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 20000 + 0 mapped (100.00% : N/A) ## 20000 + 0 primary mapped (100.00% : N/A) ## 20000 + 0 paired in sequencing ## 10000 + 0 read1 ## 10000 + 0 read2 ## 18012 + 0 properly paired (90.06% : N/A) ## 20000 + 0 with itself and mate mapped ## 0 + 0 singletons (0.00% : N/A) ## 1988 + 0 with mate mapped to a different chr ## 1988 + 0 with mate mapped to a different chr (mapQ>=5) Flag of a proper pair. samtools flag $(samtools view -f 2 aln.sam | head -1 | cut -f2) ## 0x63 99 PAIRED,PROPER_PAIR,MREVERSE,READ1 Flag of a pair (that is not a proper pair). samtools flag $(samtools view -F 2 aln.sam | head -1 | cut -f2) ## 0x61 97 PAIRED,MREVERSE,READ1","title":"Proper pair"},{"location":"#filtering-unmapped-reads","text":"Use -F 4 to filter out unmapped reads. samtools view -F 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.mapped.bam Use -f 4 to keep only unmapped reads. samtools view -f 4 -b eg/ERR188273_chrX.bam > eg/ERR188273_chrX.unmapped.bam We can use the flags subcommand to confirm that a value of four represents an unmapped read. samtools flags 4 ## 0x4 4 UNMAP","title":"Filtering unmapped reads"},{"location":"#extracting-entries-mapping-to-a-specific-loci","text":"Use samtools view and the ref:start-end syntax to extract reads mapping within a specific genomic loci; this requires a BAM index file. samtools view eg/ERR188273_chrX.bam chrX:20000-30000 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP Note that this takes into account the mapping of the entire read and not just the starting position. For example, if you specified chrX:20000-30000, a 75 bp long read that starts its mapping from position 19999 will also be returned. In addition, you can save the output as another BAM file if you want. samtools view -b eg/ERR188273_chrX.bam chrX:20000-30000 > eg/ERR188273_chrX_20000_30000.bam If you want reads mapped to a single reference (e.g. chromosome), just specify the ref and leave out the start and end values. samtools view eg/ERR188273_chrX.bam chrX | head ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECHDC>>+@::8-755-BBBFDDEHHBGGEGHEEIJIIGIJJIGEIIIJJJIIJJIGGHHHGGFFFFF@@C AS:i:0 ZS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:75 YT:Z:UP NH:i:10 ## ERR188273.5927795 385 chrX 265991 1 75M = 114048277 0 TGGGACTACAGGCGCCCGCCACCACGCCCGGCTAATTTTTTGTATTTTTAGTAGAGACGGGGTTTCACCGTGTTA =?BB??BD?FBHHBEAE@CDGG@HH=FA@GEGE;FGACCHBE6?A=ACE9)7@DCE>>5'3=338:;:>2;3?BCFFEEHHHEEGIGGHAGFBBHFBHHEHCG@<@ABG??@@?BB9GBGAFFD<DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT @@HHEGIIGAGIIIBGIIG@FECH eg/first.bam Once again, you can use flags to verify this (it also accepts hexadecimal input). samtools flags 0x0040 ## 0x40 64 READ1","title":"Extracting only the first read from paired end BAM files"},{"location":"#stats","text":"For simple statistics use samtools flagstat . samtools flagstat eg/ERR188273_chrX.bam ## 1176360 + 0 in total (QC-passed reads + QC-failed reads) ## 1160084 + 0 primary ## 16276 + 0 secondary ## 0 + 0 supplementary ## 0 + 0 duplicates ## 0 + 0 primary duplicates ## 1126961 + 0 mapped (95.80% : N/A) ## 1110685 + 0 primary mapped (95.74% : N/A) ## 1160084 + 0 paired in sequencing ## 580042 + 0 read1 ## 580042 + 0 read2 ## 1060858 + 0 properly paired (91.45% : N/A) ## 1065618 + 0 with itself and mate mapped ## 45067 + 0 singletons (3.88% : N/A) ## 0 + 0 with mate mapped to a different chr ## 0 + 0 with mate mapped to a different chr (mapQ>=5) For more stats, use samtools stats . samtools stats eg/ERR188273_chrX.bam | grep ^SN ## SN raw total sequences: 1160084 # excluding supplementary and secondary reads ## SN filtered sequences: 0 ## SN sequences: 1160084 ## SN is sorted: 1 ## SN 1st fragments: 580042 ## SN last fragments: 580042 ## SN reads mapped: 1110685 ## SN reads mapped and paired: 1065618 # paired-end technology bit set + both mates mapped ## SN reads unmapped: 49399 ## SN reads properly paired: 1060858 # proper-pair bit set ## SN reads paired: 1160084 # paired-end technology bit set ## SN reads duplicated: 0 # PCR or optical duplicate bit set ## SN reads MQ0: 905 # mapped and MQ=0 ## SN reads QC failed: 0 ## SN non-primary alignments: 16276 ## SN supplementary alignments: 0 ## SN total length: 87006300 # ignores clipping ## SN total first fragment length: 43503150 # ignores clipping ## SN total last fragment length: 43503150 # ignores clipping ## SN bases mapped: 83301375 # ignores clipping ## SN bases mapped (cigar): 83064942 # more accurate ## SN bases trimmed: 0 ## SN bases duplicated: 0 ## SN mismatches: 423271 # from NM fields ## SN error rate: 5.095663e-03 # mismatches / bases mapped (cigar) ## SN average length: 75 ## SN average first fragment length: 75 ## SN average last fragment length: 75 ## SN maximum length: 75 ## SN maximum first fragment length: 75 ## SN maximum last fragment length: 75 ## SN average quality: 36.0 ## SN insert size average: 182.7 ## SN insert size standard deviation: 176.0 ## SN inward oriented pairs: 530763 ## SN outward oriented pairs: 1042 ## SN pairs with other orientation: 1004 ## SN pairs on different chromosomes: 0 ## SN percentage of properly paired reads (%): 91.4","title":"Stats"},{"location":"#samtools-calmdfillmd","text":"The calmd or fillmd tool is useful for visualising mismatches and insertions in an alignment of a read to a reference genome. The -e argument changes identical bases between the read and reference into = . samtools view -b eg/ERR188273_chrX.bam | samtools fillmd -e - genome/chrX.fa > eg/ERR188273_chrX_fillmd.bam head eg/ERR188273_chrX_fillmd.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -b eg/ERR188273_chrX.bam ## @PG ID:samtools.1 PN:samtools PP:samtools VN:1.19.2 CL:samtools fillmd -e - genome/chrX.fa ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP ## ERR188273.4711308 329 chrX 233717 0 5S70M = 233717 0 CGGGT====================================================================== @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.14904746 99 chrX 251271 60 75M = 251317 121 =========================================================================== @@HHEGIIGAGIIIBGIIG@FECHDGHCHHHGHHFFFFFDEACC@ ## @ERR188273.14904746 ## GAAAAATGGGCCCAGGGGACCGGCGCTCAGCATACAGAGGACCCGCGCCGGCACCTGCCTCTGAGTTCCCTTAGT ## + ## @@HHEGIIGAGIIIBGIIG@FECH eg/ERR188273_chrX_rand.bam","title":"Random subsampling of BAM file"},{"location":"#count-number-of-reads","text":"Use samtools idxstats to print stats on a BAM file; this requires an index file which is created by running samtools index . The output of idxstats is a file with four tab-delimited columns: Reference name Sequence length of reference Number of mapped reads Number of unmapped reads samtools idxstats eg/ERR188273_chrX.bam ## chrX 156040895 1126961 45067 ## * 0 0 4332 We can use this with awk to calculate: The number of mapped reads by summing the third column. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3} END {print s}' ## 1126961 The number of reads, which is the sum of mapped and unmapped reads. samtools idxstats eg/ERR188273_chrX.bam | awk '{s+=$3+$4} END {print s}' ## 1176360","title":"Count number of reads"},{"location":"#obtaining-genomic-sequence","text":"Use faidx to fetch genomic sequence; coordinates are 1-based. We need to first index the reference FASTA file that was used to map the reads. samtools faidx genome/chrX.fa Now we can obtain the sequence. samtools faidx genome/chrX.fa chrX:300000-300100 ## >chrX:300000-300100 ## ctgagatcgtgccactgcactccagcctgggcgacagagcgagactccatctcaaaaaaa ## aaaaaaaaaaaaaagaTggggtctctctatgttggccaggt","title":"Obtaining genomic sequence"},{"location":"#comparing-bam-files","text":"The output from mpileup can be used to compare BAM files. The commands below generates alignments using bwa and minimap2 . len=100 n=10000 m=300 script/generate_random_seq.pl 30 1000000 1984 > test_ref.fa script/random_paired_end.pl -f test_ref.fa -l ${len} -n ${n} -m ${m} bwa index test_ref.fa 2> /dev/null bwa mem test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_bwa.bam minimap2 -ax sr test_ref.fa l${len}_n${n}_d${m}_1984_1.fq.gz l${len}_n${n}_d${m}_1984_2.fq.gz 2> /dev/null | samtools sort - -o aln_mm.bam The BAM files can be used with mpileup to compare the depths. samtools mpileup -s -f test_ref.fa aln_bwa.bam aln_mm.bam | head -20 ## [mpileup] 2 samples in 2 input files ## 1 4020 C 1 ^]. E ] 1 ^]. E ] ## 1 4021 T 1 . J ] 1 . J ] ## 1 4022 C 1 . J ] 1 . J ] ## 1 4023 T 1 . J ] 1 . J ] ## 1 4024 G 1 . J ] 1 . J ] ## 1 4025 T 1 . J ] 1 . J ] ## 1 4026 T 1 . J ] 1 . J ] ## 1 4027 A 1 . J ] 1 . J ] ## 1 4028 T 1 . J ] 1 . J ] ## 1 4029 A 1 . J ] 1 . J ] ## 1 4030 G 1 . J ] 1 . J ] ## 1 4031 C 1 . J ] 1 . J ] ## 1 4032 G 1 . J ] 1 . J ] ## 1 4033 G 1 . J ] 1 . J ] ## 1 4034 G 1 . J ] 1 . J ] ## 1 4035 A 1 . J ] 1 . J ] ## 1 4036 T 1 . J ] 1 . J ] ## 1 4037 T 1 . J ] 1 . J ] ## 1 4038 A 1 . J ] 1 . J ] ## 1 4039 C 1 . J ] 1 . J ] Another approach is to use deepTools and the bamCompare command. The bigWig output file shows the ratio of reads between b1 and b2 in 50 bp (default) windows.","title":"Comparing BAM files"},{"location":"#converting-reference-names","text":"One of the most annoying bioinformatics problems is the use of different chromosome names, e.g. chr1 vs 1, in different references even when the sequences are identical. The GRCh38 reference downloaded from Ensembl has chromosome names without the chr : >1 dna:chromosome chromosome:GRCh38:1:1:248956422:1 REF Whereas the reference names from UCSC has the chr : >chr1 AC:CM000663.2 gi:568336023 LN:248956422 rl:Chromosome M5:6aef897c3d6ff0c78aff06ac189178dd AS:GRCh38 Luckily you can change the reference names using samtools reheader but just make sure your reference sequences are actually identical. samtools view eg/ERR188273_chrX.bam | head -2 ## ERR188273.4711308 73 chrX 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 chrX 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP View header samtools view -H eg/ERR188273_chrX.bam ## @HD VN:1.0 SO:coordinate ## @SQ SN:chrX LN:156040895 ## @PG ID:hisat2 PN:hisat2 VN:2.2.0 CL:\"/Users/dtang/github/rnaseq/hisat2/../src/hisat2-2.2.0/hisat2-align-s --wrapper basic-0 --dta -p 4 -x ../raw/chrX_data/indexes/chrX_tran -1 /tmp/4195.inpipe1 -2 /tmp/4195.inpipe2\" ## @PG ID:samtools PN:samtools PP:hisat2 VN:1.19.2 CL:samtools view -H eg/ERR188273_chrX.bam Substitute header with new name. samtools view -H eg/ERR188273_chrX.bam | sed 's/SN:chrX/SN:X/' > eg/my_header Save bam file with new ref and check it out. samtools reheader eg/my_header eg/ERR188273_chrX.bam > eg/ERR188273_X.bam samtools view eg/ERR188273_X.bam | head -2 ## ERR188273.4711308 73 X 21649 0 5S70M = 21649 0 CGGGTGATCACGAGGTCAGGAGATCAAGACCATCCTGGCCAACACAGTGAAACCCCATCTCTACTAAAAATACAA @@@F=DDFFHGHBHIFFHIGGIFGEGHFHIGIGIFIIIGIGIGGDHIIGIIC@>DGHCHHHGHHFFFFFDEACC@ AS:i:-5 ZS:i:-5 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:70 YT:Z:UP NH:i:2 ## ERR188273.4711308 133 X 21649 0 * = 21649 0 CTACAGGTGCCCGCCACCATGCCCAGCTAATTTTTTTTGTATTTTTAGTAGAGATGGGGTTTCACTGTGTTGGCC CB@FDFFFHHGFHIJJJJIIIIIIIGGGIJGIIJJJJJJFFHIIIIGECHEHHGGHHFF?AACCDDDDDDDDBCD YT:Z:UP","title":"Converting reference names"},{"location":"#coverage","text":"Coverage can mean the: average depth of each covered base percentage of bases covered samtools depth and samtools mpileup can be used to indicate the depth of each covered base (and used to calculate the average depth. samtools coverage will provide both the average depth and percentage of bases covered per chromosome/reference sequence. samtools depth will return three columns: reference, position, and coverage. samtools depth -@ 4 eg/ERR188273_chrX.bam > ERR188273_depth.tsv head ERR188273_depth.tsv ## chrX 21649 1 ## chrX 21650 1 ## chrX 21651 1 ## chrX 21652 1 ## chrX 21653 1 ## chrX 21654 1 ## chrX 21655 1 ## chrX 21656 1 ## chrX 21657 1 ## chrX 21658 1 The average depth can be calculated by summing the third column and dividing by the total number of bases (be sure to use -a with samtools depth as that will output all positions including zero depth). samtools depth -@ 4 -a eg/ERR188273_chrX.bam | perl -ane '$t += $F[2]; END {$cov = $t / $.; printf \"Bases covered:\\t%.3f\\nCoverage:\\t%.3f\\n\", $., $cov}' ## Bases covered: 156040895.000 ## Coverage: 0.532 The samtools mpileup command also provides depth information (but not for reads that have a mapping quality of 0, by default) with some additional information: Sequence name 1-based coordinate Reference base (when used with -f ) Number of reads covering this position Read bases Base qualities Alignment mapping qualities (when used with -s ) samtools mpileup -f genome/chrX.fa -s eg/ERR188273_chrX.bam > ERR188273_mpileup.tsv head ERR188273_mpileup.tsv ## [mpileup] 1 samples in 1 input files ## chrX 251271 g 1 ^]. @ ] ## chrX 251272 a 1 . @ ] ## chrX 251273 a 1 . < ] ## chrX 251274 a 1 . D ] ## chrX 251275 a 1 . D ] ## chrX 251276 a 1 . D ] ## chrX 251277 t 1 . D ] ## chrX 251278 g 1 . D ] ## chrX 251279 g 1 . F ] ## chrX 251280 g 1 . B ] Note that the start of the samtools mpileup output differ from the start of the samtools depth output. This is because mpileup performs some filtering by default. In the case of this example, read pairs that are not both mapped will be ignored. To count these \u201corphan\u201d reads, use the --count-orphans argument. samtools mpileup -f genome/chrX.fa --count-orphans -s eg/ERR188273_chrX.bam > ERR188273_mpileup_orphans.tsv head ERR188273_mpileup_orphans.tsv ## [mpileup] 1 samples in 1 input files ## chrX 21649 g 0 * * * ## chrX 21650 a 1 . D ! ## chrX 21651 t 1 . F ! ## chrX 21652 c 1 . F ! ## chrX 21653 a 1 . H ! ## chrX 21654 c 1 . G ! ## chrX 21655 g 1 . H ! ## chrX 21656 a 1 . B ! ## chrX 21657 g 1 . H ! ## chrX 21658 g 1 . I ! In addition mpileup performs \u201cper-Base Alignment Quality\u201d (BAQ) by default and will adjust base quality scores. The default behaviour to to skip bases with baseQ/BAQ smaller than 13. If you are finding discrepancies between mpileup \u2019s coverage calculation with another coverage tool, you can either set --min-BQ to 0 or use --no-BAQ to disable BAQ. I have an old blog post on using mpileup . samtools coverage will provide the following coverage statistics: rname - Reference name / chromosome startpos - Start position endpos - End position (or sequence length) numreads - Number reads aligned to the region (after filtering) covbases - Number of covered bases with depth >= 1 coverage - Proportion of covered bases [0..1] meandepth - Mean depth of coverage meanbaseq - Mean base quality in covered region meanmapq - Mean mapping quality of selected reads samtools coverage eg/ERR188273_chrX.bam ## #rname startpos endpos numreads covbases coverage meandepth meanbaseq meanmapq ## chrX 1 156040895 1110685 3402037 2.18022 0.532299 36.3 59.4 The example BAM file only contains reads for chrX hence the statistics are only returned for chrX . Returning to our coverage definition at the start of this section: average depth of each covered base = meandepth percentage of bases covered = covbases The mosdepth tool can also calculate depth (and much faster than samtools depth ) per base or within a given window. The output is given in a BED file, where the fourth column indicates the coverage. mosdepth ERR188275 eg/ERR188273_chrX.bam gunzip -c ERR188273.per-base.bed.gz | head ## gzip: ERR188273.per-base.bed.gz: No such file or directory Coverage in using a 500 bp window. mosdepth -n --fast-mode --by 500 ERR188275_500 eg/ERR188273_chrX.bam gunzip -c ERR188275_500.regions.bed.gz | head ## chrX 0 500 0.00 ## chrX 500 1000 0.00 ## chrX 1000 1500 0.00 ## chrX 1500 2000 0.00 ## chrX 2000 2500 0.00 ## chrX 2500 3000 0.00 ## chrX 3000 3500 0.00 ## chrX 3500 4000 0.00 ## chrX 4000 4500 0.00 ## chrX 4500 5000 0.00","title":"Coverage"},{"location":"#stargazers-over-time","text":"","title":"Stargazers over time"}]}
\ No newline at end of file
diff --git a/sitemap.xml.gz b/sitemap.xml.gz
index e12acca..30685f5 100644
Binary files a/sitemap.xml.gz and b/sitemap.xml.gz differ