Hot answers tagged

12 votes

Difference between BWA-backtrack and BWA-MEM

TL;DR: BWA-backtrack is based on backtracking. This approach is appropriate only when the dissimilarity between the reads and the reference is low, or when you want to find all best hits or enumerate ...
11 votes

Random access on a FASTQ file

Arbitrary record access in constant time To get a random record in constant time, it is sufficient to get an arbitrary record in constant time. I have two solutions here: One with ...
  • 2,156
9 votes

Why does cutadapt remove low quality bases from the ends of reads only?

If you check the read QC statistics of an Illumina run in e.g. fastQC, you will see that at the end of the read the quality decreases. This is because of exhaustion of chemicals at the end of the run. ...
  • 3,551
9 votes

Difference between BWA-backtrack and BWA-MEM

To quote the Introduction to BWA on sourceforge: BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the ...
9 votes
Accepted

Extracting all reads from bam file which match read IDs in another file

It is still slow but grep has a -f option to take in a file samtools view inbam.bam | grep -f read_names.txt > read_locs.txt
  • 2,574
8 votes

Filtering step for read counts data

If you have your counts in a data.frame called counts, something like this might work: ...
  • 3,551
7 votes
Accepted

Random access on a FASTQ file

As wkretzsch suggested this was worthy of an actual answer, I feel the obvious solution is missing here; index the FASTQ. Index it As much as I typically hesitate ...
7 votes
Accepted

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

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ...
  • 19.4k
7 votes
Accepted

How to remove all BAM read groups from all reads (not just the header)?

You might have to manually strip those auxiliary tags off: ...
  • 19.4k
7 votes
Accepted

Why does cutadapt remove low quality bases from the ends of reads only?

I am commenting on this part: The algorithm only removes low quality bases from the end until it reaches a good quality base. If there is a bad quality base beyond that, it is not trimmed. ...
  • 6,121
6 votes

Counting repeated kmers sequences that match at least x % of reads sequence

The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution. EDIT (response to new ...
  • 6,121
6 votes
Accepted

How to check whether all BAM read contain defined read groups?

According to the man page, running samtools stats --split RG <file1.bam> should produce summary statistics separated by read group. If it doesn't produce a ...
  • 12.5k
5 votes

Random access on a FASTQ file

You could shuffle the FASTQ once and then read sequences off the top of the file as you need them: ...
  • 2,156
5 votes
Accepted

What is the right way of calculating a Phred score by hand?

Not sure if my explanation is any good but let's try... First, let's convince ourselves that converting phred to probability is the right thing to do (i.e. your Way 2). Imagine a read of length 10000 ...
  • 689
4 votes
Accepted

Why are there missing calls in a VCF file from exome sequencing?

Missing variant calls due to lack of coverage shouldn't happen in the targeted capture region and I'd think most of these would come from off-target regions where some samples had reads mapped. I'd ...
  • 393
4 votes

How to check whether all BAM read contain defined read groups?

using samjdk and invoking the function getReadGroup() getReadGroup() returns The SAMReadGroupRecord from the SAMFileHeader for this SAMRecord, or null if 1) this record has no RG tag, or 2) the ...
  • 1,473
4 votes
Accepted

Select top 100 genes ranked by variance in read counts

I assume by "reverse sort of variance" you mean "highest variance". Assuming you made a matrix out of that (set the row names to the first column and then remove it) and called it ...
  • 19.4k
4 votes

What is and how to detect a dephased read

Phasing in Illumina data means that the read you get is shifted one base before or after what the sequence really should be. They must be expecting some constant 6 bases at the end of each read, and ...
  • 1,792
4 votes

Extract reads from bam files by their @RG

Typically I use samtools for operations like this. Specifically I use samtools view with either -r or ...
  • 2,574
3 votes

While opening a Bam file with the SeqAn library I get 'seqan::FileOpenError'

Remove the getAbsolutePath call — your path is already absolute. getAbsolutePath garbles it, as you can see in the error message....
3 votes
Accepted

Scaling by linear regression against the number of reads

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress ...
  • 200
3 votes

Filtering step for read counts data

While this answers explains how to do it I want to address when and why and which thresholds to do it. Filtering the genes with low counts is usually done because the counts are not reliable it would ...
  • 4,652
3 votes

Why are there missing calls in a VCF file from exome sequencing?

I've just been generating data like this, so can tell you about why/how missing calls are created in my dataset. There are two main reasons: 1. Sequencing failure When reads don't map across the ...
  • 12.5k
3 votes

Random access on a FASTQ file

One possibility is to: reformat the data such that each record is a single line containing the read description, bases, and quality scores pad out each record to a maximum length in each field such ...
  • 31
3 votes

Reads mapped to exonic, intronic and intergenic regions

As @DevonRyan mentioned, it's very likely that those samples were degraded, which is good justification for excluding them from subsequent analysis.
3 votes

If a gene is expressed at a level of 1/1200 compared to the average gene, how is probability 50:50 that we have a read mapped to it?

I actually disagree with that.. I guess I write this down as a discussion. The expected value for an average gene is 1200. For this gene at 1/1200 expression, you expect 1 read. However, because of ...
  • 1,658
3 votes

If fastp output is not a good measure of FASTQ correctness, what is?

On a regular basis, using FastQC is quite enough as to assess the quality of FastQ-formatted data. It gives you plenty of details as can be seen in the screenshots you shared. As others have commented ...
  • 705
2 votes

Random access on a FASTQ file

I wrote a tool called sample that you can use to do random sampling without reading the entire file into memory. It can be used where GNU ...
2 votes

Random access on a FASTQ file

One of the most thorough treatments of this question (or a similar question: grabbing a random subset of reads) was given by Jared Simpson in a blog post a few years ago. http://simpsonlab.github.io/...
2 votes

Genome assembly from error-prone reads

For short reads, the typical and the most widely used solution is to correct away sequencing errors before assembly. You can correct errors with k-mer spectrum, a trie or multi-alignment. There are ...
  • 6,121

Only top scored, non community-wiki answers of a minimum length are eligible