15
votes
Accepted
Are soft-clipped bases used for variant calling in samtools + bcftools?
No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or ...
- 19.6k
15
votes
Accepted
Merge hundreds of small BAM files into a single BAM file
samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not ...
- 19.6k
13
votes
using python to write bioinformatics pipelines tutorial
Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about ...
- 788
12
votes
Convert bam file to highly compressible bam
You can just set the fields you don't need to *:
...
- 19.6k
10
votes
Accepted
How to count the number of mapped read in 100-bp window from a BAM/SAM file
one-liner
Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ...
- 3,141
9
votes
Accepted
Quantifying reads mapping to multiple loci
You almost had the correct python code already, you just need to filter out secondary alignments:
...
- 19.6k
9
votes
How to subset a BAM by a list of QNAMEs?
use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads
READ_LIST_FILE (File) Read List File containing reads that will be included or excluded ...
- 1,501
9
votes
Accepted
Extracting all reads from bam file which match read IDs in another file
It is still slow but grep has a -f option to take in a file
samtools view inbam.bam | grep -f read_names.txt > read_locs.txt
- 2,574
8
votes
Merge hundreds of small BAM files into a single BAM file
Merging sorted files is a linear operation, so any well-implemented tools that do it will do it with approximately the same efficiency. So samtools merge (use the ...
- 956
8
votes
How to get fasta alignment file from SAM/BAM file?
I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align ...
- 6,485
8
votes
Samtools sort: most efficient memory and thread settings for many samples on a cluster
I usually use samtools sort -@4 -m4g. Samtools will be faster with more threads and more RAM, but I feel the gain is not significant. I haven't done a proper ...
- 6,485
7
votes
Accepted
Read counts from BAM file
Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR_* or NM_* accession codes with your ...
- 3,571
7
votes
Accepted
variant calling on ChIP-seq style data: samtools mpileup with minimal filters
I used this in the past for ChIP-seq data and it generated SNVs:
...
- 2,169
7
votes
Accepted
What are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped?
You'll get the exact same index (the amb, ann, bwt, pac ...
- 19.6k
7
votes
Accepted
How to remove all BAM read groups from all reads (not just the header)?
You might have to manually strip those auxiliary tags off:
...
- 19.6k
7
votes
How to subset a BAM by a list of QNAMEs?
UPDATE 2021/06/28: since version 1.12, samtools now accepts option -N, which takes a file containing read names of interest.
Using samtools 1.12 or greater:
...
- 175
7
votes
Accepted
What is "block-compressed" file in samtools?
A block compression usually refers to compressing your file into a series of small blocks (with a tool like bgzip). This allows indexing in that the index can record which record lives in which block ...
- 3,301
7
votes
Accepted
How to Sort and Index a SAM file without converting it to BAM?
Not only will you save disk space by converting to BAM, but BAM files are faster to manipulate than SAM. Source: Dave Tang's SAMTools wiki.
sort supports ...
- 613
6
votes
Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?
I modified your original question: as you are extracting 4 fields, you are not outputting BAM. The answer to the modified question is: yes, you can write a C program with htslib (or with bamtools, ...
- 6,485
6
votes
Accepted
Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?
The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a ...
- 13.8k
6
votes
Accepted
samtools depth print out all positions
You might try using bedtools genomecov instead. If you provide the -d option, it reports the coverage at every position in the BAM file.
...
- 540
6
votes
Accepted
How to check whether all BAM read contain defined read groups?
According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a ...
- 13.8k
6
votes
Does picard markduplicate toggle PCR duplicate samflag
Yes, if MarkDuplicates encounters a pair that's marked as a duplicate that it considers (for whatever reason) to not be a duplicate then it will unset the duplicate mark. You can test this yourself by ...
- 19.6k
6
votes
Accepted
Is there an efficient way to check an input BAM in R?
You won’t know whether it’s a valid BAM file until you attempt reading it in full; even Devon’s method (which I recommend) only does a superficial check — i.e. it checks whether (parts of) the BAM ...
- 4,845
6
votes
Frequency of specific viral sequence in .BAM or .fastq
Most read aligners will report unaligned reads as well, which presumably will include your viral sequences. I would ask them to formally confirm that the BAM files will contain unaligned reads before ...
- 5,060
6
votes
Accepted
Convert BAM to properly paired FASTQ files
As noted in the comments, the problem is “some reads fall in the target region but their pairs fall outside it”, leading to non-trivial numbers of singleton reads coming out of ...
- 956
5
votes
Accepted
Converting a VCF into a FASTA given a reference with Python, R
You can whip up something quite easily in Python using pyfaidx, which allows you to access a FASTA file by genomic coordinates. So you'd just have to pull the coordinates from the VCF to get the ...
- 535
5
votes
variant calling on ChIP-seq style data: samtools mpileup with minimal filters
Another approach is htsbox. You can get a candidate list with:
htsbox pileup -Cvcf ref.fa -q20 -Q20 -s5 file.bam > out.vcf
Here, ...
- 6,485
5
votes
using python to write bioinformatics pipelines tutorial
BioPython has some good tools for processing reads and alignments.
http://biopython.org/DIST/docs/tutorial/Tutorial.html
There is a python library wrapping samtools so many of the samtools calls can ...
- 2,574
5
votes
Accepted
Modifying a SAM header after removing all non-primary reads
Your first command only needs a slight modification to add in -h. This will create a SAM file with a header.
...
- 13.8k
Only top scored, non community-wiki answers of a minimum length are eligible