# Question
Title: Classification (supervised learning) of expression data on pathway level
I was wondering if there is any way to apply classification algorithms (e.g random forest) on microarray data but not using the genes as predictors/features but the pathway they belong to.

The thing is that the expressions of the genes that belong to the same pathway should somehow, grouped together and transformed into a new value that will describe that pathway for each specific sample.

Let me give an example.

Assuming that we have `g1, g2, g3, g4, g5, g6, g7, g8, g9, g10` as genes and `p1, p2, p3` as pathways.

And `g1, g2, g3` belong to `p1` while `g3, g4, g5, g6`belong to `p2` and `g2, g6, g7, g8, g9, g10` belong to `p3`

One thought is to calculate the mean expression of the above gene groups per sample and then run the classification algorithm on top of these values. But I would like to hear/read your opinion on such an approach and also your suggestions on using a different method to do this.

The final goal would be the extraction of the most significant classification predictors (pathways).

# Answer
> 6 votes
Have a look at the GSVA package. It allows to convert a matrix with `genes x Samples` to a `pathways x Samples` using several methods ssgsea, gage, gsva...

Afterwards you can use that matrix as input for differential expression of pathways or classification algorithms or whatever.

---

However it depends on the input of the "pathways" you give it. Make sure that the pathways are correct. Also consider that the pathway topology might change with conditions.

---
Tags: microarray, pathway
---


# Question
Title: gffread: GFaSeqGet errors on coordinate overhang
***Disclaimer****: I had this issue posted on Tuxedo Tools users group and shared it on Twitter, but could not get an answer to this from the developers, nor find documentation of this issue online. So, I'll share my solution here below and see if some of you have a better solution to add.*

Using gffread utility to convert transcripts from GTF as output of Cufflinks to FASTA, I get some `GFaSeqGet` errors, as follows (here is the full log file):

```
$ gffread transcripts.gtf -g ../genome.fa -w transcripts.fasta
Error (GFaSeqGet): end coordinate (32519) cannot be larger than sequence length 32512
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (6327) cannot be larger than sequence length 6326
Error (GFaSeqGet): end coordinate (4651) cannot be larger than sequence length 4648
Error (GFaSeqGet): subsequence cannot be larger than 3551
Error getting subseq for CUFF.29386.2 (1..3568)!

```

That's a bit weird since I used the same genome file as the one used to generate the alignment file required by Cufflinks, and the coordinates overhang does not make a lot of sense to me.

The alignment file (SAM) was generated using Hisat2:

```
$ hisat2
    --rna-strandness RF
    --dta-cufflinks
    -x ./hisat_idx
    -q
    -1 ../1A_1.fq.gz,../1B_1.fq.gz,../2A_1.fq.gz,../2B_1.fq.gz,../3A_1.fq.gz,../3B_1.fq.gz
    -2 ../1A_2.fq.gz,../1B_2.fq.gz,../2A_2.fq.gz,../2B_2.fq.gz,../3A_2.fq.gz,../3B_2.fq.gz
    -S out.sam

```

Then sorted:

```
$ cat out.sam | sort -k 3,3 -k 4,4n -T . > out_sorted.sam

```

My question is: are only the few sequences reported in the log affected? Or is it a more general problem? The answer seem to be the second one.

Out of gffread, alongside the log file, I actually got a FASTA file with no additional error messages during its generation (537M out of 404M GTF file) - but there are around 13k transcripts which are present in the GTF but seem to be missing from the FASTA:

```
$ cat transcripts.gtf | awk '$3=="transcript"{print $0}' | wc -l
231139
$ cat transcripts.fasta | grep '^>' | wc -l
218075

```

Do you know what would be a cause of this problem?

# Answer
Basically, the source of the problem is issues with hard/soft clipped bases in short-reads RNAseq reads alignments to reference genome in post-transcriptome assembly steps.

Truncating the coordinates in the gtf file to match the reference is a solution. For this, I will share a perl solution (non-official) that will do the job.

```
perl cufftrim.pl genome.fa.fai transcripts_to_fix.gtf > transcripts_fixed.gtf

```

Where:

* `genome.fa.fai` is the fasta index of the genome file
* `transcripts_to_fix.gtf` is the transcript set out of the assembler with issues
* `transcripts_fixed.gtf` is the new transcript file with the repaired end coordinates

Content of `cufftrim.pl`:

```
#!/usr/bin/perl
use strict;
use Getopt::Std;

my $usage = q/Usage:
 cufftrim.pl genome.fa.fai transcripts.gtf > transcripts.trim.gtf
/;
umask 0002;
getopts('o:') || die($usage."\n");
my $outfile=$Getopt::Std::opt_o;
if ($outfile) {
  open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
  select(OUTF);
  }
# --
my $fai=shift(@ARGV) || die ("$usage\nSamtools fasta index required!\n");
my $gtf=shift(@ARGV) || die ("$usage\nCufflinks' transcripts file required!\n");
my %slen;

open(FAI, $fai) || die("Error opening $fai\n");
open(GTF, $gtf) || die("Error opening $gtf\n");
while (<FAI>) {
 chomp;
 next if m/^#/;
 my @t=split(/\t/);
 die("Parsing error: $t[1] length for sequence $t[0] ?\n")
    unless $t[1]>10;
 $slen{$t[0]}=int($t[1]);
}
close(FAI);

# now process the GTF file
my ($tcount, $ecount); #trim counts
while (<GTF>) {
 my @t=split(/\t/);
 my $l=$slen{$t[0]} || die("Error: contig $t[0] not found in the fasta index!\n");
 if ($t[4]>$l) {
   $t[4]=$l;
   $t[2] eq 'exon' ? ++$ecount : ++$tcount;
 }
 print join("\t", @t);
}
close(GTF);
print STDERR "$tcount transcripts, $ecount exons adjusted.\n";
# --
if ($outfile) {
 select(STDOUT);
 close(OUTF);
 }

```

Filtering the SAM file prior to assembly to exclude clipping might also be an alternative solution.
> 2 votes

---
Tags: rna-seq, sequence-alignment, software-recommendation, base-clipping
---


# Question
Title: Can a data file in VCF format be converted into FASTA?
I'm considering purchasing the 'MyGenome' product by Veritas Genetics to analyze my genome for a project. I'd like the data to be in FASTA format, but Veritas only provides VCF data. Is it possible to convert this VCF data into FASTA format?

# Answer
> 5 votes
You can try gatk function `FastaAlternateReferenceMaker`

```
java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \
   [--snpmask mask.vcf]

```

---
Tags: fasta, vcf, ngs, file-formats, variant-calling
---


# Question
Title: coexnet package for network co-expression analysis, is looking for a filelist.txt that doesn't exist
On my way to create a co-expression matrix for a microarray experiment, using the `coexnet` package I had to use a function called `getAffy()` which according to the packages' manual it :

> searches in your current or designated path file the folder with the GSE accession number and reads the filelist.txt file that contains the name of each GSM sample, in order to recognize them and join them in an unified AffyBatch object.

So in some way, this file is needed to create the AffyBatch object. The thing here is that there is no such a file downloaded from the `getInfo()` function.

Googling out the `filelist.txt` I reached only this result from GEOQuery manual, referring to that file. And as you can see there is a command in it `eList2 <- getGEOSuppFiles("GSE11675")` that is believed that it returns also the filelist.txt. I tried even with the example of the previously mentioned command, but no filelist.txt file ever fetched.

Is there any way to setup my own filelist.txt file? Does anyone know what it contains and in what form?

# Answer
Finally, I found it. `filelist.txt` file is accessible from the FTP server of the NCBI. For example for that specific microarray is stored here

So there might be a problem at software side and they don't get downloaded.
> 1 votes

---
Tags: microarray, geoquery
---


title: Identifying outlier loci in R

text:
# Question
Title: Identifying outlier loci in R
I have done some genome scan analysis and I want to identify the outlier loci. For an expected value of the FRD equal to q = 10, a list of candidate loci can be obtained by using the Benjamin-Hochberg procedure as follows:

```
q=0.1
L = length(adjusted.p.values)
W = which(sort(adjusted.p.values) < q * (1:L) / L)
candidates = order (adjusted.p.values)[w]

```

and this is part of the outliers and my output

```
head(candidates)
 [1]  3  4  9   10  13

```

But I need the output in a different type. I want to assign TRUE if a locus has been chosen as outlier and FALSE if it has not

```
  head (desired_output)

  [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE

```

Any idea how can I modify the Code?

# Answer
As Devon said, removing the `which` will give you `W` as a logical vector that can be used to index the *sorted* p-values. To get a logical vector that can be used to subset the original (unsorted) values you’ll need to go a step further:

```
candidates = is.element(1 : L, order(adjusted.p.values)[W])

```

However, I’m not sure this will work with your data: the name `adjusted.p.values` suggests that these values are *already* FDR-adjusted. So you cannot/should not perform additional correction on them. If they are already adjusted, then you can get candidate loci simply by

```
candidates = adjusted.p.values < q

```
> 3 votes

# Answer
Remove the `which()`, which is converting the boolean vector to a vector if indices to `TRUE`.
> -1 votes

---
Tags: r
---


# Question
Title: Can I export the SNPs between whole genome alignments using the command line?
I am using Mauve to align two whole genomes. I use the following command to get the alignment in xmfa format:

`/Applications/Mauve.app/Contents/MacOS/progressiveMauve ref1.fa ref2.fa --output alignment.xmfa`

I then open the xmfa file in the Mauve GUI and export the SNPs. Is there a way of doing this using a command? I have tried to follow the procedure described here but the link to snapshot releases of Mauve does not work anymore.

I could not find anything in the Mauve command line help:

`/Applications/Mauve.app/Contents/MacOS/progressiveMauve --help`

nor on the Mauve website.

I tried the procedure described by Konrad Rudolph. I wrote:

MAUVE\_DIR=/Applications/Mauve.app/Contents/Java

in the terminal. Is this correct? Don't I need an export and quotation marks to create a bash variable? I tried adding an export and using quotation marks for the file but the error below remained the same.

followed by:

`export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"`

followed by:

java org.gel.mauve.analysis.SnpExporter -f alignment.xmfa -o aln.snps

and it seems to work.

# Answer
> 2 votes
It works for me with the download from download page. The download archive contains `Mauve.jar` so the command line in the answer you link to should work.

However, you might need to include all the jar files (not just `Mauve.jar`) in the classpath for this to command work, though. For instance, the following works for me:

```
MAUVE_DIR=‹the-dir-containing-Mauve.jar›
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
java org.gel.mauve.analysis.SnpExporter -f aln.xmfa -o aln.snps

```

… there’s probably a smarter way to make this work but I don’t use Java enough to know how.

---
Tags: sequence-alignment, snp, variation, mauve
---


# Question
Title: How to use C++ htslib to read VCF contig name and size?
A typical VCF file has:

```
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>

```

I would like to use `htslib` in C++ to read it. My attempt:

```
htsFile *fp = bcf_open("my.vcf", "r");
bcf_hdr_t *hdr = bcf_hdr_read(fp);

```

In https://github.com/samtools/htslib/blob/develop/htslib/vcf.h, I'm not able to find a function that can do that for me.

How to read `chr1` and `248956422` in C++?

# Answer
```
#include "htslib/vcf.h"

int main(int argc, char *argv[])
{
    htsFile *fp;
    bcf_hdr_t *hdr;
    bcf_idpair_t *ctg;
    int i;
    if (argc == 1) {
        fprintf(stderr, "Usage: print-ctg <in.vcf>\n");
        return 1;
    }
    fp = vcf_open(argv[1], "r");
    hdr = vcf_hdr_read(fp);
    ctg = hdr->id[BCF_DT_CTG];
    for (i = 0; i < hdr->n[BCF_DT_CTG]; ++i)
        printf("%s\t%d\n", ctg[i].key, ctg[i].val->info[0]);
    bcf_hdr_destroy(hdr);
    vcf_close(fp);
    return 0;
}

```

On stability: this use has been in htslib forever. In general, functions/structs/variables in the public headers are meant to be stable. However, there is no guarantee that future versions will always keep the same APIs.
> 4 votes

---
Tags: c++, htslib
---


# Question
Title: How to make chromosome color maps for bed ranges
I have genomic .bed file data of 4 different types; type A,B,C,D. These are some examples-

Type A:

```
1   101380000   101710000   A
1   110085000   110320000   A

```

Type B:

```
1   100930000   101335000   B
1   10430000    10560000    B

```

Type C:

```
1   101775000   101825000   C
1   103570000   104070000   C

```

Type D:

```
1   153325000   153475000   D
1   156355000   156415000   D

```

This data is for chr1, but I such data for all the chromosomes. I want to draw a chromosome map (possibly using the Gviz package?) which can assign different colors to the 4 types of genomic ranges and make me an idiogram like plot for the chromosomes I wish to see, so I know where my ranges physically are, on the chromosome.

Can this be done using the **Gviz** package in R? If not, what is the best way and code to do this?

# Answer
You can use karyoploteR to plot your regions on an ideogram quite easily.

**Disclosure:** I'm one of the authors of the tool.

For this example I'll start creating 4 GRanges with random regions, but with your data you should load your bed into R and split it into four GRanges.

```
library(karyoploteR)

typeA <- createRandomRegions(nregions = 100, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 20, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 1000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 100, length.mean = 2e6, length.sd = 5e6)

colors <- c(A="#FFD700", B="#00BFFF", C="#FF7F50", D="#90EE90")

```

And we can plot them. We have to start creating an empty ideogram with `plotKaryotype` (we are working with human here, but you can use any genome you want) and then plot the regions in the GRanges with `kpPlotRegions`. As a bonus, if two regions overlap it will plot them one above the other.

```
kp <- plotKaryotype(genome = "hg19")

kpPlotRegions(kp, data=typeA, r0=0, r1=0.2, col=colors["A"])
kpPlotRegions(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"])
kpPlotRegions(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"])
kpPlotRegions(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"])
legend(x = "right", legend = names(colors), fill = colors, border=NA)

```

And you'll get something similar to this

In this example we have fairly large regions and not a lot of them. As @heathobrien pointed out, plotting the actual regions in the whole genome plot might not be the best visualization option. You can either use `zoom` to create different plots of smaller parts of the genome or change your visualization to show the density of regions along the genome with `kpPlotDensity`.

Plotting the density and tweaking a bit the visualization (and changing the example to have 10 times more regions) with this codei

```
typeA <- createRandomRegions(nregions = 1000, length.mean = 5e6, length.sd = 3e6)
typeB <- createRandomRegions(nregions = 200, length.mean = 15e6, length.sd = 10e6)
typeC <- createRandomRegions(nregions = 10000, length.mean = 5e5, length.sd = 2e5)
typeD <- createRandomRegions(nregions = 1000, length.mean = 2e6, length.sd = 5e6)

kp <- plotKaryotype(genome = "hg19", plot.type = 4,
                    ideogram.plotter = NULL, labels.plotter = NULL, main="Density")
kpAddCytobandsAsLine(kp)
kpAddChromosomeNames(kp, srt=45)

kpPlotDensity(kp, data=typeA, r0=0, r1=0.2, col=colors["A"], window.size     = 10e6)
kpAddLabels(kp, labels = "Type A", r0=0, r1=0.2)
kpPlotDensity(kp, data=typeB, r0=0.25, r1=0.45, col=colors["B"], window.size = 10e6)
kpAddLabels(kp, labels = "Type B", r0=0.25, r1=0.45)
kpPlotDensity(kp, data=typeC, r0=0.5, r1=0.7, col=colors["C"], window.size = 10e6)
kpAddLabels(kp, labels = "Type C", r0=0.5, r1=0.7)
kpPlotDensity(kp, data=typeD, r0=0.75, r1=0.95, col=colors["D"], window.size = 10e6)
kpAddLabels(kp, labels = "Type D", r0=0.75, r1=0.95)

```

it would look something like this

The plots could be tweaked a lot changing various parameters and you could add other data on top of it if needed. You can find more information on the Bioconductor landing page or at the karyoploteR's Tutorial and Examples page.
> 12 votes

---
Tags: r, gviz
---


# Question
Title: How to fix 5'UTR annotation with RNA-seq data?
I am studying the transcriptome of Arabidopsis. Interestingly, the 5' UTR of the latest annotation is usually too long. Here is an example (this gene is right to left).  You can see both EST (orange bars) and the RNA-seq reads cover a much smaller region of the annotated 5' UTR. Unfortunately, many 5'UTR of the genes have the same problem. Is there a way to fix that? I was wondering if any assembly tools can help me to get the correct 5'UTR ranges. Thanks!!!

# Answer
1. Genes, particularly transcribed genes in Eukaryotes, can have alternate start and stop sites (and splicing). Alternate transcriptional start sites may be found in different tissues or at different developmental periods due to alternate transcription factor binding sites being used. The only reliable way to determine the transcriptional start and stop sites is by read coverage.
2. Genome annotations are generally the product of bioinformatic pipelines and are mostly predictions based upon computational models. Different gene-finders may call different boundaries; for example, in eukaryotic genomes, an exon may be assigned to the wrong up/down-stream gene. Arabadopsis, as model organism, has had it's gene boundaries considerably refined/corrected/validated, but even so, there could still exist errors in the annotation or, more likely, the boundaries indicated are the most common start/stop sites (see point 1 above).

In short, your UTRs are whatever your experimental data (reads) show it to be. Don't worry too much about it not matching up to the standard annotation boundaries.
> 1 votes

---
Tags: rna-seq, sequence-annotation
---


title: Random access on a FASTQ file

text:
# Question
Title: Random access on a FASTQ file
I would like to select a random record from a large set of `n` unaligned sequencing reads in `log(n)` time complexity (big O notation) or less. A record is defined as the equivalent of four lines in FASTQ format. The records do not fit in RAM and would need to be stored on disk. Ideally, I would like to store the reads in a compressed format.

I would prefer a solution that does not require any extra files such as for example a reference genome.

The title of this question mentions a FASTQ only because FASTQ is a common format for storing unaligned reads on disk. I am happy with answers that require a single limited transformation of the data to another file format in time complexity order `n`.

## Update

A clarification: I want the random record to be selected with probability `1/n`.

# Answer
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 to jump to a solution that requires a script or framework (as opposed to just unix command line tools), there is sadly no `samtools fqidx` (perhaps there should be), and existing answers suggest a lot of munging. Whilst they probably work, some appear cumbersome and have many steps in which you could make a mistake.

So, to keep things simple - a quick and dirty alternative approach might be to use `biopython`, seeing as it already has this functionality implemented to do this, and if installed, is trivial to use:

```
from Bio import SeqIO
fq = SeqIO.index("myfastq.fq", "fastq")

```

Once you've acquired an index, you'll get fast random access for any read:

```
# Random access by read name (I like owls)
record = fq["HOOT"]
record
#> SeqRecord(seq=Seq('ACGTACGT', SingleLetterAlphabet()), id='HOOT', name='HOOT', description='HOOT', dbxrefs=[])

# We can get the sequence
record.seq
#> Seq('ACGTACGT', SingleLetterAlphabet())

# and qualities
record.letter_annotations
#> {'phred_quality': [41, 41, 41, 41, 41, 41, 41, 41]}

```

If you want to select arbitrary random records, you could use something like `randrange` to select between 0 and the length of the references list.

```
from random import randint

# Coerce keys to a list, as `dictionary-keyiterator` has no
#   __getitem__ attribute that would allow integer indices
# Note also that this doesn't necessarily guarantee a sorted order
#  but I guess that doesn't matter if you just want random records
key_list = list(fq.keys())

# Select a random key
# Note we use len(key_list)-1 as randint endpoints are inclusive
random_readname = key_list[ randint(0, len(key_list)-1) ]

# Get your record
rand_record = fq.get( random_readname )

```

If you want multiple records, you'll probably want `sample` (to avoid replacement) instead:

```
from random import sample
N = 100
random_indices = sample(xrange(len(key_list)), N)

for key_i in random_indices:
    random_readname = key_list[ key_i ]
    rand_record = fq.get( random_readname )
    # ...

```

For what it's worth, I think `biopython` holds this index in RAM, so if your file is absolutely massive, you might need to be more clever. If that's the case, you could iterate through the `FASTQ` once, and output the readname, file offset and length - akin to a `FASTA` `fai`.

**Update** See `samtools fqidx`
> 7 votes

# Answer
# 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 `tabix` and one with `grabix`. I think the `grabix` solution is more elegant, but I am keeping the `tabix` solution below because `tabix` is a more mature tool than `grabix`.

Thanks to user172818 for suggesting `grabix`.

### Update

This answer previously stated that `tabix` and `grabix` perform lookups in `log(n)` time. After taking a closer look at the grabix source code and the tabix paper I am now convinced that lookups are independent of `n` in complexity. However, both tools use an index that scales in size proportionally to `n`. So, the loading of the index is order `n`. However, if we consider the loading of the index as "...a single limited transformation of the data to another file format...", then I think this answer is still a valid one. If more than one record is to be retrieved, then the index needs to be stored in memory, perhaps with a framework such as pysam or htslib.

## Using `grabix`

1. Compress with `bgzip`.
2. Index the file and perform lookups with `grabix`

In bash:

```
gzip -dc input.fastq.gz | bgzip -c > output.fastq.gz

grabix index output.fastq.gz

# retrieve 5-th record (1-based) in log(n) time
# requires some math to convert indices (4*4 + 1, 4*4 + 4) = (17, 20)
grabix grab output.fastq.gz 17 20

# Count the number of records for part two of this question
export N_LINES=$(gzip -dc output.fastq.gz | wc -l)

```

## Using `tabix`

The tabix code is more complicated and relies on the iffy assumption that `\t` is an acceptable character for replacement of `\n` in a FASTQ record. If you are happy with a file format that is close to but not exactly FASTQ, then you could do the following:

1. Paste each record into a single line.
2. Add a dummy chromosome and line number as the first and second column.
3. Compress with `bgzip`.
4. Index the file and perform lookups with `tabix`

Note that we need to remove leading spaces introduced by `nl` and we need to introduce a dummy chromosome column to keep `tabix` happy:

```
gzip -dc input.fastq.gz | paste - - - - | nl | sed 's/^ *//' | sed 's/^/dummy\t/' | bgzip -c > output.fastq.gz
tabix -s 1 -b 2 -e 2 output.fastq.gz

# now retrieve the 5th record (1-based) in log(n) time
tabix output.fastq.gz dummy:5-5

# This command will retrieve the 5th record and convert it record back into FASTQ format
tabix output.fastq.gz dummy:5-5 | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'

# Count the number of records for part two of this question
export N_RECORDS=$(gzip -dc output.fastq.gz | wc -l)

```

# Random record in constant time

Now that we have a way of retrieving an arbitrary record in `log(n)` time, retrieving a random record is simply a matter of getting a good random number generator and sampling. Here is some example code to do this in python:

## Using grabix

```
# random_read.py
import os
import random

n_records = int(os.environ["N_LINES"]) // 4
rand_record_start = random.randrange(0, n_records) * 4 + 1
rand_record_end = rand_record_start + 3
os.system("grabix grab output.fastq.gz {0} {1}".format(rand_record_start, rand_record_end))

```

## Using tabix

```
# random_read.py
import os
import random

n_records = int(os.environ["N_RECORDS"])
rand_record_index = random.randrange(0, n_records) + 1
# super ugly, but works...
os.system(
    "tabix output.fastq.gz dummy:{0}-{0} | perl -pe 's/^dummy\t\d+\t//' | tr '\t' '\n'".format(
        rand_record_index)
)

```

And this works for me:

```
python3.5 random_read.py

```

# Disclaimer

Please note that `os.system` calls a system shell and is vulnerable to shell injection vulnerabilities. If you are writing production code, then you probably want to take extra precautions.

Thanks to Chris\_Rands for raising this issue.
> 11 votes

# Answer
You could shuffle the FASTQ once and then read sequences off the top of the file as you need them:

```
gzip -dc input.fastq.gz | paste - - - - | shuf | tr '\t' '\n'| gzip -c > output.fastq.gz

```

I would recommend pigz as a replacement for gzip in the compression step if you have it available.

The downside of this approach is that you only get `n` reads before you need to run the shuffle again, and apparently `shuf` holds all data in RAM, so it would die with an out of memory error if the FASTQ file does not fit into RAM as is specified in the question.

Using `sort -R` is complexity `n log(n)` and uses temporary files, so it should not run out of memory:

```
gzip -dc input.fastq.gz | paste - - - - | nl | sort -R | perl -pe 's/\s*\d+\t//' | tr '\t' '\n'| gzip -c > output.fastq.gz

```

The `nl` and `perl` commands are necessary to make sure that identical records are not sorted next to each other.
> 5 votes

# Answer
One possibility is to:

1. reformat the data such that each record is a single line containing the read description, bases, and quality scores
2. pad out each record to a maximum length in each field such that every record in the file is the same number of bytes
3. the total number of records can now be calculated as file size / record size
4. choose a random record number between 0 and the total number of records
5. binary search over the reformatted file until you obtain your read

This would get you the log(n) lookup time you want.

Of course, the data wouldn't be compressed. You could 2-bit encode the bases and quantize the quals to save some space, but that'd be lossy and is perhaps not what you're looking for. Alternatively, you could block-gzip the reformatted data and keep a record of how many blocks are in the file and how many reads are in each block (since the filesize will no longer reflect the number of records in the file). Then to obtain a specific read, you'd calculate the block number it'll appear in, decompress the block, and return the appropriate read.
> 3 votes

# Answer
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 `shuf` fails for lack of sufficient memory.

It requires two passes through the file to do a random sample, but the second pass is generally fast(er) as it uses `mmap` routines to do cached reads.

If you do repeated samples, the repeated samples are also `mmap`-ed (cached) and will run quickly.

You might use it on a FASTQ file like so:

```
$ sample -k 1234 -l 4 in.fq > out.fq

```

It parses the input file into records by every four newline characters (such as the format of a FASTQ file), reading line offset positions into memory. So the memory overhead is relatively very low.

It then applies reservoir sampling on those line offsets to write out a random sample (say, `1234` records in this example) to standard output.
> 2 votes

# Answer
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/2015/05/19/io-performance/

If you just want to grab a single random read, Jared's benchmarks suggest that seeking to a random position in the file and then retrieving the next complete read should be the most performant option.

However, more generally, if you want to take a random subset of reads, there are many factors that affect performance.
> 2 votes

# Answer
To emphasise the issue (as outlined in the `1/n` update), consider an input file with one record that is 5 million bases long, and one records that is 100 bases long. You want an equal probability of selecting either of these two records. Any random seek methods will overwhelmingly pick out the long record.

I expect that indexing the locations of record starts is really the only workable option here, particularly if multi-line records are possible. Create an index file containing the start locations of each record (as identically sized integers, e.g. 64-bit), then sample from the index file (which is an identical length for each record) to fetch the start location. I'd envisage that this file would *only* contain the start locations; any additional metadata (including sequence name) would require seeking in the original file.

Once the indexing is done, the file can be compressed with bgzip, with specific offsets retrieved using the `-b` and `-s` options. However, I expect that compression would not be particularly efficient if multiple random records were desired.
> 1 votes

# Answer
Here's another approach that doesn't require any indexing, using BEDOPS `bedextract` to do a `log(n)` sample on a sorted BED file. Your sample contains random records with equal probability `1/n`.

This approach requires a single `O(n)` pass through the file to transform it to a BED file:

```
$ cat records.fastq | paste - - - - | awk '{ print "chrZ\t"s"\t"(s+1)"$0 }' > records.bed

```

Store the intervals in a separate file:

```
$ cut -f1-3 records.bed > intervals.bed

```

To do a random sample of `k` elements, shuffle the intervals file and preserve the order of shuffled elements.

You can do this with the `sample` tool I outlined earlier:

```
$ sample -k ${K} -s intervals.bed > intervals-sample.bed

```

Or you can `shuf` and `sort-bed` to do the same thing:

```
$ shuf -n ${K} intervals.bed | sort-bed - > intervals-sample.bed

```

There's an `O(klog(k))` cost here, but if `k <<< n`, i.e., you're working with whole-genome scale input, this cost is amortized over the `log(n)` search performance.

Next, use `bedextract` to do a binary search on the records, and delinearize to get back to FASTQ:

```
$ bedextract records.bed intervals-sample.bed | cut -f4 | tr '\t' '\n' > sample.fq

```

With Unix I/O streams, this can be done in one pass:

```
$ sample -k ${K} -s intervals.bed | bedextract records.bed - | cut -f4 | tr '\t' '\n' > sample.fq

```

By baking the sort order into `records.bed`, you're guaranteed the ability to do a binary search, which is `log(n)`.

*Note:* Further, by linearizing the FASTQ input to a BED file and querying on BED intervals, you have equal probability of picking any one interval (interval == FQ record). You can draw an unbiased sample without the hassle of creating and storing a separate index.
> 1 votes

# Answer
I wrote a tool called strandex that matches entire FASTQ records starting at a random offset inside a decompressed file (since non-block compressed gzip files are not seek-able). It was a response to this comment about the `screed` FASTQ parser in khmer, and was initially just a proof of concept but it's been useful for me and a few other people. The general idea is to:

1. seek to a random offset in the file
2. read a chunk of the file
3. test if the regex pattern `@.+[\n\r]+.+[\n\r]+\+.*?[\n\r].+[\n\r]` matches
4. if so, extract the matching FASTQ record using regex capture groups
5. if no match, goto step 2

The CLI script and python library can be installed with `pip install strandex`, and by default the script samples `-n` reads using file offsets starting at a random number (reproducable by setting `-s` seed), then moves a number of bytes to uniformly sample the entire file, determined by the file size. In this way the process is not truly random, but random enough, and avoids sampling any part of the file too much.

If you specify `-n` reads less than the total number in the file you get *down*sampling, and if `-n` is greater than the total number of reads you get *up*sampling.

The advantage of this method over others is - it's a one-pass method, and since it's using the python C internals regex engine it's quite fast.
> 1 votes

# Answer
I believe some of these answers misinterpret the original question. The OP (I believe) doesn't wish to find particular sequences, so indexing is unnecessary. I believe the OP just wants a random sampling of the dataset.

1. If the number of reads you need in small, sure you can seek on the file and return the next complete record. Remember the read IDs to sample without replacement. Jumping around on disk is slow, however, and you'll be reading 4k or whatever your disk buffer size is every time.
2. If the number of reads you need is large, read the entire file and decide for each read whether or not to sample it.
3. The way I do it (for Illumina): estimate the number of reads using the file size and read length and header length, decide how many reads I want, skip the first bunch of reads (due to edge effect, e.g. 20-100k; you can use seek or just read from start and discard), and read the next batch sequentially. I suppose it depends on what you're using the reads for; i.e. if you want random reads from your sample for biological analysis, or you want random reads from your run for QC analysis. This method is more appropriate for the former, but works well enough for the latter too.
> 1 votes

# Answer
I suppose I can provide code for this if needed, but keep in mind that it'd probably have to be done in C.

One could make this easier by using a bgzf compressed fastq file. Yes, BAM uses that already, but since Illumina's software now defaults to producing that there should be a reasonable amount of it already available.

1. Count the number of blocks in the file.
2. Randomly select one of these blocks.
3. Perform reservoir selection on the entries in that block.

That won't be strictly O(log N), but it's lower on memory and IO requirements than parsing through the whole file and doing shuffling or otherwise completely munging the data into an otherwise not-terribly-useful format. This also has the benefit of allowing selecting more than one random entry a bit faster, since you can memoize things.
> 0 votes

# Answer
A quick and dirty solution could be to convert the FASTQ file to two FASTAs, one of the them storing the bases and the other one storing the qualities, and then use standard methods for random access for FASTA.

The only complication is that base qualities can contain `>`. However, this is a problem only when it is the first character of a line and we can fix it by prepending, e.g., `I` to each quality sequence.

If also other people think that this could be a sufficient solution, I will expand the answer and add all the commands.
> 0 votes

---
Tags: sam, fastq, reads, benchmarking
---


# Question
Title: What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files?
I used to work with publicly available genomic references, where basic statistics are usually available and if they are not, you have to compute them only once so there is no reason to worry about performance.

Recently I started sequencing project of couple of different species with mid-sized genomes (~Gbp) and during testing of different assembly pipelines I had compute number of unknown nucleotides many times in both raw reads (in fastq) and assembly scaffolds (in fasta), therefore I thought that I would like to optimize the computation.

* For me it is reasonable to expect 4-line formatted fastq files, but general solution is still prefered
* It would be nice if solution would work on gzipped files as well

Q : What is the fastest way (performance-wise) to compute the number of unknown nucleotides (Ns) in fasta and fastq files?

# Answer
For FASTQ:

```
seqtk fqchk in.fq | head -2

```

It gives you percentage of "N" bases, not the exact count, though.

For FASTA:

```
seqtk comp in.fa | awk '{x+=$9}END{print x}'

```

This command line also works with FASTQ, but it will be slower as awk is slow.

EDIT: ok, based on @BaCH's reminder, here we go (you need kseq.h to compile):

```
// to compile: gcc -O2 -o count-N this-prog.c -lz
#include <zlib.h>
#include <stdio.h>
#include <stdint.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

unsigned char dna5tbl[256] = {
    0, 1, 2, 3,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 5, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
    4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
};

int main(int argc, char *argv[]) {
    long i, n_n = 0, n_acgt = 0, n_gap = 0;
    gzFile fp;
    kseq_t *seq;
    if (argc == 1) {
        fprintf(stderr, "Usage: count-N <in.fa>\n");
        return 1;
    }
    if ((fp = gzopen(argv[1], "r")) == 0) {
        fprintf(stderr, "ERROR: fail to open the input file\n");
        return 1;
    }
    seq = kseq_init(fp);
    while (kseq_read(seq) >= 0) {
        for (i = 0; i < seq->seq.l; ++i) {
            int c = dna5tbl[(unsigned char)seq->seq.s[i]];
            if (c < 4) ++n_acgt;
            else if (c == 4) ++n_n;
            else ++n_gap;
        }
    }
    kseq_destroy(seq);
    gzclose(fp);
    printf("%ld\t%ld\t%ld\n", n_acgt, n_n, n_gap);
    return 0;
}

```

It works for both FASTA/Q and gzip'ed FASTA/Q. The following uses SeqAn:

```
#include <seqan/seq_io.h>

using namespace seqan;

int main(int argc, char *argv[]) {
    if (argc == 1) {
        std::cerr << "Usage: count-N <in.fastq>" << std::endl;
        return 1;
    }
    std::ios::sync_with_stdio(false);
    CharString id;
    Dna5String seq;
    SeqFileIn seqFileIn(argv[1]);
    long i, n_n = 0, n_acgt = 0;
    while (!atEnd(seqFileIn)) {
        readRecord(id, seq, seqFileIn);
        for (i = beginPosition(seq); i < endPosition(seq); ++i)
            if (seq[i] < 4) ++n_acgt;
            else ++n_n;
    }
    std::cout << n_acgt << '\t' << n_n << std::endl;
    return 0;
}

```

On a FASTQ with 4-million 150bp reads:

* The C version: ~0.74 sec
* The C++ version: ~2.15 sec
* An older C version without a lookup table (see the previous edit): ~2.65 sec
> 26 votes

# Answer
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:

1. R with the `ShortRead` package (even if not the fastest, certainly a super convenient method).
2. A pipeline of `grep -v "^>" | tr -cd A | wc -c`
3. A pipeline of `grep -v "^>" | grep -io A | wc -l`
4. A pipeline of `seqtk comp | awk`
5. Something custom in C/C++

The R code can't actually count N records for some reason when I try it, so I counted `A` for everything.

I'll exclude all python solutions, because there's no conceivable universe in which python is fast. For the C/C++ solution I just used what @user172818 posted (everything I tried writing took the same amount of time).

Repeating 100 times and taking the average (yes, one should take the median):

1. 1.7 seconds
2. 0.65 seconds
3. 15 seconds
4. 1.2 seconds
5. 0.48 seconds

Unsurprisingly, anything in straight C or C++ is going to win. `grep` with `tr` is surprisingly good, which surprises me since even though `grep` is very very optimized I still expected it to have more overhead. Piping to `grep -io` is a performance killer. The R solution is surprisingly good given the typical R overhead.

**Update1:** As suggested in the comments

6. `sum(x.count('A') for x in open('seqs.fa','r') if x[0] != '>')` in python

This yields a time of

6. 1.6 seconds (I'm at work now, so I've adjusted the time as best I can to account for the different computers)

**Update2:** More benchmarks from the comments

7. `awk -F A '!/^>"/ {cnt+=NF-1}END{print cnt}' foo.fa`
8. `perl -ne 'if(!/^>/){$count += tr/A//} END{print "$count\n"}' foo.fa`

These yield times of:

7. 5.2 seconds
8. 1.18 seconds

The awk version is the best hack I could come up with, there are likely better ways.

**Update 3**: Right, so regarding fastq files, I'm running this on a 3.3GB gzipped file with 10 repetitions (this takes a bit of time to run), so I'm not going to initially limit tests to commands that I can trivially modify to handle compressed files (after all, uncompressed fastq files are an abomination).

1. `seqtk fqchk`
2. `sed -n '2~4p' <(zcat foo.fastq.gz) | grep -io A | wc -l`
3. `sed -n '1d;N;N;N;P;d' <(zcat Undetermined_S0_R1_001.fastq.gz) | tr -cd A | wc -c`
4. The C example from user1728181 (the comparison to C++ is already provided there, so no need to include it).
5. `bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' foo.fastq.gz`

The average times are:

1. 1 minute 44 seconds
2. 2 minutes 6 seconds
3. 1 minute 52 seconds
4. 1 minute 15 seconds
5. 3 minutes 8 seconds

N.B., these will often not handle fastq files with entries spanning more than 4 lines. However, such files are the dodo birds of bioinformatics, so that's not practically important. Also, I couldn't get `ShortRead` to work with the compressed fastq file. I suspect there's a way to fix that.
> 18 votes

# Answer
I think that this should be pretty fast:

FASTA:

```
grep -v "^>" seqs.fa | tr -cd N | wc -c

```

FASTQ:

```
sed -n '1d;N;N;N;P;d' seqs.fq | tr -cd N | wc -c

```

See this answer on SO about how to count characters in BASH using different approaches.
> 10 votes

# Answer
**FASTQ**

As it was pointed out, fastq can be complicated. But in a simple case when you have four lines per record, one possible solution in bash is:

```
sed -n '2~4p' seqs.fastq | grep -io N | wc -l

```
* `sed -n '2~4p'` will print every fourth line
* `grep -o N` will output a line with `N` for every matching symbol
* `wc -l` will count the lines

I suspect this python approach will work faster:

```
cat seqs.fastq | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin))"

```

**FASTA**

Coreutils:

```
grep -v "^>" seqs.fasta | grep -io N | wc -l

```

Python:

```
cat seqs.fasta | python3 -c "import sys; print(sum(line.upper().count('N') for line in sys.stdin if not line.startswith('>')))"

```
> 8 votes

# Answer
## Using bioawk

With bioawk (counting "A" in the *C. elegans* genome, because there seem to be no "N" in this file), on my computer:

```
$ time bioawk -c fastx '{n+=gsub(/A/, "", $seq)} END {print n}' genome.fa
32371810

real    0m1.645s
user    0m1.548s
sys     0m0.088s

```

bioawk is an extension of awk with convenient parsing options. For instance `-c fastx` makes the sequence available as `$seq` in fasta and fastq format (**including gzipped versions**), so the above command **should also handle fastq format robustly**.

The `gsub` command is a normal awk function. It returns the number of substituted characters (got it from https://unix.stackexchange.com/a/169575/55127). I welcome comments about how to count inside awk in a less hackish way.

On this fasta example it is almost 3 times slower than the `grep`, `tr`, `wc` pipeline:

```
$ time grep -v "^>" genome.fa | tr -cd "A" | wc -c
32371810

real    0m0.609s
user    0m0.744s
sys     0m0.184s

```

(I repeated the timings, and the above values seem representative)

## Using python

### readfq

I tried the python readfq retrieved from the repository mentioned in this answer: https://bioinformatics.stackexchange.com/a/365/292

```
$ time python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('A') for _, seq, _ in readfq(stdin)))" < genome.fa
32371810

real    0m0.976s
user    0m0.876s
sys     0m0.100s

```

### "pure python"

Comparing with something adapted from this answer: https://bioinformatics.stackexchange.com/a/363/292

```
$ time python -c "from sys import stdin; print sum((line.count('A') for line in stdin if not line.startswith('>')))" < genome.fa
32371810

real    0m1.143s
user    0m1.100s
sys     0m0.040s

```

(It is slower with python 3.6 than with python 2.7 on my computer)

### pyGATB

I just learned (08/06/2017) that GATB includes a fasta/fastq parser and has recently released a python API. I tried to use it yesterday to test another answer to the present question and found a bug. This bug is now fixed, so here is a pyGATB-based answer:

```
$ time python3 -c "from gatb import Bank; print(sum((seq.sequence.count(b\"A\") for seq in Bank(\"genome.fa\"))))"
32371810

real    0m0.663s
user    0m0.568s
sys     0m0.092s

```

(You can also do `sequence.decode("utf-8").count("A")` but this seems a little slower.)

Although I used python3.6 here (pyGATB seems python3-only), this is faster than the other two python approaches (for which the reported timings are obtained with python 2.7). This is even almost as fast as the `grep`, `tr`, `wc` pipeline.

### Biopython

And, to have even more comparisons, here is a solution using `SeqIO.parse` from Biopython (with python2.7):

```
$ time python -c "from Bio import SeqIO; print(sum((rec.seq.count(\"A\") for rec in SeqIO.parse(\"genome.fa\", \"fasta\"))))"
32371810

real    0m1.632s
user    0m1.532s
sys     0m0.096s

```

This is a bit slower than the "pure python" solution, but perhaps more robust.

There seems to be a slight improvement with @peterjc's suggestion to use the lower level `SimpleFastaParser`:

```
$ time python -c "from Bio.SeqIO.FastaIO import SimpleFastaParser; print(sum(seq.count('A') for title, seq in SimpleFastaParser(open('genome.fa'))))"
32371810

real    0m1.618s
user    0m1.500s
sys     0m0.116s

```

(I did a series of timings and tried to take one that seemed representative, but there's a lot of overlap with the higher-level parser's timings.)

### scikit-bio

I read about scikit-bio today (23/06/2017), which has a sequence reader.

```
$ time python3 -c  "import skbio; print(sum(seq.count('A') for seq in skbio.io.read('genome.fa', format='fasta', verify=False)))"
32371810

real    0m3.643s
user    0m3.440s
sys     0m1.228s

```

### pyfastx

I read about pyfastx today (27/08/2020), which has a sequence reader.

```
$ time python3 -c "from pyfastx import Fasta; print(sum(seq.count('A') for (_, seq) in Fasta('genome.fa', build_index=False)))"
32371810

real    0m0.781s
user    0m0.740s
sys     0m0.040s

```

(The timing is done using python 3.6, on the same machine as previous timings from 3 years ago, so it should be comparable.)

### `fastx_read` from mappy

I realize (27/08/2020) that I hadn't added a benchmark for `fastx_read` from mappy which I regularly use since quite some time already:

```
$ time python3 -c "from mappy import fastx_read; print(sum(seq.count('A') for (_, seq, _) in fastx_read('genome.fa')))"
32371810

real    0m0.731s
user    0m0.648s
sys     0m0.080s

```

## Benchmarks on a fastq.gz file

I tested some of the above solutions (or adaptations thereof), counting "N" in the same file that was used here: https://bioinformatics.stackexchange.com/a/400/292

### bioawk

```
$ time bioawk -c fastx '{n+=gsub(/N/, "", $seq)} END {print n}' SRR077487_2.filt.fastq.gz
306072

real    1m9.686s
user    1m9.376s
sys     0m0.304s

```

### pigz + readfq python module

readfq doesn't complain and is very fast when I pass directly the compressed fastq, but returns something wrong, so don't forget to manually take care of the decompression.

Here I tried with `pigz`:

```
$ time pigz -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))"
306072

real    0m52.347s
user    1m40.716s
sys     0m8.604s

```

The computer has 16 cores, but I suspect the limiting factor for `pigz` is reading from the disk: the processors are very far from running full speed.

And with `gzip`:

```
$ time gzip -dc SRR077487_2.filt.fastq.gz | python -c "from sys import stdin; from readfq import readfq; print sum((seq.count('N') for _, seq, _ in readfq(stdin)))"
306072

real    0m49.448s
user    1m31.984s
sys     0m2.312s

```

Here gzip and python both used a full processor. Resource usage balance was slightly better. I work on a desktop computer. I suppose a computing server would take advantage of `pigz` better.

### pyGATB

```
$ time python3 -c "from gatb import Bank; print(sum((seq.sequence.count(b\"N\") for seq in Bank(\"SRR077487_2.filt.fastq.gz\"))))"
306072

real    0m46.784s
user    0m46.404s
sys     0m0.296s

```

### biopython

```
$ time python -c "from gzip import open as gzopen; from Bio.SeqIO.QualityIO import FastqGeneralIterator; print(sum(seq.count('N') for (_, seq, _) in FastqGeneralIterator(gzopen('SRR077487_2.filt.fastq.gz'))))"
306072

real    3m18.103s
user    3m17.676s
sys     0m0.428s

```

### pyfastx (added 27/08/2020)

```
$ time python3 -c "from pyfastx import Fastq; print(sum(seq.count('N') for (_, seq, _) in Fastq('SRR077487_2.filt.fastq.gz', build_index=False)))"
306072

real    0m51.140s
user    0m50.784s
sys     0m0.352s

```

### `fastx_read` from mappy (added 27/08/2020)

```
$ time python3 -c "from mappy import fastx_read; print(sum(seq.count('N') for (_, seq, _) in fastx_read('SRR077487_2.filt.fastq.gz')))"
306072

real    0m52.336s
user    0m52.024s
sys     0m0.300s

```

## Conclusion

pyGATB and bioawk both handle transparently compression (gzipped or not) and format differences (fasta or fastq). pyGATB is quite a new tool, but seems more efficient compared to the other python modules I tested.

(Update 27/08/2020) pyfastx and mappy seem quite convenient too, and are only slightly slower than pyGATB.
> 6 votes

# Answer
### Decompression of gzipped FASTQ is the main issue

If we take real world gzipped FASTQ (which as the OP suggested would be beneficial) rather than trivial FASTA as the starting point then the real issue is actually decompressing the file not counting the Ns and in this case the C program `count-N` is no longer the fastest solution.

Additionally it would be good to use a specific file for benchmarking which actually has Ns, because you'll get some quite interesting execution time differences with some methods counting the more frequently occurring As rather than Ns.

One such file is:

http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence\_read/SRR077487\_2.filt.fastq.gz

It's also worth checking the various solutions return the correct answer there should be 306072 Ns in the above file.

Next note that decompression of this file redirected to `/dev/null` is slower with `zcat` and `gzip` (which are both gzip 1.6 on my system) than say a parallel implementation of gzip like Mark Adler's `pigz`, which appears to use 4 threads for decompression. All timings represent an average of 10 runs reporting real (wall clock time).

```
time pigz -dc SRR077487_2.filt.fastq.gz > /dev/null

real    0m29.0132s

time gzip -dc SRR077487_2.filt.fastq.gz > /dev/null

real    0m40.6996s

```

There is an ~11.7 second difference between the two. Next if I then try to benchmark a one-liner which performs on FASTQ and gives the correct answer (Note I've yet to encounter FASTQ which is not 4 line, and seriously who generates these files!)

```
time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c
306072

real    0m34.793s

```

As you can see the counting adds ~5.8 second to the total run time versus `pigz` based decompression. Additionally this time delta is higher when using `gzip` ~6.7 seconds above gzip decompression alone.

```
time gzip -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | tr -cd N | wc -c
306072

real    0m44.399s

```

The `pigz awk tr wc` based solution is however ~4.5 seconds faster than the `count-N` based C code solution:

```
time count-N SRR077487_2.filt.fastq.gz
2385855128      306072  0

real    0m39.266s

```

This difference appears to be robust to re-running as many times as I like. I expect if you could use pthread in the C based solution or alter it to take the standard out from pigz it would also show an increase in performance.

Benchmarking another alternative `pigz grep` variant appears to take more or less the same time as the `tr` based variant:

```
time pigz -dc SRR077487_2.filt.fastq.gz | awk 'NR%4==2{print $1}' | grep -o N | wc -l
306072

real    0m34.869s

```

Note that the `seqtk` based solution discussed above is noticeably slower,

```
time seqtk comp  SRR077487_2.filt.fastq.gz | awk '{x+=$9}END{print x}'
306072

real    1m42.062s

```

However it's worth noting `seqtk comp` is doing a bit more than the other solutions.
> 6 votes

# Answer
Honestly, the easiest way (especially for FASTQ) is probably to use a dedicated parsing library, such as R/Bioconductor:

```
suppressMessages(library(ShortRead))

seq = readFasta(commandArgs(TRUE)[1]) # or readFastq
cat(colSums(alphabetFrequency(sread(seq))[, 'N', drop = FALSE]), '\n')

```

This may not be the fastest, but it’s *pretty fast*, since the relevant methods are highly optimised. The biggest overhead is actually probably from loading `ShortRead`, which is an unjustifiable slog.
> 5 votes

# Answer
If it is raw speed you're after, then writing an own little C/C++ program is probably what you need to do.

Fortunately, the worst part (a fast and reliable parser) has already been tackled: the readfq from Heng Li is probably the fastest FASTA/FASTQ parser around.

And it's easy to use, the example on GitHub can easily be expanded to do what you need. Just add in a parameter parser (for the filename you want to analyse), program a simple 'N'-counter loop and have the results printed to stdout. Done.

Else, if it is simplicity you are after, see Konrad's answer for an R based solution. Or a very simple script using Biopython.
> 5 votes

# Answer
For FASTA files, I've implemented a relatively efficient method in pyfaidx v0.4.9.1. This post made me realize that my previous code was quite slow and easy to replace:

```
$ pip install pyfaidx
$ time faidx -i nucleotide ~/Downloads/hg38.fa
name    start   end A   T   C   G   N   others
chr1    1   248956422   67070277    67244164    48055043    48111528    18475410
chr10   1   133797422   38875926    39027555    27639505    27719976    534460
chr11   1   135086622   39286730    39361954    27903257    27981801    552880
chr11_KI270721v1_random 1   100316  18375   19887   31042   31012   0
...
chrX    1   156040895   46754807    46916701    30523780    30697741    1147866
chrY    1   57227415    7886192 7956168 5285789 5286894 30812372
chrY_KI270740v1_random  1   37240   8274    12978   7232    8756    0

real    2m28.251s
user    2m15.548s
sys 0m9.951s

```
> 3 votes

# Answer
Here are perl 1-liners for fasta and fastq files that are fast enough if you just want totals:

`$ perl -ne 'BEGIN { my $n=0 }; /^>/ or $n += s/N//gi; END { print "N=$n\n" }' contigs.fasta`

For each line of a FASTA file, if not header (/^\>/) then it's a sequence line. Substitute Ns, which returns the number of substitutions made, which is added to the total.

`$ zcat reads.fastq | perl -ne 'BEGIN { my $n=0 }; /^[ATCGN]+$/i and $n += s/N//gi; END { print "N=$n\n" }'`

For each sequence line of a FASTQ file, matched by /^\[ATCGN\]+$/i (i.e. all NTs), then count substitutions made in the same way.

You may use pigz instead of zcat (or gunzip -c) if you have pigz installed.

An uncompressed fastq file of 14M 250bp Illumina reads completed in 40sec.
> 0 votes

---
Tags: fasta, fastq, benchmarking
---


# Question
Title: Why can using more reads lead to a lower quality assembly?
I am experimenting with adding additional reads to the input files I'm giving SOAPdenovo2, and there comes a point where a good contig I've been watching actually stops showing up. Does anyone have a quick answer to why that might happen?

# Answer
> 4 votes
The more reads you add the more errors you add into the assembly. This is because additional duplicate reads don't add nodes/edges to the de Bruijn graph, but those with errors do. By preselecting those reads aligning to your desired contigs, you may be further exacerbating this effect and further hindering SOAPdenovo2's ability to do k-mer correction (after all, the observed k-mer frequency of errors would be higher). This then results in additional incorrect paths through the graph that the assembler has to traverse.

# Answer
> 0 votes
I think the reason that no one has submitted an actual answer yet is because this question is really dependent on the size of the genome you are sequencing, its repetitive content, the level of heterozygosity, the amount of contamination, the type of data (illumina, long illumina, long reads), the library type, library peculariaties, PCR error, the heterozygosity of the contig in question, the complexity of the contig in question (high information content? repeats? AT/GC%?), and probably some other things that I haven't thought of.

Like @user172818 said, all of those things ON TOP OF all of the quirks of different assemblers make this a hard question to answer for you without a full diagnosis of your data and the genome of the thing you are sequencing/assembling.

# Answer
> 0 votes
Yes, this is commonly observed and both the previous answers (Devon Ryan and conchoecia) are spot on so I won't reiterate their explanations on why.

Try bfc (bloom-filter correction) to "fix" likely erroneous base-calls (i.e. uncharacteristically rare kmers) in the reads prior to assembly. This will not only help your assembly, but it shall require less RAM and time.

I also recommend trying running your assembler with multiple kmer parameters, as some values will work better than others for a particular dataset. You can simply choose the best assembly or go through a more complicated effort to merge them, depending on your circumstances.

---
Tags: assembly
---


# Question
Title: Can I use my computer while MinION is sequencing without negatively affecting the run?
I have a Mac Book Pro and I am about to start a sequencing run on the MinION. MinKNOW is going to be running for 2 days. Can I use my computer during sequencing? Can I browse the web and use Excel, etc? I am not speaking about running a parallelised genome assembly on all my cores using my entire RAM but more mundane activities, like checking my emails or taking notes. Furthermore, can I lock my screen using Control + Shift + OFF without negatively affecting the sequencing?

# Answer
Yes, that should be fine, especially if basecalling is disabled. Locking a screen should be fine as well.

I prefer to do basecalling after a run anyway, because there's a chance that it could chew up too much processor time and result in lots of skipped reads.
> 4 votes

# Answer
Since biological samples and reagents cost time and money, I wouldn't push it too much. Checking your email would be fine, but I wouldn't be browsing the web (besides gmail/calendar) in case a plug-in (e.g. flash), bad javascript, etc. causes any trouble.

I'd disable any auto-updates too.

Locking the screen would be fine.

Personally, I'd buy a cheap chromebook and let the Mac work.
> 1 votes

---
Tags: nanopore, ngs, minion, minknow
---


# Question
Title: Why does cutadapt remove low quality bases from the ends of reads only?
I use cutadapt to remove low quality bases from my Illumina reads. 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. Why? Why doesn't the algorithm remove all low quality bases? It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

# Answer
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.

According to its user guide, cutadap is designed this way: it trims off bases from the 3'-end until it sees a base with quality higher than a threshold. This is not a good algorithm. For example, when we see a quality string `30,30,30,3,3,3,3,3,20,3,3`, we would prefer to trim it down to `30,30,30` because a `3,3,3,3,3,20` tail is still not useful.

Some tools such as trimmomatic uses a sliding window to avoid this issue. In my view, a better algorithm is the so-called the "modified Mott algorithm" used by phred 20 years ago. seqtk among others implements this algorithm. Note that the Mott algorithm always trims from both ends, which is sometimes not preferred. BWA implements a variant of this algorithm, trimming from the 3'-end only.

In practice, though, different quality trimming algorithms probably work equally well because it is relatively rare to see a high-quality base in the middle of a low-quality tail. In addition, modern aligners can well handle low-quality tails; assemblers and error correctors can correct through such tails, too. It is usually not necessary to apply quality trimming.

> It could replace low quality bases in the middle of a read by an N, for example. Or is it that the low quality base is still most probably correct and hence one does not want to lose the base information for mapping?

With a low-quality base, you might have an error/mismatch, but with an N, you always have a mismatch. Masking low-quality bases to N is not as good.

EDIT: responding to the following comment from OP:

> What if the base is low-quality and then it matches the reference by mistake? Wouldn't it be better if it did not match and was penalized? I guess this should be modeled mathematically.

If a read could be mismapped due to a sequencing error, its true location is often in the genome. In this case, a competent mapper will give the alignment a low mapping quality (mapQ). A quality-aware mapper such as novoalign will penalize mapQ even more. In comparison, if you hard mask this sequencing error to "N", you will get a mapQ=0. You can see the difference between the two approaches mostly comes from mapQ.

In modern Illumina data, it is quite often to see a Q8 base in the middle of high-quality bases. \>80% of them (in theory) are still correct. My *hunch* is that hard masking all of them would lead to considerable data loss.
> 7 votes

# Answer
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. This is a general trend seen in all runs, therefore you can remove these low quality bases from the end of your run. If you have incidentally a bad quality base in the middle of a read this is not the same general trend but is sporadic. Which is hard to automatically remove with tools such as cutadapt. If you really want an N for each low quality bases within the read, you can make a simple program for that, but why would you need that? Most aligners can handle a few bad quality bases within a read, and for variant calling you'll need many reads (so the incidental error base will be overruled by the good ones).
> 9 votes

# Answer
For Illumina (and 454) reads, the quality decreases with read length. It's not linear and is run/library-dependent. It has less to do with exhaustion of reagents and more to do with the strands on a spot being out of phase due to incomplete/missed base incorporation during sequencing.

It is common practice to trim off the low quality 3' ends as the first step in QC. Most analyses and most can deal with a few low errors in the middle of the read and what you do about such sequencing errors may vary depending on what type and use of the data. But in general, nobody wants the garbage at the end of the reads.

P.S. Even for the same library prep, the fwd/rev runs will have different quality curves. P.P.S. Sometimes the first couple bases are also poor too.
> 1 votes

---
Tags: reads, quality-control, trimming, qc, cutadapt
---


title: Pooling data in metagenome assembly

text:
# Question
Title: Pooling data in metagenome assembly
I have 12 human gut microbiome WGS Nextseq reads (151 bp paired end). What will be an effective strategy to assemble a metagenome?

Let us say I have already filtered the fastq for quality, adapter sequence and host contamination (human, in this case).

1) Should I concatenate all the R1 reads as one single R1 read and one Single R2 read?

```
cat Sample[1..12].R1 > Single_R1.fastq
cat Sample[1..12].R2 > Single_R2.fastq

```

and then use Diginorm to normalize the Single\_R1.fastq and Single\_R2.fastq. Subsequently, feed these fastq files into any metagenome assembler such as Megahit, MetaSPAdes?

Normalize the output by using CD-HIT or similar tool to remove duplicates and filter by contig length.

OR

2) Perform metagenome assembly for each of the samples individually after applying filtering, removing adapters and host contamination.

```
R1=(*_R1_001.filtered.fastq)
R2=(*_R2_001.filtered.fastq)
for ((i=0;i<=${#R1[@]};i++)); do
  /bin/metagenome-assembler -1 "${R1[i]}" -2 "${R2[i]}" -o ${R1[i]%.*}.contigs.fa;
done

```

Followed by combining all the contigs.fa into one mega\_contigs.fa

```
cat *.contigs.fa > Mega_contigs.fa

```

and use CD-HIT or similar tool to remove duplicates.

# Answer
> 2 votes
With metagenomes, there's several strategies and no one-size-fits-all. You can spend a lot of time playing with the data to try to get the "best assembly." But it depends on what you're trying to do, e.g. whether or not you are mostly interested in the highest abundant species or are interested in the less abundant ones too.

But to propose a strategy that seems reasonable and hasn't been mentioned yet, my initial thoughts are that since the gut microbiome isn't complex, you can likely get decent assemblies from each run separately. Then you can run a contig clustering algorithm (e.g. MetaBAT) to bin the contigs into genomes. Next compare the bins from the various assemblies (i.e. cluster the clusters).

It may or may not be worth the trouble to try to combine the related bins' contigs into a consensus (pan)genome sequence. If you want to, methods that have been used successfully include shredding the contigs into overlapping chunks and assembling with an OLC assembler. But if you're interested in KO:read-counts (for example), a consensus genome sequence isn't really required. Cross-mapping the unassembled reads can be useful too (i.e. a sample may have a low-abundance OTU that doesn't produce good contigs, but it can be associated with an OTU that is better assembled in another sample).

If you provide more detail on your goal, I can comment further.

# Answer
> 1 votes
If your goal is to bin the resulting contigs into genomes, then you should do option #1, pooling the reads and assembling into one set of contigs.

# Answer
> 1 votes
I also think pooling is the better option, followed by partitioning by coverage/taxonomy of contigs.

Maybe check out BlobTools, which helps with filtering read-pairs by taxonomy of contigs they contribute to and also does nice visualisations of assemblies.

Workflow B seems to be what you want.

*Disclaimer: I am the developer of this tool.*

---
Tags: assembly, metagenome
---


# Question
Title: Does the DNA or RNA of some ethnic groups map better than others' to the human reference genome sequence?
I believe that the human genome reference took DNA samples from different people and that some natural variation is included in extra contigs.

However, the human reference genome comes from a limited number of people.

Can the limited number of donors and the linear representation of the genome (as opposed to a graph representation that could potentially include variation) bias an analysis in favour of one ethnic group over another? Are there studies that have compared the alignment rates/mappabilities of reads of people coming from different geographic regions? Do some ethnic groups' reads align better to the reference?

# Answer
For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.

In studies where they compare different ethnic groups (e.g., 1000 G) they rather not discuss such matters, but act like the analysis is 100% perfect.

There is a nice comment in Genome Biology about this issue.
> 6 votes

---
Tags: sequence-alignment, genome, human-genome, mappability
---


title: Searching for protein motifs with MEME

text:
# Question
Title: Searching for protein motifs with MEME
I am trying to look for motifs in a set of proteins returned by a mass spec experiment using MEME. The experiment was carried out in the yeast Candida albicans.

Does anyone know if it is worth calculating either a custom background model or a position specific prior using the non-hit proteins? I don't generally bother for nucleotide motif searching, but I've never done protein motif searching before.

Also advice of whether to use the ZOOPS or ANR models appreciated: is there any reason not to use ANR other than run time?

# Answer
> 2 votes
Some answers by the author Charles Grant via the MEME Suite Google group may be useful to you.

Via: https://groups.google.com/forum/#!topic/meme-suite/yNascbE8Tig/discussion

> Providing your own higher-order background model can greatly improve MEME's ability to discover motifs...
>
> The command-line utility `fasta-get-markov`, included in the MEME Suite download, is used to generate custom background Markov models. The input to `fasta-get-markov` is a FASTA file containing "background" sequences. Ideally, these "background" sequences will be different from the sequences you are analyzing with MEME, but as similar in nature as possible. For example, if you wanted to discover motifs in certain intergenic regions, you might use other sequence data from other intergenic regions to generate the background. The larger the set of "background" sequences is, the better the results will be...
>
> Typically, you should not specify an order larger than 3 for DNA sequences, or larger than 2 for protein sequences. However, if your input sequences contain higher-order non-random effects that are getting in the way of motif finding, you can follow the following "rules of thumb":
>
> * Use a background model at least four orders less than the shortest motifs you are looking for. So, if you want to find motifs as short as six, I wouldn't use a model higher than order two.
> * For an accurate model of order N, you need to use a FASTA file as input to `fasta-get-markov` with at least 10 times 4\**(N+1) DNA characters in it.
>
> So,
>
> * order-3 requires 2560 characters
> * order 4 requires 10240 characters
> * order 5 requires 40960 characters etc.

I'm not certain if that input requirement scales differently for proteins (i.e., 20^(N+1)). That may be worth following up on the Google group.

And via: https://groups.google.com/forum/#!topic/meme-suite/UiDhdY6WlK4/discussion

> You should choose the motif occurrence model based on your understanding of the experimental protocol that generated the sequence data and the underlying biology. MEME is trying to find short sequences that are statistically over-represented in your sequence data. To do this, it has to assume a model for how many occurrences of a motif there will be in each sequence. The nature of your experiment should be the basis for the model you choose. For example, data from a ChIP-Seq experiment will typically contain exactly one occurrence of a motif per sequence. On the other hand, sequences drawn from a database of upstream regions of co-regulated genes might contain any number of motif occurrences, including 0.

Via: https://groups.google.com/d/topic/meme-suite/QJiJsy1QxYk/discussion

> The MEME algorithm run time is cubic with respect to the number of input sequences, therefore, it is unsuitable for OOPS (only one per sequence) analyses that have more than 1000 sequences. For ZOOPS (zero or one per sequence) and ANR (any number of repeats) analyses, having more than 1000 sites are also intractable in terms of run time. For example, on a 3 GHz CPU, initialization of the p-value table takes about 10 min for 1000 sites/sequences in OOPS, about a day for 5000 sites, and about 1 week for 10,000 sites. Then the quadratic run time with respect to the total number of characters kicks in next.

Thus it seems like the ZOOPS and ANR models have similar runtime complexity, and use of the MPI-based parallel `meme_p` binary may be helpful, depending on your analysis:

> The parallel version of MEME scales up to about 128 processors. Please see http://www.sdsc.edu/~tbailey/papers/cabios96.pdf for a discussion of the parallel program. You must bear in mind however, that doubling the total number of sequences input to MEME means that you will need 8 times more processors for the job to finish in the same amount of time.

---
Tags: proteins, motifs
---


# Question
Title: Can pair end reads with high heterozygosity be used to polish PacBio assembly?
I have used canu (correct) and smartdenovo to assembly PacBio reads. Next I am going to polish my assembly using illumina pair end reads by Pilon. However, I found there is high heterozygosity in my pair end reads when I used jellyfish to do k-mer analysis. There are two peaks in the k-mer frequencies distribution plot. The size of my genome is about 190M, and it is a diploid organism.
What I am doubting is that can I use this pair end data to polish my assembly?

If not, I have to do sequencing again.

# Answer
Resequencing your sample is not going to make it any less heterozygous. I am assuming that you are working with a diploid organism? Please update your question and specify the genome size and ploidy of your organism (the type of organism will help us update our answers, too).

The answer, yes, you can use Pilon + short reads to correct the PacBio assembly even with a highly heterozygous sample. To do so you should use the `--diploid` option in Pilon to let it know that there will be heterozygous sites. The way that it corrects heterozygous SNPs and heterozygous indels is different, and it seems to perform best with SNPs and sometimes does not correctly identify heterozygous indels (source).

Pilon itself doesn't handle phasing, but just corrects with whatever reads are in the `bam` file that it is provided with. The end result of providing a `bam` with reads from both haplotypes is that the corrected sequence will be a mosaic of the dominant allele at the heterozygous site as was observed in the short read data. This may or may not select for some alleles based on sequencing coverage bias of Illumina sequencers.
> 3 votes

---
Tags: pacbio
---


title: Installing Rostlab's PredictProtein

text:
# Question
Title: Installing Rostlab's PredictProtein
I am trying to install Rostlab's PredictProtein in a Virtual Machine as detailed here https://rostlab.org/owiki/index.php/How\_to\_Run\_The\_PredictProtein\_Machine\_Image\_with\_VirtualBox:

> 1. Download PredictProtein Machine Image to your local drive. You will need the Bio-Linux-based 'VMWARE/VMDK' format version.
>     * Choose the Debian-based version in case you do not need a graphical user interface. This image is only half the size of the Bio-Linux-based image.
> 2. Download a database image to your local drive. We recommend the 'raw' format version.
> 3. Create New Virtual Machine: Machine -\> New
> 4. Name your virtual machine as 'PredictProtein'
> 5. Choose Linux and Debian as operation system
> 6. Locate and select the VMware .vmdk file as your boot disk
> 7. Give at least 2GB of RAM to your virtual machine
> 8. Click 'Finish'
> 9. Boot your new virtual machine.
> 10. If the message "The disk drive for /mnt/local-storage is not ready yet or not present" appears, press 's'.

The problem comes at this step:

> 11. Open a terminal and enter:
>
> ```
> sudo mkfs.ext4 -F -m 0 -L local-storage /dev/sdb
> sudo mount -a
> sudo mkdir /mnt/local-storage/rostlab-data && cd /mnt/local-storage/rostlab-data
> sudo mount -a
> sudo mount /mnt/rostlab-data
> sudo tar -xvJf /mnt/rostlab-data/rostlab-data.txz
> sudo umount /mnt/rostlab-data && cd
>
> ```

When I run `sudo mount /mnt/rostlab-data`, I get:

```
mount: special device LABEL=rostlab-data does not exist

```

How should I proceed?

Additional notes -

The output for

```
cat /etc/fstab

```

is

```
# /etc/fstab: static file system information.
#
# Use 'blkid -o value -s UUID' to print the universally unique identifier
# for a device; this may be used with UUID= as a more robust way to name
# devices that works even if disks are added and removed. See fstab(5).
#
# <file system> <mount point>   <type>  <options>       <dump>  <pass>
proc            /proc           proc    nodev,noexec,nosuid 0   0
# / was on /dev/sda1 during installation
UUID=9a9418ff-48f9-421a-b25a-82161839351d / ext4 errors=remount-ro 0 1
# swap was on /dev/sda5 during installation
UUID=a96757f4-8cc7-410f-95ee-b5ae642eb9ed none swap sw  0       0

LABEL=local-storage /mnt/local-storage auto defaults,nofail 0   2
LABEL=rostlab-data /mnt/rostlab-data auto defaults,noauto,ro 0  0

/mnt/local-storage/rostlab-data /usr/share/rostlab-data none bind,nofail 0 0

```

The output for

```
blkid

```

is

```
/dev/ramzswap0: TYPE="swap"
/dev/sda1: UUID="9a9418ff-48f9-421a-b25a-82161839351d" TYPE="ext4"
/dev/sda5: UUID="a96757f4-8cc7-410f-95ee-b5ae642eb9ed" TYPE="swap"
/dev/sdb: LABEL="local-storage" UUID="df841028-8dd1-4cd2-b1bb-6ae1f6887510" TYPE="ext4"

```

# Answer
It looks like the image you downloaded is missing the partition with the label `rostlab-data`. Your `blkid` output shows that you have the `local-storage` partition, but not the `rostlab-data` which is automatically added to your `/etc/fstab`.

I am afraid you'll have to contact the developers and let them know. Sounds like a bug.
> 0 votes

---
Tags: proteins, software-installation, linux
---


title: ROC analysis, CAFA 2 experiment

text:
# Question
Title: ROC analysis, CAFA 2 experiment
My question is about ROC curves used in CAFA2 experiment. In this paper they used the ROC analysis for the term-centric evaluation. In order to perform ROC curve analysis we should have a continuous variable and a classifier (categorical) variable. In CAFA2 experiment, which variables did they consider as a categorical and continuous? Thank you.

# Answer
The paper doesn't seem to have been written for general readers. Unfortunately I don't understand the equations enough to be able to explain what the categorical and continuous variables are. Just in case it helps, here are the equations from the paper showing the precision/recall calculations:

$$\begin{array}{@{}rcl@{}} \text{pr}(\tau) &=& \frac{1}{m(\tau)}\sum\limits\_{i=1}^{m(\tau)} \frac{{\sum\nolimits}\_{f} \unicode{x1D7D9}\left(f \in P\_{i}(\tau) \wedge f \in T\_{i}\right)}{\sum\_{f} \unicode{x1D7D9}\left(f \in P\_{i}(\tau) \right)},\\ \text{rc}(\tau) &=& \frac{1}{n\_{e}}\sum\limits\_{i=1}^{n\_{e}} \frac{{\sum\nolimits}\_{f} \unicode{x1D7D9}\left(f \in {P}\_{i}(\tau) \wedge f \in T\_{i}\right)}{{\sum\nolimits}\_{f} \unicode{x1D7D9}\left(f \in T\_{i} \right)}, \\ F\_{\max} &=& \max\_{\tau} \left\\{ \frac{2\cdot \text{pr}(\tau)\cdot \text{rc}(\tau)}{\text{pr}(\tau) + \text{rc}(\tau)} \right\\}, \end{array} $$

> where P i (τ) denotes the set of terms that have predicted scores greater than or equal to τ for a protein sequence i, T i denotes the corresponding ground-truth set of terms for that sequence, m(τ) is the number of sequences with at least one predicted score greater than or equal to τ, $\unicode{x1D7D9}\left (\cdot \right)$ is an indicator function, and n e is the number of targets used in a particular mode of evaluation.
> 2 votes

# Answer
From a quick reading of the paper it looks like they are predicting biological function of proteins. They frame this as a classification problem, where they try to classify a protein as pertaining to a set of terms (Abnormality of the lymph nodes, Abnormality of skin physiology, etc). It looks like they predict for each protein whether it falls into this category and create a ROC curve to plot the false positive and true positive rates of prediction.

In this case, the continuous variable would be your predictions (Abnormality of the lymph nodes: predicted 0.65 probability, Abnormality of skin physiology: predicted 0.25 probability). The categorical variable would then be the ground truth that you are trying to predict (Abnormality of the lymph nodes: known 1, Abnormality of skin physiology: known 0). Based on these two variables you can make your ROC curve.
> 2 votes

---
Tags: proteins
---


# Question
Title: Have DNA motifs 6-12bp long, trying to get conservation scores
I have about 200 short nucleotide motifs (6-12 bp in length) from the human genome, and I'm trying to see how conserved they are across vertebrates.

I was thinking that I'd need to make a bed file for each motif that lists all of its occurrences in the human genome. From there, I could map the beds to a bigwig files of PhastCons scores (essentially doing the reverse of what the PhastCons software was designed to do). Does that sound like the best approach?

I'm getting stuck at the step of going from motifs to bed files. I've tried using BLAST to find all occurrences of motifs, but their short length is causing issues.
I've tried messing with the e-value threshold, word size, and filter parameters, but I still don't get any hits.

Is there a work-around for this issue, or should I just rethink my entire approach?

# Answer
> 7 votes
**In case you have only ACGT in your motifs**

The short motifs make it sound as if you are in the business of looking for a kmer counter. You can choose to use existing software or build your own.

1. Using existing software might be your easiest path. An older post from 2014 will probably give you a first idea what's out there: http://homolog.us/blogs/blog/2014/04/07/kmer-counting-a-2014-recap/ . Note that a couple of algorithms mentioned there got successors, so it is worthwhile digging a bit around. The small kmer size will make most of them usable for your needs.
2. As the maximum size of your kmers is comparatively small (12 nt need 24 bits, i.e., max 16.7 million entries in your kmer table), you should be able to easily roll your own kmer counting in about any language you like and on about any of nowadays computer. The pseudocode section on the Wikipedia entry for kmers will give you first pointers for that. Might be a bit more work, but maybe more flexible depending on your needs.

**In case you have IUPAC bases (N, W, etc.) in your motifs**

I don't know any any pre-existing software doing what you need. I could imagine that the short motifs make using regular expressions doable for this kind of search, but I may be wrong. Testing this should be easy in a simple script as all major programming languages have modules or libraries for REs. Even if it should take a couple of hours to run on your data set, that would be good enough for a one-off calculation.

# Answer
> 3 votes
To scan motifs in a genome (or database) I would use FIMO which will give you the exact locations of these motifs in your genome.

Once you have the locations, you can use a `phastCons` bigiwig from UCSC to calculate the basewise conservation scores. However, please remember that `phastCons` scores are smoothed across windows and it might not be the best metric if you are trying to compare the conservation levels at your motif matching sites as compared to the sequences flanking them.

I wrote a package a while back to do this, including doing de-novo motif discovery. However, it might be an overkill for your use case.

# Answer
> 1 votes
Where transcription factor binding sites are concerned, in our papers we have used per-base phyloP conservation data instead of smoothed phastCons scores.

We use BEDOPS `bedmap` to map scores over multiple, padded binding sites (BED-formatted) for a given motif model. The phyloP scores are WIG files obtained from UCSC goldenpath and converted to BED via `wig2bed`.

The resulting matrix of binding site intervals and their per-base scores can be turned into a ranked heatmap or aggregated to determine average per-base conservation for a motif model.

While TF binding sites innately show high information content and therefore high conservation, matrices can be further sorted by score-maps of ChIP-seq or DNaseI-seq tag density.

This can help filter for interesting low-noise, high-signal patterns and is useful when the matrix is windowed and the window includes low-information regions with a few high-information residues outside of the main motif — as is the case with CTCF, for instance.

---
Tags: k-mer, motifs
---


# Question
Title: Stranded vs. unstranded library preparation protocols in RNAseq
I've been reading this paper lately:

Sailfish Enables Alignment-Free Isoform Quantification from RNA-Seq Reads Using Lightweight Algorithms

I don't really understand the second paragraph under Quantification in the Online Methods section. I quote the author:

> In an unstranded protocol, sequencing reads, and hence the k-mers they contain, may originate from transcripts in either the forward or reverse direction.

As far as I understand we are sequencing cDNA here and it is a single stranded molecule so all the reads will have the same direction (unlike the DNA where the reverse complement strand is also there), so I don't see the author's point. Can someone clarify please? Thanks.

# Answer
Short answer: cDNA is double-stranded.

I'll try to explain the concepts surrounding this with a nanopore read, where the cDNA component of the read (from the mouse genome) is indicated in a light grey, and the adapter components are either coloured or dark grey:

The cDNA library from this RNA was prepared by adding a polyT adapter (binding to the polyA sequence of the RNA), and carrying out reverse transcription to the start of the transcript. A strand-switching process was then used to prime and extend along the opposite strand of the DNA, creating a **double-stranded** full-length cDNA template. This template was then amplified using PCR to increase cDNA concentration and add in barcoding sequences.

You may notice that this particular read has a polyT head, rather than a polyA tail, attached to the `VNP` primer in the forward orientation. I can tell by this feature that this read comes from the reverse-complement strand of a DNA transcript, and so I know that the actual RNA transcript sequence is the reverse-complement of this read. Both the forward and reverse-complement sequences could potentially come through in the sequencer, because the adapters for sequencing are attached to both sides of the cDNA template, but the transcript direction can still be determined by looking at the orientation of the adapter sequences. When mapping this sequence to the genome, only a single strand (and its associated set of kmers) needs to be considered.

When considering what happens after mapping to the genome, there are four possible match situations:

1. Adapter is in reverse orientation, genome is in reverse orientation
2. Adapter is in forward orientation, genome is in reverse orientation
3. Adapter is in reverse orientation, genome is in forward orientation
4. Adapter is in forward orientation, genome is in forward orientation

In the case of the library preparation above, situations 2 and 3 correspond to a transcript that is encoded on the forward strand of the genome, whereas situations 1 and 4 correspond to a reverse-encoded transcript. Different library preparation methods could have a different way of encoding transcript direction. See here for more details on how these match situations can be used to produce strand-specific BAM files.

If I were to use an unstranded protocol (e.g. generating a cDNA library through other means, and attaching adapters randomly to the forward or reverse transcript strand), then the adapter direction wouldn't give me information about which strand this sequence came from. I would need to consider both strands when mapping (or more correctly, when identifying the native direction of the originating transcript), because the adapter sequence has been decoupled from the transcript direction.

On Illumina systems, this strand direction calculation is performed via the chemistry on the flow cell, and the adapters are \[usually\] not sequenced by the sequencer at all. Adapters from one direction bind to complementary sequences on the flow cell and are extended to produce the first read, and adapters from the other direction bind to produce the second read. Even with a stranded protocol, adapters from both directions are bound (because the cDNA template is double-stranded), but the binding order (as first or second read) can be used to determine the native direction of the originating transcript.
> 8 votes

# Answer
There are a number of library preparation methods out there. Some, like the popular dUTP-based methods, result in libraries where the orientation of a read's alignment confers information regarding the strand of the transcript from whence it came. Other protocols do more PCR cycles and don't use dUTP to allow distinguishing strands. These are then termed "unstranded" libraries and one can no longer determine the strand/direction of the originating transcript from the orientation of a read's alignment. In general, a transcript will then have reads pointing in both directions in more or less equal proportions.

You tend to not see unstranded libraries much any more. Early on they were significantly more expensive, but that's not really the case any more. There also used to be two types of stranded protocols, whereas you mostly see dUTP-based protocols these days.

BTW, you'll see these same settings in most RNAseq programs. The equivalent to an unstranded protocol setting in featureCounts is `-s 0`, in tophat2 it was `--library-type fr-unstranded` (no one should use tophat2 any more, it's just an example) and so on in other programs.
> 5 votes

---
Tags: rna-seq, k-mer
---


# Question
Title: How can I use annotations to remove variants not relevant to cancer risk?
I currently have ~180 whole germlines and around 10M SNPs/indels. I would like to build a predictive model using Machine Learning (ML) techniques to predict cancer risk according to these germline variants. The thing is, most of these 10M variants are not relevant to cancer, so first I have to annotate them and remove non-relevant SNPs/indels.

Which tools/measurements do you think are important to use in order to filter out irrelevant variants and keep only, let's say, around 10,000 SNPs/indels?

I guess the most obvious is CADD score, so I could keep those SNPs/indels with PHRED \>= 10. Which one would you suggest taking into account that downstream analysis would be ML algorithms that take as features variants that **really** contribute to cancer?

# Answer
> 5 votes
While your question is specific to cancerous germline mutations, I'd suggest you look at the COSMIC database of somatic mutations to include in your analysis.

There are other factors to include in this kind of analysis you're suggesting, such as predictive deleterious effects (PolyPhen for example can perform such predictions).

If you have 10M variants/indels, for the same form of cancer, then look for common variants, or maybe the frequency of variants identified across exons of a gene.

What you've asked is highly broad and we'd need more specific details about your cohort to suggest more specific approaches.

# Answer
> 3 votes
I don't think it will be possible to do what you ask, right now with current knowledge. Selecting variants relevant to cancer risk is still an open problem and usually requires quite a lot of human intervention.

You can use different measures of population frequency to filter common variants with the assumption that frequent variants won't be pathogenic.

You can use multiple in-silico predictors of pathogenicity to try to remove "benign" or "tolerated" variants, but predictors are far from perfect and you have to take into account that many of them are essentially measuring the same, so when combining their outputs you may be biasing your analysis towards nucleotide conservation, etc...

You can use in-silico predictors of splicing effects. These are even farther from perfect and return very noisy data.

You can take into account where the variants lie in the gene. In many genes, variants in the last few nucleotides won't be pathogenic, but may be important in others. In some genes only variants in some domains have an effect, while not in others.

You have to take into account if only one or both alleles of certain genes are affected. For that you may need phasing information. And obviously the list of genes.

Some genes have an impact on cancer risk when they are "deactivated" (or broken somehow) (tumor supressors) while others have an effect when their function is improved (oncogenes).

And for many many genes you'll simply not know if they have any effect in cancer risk at all.

You can get data from different databases and annotate your genomes to remove many variants but in the end you won't have very high quality data to enter to your ML algorithms. Maybe this will work for your project, but take into account that this "first step" of data cleaning is not easy or obvious with current biological knowledge.

---
Tags: ngs, sequence-annotation, machine-learning, cancer
---


# Question
Title: How to get the count of each kmer past 255 using khmer
I have a Fastq file and I want to get the exact count of each possible kmer from this file.

On a previous post called How to use khmer to count k-mers? Daniel Standage proposed a custom script based on khmer methods that you can see down below :

```
>>> outfile = open('outfile.txt', 'w')
>>> seenkmers = set()  # Consumes a lot of memory for large input!!!
>>> for read in khmer.ReadParser('reads.fq.gz'):
...     for kmer in counts.get_kmers(read.sequence):
...         if kmer not in seenkmers:
...             print(kmer, counttable.get(kmer), file=outfile)
...             seenkmers.add(kmer)

```

The code works if I switch `counttable.get(kmer)` with `counts.get(kmer)` but the main problem is that when I try on a customized fastq file with 300 occurrences of one kmer, the count will still be stuck at 255 and can't go further, and for some reason I can't pass through this limit.

Is this coming from the minimum table size ?

Does anyone know if there is a parameter that can work on a custom script to disable this limit ?

# Answer
255 is the default maximum size of a `Counttable` in khmer. You want to do the following:

```
import khmer
counts = khmer.Counttable(31, 1e7, 1)  # Feel free to change the parameters
counts.set_use_bigcount(True)

```

The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (at some point you have to increase the size of the number used in the hash table).
> 6 votes

---
Tags: python, fastq, k-mer, khmer
---