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,
135
questions
1
vote
1
answer
6k
views
Samtools Index: Chromosome Blocks not Continuous
I am working with short-read whole-genome sequences from the NCBI's SRA. I have aligned and sorted all of my short-read sequences and am attempting to index each sequence into .bai format using ...
- 169
2
votes
1
answer
1k
views
How to extract unmatched reads using bwa and samtools?
I have a single read (NOT paired) that I need to pass through the workflow described in Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying ...
1
vote
1
answer
152
views
How can I get unmatched reads for defective genomes analysis using bwa and samtools?
I am trying to follow the Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.
According to ...
2
votes
4
answers
385
views
How to make work programs from the $PATH?
I am trying to analyze my RNAseq reads for defective genomes and I use this program (http://www.di-tector.cyame.eu/) that is suggested by Beauclair et al (https://pubmed.ncbi.nlm.nih.gov/30012569/) ...
4
votes
1
answer
241
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 ...
- 2,166
2
votes
1
answer
70
views
Phasing problem using "samtools phase"
so I'm having a problem using samtools phase. I'm trying to phase my .bam archives that where generated by aligning my samples with my reference genome. Following the documentation I tried the command:...
- 21
2
votes
3
answers
508
views
Accessing .bam/.cram files from AWS S3?
What's the best/easiest way for accessing .bam/.cram files from S3?
I have .bam/.bai .cram/.crai and .bed + .gff files that I want to make available easily to others so they can browse in IGV and ...
- 2,274
1
vote
2
answers
4k
views
merging two/or more bed file into one bed file
I am trying to merge two bed files (more in future) to one. my bed files are something like :
.
I need to merge them in a way to have the shared chromosome location.
Is there a way to do that ?
- 95
2
votes
1
answer
469
views
Cuffmerge: EOF marker is absent. The input is probably truncated
I need to merge my all transcripts.gtf from cufflinks output follow this command line :
cuffmerge -o merged_gtf_output -p 15 -s ref.fasta -g anot.gtf assembly.txt
...
- 151
2
votes
2
answers
1k
views
Marking optical or PCR duplicates with picard vs. samtools flagstat
I am trying to identify the best library prep method for noninvasive prenatal test samples, to be sequenced with Illumina Novaseq.
One metric that I am evaluating is the number/ percentage of PCR or ...
- 273
1
vote
0
answers
44
views
Indel quality from mpileup
Could someone please confirm that I understand this correctly? I have the res and qual columns from mpileup and I would like to match them to get the qual per base. It seems that the indel initiation (...
- 559
2
votes
1
answer
83
views
Why does samtools mpileup sometimes include ref bases (other than ',' or '.')?
This is my first post here. I can think of no way of giving an easily reproducible example, as per the stack ethos. SO apologies in advance and any feedback on question format appreciated.
I have ...
- 559
2
votes
2
answers
2k
views
How to filter a SAM file by a bed file?
I have a large BAM file and I want to filter it by just what's defined in a bed file, and write it to a new file.
Is it possible to do this with samtools, and if so, how? Is there a better tool I ...
- 161
1
vote
1
answer
128
views
samtools / bamUtil | Meaning of <wildcard> as Reference Name
I have been working with several tools on bam files recently and I am not sure how I should interpret some of the outputs.
samtools idxstats mapped.sorted.bam
<...
- 352
2
votes
1
answer
217
views
How to obtain file offsets corresponding to range query in block-gzipped file
How can we retrieve the file offsets of the zipped-chunks in a block-gzipped file that contain records for a region/window of interest (e.g. Chr01:100-Chr02:145)?
It seems like the CLI of ...
- 319
3
votes
1
answer
566
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,111
0
votes
1
answer
288
views
Single-cell sequencing dataset has too many barcodes
I am analyzing a single-cell sequencing dataset from the website 10xgenomics, with 2000 cells. It is a BAM file and I am trying to obtain the individual cells per sample. I used the command
...
0
votes
1
answer
83
views
Obtaining HGDP project data in fasta format
I need to obtain sample data from modern humans in fasta format. I just need some megabytes of data from every individual. I actually use a script that obtains the cram file from here (ftp.1000genomes....
- 349
3
votes
1
answer
408
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
1
answer
292
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.
4
votes
2
answers
490
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 ...
- 744
2
votes
2
answers
884
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 ...
- 175
1
vote
4
answers
76
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
...
- 105
2
votes
1
answer
261
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. ...
- 21
5
votes
2
answers
393
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:
...
- 585
4
votes
1
answer
750
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 ...
- 41
6
votes
1
answer
2k
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 ...
- 4,815
3
votes
2
answers
4k
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 ...
- 273
2
votes
3
answers
1k
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 ...
- 109
3
votes
2
answers
304
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 ...
- 836
0
votes
1
answer
2k
views
samtools view: writing to standard output failed: Broken pipe
I have the following alignment script which uses BWA:
...
- 453
5
votes
1
answer
356
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 ...
- 5,010
3
votes
1
answer
2k
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 ...
- 836
5
votes
3
answers
4k
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 ...
- 357
1
vote
2
answers
91
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 ...
- 836
5
votes
0
answers
203
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 ...
- 153
5
votes
2
answers
184
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,010
3
votes
2
answers
2k
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 ...
- 560
4
votes
2
answers
5k
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:
...
- 631
1
vote
0
answers
307
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 <...
- 621
3
votes
1
answer
551
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
5
answers
6k
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:
...
- 1,417
4
votes
1
answer
179
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.
...
- 2,719
9
votes
4
answers
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 ...
- 631
4
votes
1
answer
1k
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 ...
Haley
11
votes
1
answer
277
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 ...
- 4,815
5
votes
2
answers
6k
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?
10
votes
1
answer
198
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,010
5
votes
3
answers
419
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.
...
- 1,373
5
votes
1
answer
284
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 ...
- 53