Episode #125 of the Stack Overflow podcast is here. We talk Tilde Club and mechanical keyboards. Listen now

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
1
vote
2answers
33 views

bedtools intersect on very large .bam file - -sorted confusion

I need to do a bedtools intersect operation on a very large .bam file. When I use the standard bedtools intersect operation, the process consumes all the memory on my system. From the bedtools ...
0
votes
3answers
44 views

Find indels between two short sequences

I have two sequences, say AAAGCTCGAGG and AAAGCGAGG. I need a convenient tool which shows me insertions or deletions between these, i.e. in this case something like ...
1
vote
1answer
24 views

How to take into account alternative bwa mem mapping when computing coverage

when mapping short reads with bwa mem, if a read has alternative mapping positions they are reported by bwa mem in the X0 and XA tags. Now, let's say I want to compute the coverage of my bam file. ...
4
votes
2answers
101 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: ...
2
votes
0answers
67 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
205 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 ...
2
votes
2answers
176 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
184 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 ...
0
votes
0answers
49 views

Reconstructing reference FASTA from CRAM file

I would like to reconstruct (or download) the original FASTA reference file that was used to create a CRAM file in my possession. I realize that I can download the contigs specified in the CRAM header ...
2
votes
2answers
145 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
230 views

samtools view: writing to standard output failed: Broken pipe

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

Running htseq-count over BAM files

I am trying to derive an expression matrix from BAM files using htseq-count on the server. These are bulk RNASeq BAM's by the way. I have read the ...
4
votes
3answers
607 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
63 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 ...
3
votes
0answers
65 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
65 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
464 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 ...
2
votes
1answer
1k 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
114 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
123 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: ...
2
votes
5answers
803 views

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

I have some problems with creating a .bed file for hg19 which I will be able to load on IGV. The fasta file contains rows of this form: ...
4
votes
1answer
46 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. ...
6
votes
4answers
619 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 ...
3
votes
0answers
124 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
86 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
2k 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
55 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
115 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. ...
4
votes
0answers
79 views

Samtools/bcftools calling indels in noisy reads

I have very noisy reads and am trying to call SNPs/indels. I'm runinig into some trouble when truing to use the samtools mpileup | bcftools call combination. It ...
2
votes
1answer
70 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
104 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
479 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 ...
6
votes
2answers
238 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
2answers
404 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
322 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
123 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
174 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
560 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
170 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
144 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 ...
4
votes
6answers
2k 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
358 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
110 views
5
votes
1answer
628 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
695 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 ...
4
votes
1answer
738 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
327 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). ...
10
votes
3answers
3k 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
89 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:...