Devon Ryan
Moderator
  • Member for 4 years, 8 months
  • Last seen this week
What is the fastest way to get the reverse complement of a DNA sequence in python?
Accepted answer
23 votes

I don't know if it's the fastest, but the following provides an approximately 10x speed up over your functions: import string tab = string.maketrans("ACTG", "TGAC") def reverse_complement_table(seq):...

View answer
Convert a BAM file from one reference to another?
Accepted answer
23 votes

You're the second person I have ever seen using NCBI "chromosome names" (they're more like supercontig IDs). Normally I would point you to a resource providing mappings between chromosome names, but ...

View answer
What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
18 votes

5 hours and no benchmarks posted? I am sorely disappointed. I'll restrict the comparison to just be fasta files, since fastq will end up being the same. So far, the contenders are: R with the ...

View answer
Why would someone use a CRAM instead of a BAM?
Accepted answer
17 votes

Whenever you want to save space (this can be a substantial savings). Until quite recently (samtools/htslib 1.7), only CRAM supported long CIGAR strings. If you need to guarantee that any random ...

View answer
How exactly is "effective length" used in FPKM calculated?
15 votes

The effective length is $\tilde{l}_i = l_i - \mu + 1$ (note the R code at the bottom of Harold's blog post), which in the case of $\mu < l_i$ should be 1. Ideally, you'd use the mean fragment ...

View answer
Are soft-clipped bases used for variant calling in samtools + bcftools?
Accepted answer
15 votes

No, samtools (and therefore bcftools) does not use soft-clipped bases. You can quickly confirm this by using either samtools depth or samtools mpileup to look at a region with a soft-clipped alignment....

View answer
Understanding DESeq2 design, contrast and results
14 votes

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest: dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment) The above would give ...

View answer
Merge hundreds of small BAM files into a single BAM file
Accepted answer
14 votes

samtools merge merged.bam *.bam is efficient enough since the input files are sorted. You can get a bit faster with sambamba and/or biobambam, but they're not typically already installed and IO ...

View answer
Can open-source software be peer-reviewed and published?
Accepted answer
13 votes

There are whole journals based primarily around publishing open source tools. The primary example of that is "Bioinformatics", where a lot of the open-source tools are published. We've also had luck ...

View answer
Using the t-SNE algorithm on microarray data + an error bonus
Accepted answer
13 votes

Converting your data.frame to a matrix (and then removing the data.frame) will often free up enough memory that you won't run into this. Note that a matrix is more memory efficient than a data.frame ...

View answer
What Ensembl genome version should I use for alignments? (e.g. toplevel.fa vs. primary_assembly.fa)
13 votes

There's rarely a good reason to use a hard-masked genome (sometimes for blast, but that's it). For that reason, we use soft-masked genomes, which only have the benefit of showing roughly where repeats ...

View answer
Definition of "seed" in sequence alignment
Accepted answer
12 votes

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the ...

View answer
Convert bam file to highly compressible bam
12 votes

You can just set the fields you don't need to *: samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | ...

View answer
How can I extract normalized read count values from DESeq2 results?
Accepted answer
12 votes

The normalized counts themselves can be accessed with counts(dds, normalized=T). Now as to what the baseMean actually means, that will depend upon whether an "expanded model matrix" is in use or not. ...

View answer
Difference between CPM and TPM and which one for downstream analysis?
Accepted answer
10 votes

You can find the various equations in this oft-cited blog post from Harold Pimentel. CPM is basically depth-normalized counts, whereas TPM is length-normalized (and then normalized by the length-...

View answer
Determine if a gene is mitochondrial or not
Accepted answer
10 votes

You can't determine if a gene is mitochondrial purely from its name. The simplest route would be to use the wormbase biomart and either download the gene names and their associated chromosomes or use ...

View answer
How to detect a mutation and predict its consequence?
10 votes

I'll follow up to the great answer from Kamil S Jaron: Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar ...

View answer
What methods are available to find a cutoff value for non-expressed genes in RNA-seq?
10 votes

A common method is to use zFPKMs, which you can find implemented in R here. Having said that, there's an inherent problem in declaring a difference between two things on either side of a given ...

View answer
Edit FASTA header using sed
9 votes

awk 'BEGIN{FS="::"}{if($1~">"){printf(">%s_%s\n",$2,substr($1,2))}else{print $0}}' input_file.fa > output.fa An explanation of the non-obvious bits: BEGIN{FS="::"} split columns using :: as ...

View answer
What is a simple command line tool for doing Needleman-Wunsch pair-wise alignment on the command line
9 votes

Install biopython with conda and use the following script: #!/usr/bin/env python import sys from Bio import Align aligner = Align.PairwiseAligner() aligner.mode = "local" alignments = aligner.align(...

View answer
Running SnakeMake on cluster
Accepted answer
9 votes

You can omit --nodes, you need the following: #!/bin/bash #SBATCH --ntasks-per-node=1 #SBATCH -c threads_from_snakemake #SBATCH -p some_partition Commands go here For slurm, you may want to modify ...

View answer
How to validate that BAMs have been downloaded correctly?
Accepted answer
9 votes

samtools quickcheck is all you need. From the manual: Quickly check that input files appear to be intact. Checks that beginning of the file contains a valid header (all formats) containing at least ...

View answer
R package development: How does one automatically install Bioconductor packages upon package installation?
Accepted answer
9 votes

There's a trick to this where one needs to add biocViews: to the package Description. That's the only solution I've ever seen to allowing automatic installation of bioconductor dependencies. If you ...

View answer
How to count reads in bam per bed interval with bedtools
Accepted answer
9 votes

The order of -a and -b switched at some point. You want: bedtools coverage -a All_peaks.bed -b file.bam > file.cov.txt For reference, this is the end of the help output in version 2.25: Default ...

View answer
Extract nanopore read ID & start times from fastq file
Accepted answer
9 votes

awk '{if(NR%4==1) print $1, $5}' file.fastq | sed -e "s/ start_time=/, /" -e "s/^@//" The awk command gets the first of every 4 lines, printing the first and fifth "word". sed is then used to strip ...

View answer
When to account for the blacklisted genomic regions in ChIP-seq data analyses?
9 votes

Aside: Cross-correlation is largely meaningless, regardless of what some of the ENCODE folks might argue. When we process our DEEP samples we don't even look at that value. Regardless, if you're ...

View answer
Quantifying reads mapping to multiple loci
Accepted answer
9 votes

You almost had the correct python code already, you just need to filter out secondary alignments: def get_reads_hist(bam): bam = pysam.AlignmentFile(bam, 'rb') counts = Counter() for ...

View answer
Bootstrap analysis - why it is called bootstrap?
Accepted answer
9 votes

The method itself has nothing really to do with boots or straps, just as the jack-knife method also has nothing to do with knives. In the case of bootstrapping, the goal is to determine the accuracy ...

View answer
Why does the SARS-Cov2 genome has letter t
8 votes

We sequence and therefore typically report assemblies as DNA sequences, even if they're actually RNA.

View answer
How can I count the number of reads that support a variant in a bam file?
Accepted answer
8 votes

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, ...

View answer
1
2 3 4 5
15