Skip to main content

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
0 answers
16 views

pysam ignoring some reads?

I want to modify TEtools to be able to use pysam so it can accept bam. My version returns less read than the TE tools version. Can someone explain to me why my code returns fewer reads? From my ...
RomainL.'s user avatar
  • 165
4 votes
1 answer
122 views

tabix errors when accessing 1000 Genomes data: "[E::bgzf_read] bgzf_read_block error -1 after 50219 of 52392 bytes" and Could not load .tbi/.csi index

I am trying to access 1000 Genomes (1KG) data using tabix as per the 1KG tabix documentation "How do I get a genomic region sub-section of your files?" ...
N Brouwer's user avatar
  • 141
4 votes
2 answers
155 views

Running bwa mem and samtools in one script on a linux cluster

I am trying to write a script that will run this line from https://github.com/kmceres/Mbovis_pangenome/blob/main/Parameter_file.md: ...
Lucy Kelly's user avatar
1 vote
2 answers
34 views

How to Call SNPs from Assembled Fasta Files?

I have experience in obtaining a variant calling file from NGS sequencing files using the Samtool/bedtools workflow. However, when we receive pre-assembled sequences, it becomes challenging to apply ...
Gerald Vasquez Aleman's user avatar
1 vote
1 answer
83 views

How to subset BAM file based on read length range (120-180) bp?

Hi I'm wondering how I can subset a BAM file based on read length. Precisely I want just read lengths of 120 to 180 bp reads in the new BAM file. I'm trying several way, but none of them giving the ...
Deb's user avatar
  • 319
1 vote
1 answer
35 views

Problem while mapping reads to mtDNA (SortSam)

I am trying to map MiSeq reads to a reference genome and extract mutations using MToolBox, which implements gsnap, GATK, Picard, and other tools. When running the tool with example data, there were no ...
uri's user avatar
  • 85
4 votes
1 answer
62 views

Different line length in sequence 'chrY'

I just downloaded a reference genome (using wget), and attempted to use it with samtools view. However, I received the following error: ...
Wouter De Coster's user avatar
2 votes
1 answer
83 views

Pysam - fetch reads within region

I'm using pysam (and also samtools) to find all reads in BAM files that are within a specific region py_bamfile.fetch('4',42266768,42268410) and the following read ...
Tien's user avatar
  • 53
2 votes
1 answer
89 views

Aligning FASTQs to FASTA reference

I'd like to align some FASTQs to an average mtDNA FASTA file that I have downloaded so I can have the human mtDNA isolated from those FASTQs. For that, I used bowtie2. Can I expect that after running ...
dyxcvi's user avatar
  • 151
1 vote
1 answer
36 views

Discrepancy in Depth of Coverage Estimation between Pileup and samtools depth

I'm encountering a discrepancy in the depth of coverage assessment for specific coordinates when comparing results obtained from pileup files and samtools depth. I've processed pileup files and ...
Irina's user avatar
  • 11
1 vote
1 answer
114 views

Get a certain gene sequence from bam/vcf and reference

I need to get a fasta sequence of a certain gene for a certain worm strain that is different from reference. I have a reference genome, BAM for the strain of interest, and coordinates of the gene. I ...
user98747's user avatar
2 votes
1 answer
39 views

I am trying to create a subset of 10k variants from 25-30 unmapped contigs of a g.vcf file including the header

My objective is to take a g.vcf.gz file and from 25-30 unmapped contigs with titles like "NW_020192317.1", I want to make a subset of ~10k variants from each of the unmapped contigs and make ...
Lauren Sabo's user avatar
3 votes
1 answer
243 views

low mapping percentage in bwa mem

After aligning paired-end reads to a reference genome, I am getting low percentage: ...
Moon's user avatar
  • 103
3 votes
1 answer
130 views

How to reduce BAM

Consider a sorted, indexed, only mapped reads containing BAM file. Is there a way to get a sub BAM based on read line numbers? This can be done iteratively using a counter but its too slow. Is there a ...
papabiceps's user avatar
1 vote
1 answer
67 views

Create bam from separate sam and header

If I have multiple SAM files without header, and a header saved in a separate file, how do I convert that to a BAM file? One solution is to create intermediate files: ...
Alexlok's user avatar
  • 394
1 vote
0 answers
203 views

Mapping statistics from the bam files

I would like to find the mapping statistics from the sorted bam files. Samtools flagstat gives the output only for a single file. What is the easiest way to find the mapping statistics for all ...
Kam's user avatar
  • 11
3 votes
1 answer
59 views

regex: samtools command to "refine" PacBio IsoSeq data?

I am working with preprocessed data (IsoSeq PacBio) and I cannot understand one of the steps: samtools view -b -e [rq] >= 0.9 flnc.bam What exactly was ...
Aleksandra Greshnova's user avatar
2 votes
2 answers
302 views

Problem with samtools command with a final hyphen

...
ii4 unsafe's user avatar
3 votes
1 answer
831 views

Nextflow Error: failed to read header from "-"

I am trying to run my nextflow pipeline, and have gotten this error: samtools sort: failed to read header from "-" I'm not sure why this error is ...
anne's user avatar
  • 65
1 vote
0 answers
48 views

How do I filter reads from a bamfile based on 5' and 3' ends?

I have data in paired-end as well as single-end format for the same sample.I would like to split those two files into two per each, where one contains sequences that have 5' end and the other one has ...
Pikabu's user avatar
  • 11
3 votes
1 answer
244 views

CRISPR genome wide screen analysis using MAGeCK-- Loss of reads/sgRNA using MAGeCK count function?

Sorry if this question seems basic. I am trying to run the MAGeCK pipeline to analyze CRISPR knockout screen data produced by someone in my lab around 5 years ago. I was given the data as bam files ...
Cassie Bishop's user avatar
4 votes
0 answers
138 views

Write a bash script to run gatk, fix errors with input, and rerun until completion

I have a bam file that I want to run through GATK's SplitNCigarReads tool. Because of the way the bam file was generated, the program will often fail, with an error message stating: ...
kylep's user avatar
  • 41
2 votes
0 answers
150 views

Difference between pileup and mpileup in VarScan and samtools

I'm going to call variants with VarScan from a pileup files created with samtools. I realized that there are in general two major possibilities to call variants, pileup as well as mpileup. The VarScan ...
Anti's user avatar
  • 153
1 vote
0 answers
149 views

Trinity de novo transcriptome assembly: samtools view failed to add PG line to the header

I have assembled the transcriptome a plant species using Trinity. Here is the command: ...
Faroll's user avatar
  • 11
2 votes
1 answer
1k views

how to retain reads with low mapping quality (MAPQ) scores when using samtools view -q

I have been using the -q option of samtools view to filter out reads whose mapping quality (MAPQ) scores are below a given ...
Patrick Bastedo's user avatar
2 votes
1 answer
256 views

Get number of reads with a single, (almost) exact match to the full length of a reference sequence

I can't find an answer to this in previous questions, so hoping someone can help me now. We have nanopore sequenced a PCR product, and I've filtered our reads to +/- 10 bp of the expected product ...
Nat R McB's user avatar
1 vote
1 answer
400 views

Different line length in fasta file

I am currently using VEP for variant annotation. I am facing an error as below: [E::fai_build_core] Different line length in sequence 'Pn9' I understand there is ...
Sowmya Pulapet's user avatar
3 votes
1 answer
2k views

How to remove the unpaired reads in sam/bam files?

I have sam and bam files for the chimeric reads, which come from two different parts of the genome (For example, the first half of the read from part of Chromosome 1 and the second half of the read ...
Wang Ming's user avatar
  • 101
0 votes
1 answer
407 views

Can I use samtools addreplacerg to replace multiple RG entries at the same time?

I have a bam file that contains two @RG lines: ...
Greg Dougherty's user avatar
2 votes
2 answers
521 views

Checking that two multiline FASTAs are identical (allowing for different order and lower/upper case)

Is there a good CLI tool that can diff two multiline fasta files, ignoring order, and comparing sequences with identical names with each other pairwise. Bonus points if lower/upper case is tolerated (...
Cornelius Roemer's user avatar
0 votes
2 answers
57 views

BAM files with no RNAME and POS, how to map contents to SNPs?

I have a set of 4 .bam files containing the exome of an individual, around 400 MB each. I used samtools to generate a 2.4 GB .sam file out of one of the .bam files, and I found it contains lines with ...
Fernando D'Andrea's user avatar
2 votes
0 answers
402 views

How to extract reads that map exclusively to a single site with 1 or zero mismatches from BAM files

I generated BAM files (sorted by coordinates) by aligning human RNA reads against the human reference genome using BWA MEM and ...
plicht's user avatar
  • 21
3 votes
2 answers
209 views

Detect mutation context in a read of a sam file

After sequencing (Illumina) of some DNA, I generated a sam file through alignment of a fastq file (using Bowtie2). Instead of using variant calling programs, I want to know the specific variances, ...
Annelisa's user avatar
1 vote
2 answers
1k views

Potential side effects of replacing read group tags in BAM file

I have a set of BAM files where the read group tags have some (default?) values, i.e.: @RG ID:RG0 LB:LB0 PU:PU0 SM:SM0 This creates issues in my downstream ...
gc5's user avatar
  • 1,813
4 votes
1 answer
489 views

samtools view command not found error

When I tried to use samtools to split a bam file based on different chromosomes, I used this command: samtools view input.bam -b chr21 | chr21.bam However, I get ...
Scott XU's user avatar
  • 145
4 votes
1 answer
658 views

BWA mem | samtools view: Intermittent parsing error

Update The issue was that bwa was running out of memory and failing, but that error wasn't floating to the top (see @Steve's answer, below). I was getting an error ...
Mark Ebbert's user avatar
  • 1,404
2 votes
1 answer
476 views

Samtools sort: most efficient approach to sort a single sample after aligning split .fastq files

Related to my other question (Samtools sort: most efficient memory and thread settings for many samples on a cluster), we need to optimize samtools sort as we ...
Mark Ebbert's user avatar
  • 1,404
5 votes
2 answers
5k views

Samtools sort: most efficient memory and thread settings for many samples on a cluster

We're preparing to analyze thousands of .bam files beginning with re-alignment, sorting, etc. Complicated question: Has anyone investigated the optimal thread and ...
Mark Ebbert's user avatar
  • 1,404
0 votes
1 answer
4k views

How to Sort and Index a SAM file without converting it to BAM?

Generally, I use samtools to sort and index BAM files. Samtools works also with BGZF ...
alec_djinn's user avatar
3 votes
0 answers
635 views

How to analyze an IGV alignment

I'm working on a project where I am analyzing the performance of an alignment workflow. My goal is to find regions in the resulting BAM file where there are outstanding discrepancies or anything that ...
pythonbeginner44's user avatar
1 vote
1 answer
58 views

Is there a sam flag for all&none of the reads?

a. Is there a SAM flag that specifies all reads? b. Is there a SAM flag that specifies none of the reads? So that if I run samtools view -f (a) -F (b) the result will be all reads of the file (as if ...
Sam's user avatar
  • 149
0 votes
1 answer
231 views

Convert VCF file to mpileup.txt

I am working on an iterative analysis that uses orthologous pipelines that require mpileup.txt files as input for a visualization step. This requires me to convert VCF files to mpileup.txt. This ...
Andrew Judell's user avatar
1 vote
2 answers
563 views

How to extract all sequences mapped to a transcript from Kallisto output

This question has also been asked on Biostars I ran Kallisto with the --pseudobam option. How do I extract all the short reads that are mapped to a single transcript (e.g. ENST00000367969.8)? As a ...
Dr. Who's user avatar
  • 21
1 vote
1 answer
833 views

How to convert multiple single-end bam files to fastq using samtools

Hi I am trying to convert bam files generated from Ion Torrent Proton sequencing to fastq format so that I can upload them to ...
Justin1609's user avatar
2 votes
1 answer
70 views

Counting the co-occurrence of variants A and B in an aligned sequencing read

I need some help getting started on this project. To simplify we want to be able to quantify the occurrence of 3 variants on each sequencing read in an alignment file as a proxy measurement for ...
Bioreeb's user avatar
  • 21
0 votes
1 answer
601 views

Is it normal to have a smaller bam file after merging bam files?

I have 2 bam files that belong to the same sample and I merged them with samtools merge. And after merging, I realized that the merged version is a bit smaller than ...
Baran Aldemir's user avatar
1 vote
2 answers
436 views

Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?

I am using htslib/sam.h to write a C++ program. As part of this program, I must extract reads occurring on specific scaffolds in specific regions from a bam file. Essentially, I want to perform the ...
annabelperry's user avatar
4 votes
4 answers
2k 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++....
annabelperry's user avatar
1 vote
0 answers
884 views

Conversion of SAM to BAM files

I am very new to micro RNA analysis. I have been using H. sapiens, GRCh38 + major index as given in the Bowtie Website to align with my trimmed FASTQ file . The command I am using to make very ...
Aranyak Goswami's user avatar
2 votes
1 answer
46 views

Python package or CLI tool outputting mutations of sequence with respect to reference

I'm new to bioinformatics, and I'm looking for a tool that takes two FASTA files, one containing a reference, let's call it A.fasta, the other another sequence, let'...
Cornelius Roemer's user avatar