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 ...
user avatar
  • 19.3k
14 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 ...
user avatar
  • 19.3k
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 ...
user avatar
  • 758
12 votes

Convert bam file to highly compressible bam

You can just set the fields you don't need to *: ...
user avatar
  • 19.3k
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: ...
user avatar
  • 19.3k
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 ...
user avatar
  • 1,473
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
user avatar
  • 2,566
9 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 ...
user avatar
  • 3,061
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 ...
user avatar
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 ...
user avatar
  • 5,645
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 ...
user avatar
  • 3,521
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: ...
user avatar
  • 2,099
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 ...
user avatar
  • 3,201
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, ...
user avatar
  • 5,645
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 ...
user avatar
  • 11.5k
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. ...
user avatar
6 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 ...
user avatar
  • 19.3k
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 ...
user avatar
  • 11.5k
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 ...
user avatar
  • 19.3k
6 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: ...
user avatar
  • 19.3k
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 ...
user avatar
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 ...
user avatar
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 ...
user avatar
6 votes
Accepted

Marking optical or PCR duplicates with picard vs. samtools flagstat

There are duplicates, in this line: 1636809 + 0 duplicates, gives 1636809/26595942 = 0.06154356 According to samtools documentation for flagstat: Provides counts for each of 13 categories based ...
user avatar
  • 1,638
6 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 ...
user avatar
  • 443
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 ...
user avatar
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, ...
user avatar
  • 5,645
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 ...
user avatar
  • 2,566
5 votes
Accepted

How to subset a BAM by a list of QNAMEs?

(1) It won't be super fast but you can provide grep with a file of QNAMES. samtools view file.bam | grep -f 'qnames.txt > subset.sam where qnames.txt has <...
user avatar
  • 2,566
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. ...
user avatar
  • 11.5k

Only top scored, non community-wiki answers of a minimum length are eligible