Skip to main content
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 ...
Karel Břinda's user avatar
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 ...
winni2k's user avatar
  • 2,286
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. ...
benn's user avatar
  • 3,591
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 ...
Michael Hall's user avatar
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
Bioathlete's user avatar
  • 2,574
8 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 ...
Sam Nicholls's user avatar
8 votes

Filtering step for read counts data

If you have your counts in a data.frame called counts, something like this might work: ...
benn's user avatar
  • 3,591
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 ...
Devon Ryan's user avatar
  • 19.7k
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: ...
Devon Ryan's user avatar
  • 19.7k
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. ...
user172818's user avatar
  • 6,545
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 ...
user172818's user avatar
  • 6,545
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 ...
gringer's user avatar
  • 14.4k
6 votes
Accepted

How many false positive duplicates are marked using just the position of first unclipped base?

There's a tool in Picard (full disclaimer: I wrote it...) called CollectIndependentReplicateMetrics which performs the test you are looking for, I think. It measures the percent of duplication that is ...
Yossi Farjoun's user avatar
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: ...
winni2k's user avatar
  • 2,286
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 ...
dariober's user avatar
  • 699
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 ...
Vivek's user avatar
  • 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 ...
Pierre's user avatar
  • 1,536
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 ...
Devon Ryan's user avatar
  • 19.7k
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 ...
swbarnes2's user avatar
  • 1,941
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 ...
Bioathlete's user avatar
  • 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....
Konrad Rudolph's user avatar
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 ...
DCZ's user avatar
  • 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 ...
llrs's user avatar
  • 4,713
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 ...
gringer's user avatar
  • 14.4k
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 ...
kvg's user avatar
  • 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.
Daniel Standage's user avatar
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 ...
StupidWolf's user avatar
  • 1,688
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 ...
JRodrigoF's user avatar
  • 827
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 ...
Alex Reynolds's user avatar
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/...
Daniel Standage's user avatar

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