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
2
votes
1answer
34 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 ...
0
votes
1answer
22 views

How to convert a Pileup file to VCF format with Hg19 alignment

I would like to know how to convert a pileup file for an adna y chromosome to vcf with Hg19 alignment.
2
votes
2answers
88 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 ...
1
vote
2answers
61 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
4answers
58 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
28 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. ...
5
votes
2answers
121 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
83 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
293 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
324 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
254 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
52 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
156 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
333 views

samtools view: writing to standard output failed: Broken pipe

I have the following alignment script which uses BWA: ...
4
votes
1answer
81 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
343 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
737 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
65 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
70 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
69 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
561 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
2k 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
127 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
155 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
988 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
48 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
699 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
147 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
97 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
59 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
140 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
87 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
74 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
124 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
539 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
259 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
408 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
358 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
136 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
195 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
616 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
182 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
155 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 ...
5
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
381 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
115 views
5
votes
1answer
651 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
778 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
795 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 ...