Questions tagged [read-mapping]

When dealing with associating sequencing reads to a matching position in a set of reference sequences (typically a genome).

Filter by
Sorted by
Tagged with
15 votes
2 answers
13k views

Meaning of BWA-MEM MAPQ scores

Does anyone know what the MAPQ values produced by BWA-MEM mean? I'm looking for something similar to what Keith Bradnam ...
ijoseph's user avatar
  • 253
15 votes
1 answer
4k views

Why is bwa-mem the standard algorithm when using bwa?

The industry standard for aligning short reads seems to be bwa-mem. However, in my tests I have seen that using bwa backtrack (bwa-aln + bwa-sampe + bwa-samse) performs better. It is slightly slower, ...
terdon's user avatar
  • 9,662
14 votes
2 answers
4k views

How does the BWA-MEM algorithm assign its mapping qualities?

Is there any resource (paper, blogpost, Github gist, etc.) describing the BWA-MEM algorithm for assigning mapping qualities? I vaguely remember that I have somewhere seen a formula for SE reads, which ...
Karel Břinda's user avatar
13 votes
4 answers
620 views

Library for computing BWT-based alignments

I am writing a software tool to which I would like to add the ability to compute alignments using the efficient Burrows-Wheeler Transform (BWT) approach made popular by tools such as BWA and Bowtie. ...
Daniel Standage's user avatar
12 votes
2 answers
3k views

Difference between BWA-backtrack and BWA-MEM

Many of my colleagues recommend I use BWA-MEM instead of regular old BWA. The problem is I don't understand why and reading the BWA man page doesn't seem to help the matter. What is the difference ...
David Ross's user avatar
9 votes
3 answers
849 views

Visualisation of long read RNA-Seq splicing

I have a dataset of Oxford Nanopore cDNA reads. Many of my reads are full-length or close to full-length transcripts, and I and am interested in examining alternative splicing. For this, I would like ...
Scott Gigante's user avatar
8 votes
2 answers
3k views

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

This is a problem I have to solve frequently, and I'd be interested in knowing what other methods people use to solve the same problem. About twice a year, I get asked to determine variants from ...
gringer's user avatar
  • 13.8k
8 votes
2 answers
1k views

What is mate rescue in bwa mem?

BWA mem has the -S and -P tags for skip mate rescue and skip pairing; mate rescue performed unless -S also in use. What do ...
juniper-'s user avatar
  • 890
8 votes
1 answer
4k views

What are "split reads" and "intron clusters?"

I'm working through the example data set for LeafCutter and the documentation mentions "split reads": This will cluster together the introns fond [sic] in the junc files listed in test_juncfiles....
CelineDion's user avatar
8 votes
1 answer
142 views

Highly heterozygous reads mapping

I have short (67bp) Hi-C reads from a highly heterozygous organism (~15% between-haplotype SNP divergence) and I have both reference haplotypes. I wanted to try and compare different haplotyping ...
cmdoret's user avatar
  • 595
7 votes
2 answers
3k views

Definition of "seed" in sequence alignment

I would like to know what is meant by "seed" for various sequence aligners. How is it important?
user3138373's user avatar
7 votes
2 answers
2k views

What is the correct way to map Hi-C data with bwa mem?

Library Prep I have a Hi-C library prepped using an enzyme that cuts at GATC, so it leaves GATCGATC as the junction sequences. ...
conchoecia's user avatar
  • 3,141
7 votes
1 answer
1k views

How to trim adapter sequences from GSE65360 in order to map the reads?

I am trying to map single-cell chromatin accessibility data from doi:10.1038/nature14590, obtained using scATAC-seq, to the reference genome. Example paired-end reads are here. The reads contain ...
firefly's user avatar
  • 73
7 votes
1 answer
103 views

How to interpret contig-alignment.psa produced by velvet

I'm using velvet to align given reads of RNA to given CDSs (i.e. coding areas and genes) of an organism, so I can generate gene-expression profiles. But after using ...
hhoomn's user avatar
  • 325
6 votes
4 answers
6k views

Why total RNA-seq usually yields low mapping rate?

Maybe this is a silly question, but I'm really wondering why we usually get low mapping rates if we map total RNA-seq, but not poly(A)-enriched (in particular, for human, mouse, and zebrafish datasets)...
kaka01's user avatar
  • 111
6 votes
3 answers
342 views

Filter Trinity transcriptome based on RNASeq reads

I have recently generated a genome-guided transcriptome with Trinity, and would like to apply an additional filter to exclude transcripts that don't have good support from the RNASeq reads. This is ...
gringer's user avatar
  • 13.8k
6 votes
1 answer
159 views

What is the difference between SAM mapping quality and Blast E-value?

Blast reports E-values, but short-read mappers report mapping qualities. Are they the same thing? Can they be converted to each other? If not, why blast doesn't report mapping quality while short-read ...
medbe's user avatar
  • 807
6 votes
2 answers
206 views

How to estimate whether a long-read is meaningful sequence?

The setup Imagine that I work on an organism without a reference genome, and that the closest reference genome I can get is quite diverged. E.g. ~10% diverged in terms of SNVs when measured with ...
roblanf's user avatar
  • 952
5 votes
1 answer
2k views

What is local realignment and what is the problem it solves?

I am trying to understand local realignment but I could not get a clear idea of what is the problem solved by it. For example, reading Homer and Nelson (2010): Because alignment algorithms map ...
gc5's user avatar
  • 1,783
5 votes
2 answers
191 views

RapMap: reference transcriptome for simulated reads

I am very new to bioinformatics and trying to repeat the benchmark in the RapMap paper with an experimental tool working in a similar but different fashion. In the paper (taken from their github) ...
JMC's user avatar
  • 153
5 votes
2 answers
193 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 ...
Daniel Standage's user avatar
5 votes
1 answer
469 views

Why is SMALT better for microbial genomics than other mappers?

SMALT seems to be one of the most used read mappers for bacterial data, see, e.g., this query. I do not say that it is not a great mapper, but I cannot easily see what are its main strengths compared ...
Karel Břinda's user avatar
5 votes
1 answer
392 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 ...
Daniel Standage's user avatar
4 votes
3 answers
530 views

what percentage of the human genome is MAPQ=0?

When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, what percentage of the non-N fraction of the human genome is MAPQ=0? This is, what part corresponds to ...
719016's user avatar
  • 2,324
4 votes
3 answers
106 views

mapping nucleotide sequences

I am trying to do mapping with multiple sequences against a reference genome. I tried the following code: The outputfile is .fastq.sam. I need only .sam. How do I get it? ...
hgf's user avatar
  • 43
4 votes
1 answer
218 views

Kallisto error: index input file could not be opened!

I am utilizing Kallisto in Anaconda/miniconda for RNA sequencing. I have successfully made ...
Tristan's user avatar
  • 43
4 votes
1 answer
5k views

Hisat2 : which option should mention for strand specific library read

I was trying hisat2 I get confused by the strand options... The question has been already asked on github here but didn't get any satisfactory answer. I'm surprise how much this strand information ...
juke34's user avatar
  • 301
4 votes
1 answer
245 views

Problem: "Pair-end" reads scRNA seq data (Drop-seq)

In case of Drop-seq, we have paired end data. Read 1: Cell code + UMI (unique molecule identifier) Read 2: The transcript information But I have a problem/doubt with the sample I am working on. ...
jayesh kumar's user avatar
4 votes
1 answer
106 views

Reference genome for allele specific expression

We are trying to sort out a pipeline for doing allele specific expression. Our plan is to call SNPs from RNA-seq data and combine with known SNP annotations. A well known problem in ASE is reference ...
Ian Sudbery's user avatar
  • 3,301
4 votes
1 answer
2k views

Can I run STAR without an annotation file?

I wish to use Rascaf to scaffold a fragmented draft genome. For this, I need to provide a BAM file of aligned RNA-seq reads and the draft genome. So, I indexed the draft genome with STAR like this: <...
Biomagician's user avatar
  • 2,459
4 votes
1 answer
79 views

Best way to detect long insertions in bisulfite sequencing data?

I am interested in identifying indels in whole genome bisulfite sequencing data (76bp paired end). Currently, I do this by setting the -rfg and ...
Ben D.'s user avatar
  • 397
4 votes
2 answers
713 views

Problems to Extract uniquely mapping reads from BWA MEM alignment

I did a mapping of genomic paired-end reads to a reference assembly using bwa mem. I need to extract the reads that mapped only once to my reference. For that, I have tried to follow the method ...
Marvin's user avatar
  • 41
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 ...
Lot_to_learn's user avatar
3 votes
3 answers
85 views

Which reference to use for read mapping for popular model organisms

What is the "best" assembly for the popular model organisms: human (GRCh37 and GRCh38 are obvious, I'd pick whatever bwakit uses) mouse (GRCm37/GRCm38, OK) but what about non-human/mouse ones? ...
Manuel's user avatar
  • 588
3 votes
1 answer
363 views

How to check what software has been used to generate bam file?

I have several bam files from unknown origin. I'd like to know what software was used to obtain those bam files, e.g. bwa-mem. Is it possible to somehow check it?
Adamm's user avatar
  • 206
3 votes
2 answers
211 views

Contamination on genome assembly

I had a question for the community. I have a genome of a new species that has been sequenced via 150pb Illumina paired-end. To verify the quality of the assembly I used the ...
chippycentra's user avatar
3 votes
2 answers
492 views

How to map reads shorter than 32bp with minimap2?

I mapped fasta sequences to a reference. The fasta sequences range in length from 17bp up to several hundred bp. I used minimap2 with the following command: ...
conchoecia's user avatar
  • 3,141
3 votes
1 answer
2k views

Counts obtained by featureCounts seem much less than observed coverage

I have surprisingly low counts when running featureCounts on some (single-end) RNA-seq data mapped on C. elegans genome using ...
bli's user avatar
  • 3,100
3 votes
1 answer
3k views

Strange per sequence GC content results

I would be grateful if someone could take a quick look at these FASTQC results. This is rna-seq paired-end data. From the FASTQC manual, an unusual distribution seems to be suggestive of contamination ...
beginner's user avatar
  • 631
3 votes
2 answers
73 views

Reference genome (bwa -mem)

I want to do mapping with bwa mem. I downloaded a reference genome from ncbi database which is (Danaus plexippus). I got one ncbi_database.zip file which contain 2 ...
Moon's user avatar
  • 103
3 votes
1 answer
53 views

How do I generate a variant list using Illumina reads from a Salmonella genome?

I am planning on performing phylogenetic analysis of Salmonella specimens using WGS data from PulseNet and GenomeTrakr, for the purpose of public health surveillance & to provide context for ...
Erik D's user avatar
  • 31
3 votes
1 answer
50 views

RNASeq read coverage in protein space?

I used samtools depth to find the per base-pair read coverage over a number of isoform contigs from my Trinity assembly. I have also conducted a multiple sequence ...
CephBirk's user avatar
  • 151
3 votes
0 answers
76 views

Patient-sample mapping in GSE72056 dataset

I want to use the single-cell data from the following expression profiling which concerns the RNA-seq of metastatic melanoma: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72056 The data was ...
hadi's user avatar
  • 31
2 votes
4 answers
522 views

Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?

When we get raw sequenced reads in FASTQ format, we usually need to do some pre-processing (QC, demultiplexing, etc.) prior to sequence alignment. The only way to register the results of this ...
okartal's user avatar
  • 137
2 votes
2 answers
2k views

Error with BWA Mem input having multiple fastq files using cat and process substitution

Running BWA Mem on a number of paired end fastq files using process substitution on the inputs results in this error: ...
111's user avatar
  • 195
2 votes
1 answer
636 views

Understanding PAF format

Does anyone have any experience with PAF files (minimap2 output)? I don’t fully understand what my output means. I’ve looked at the minimap2 explanation (https://github.com/lh3/miniasm/blob/master/PAF....
rimo's user avatar
  • 923
2 votes
1 answer
525 views

What do the symbols mean in minimap2's gap cost equation?

Heng Li's paper on minimap2 is understandable to me right up until the moment that it gives this definition of gap cost: 2.2.1 Alignment with 2-piece affine gap cost Minimap2 performs DP-based ...
Mark Amery's user avatar
2 votes
2 answers
218 views

bedtools: get name column of alignment file A

I'm using the coverage function of bedtools to check whether a certain set of regions (file B) has any overlap with known ...
Flagon13's user avatar
  • 105
2 votes
1 answer
55 views

What tool I can use to map short-reads sequences to reference genome and get specific mapped size

I have about 300 90-bp sequences which I would like to map to a reference genome to make it longer to a 300bp sequence, wherein the 90bp is at the middle. Anyone knows what bioinformatics tool I can ...
Umi Rayearth's user avatar
2 votes
1 answer
634 views

Oxford Nanopore mapping Quality and sequencing error

I'm analyzing some Oxford Nanopore sequencing reads, with several articles stating the sequencing error is between 5% and 15%. Since its 1D and not 1D^2 sequencing i would assume the sequencing error ...
Rasmus Henriksen's user avatar