Questions tagged [samtools]

Questions specific to interacting with and post-processing sequence alignments using the SAMtools package. For questions specifically about formats SAM, BAM or CRAM use tags sam or bam,

Filter by
Sorted by
Tagged with
15
votes
2answers
7k views

Merge hundreds of small BAM files into a single BAM file

I am working with over a million (long) reads, and aligning them to a large genome. I am considering running my alignment jobs in parallel, distributing horizontally across hundreds of nodes rather ...
14
votes
1answer
1k views

Are soft-clipped bases used for variant calling in samtools + bcftools?

If there are soft clipped base pairs specified in the CIGAR string for a read in a SAM/BAM file, will these be used for variant calling in a samtools + bcftools workflow? The GATK HaplotypeCaller, ...
11
votes
3answers
6k views

Converting a VCF into a FASTA given a reference with Python, R

I am interested in converting a VCF file into a FASTA file given a reference sequence with Python or R. Samtools/BCFtools (Heng Li) provides a Perl script ...
11
votes
1answer
187 views

Can I create a CRAM file with a relative reference path?

I’m trying to create a CRAM file that stores its path to the FASTA reference as a relative path, rather than an absolute path, so that I can move the files around. Unfortunately I can’t get this to ...
10
votes
1answer
821 views

Quantifying reads mapping to multiple loci

I have been using STAR for our RNA-Seq samples. The final.out log file reports percentage of uniquely mapped reads along with percentage of reads that map to ...
9
votes
2answers
4k views

samtools depth print out all positions

I am trying to use samtools depth (v1.4) with the -a option and a bed file listing the human chromosomes chr1-chr22, chrX, chrY, and chrM to print out the coverage ...
9
votes
1answer
115 views

Behavior of `--reference` flag with samtools sort

I'm about to run samtools sort (version 1.6) and I'm reviewing the available configuration options. ...
8
votes
2answers
430 views

Convert bam file to highly compressible bam

I have a large collection of bam files and I want to post-process each of them into another bam where I can make queries about: the reads position and pair-endness, insert sizes, MAPQ and other ...
8
votes
1answer
1k views

What are all the reference files produced by bwa index, and are these dependent upon whether the reference is zipped?

I have indexed a gzipped reference with bwa: bwa index reference.fa.gz, which produces a series of other files ...
7
votes
7answers
6k views

How to subset a BAM by a list of QNAMEs?

I have a text file 'qnames.txt' with QNAMEs in the following format: EXAMPLE:QNAME1 EXAMPLE:QNAME2 EXAMPLE:QNAME3 EXAMPLE:QNAME4 EXAMPLE:QNAME5 I would like to ...
7
votes
3answers
618 views

Is there a way to retrieve several SAM fields faster than `samtools view | cut -f`?

I am constructing a bit of software which pipes the outputs of the bam file via samtools view into a script for parsing. My goal is to (somehow) make this process ...
7
votes
3answers
3k views

How to check whether all BAM read contain defined read groups?

I'm trying to investigate whether there are errors within my BAM. After looking at the BAM header to see whether the read groups exists (using samtools, i.e. ...
7
votes
2answers
277 views

variant calling on ChIP-seq style data: samtools mpileup with minimal filters

I am running samtools mpileup (v1.4) on a bam file with very choppy coverage (ChIP-seq style data). I want to get a first-pass list of positions with SNVs and their frequency as reported by the read ...
7
votes
4answers
2k views

Double-counting coverage of overlapped read pairs

EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region. I used samtools depth to calculate the depth of ...
7
votes
2answers
467 views

How does htslib/samtools access optional BAM fields?

I am confused regarding how various software libraries deal with optional fields in a BAM: Based upon the BAM specification, there are 11 mandatory fields to a BAM: ...
7
votes
1answer
243 views

scanBam from Rsamtools is not importing one of my reads into R

I have this read in my BAM file. It maps on chromosome 1. I open this BAM file in IGV, and I can see the alignment on chromosome 1. But when I open this file in R with Rsamtools: ...
7
votes
1answer
145 views

How can I speed up INDEL calling/correction on BAM files?

The samtools mpileup command has quite a neat feature that it is able to correct mapping errors associated with misalignment of INDELs. By default, the ...
7
votes
1answer
357 views

samtools mpileup empty when filtering out flags

I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem. I am heavily filtering reads that are not paired and of a certain orientation. Applying the filtering, I ...
6
votes
5answers
1k views

Subset smaller BAM to contain several thousand rows from multiple chromosomes

There are many cases whereby I would like to subset a BAM to create a small file in order to work with (e.g. algorithmic testing, debugging, etc.) Normally I do the following, which will subset the ...
6
votes
1answer
1k views

Read counts from BAM file

I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for ...
6
votes
1answer
1k views

Convert BAM to properly paired FASTQ files

I thought I had figured this one out. But apparently not. I need to convert a BAM file of paired-end alignments to two FASTQ files of paired reads to realign them, with a twist: I only want reads ...
6
votes
1answer
1k views

How to remove all BAM read groups from all reads (not just the header)?

I have problem with one my BAMs---it appears to have invalid read groups. Normally when I have such a problem, I remove all the read groups from the BAM header as follows: ...
6
votes
1answer
557 views

How to build a BAM header file with htslib in C++?

I'd like to use C++ to generate a new BAM file programmatically. This is an example how to use htslib to generate a new BCF file on the fly. https://github.com/samtools/htslib/blob/develop/test/...
6
votes
1answer
550 views

low-memory high-speed replacement for Picard MarkDuplicates

I am running Picard MarkDuplicates with the following parameters below. On the file described, it takes about 41.6Gb of RAM memory and about 20-25 minutes to compute (only uses 1 core AFAICS). ...
5
votes
3answers
3k views

How to get fasta alignment file from SAM/BAM file?

I am not talking about consensus sequence, I know how to get consensus sequence using mpileup in samtools/bcftools. As I understand , SAM/BAM files are basically sequence alignment format so it's ...
5
votes
3answers
272 views

Is there an efficient way to check an input BAM in R?

I'm writing a function in R for an R package which takes as input a BAM. ...
5
votes
2answers
255 views

samtools mpileup skipping read

I run the command: samtools mpileup -O -s -q20-B -Q20 -f hg19.fa -r chr1:569929-569931 myFile.bam and get: ...
5
votes
1answer
155 views

Does picard markduplicate toggle PCR duplicate samflag

I have a RNA-seq bam file and there are few reads that are puzzling me. According to the bam header, this bam file is sorted by coordinate, created using tophat and markduplicate step is not done. ...
5
votes
2answers
5k views

Difference between samtools mark duplicates and samtools remove duplicates?

What is difference between samtools mark duplicates and remove duplicates ? Is it necessary to mark duplicates before removing duplicates with samtools?
5
votes
2answers
130 views

Control width of samtools tview "snapshot" when redirecting

I have a large number of loci I would like to examine manually with samtools tview. Rather than typing or copy-n-pasting dozens of coordinates, I was hoping to ...
5
votes
1answer
157 views

Samtools/bcftools calling indels in noisy reads

I have very noisy nanopore reads and am trying to call SNPs/indels. I'm runinig into some trouble when truing to use the samtools mpileup | bcftools call ...
4
votes
5answers
3k views

using python to write bioinformatics pipelines tutorial

I was wondering if there is a tutorial or a small code snippet to understand how to write bioinformatics pipeline using python, for example use a aligner (say hisat) get the output and process it ...
4
votes
1answer
89 views

What is "block-compressed" file in samtools?

SAMtools returned an error message for my gzipped genome FASTA: Indexed block-compressed FASTA file cannot be handled The source code for the message is here. ...
4
votes
2answers
127 views

sum of products of mapq and mapped bases for each read in a from a BAM file

Given a BAM file, I'd like to calculate the sum, over all reads, of the mapping quality and the number of mapped bases (i.e. number of M's in the CIGAR string). For example, given two reads like this:...
4
votes
1answer
156 views

Is there a tool that can perform a read-group-aware mpileup from a single file?

I would like to perform a samtools mpileup from a single file that contains thousands of read groups with different SM tags. I could split the bam by read group ...
4
votes
1answer
385 views

Error given while trying to index a BAM file with Samtools Index - NO COOR?

I am currently working on my own Metagenomic pipeline, utilizing Bowtie 2 to map. Bowtie 2 outputs a SAM file, which I convert to a .BAM and sort it using Samtools. When I try to utilize Samtools to ...
4
votes
1answer
209 views

Displaying soft-clipped nucleotides in samtools tview

There are nicer genomics visualization tools available, but the samtools tview command is almost always my go-to for a quick first look at read alignments. I just ...
4
votes
0answers
137 views

convert supplementary reads to primary in sam or bam

I have a problem with a tool, that possibly ignores supplementary reads. I want to find out a bit how this tool works by converting all supplementary reads to primary reads, and then change the names ...
3
votes
5answers
4k views

How to create a .bed file from .fasta? [duplicate]

I have some problems in creating a .bed file for hg19, so I will be able to visualize the .bed file in IGV. The .fasta file contains rows of this form: ...
3
votes
2answers
260 views

Frequency of specific viral sequence in .BAM or .fastq

I was wondering if anyone had any experience 'counting' the frequency of a particular viral sequence in an individual's fastq sequence. So basically, counting the number of occurrences of a particular ...
3
votes
1answer
4k views

Extracting all reads from bam file which match read IDs in another file

I have a long list of read IDs of interest to me in a file called read_names.txt. it is simply in the format: ...
3
votes
1answer
243 views

pysam or piping samtools view to a python script

Is it faster to use the PySam package to run a python script on a bam for read in samfile.fetch('chr1', 100, 120): print read compared to using a pipe and ...
3
votes
2answers
293 views

What is a samtools mpileup reference skip?

The samtools documentation for mpileup states: At this column, a dot stands for a match to the reference base on the forward strand, a comma for a match on the reverse strand, a '>' or '<' for a ...
3
votes
3answers
83 views

Unable to open .bam file in C++ using SeqAn due to 'seqan::UnknownExtensionError'

I am trying to open .bam files in C++ to extract reads occurring at specific scaffolds and loci. I essentially want to call "samtools view sample.bam -o sample.sam scaffold:pos-pos" from C++....
3
votes
1answer
273 views

What is the standard way to measure contig sequence lengths in a BAM?

What is the standard way to measure contig sequence lengths in a BAM? My understanding is that the community would use samtools idxstats to compute this ...
3
votes
1answer
2k views

sorting BAM file error using samtools

I have few bam files and would like to get read counts using samtools idxstats [Data is aligned to hg19 transcriptome]. To use that command I need a sorted bam ...
3
votes
1answer
360 views

How to interpret methylation calls from Bismark on opposite strands?

I'm looking at Reduced representation bisulfite sequencing (RRBS) data from ENCODE, and to align the FASTQ files I've used Bismark with Bowtie 1. When I load the resulting BAM file into ...
3
votes
1answer
295 views

How to generate two fasta files from samtools phase output?

samtools phase is an easy-to-use tool to phase variants from a bam file and a reference. For a sequence that has two haplotypes, like this: ...
3
votes
2answers
171 views

Filter bam using SNP list in bed format with minimum mapping quality and base quality

I have a bam file and a bed file that defines a list of SNPs. I would like to filter the bam file to contain only those reads with a minimum mapping quality that overlap at least one SNP with a ...
3
votes
2answers
359 views

Parse BAM by insertion size and get genomic coordinates

I would like to find out where (in genomic coordinates) large insertions are found within a given BAM file, file.bam. (In terms of genomic coordinates, I just mean ...