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
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: ...
3
votes
1answer
446 views

How to get a MSA fasta from BAM/SAM?

I am working with paired end NGS miseq data from a viral genome with multiple timepoints. I have filtered and trimmed this data for quality and adapter sequences. I have then merged the filtered ...
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 ...
3
votes
2answers
3k views

Error while calling bcftools mpileup - Failed to open -: unknown file type

I have sequenced a bacterial genome with a GridIon from ONT. Basically what I want to check is whether or not trimming 50 bps at the beginning of the reads will improve alignment against the reference ...
2
votes
3answers
818 views

How to reverse complement the DNA sequences for given inverse/reverse coordinates?

I have the series of coordinates in id.txt file, whose coordinates sequences are in genome.fasta file. The coordinates of id.txt ...
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 ...
0
votes
1answer
1k views

samtools view: writing to standard output failed: Broken pipe

I have the following alignment script which uses BWA: ...
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 ...
3
votes
1answer
1k views

Running htseq-count over BAM files

I am trying to derive an expression matrix from BAM files using htseq-count. These are bulk RNASeq BAM's by the way. I have read the ...
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 ...
1
vote
2answers
85 views

Using preprocessing/alignment functions on the server

I am new to bash and the processes behind cluster computing in general and need some help with understanding some basics. After looking all over the internet and this forum (+ askUbuntu) I found ...
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 ...
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 ...
2
votes
2answers
1k views

How to count the number of mapped read in 100-bp window from a BAM/SAM file

Although I know how to get total number of mapped read using samtools flagstat (samtools flagstat file_sorted.bam) but I want to count total number of mapped read ...
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: ...
1
vote
0answers
241 views

bedtools single-nucleotide coverage in BED-specified regions for multiple BAMs

I have a BED6 (BED + name, score, strand info) file that defines some regions of interest. I also have a set of BAMs corresponding to different samples. I would like to obtain output similar to <...
3
votes
1answer
430 views

Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?

After alignment using BWA, I have removed the dupliment using the samtools(Version: 1.9). My procedure is as follows: ...
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: ...
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. ...
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 ...
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 ...
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 ...
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?
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. ...
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
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 ...
2
votes
1answer
99 views

SAM format: Does the BAM "Integer or numeric array" field no longer exist? Why?

The SAM file format specification makes clear that optional fields must follow the format TAG:TYPE:VALUE ...
1
vote
1answer
321 views

running a samtools command for multiple bam files from 1000 genomes project

I want to compute the depth of coverage only for specific intervals in phase 3, 1000 genomes project.I have not worked with 1000 genomes project before, so a bit unfamiliar with it. I do not want to ...
2
votes
0answers
1k views

Which R package to use for differential analysis with TPM values?

I'm using hisat2, stringtie tools for the RNA-Seq analysis. After stringtie using ballgown I get FPKM and TPM values for every gene. I have seen that edgeR, Deseq2 can be used for Counts data. I ...
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: ...
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 ...
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 ...
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 ...
2
votes
1answer
916 views

Modifying a SAM header after removing all non-primary reads

I subset a BAM to only include primary reads using the following samtools commands: samtools view -F 256 input.bam > input.primaryOnly.sam Now, in order to ...
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
3answers
231 views

Parsing SAM/BAM files for additional information

I used BWA-MEM to alignment and I would like to gather the some informations like total % of match, mismatch, insert/delete etc. I am wondering if there is any existing tools that produces this ...
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 ...
1
vote
2answers
516 views

Mapping Information from SAM/BAM file

I mapped raw illumina reads to longer pacbio reads and I would like to know the following information from my mapping file (SAM/BAM) How many PacBio reads are mapped to at least one illumina reads A ...
0
votes
1answer
176 views
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: ...
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 ...
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
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). ...
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 ...
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:...
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. ...
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. ...
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/...
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 ...