SAMBAMBA-high performance parallel tool for working with SAM and BAM files

 https://github.com/biod/sambamba

 

SAMBAMBA-high performance highly parallel robust and fast tool for working with SAM and BAM files 

 

function: view, index, sort, markdup, and depth.

flagstat : 1.4x faster than samtools. 

Index : similar. 

Markdup : ~ 6x faster 

View :      ~4x faster

Sort  : Sambamba has been beaten, though sambamba is notably up to 2x faster than samtools on large RAM machines (120GB+).

[bam_rmdup_core] inconsistent BAM file for pair

 https://github.com/samtools/samtools/issues/359

rmdup with no option (neither -s nor -S) expects alignment pairs with alternating +/- TLEN values. It is a special-purpose tool for removing PCR duplicates which conform to that pattern.

The warning is issued for consecutive positive TLEN values. A simple contrived example:

cat rmdupExample.sam
@SQ SN:c2   LN:9
s1  67  c2  1   0   9M  =   1   30  CTAATAATC   XXXXXXXXX   RG:Z:1st
s1  67  c2  1   0   9M  =   1   -30 CTAATAATC   XXXXXXXXX   RG:Z:1st
s1  259 c2  1   0   9M  =   1   30  CTAATAATC   YXXXXXXXX   RG:Z:2nd
s1  131 c2  1   0   9M  =   1   30  CTAATAATC   XXXXXXXXX   RG:Z:3rd
s1  131 c2  1   0   9M  =   1   -30 CTAATAATC   XXXXXXXXX   RG:Z:3rd

./samtools view -b -o /tmp/rmdupExample.bam /tmp/rmdupExample.sam
../samtools-0.1.19/samtools rmdup /tmp/rmdupExample.bam /tmp/test_input_1_rmdup.bam
[bam_rmdup_core] processing reference c2...
[bam_rmdup_core] inconsistent BAM file for pair 's1'. Continue anyway.
[bam_rmdup_core] 2 / 3 = 0.6667 in library '    '
./samtools view /tmp/test_input_1_rmdup.bam
s1  67  c2  1   0   9M  =   1   -30 CTAATAATC   XXXXXXXXX   RG:Z:1st
s1  259 c2  1   0   9M  =   1   30  CTAATAATC   YXXXXXXXX   RG:Z:2nd

changing chromosome notation in .BAM file

http://seqanswers.com/forums/showthread.php?t=22504


https://josephcckuo.wordpress.com/2016/11/17/modify-chromosome-notation-in-bam-file/


change chromosome notation "1" to "chr1"

for file in *.bam; do filename=`echo $file | cut -d "." -f 1`; samtools view -H $file | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - $file > ${filename}_chr.bam; done


for loop

for file in *.bam; do filename=`echo $file | cut -d "." -f 1`; samtools view -H $file | sed -e 's/SN:\([0-9XY]\)/SN:chr\1/' -e 's/SN:MT/SN:chrM/' | samtools reheader - $file > ${filename}_chr.bam; done

Extract Reads From a Bam File That Fall Within A Given Region

https://www.biostars.org/p/48719/


# -h : include header
# file should be indexed.bam

samtools view -h input.indexed.bam "Chr1:10000-20000" > output.sam

samtools rmdup vs Picard MarkDuplicates

samtools rmdup PE
http://seqanswers.com/forums/showthread.php?t=5959

If you have one pair of reads where read 1 starts at position 100, and the other end starts at position 200, and a second pair of reads where read 1 starts at position 100, and read 2 starts at position 250, those came from different fragments of DNA. You can tell because the read 2 start is different, even though the read 1 start is the same.

When treating the reads as paired end, none of those reads should be deleted as PCR duplicates.

However, if you ran rmdup -S, the software will not check to see if read 2 has a different start coordinate, so one of those read 1 reads will be treated as a duplicate, and deleted.


Question: Samtools Dedup Documentation
 https://www.biostars.org/p/55111/

rmdup for PE reads is pretty straightforward. It looks for identical external coordinates, meaning it only looks at the 5' start coordinates of the FR orientation pair-reads. Then it takes the pair with the highest mapping quality.
For SE reads, I've read that samtools also only looks for identical 5' start coordinates, not both start and end coordinates. I think the idea is that sequencers usually fall in quality towards the 3'. After mapping, duplicate reads have higher chance of mapping differentially towards the 3' end. So it only looks at the adapter trimmed 5' start for duplicates.


Question: Picard MarkDuplicates and SamTools rmdup algorithm documentation
https://www.biostars.org/p/105291/
SamTools rmdup 'only' compare two reads on chrom and pos (which could be wrong if two reads come from two different libraries) and **removes** reads from the BAM: information is lost.

picard set the sam flag 1024 but do not delete the reads. two pairs of reads are compared , as far as I know, using the chrom, the pos, the group-id (sample...) + (flowcell , lane, X,Y for optical dups) (,and the cigar string ?).


Picard MarkDuplicates vs samtools rmdup for variant calling with GATK
https://gatkforums.broadinstitute.org/gatk/discussion/6793/picard-markduplicates-vs-samtools-rmdup-for-variant-calling-with-gatk 
 - this post explains Picard's duplicate marking tools