text:
Pipeline (Pypiper) returns error, but individual command works fine
I am running a script file on a cluster with slurm through sbatch. Locally it runs perfectly well on 2 different linux versions (linux mint, ubuntu 12.04), however, it returns -9 error code while running on a cluster which is CentOS 7.3. When we are using Pypiper and run the whole .py script file, it fails when it starts running trimmomatic program
http://www.usadellab.org/cms/?page=trimmomatic
However, when running just trimmomatic program step by itself it does not fail and finishes successfully. We have a suspicion that something may be wrong with how Pypiper runs on the cluster, however, feeling rather stuck.
Any suggestions would be greatly appreciated.
Here is the error that is printed out:
ILLUMINACLIP: Using 2 prefix pairs, 22 forward/reverse sequences, 0
forward only sequences, 0 reverse only sequences
</pre>
Process 109391 returned: (-9). Elapsed: 0:00:51. Peak memory: (Process:
5.211GB; Pipeline: 0.0GB)

### Pipeline failed at:  (12-14 15:50:47) elapsed:0:23:09 _TIME_

Total time: 0:23:09

Change status from running to failed
Traceback (most recent call last):
  File "...../salmon_rna_seq.py", line 344, in <module>
sys.exit(main())
File "...../salmon_rna_seq.py", line 304, in main
pipeline.process()
File "...../salmon_rna_seq.py", line 277, in process
self.trim_reads()
File "...../salmon_rna_seq.py", line 195, in trim_reads
self.pipe_manager.run(cmd, self.sample_files['r1']['fastq']['trimmed_paired'])
File "...../pypiper.py", line 508, in run
process_return_code, local_maxmem = self.callprint(cmd, shell, nofail, container)  # Run command
File "...../pypiper.py", line 691, in callprint
self.fail_pipeline(e)
File "...../pypiper.py", line 1111, in fail_pipeline
raise e
OSError: Subprocess returned nonzero result.

Pypiper terminating spawned child process 108475...
child process terminated
slurmstepd: error: Exceeded step memory limit at some point.

The last line may be informative. Although, we are still not sure what we could possibly do.

Slurm is killing the process because it's using more memory per core than the default allocation (that's what slurmstepd: error: Exceeded step memory limit at some point means). I've never used pyPiper (like seemingly most people, I use snakemake), but presuming that you have it running either srun or sbatch at some point then specify something like --mem-per-cpu 10G.

Actually I tried that and it still does not work. Will try to write to the cluster support

You might need a larger value than 10 gigs. You can use scontrol to see the exact threshold getting crossed that's causing slurm to kill the job.

Awesome. We increased it to 20Gb and it run successfully.


text:
factoextra: Error in if (xlab == "x") xlab <- "x value" : argument is of length zero
I just followed the tutorial trying to do clustering on my own data, but I can't visualize it because fviz_cluster function is giving me an error:

Error in if (xlab == "x") xlab <- "x value" : argument is of length zero

Here is the code I am running:
library("fpc")
db <- fpc::dbscan(rtsne$Y, eps = 0.15, MinPts = 5)
# Plot DBSCAN results
library("factoextra")
fviz_cluster(db, data = rtsne$Y, stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
         geom = "point",palette = "jco", ggtheme = theme_classic())

How can I make fviz_cluster function working?

Update

Filed a bug: https://github.com/kassambara/factoextra/issues/62

I actually found an answer just accidentally. It is very unfortunate that factoextra documentation does not explicitly say that data parameter should be of data.frame type. In my case rtsne$Y was a matrix. The error that factoextra is giving is very uninformative.
Here is the code that ultimately worked:
fviz_cluster(db,
----->  data = as.data.frame(rtsne$Y),
        stand = FALSE,
        ellipse = FALSE, show.clust.cent = FALSE,
        geom = "point",palette = "jco", ggtheme = theme_classic())

I am using fviz_cluster as well, but for PCA clustering. I have a separate question that I've asked but haven't gotten the perfect answer, when you put geom = "label", where does it take the labels from? My guess would be the rownames, but the rownames in my data.frame are non unique. What do I do?

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Transform and feed data into DESeq2 with DESeqDataSetFromMatrix
There is a normalized expression matrix. I split it into two and want to do DE on the two cells' subsets. I am having trouble transforming it into the format that DESeq2 would accept. The function that I would think I need to use is the following:
dds <- DESeqDataSetFromMatrix(countData = cts,
                          colData = coldata,
                          design= ~ batch + condition)

It would be perfect if I could somehow feed two matrixes separately, but I can not find it. Any suggestions would be greatly appreciated.

Update

I tried the following:
cells <- c(cell_names, cell_names_cyt)
groups <- c(rep("NE", times=length(cell_names)), rep("EX",
times=length(cell_names_cyt)))

colData <- data.frame(
  Cells <- cells,
  Groups <- groups
)

dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

But getting the error:

Error in DESeqDataSet(se, design = design, ignoreRank) :
    some values in assay are not integers

Of course, they are floats since I already normalized the matrices.

Update

I did the following since I am using Seurat:
df <- as.data.frame(as.matrix(seurat_object@raw.data))
df <- df[,colnames(df) %in% cells]

But then the following:
dds <- DESeqDataSetFromMatrix(countData = df,
                          colData = colData,
                          design= ~ Groups)

Gave me an error:

Error in DESeqDataSet(se, design = design, ignoreRank) : all variables in design formula must be columns in colData

Don't input normalized counts, use raw counts.

Yeah, I used raw counts retrieved from Seurat, but then got another error.

Error in DESeqDataSet(se, design = design, ignoreRank) : all variables
  in design formula must be columns in colData

I think that means just what it says...you don't have a "Groups" column in coldata.  Did you spotcheck your ColData data frame to make sure it looks the way it should?
> cells <- c("1","2","3")
> groups <- c("a", "b", "c")
> test <- data.frame(Cells <- cells, Groups <- groups)
> test
  Cells....cells Groups....groups
1              1                a
2              2                b
3              3                c
> colnames(test)
[1] "Cells....cells"   "Groups....groups"

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
How can I use UpSet plot for expression/RNA-seq data
I have a data.frame of 27 000 x 17 and I want to feed it into UpSetR program, but getting the error:

Error in start_col:end_col : argument of length 0

My data.frame looks like that:

And the upSetR code that I am using:
upset(df, nsets = 4, number.angles = 30, point.size = 3.5, line.size = 2,
mainbar.y.label = "Gene Intersections", sets.x.label = "Genes Per Group",
text.scale = c(1.3, 1.3, 1, 1, 2, 0.75))

What am I doing wrong here? Any suggestions would be greatly appreciated.

Update

I found here a similar issue:
https://plus.google.com/109635504520741246387/posts/Etujrp6sU1Q
They say that the input should be binary. But then, how can I calculate this binary matrix from the rlog-data? I am not sure...

UpSet plots are in many ways similar to Venn diagrams, in that the represent some sort of metric for group overlap. What I expect you want to do is get the number of DE genes in various comparisons and put those numbers into an UpSet plot.

upsetR takes as input a data.table in the long format. I found this example/code very helpful when I first started using upsetR.

I need to get the binary matrix to be able to feed it into the UpSetR. Right now I have a matrix where rows are genes and columns are groups. Somehow for each gene I need to set a threshold value to convert it to binary ‘expressed/not expressed’ probably. Groups are vectors of numbers, so I need to find their intersection somehow. I do not know which software or R package do that: I would not want to do it manually.

I am pretty confused but the problem is rather definite: I need to represent groups, their similarity, as an UpSet plot.

I think, it would be good to rephrase your question to match your ultimate goal something like "How can I use UpSet plot for expression/RNA-seq?". Say what is you sample and why you want to do the UpSet plot.

You need to generate a list of gene ID vectors to feed into upset. In the example below, I am selecting observations (genes) which have an level above the median level. You could also choose only genes above the 3rd quartile, etc. You then need to convert these from a list of character vectors to a presence/absence matrix. UpSetR has a fromList function that does just this.
In future, be more specific about what you want to achieve and supply example data for people to work with.
library("UpSetR")

tmp <- mtcars

tmp$car <- rownames(mtcars)
list <- lapply(colnames(mtcars), function(x) {
  tmp[tmp[[x]] > median(tmp[[x]]), "car"]
})
names(list) <- colnames(mtcars)

upset(fromList(list))

Created on 2019-10-17 by the reprex package (v0.3.0)

I wrote a script to create the binary data frame with my differential expression data. If it is of any help to anyone:

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
RCurl issue when installing SCnorm
I am working on the Rivanna cluster that is of SLURM-type. Because of that I do not have write permissions and need to install the newest R version myself locally, and use it. For it to be successfully built from source, I needed to also build from source libcurl4-openssl-dev, so I downloaded tar.gz file from here:
https://packages.debian.org/ru/jessie/libcurl4-openssl-dev
and did:
mkdir build
./configure --prefix=/scratch/nv4e/curl/build
make
make install

Then, in the same way I built from source R 3.5.1, the latest current version. It was compiled without any errors, but with some warnings. Then, I submitted the R script that uses SCnorm package:
https://bioconductor.org/packages/release/bioc/html/SCnorm.html
Inside of the script I am having:
source("https://bioconductor.org/biocLite.R")
biocLite("SCnorm")
biocLite("SingleCellExperiment")
library(SCnorm)
library(SingleCellExperiment)

Of course, I am submitting it through .slurm file where I have the main command:
Rscript --vanilla scnorm.R

It fails, apparently because SCnorm requires RCurl, but it can not get one, even though in the log file it is trying to access it correctly:

The libcurl library was built from source in the /scratch/nv4e/curl/build folder. The only suspicion that I have here is that maybe its version is wrong and SCnorm needs another libcurl version? Any suggestions would be greatly appreciated.

Update

I tried reinstalling libcurl library from the official website, using the latest version:
https://curl.haxx.se/
But it did not fix the issue. I have also tried using libcurl from Anaconda but it gave the same error. So, neither 7.61.0 (latest), nor 7.29.0 (Anaconda), nor 7.38.0 (Debian link) are working.

This is a bug affecting the curl compile, it's been reported here and also here, I suggest trying to compile from the newer release of the libcurl source rather than the .deb you linked.

I tried reinstalling libcurl from here: https://curl.haxx.se/, and it did not help. I see that the R is using the newly reinstalled version, and giving the same error.

Are you sure it's not hitting the old system version at compile time?

Yes, sure. Cluster support fixed it. The issue was related to gcc version, but I am not completely sure what they did.

Glad it got fixed, although if you find out exactly it would be good to post and answer since this will help others with the same issue.


text:
PCA plot shows big difference but not many differentially expressed genes are found
I got a PCA plot of bulk RNA-seq experiment that looks the following way:

It was generated by the following code:
pcaData <- plotPCA(rld_sva, intgroup=c("Group"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed()

First sva correction was run to correct for batch effects and then rlog transformed values were plugged in to plotPCA function.
The first issue that catches the eye is that 100% of the variance is explained by just 2 dimensions. I am not sure what can one say about the data in this case. The second issue is that I get only around 20 differentially expressed genes by using DESeq2 analysis (log2FoldChange > 1, p_adj < 0.05). I know that we can not directly state that if there is a large difference on PCA there will be present a plenty of differentially expressed genes, but why is it not the case? Simple logic tells me that pca shows the difference between the samples in their gene expression, so I would expect seeing a plenty of differentially expressed genes.

Here each point is a sample or a cell? As you mention explaining 100% of variance in two dimensions would mean that an error happened somewhere. What have you done to find the root of this problem? And what is your question? If it is possible, if it is normal to happen something like this or something else?

Sample. The question is why so much difference in pca does not lead to many differentially expressed genes. Second, minor question, is why the data is 2D, possible reasons for that. Sorry, I made a mistake, it is bulk RNA seq

Your two IgG replicates are separated a lot in the PCA plot, compared with the two other  replicates, which likely means that there is too much disagreement within the IgG group in order for DESeq to arrive at low enough P-values for most of the differentially expressed genes. Adding one more IgG replicate would probably help a lot.

You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points.
On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.


text:
Running python script from terminal: umap issue
Some time ago I asked the question:
RunUMAP in Seurat not working: module 'umap' has no attribute 'UMAP'
It was solved by following the answer in it, however, it is still not working when I run a python script from terminal, giving the same error:
python python/umap.py -path
"/Users/nikitavlasenko/python_scripts/output/" -data
"/Users/nikitavlasenko/python_scripts/data_files/WT_SCID_CCA.csv"

Traceback (most recent call last):
    File "python/umap.py", line 26, in
      result = umap.UMAP().fit_transform(data.T)
  AttributeError: module 'umap' has no attribute 'UMAP'

From inside spyder the same script runs perfectly, but somehow when running it from the terminal it is not successful. I am new to python, so maybe I am not pointing correctly to umap library? Maybe it should be imported in a different way? Here is the whole script:
import umap
import numpy as np

def getopts(argv):
opts = {}  # Empty dictionary to store key-value pairs.
while argv:  # While there are arguments left to parse...
    if argv[0][0] == '-':  # Found a "-name value" pair.
        opts[argv[0]] = argv[1]  # Add key and value to the dictionary.
    argv = argv[1:]
return opts

if __name__ == '__main__':
  from sys import argv
  myargs = getopts(argv)
  print(myargs)
  path_to_save = myargs['-path']
  path_to_data = myargs['-data']
  data = np.loadtxt(open(path_to_data, "rb"), delimiter=",", skiprows=1)
  result = umap.UMAP().fit_transform(data.T)
  np.savetxt(path_to_save + "/umap.csv", result, delimiter=",")

Python version that is used by spyder and in the terminal seems to be the same: 3.6. I am using Mac OSX.

Add a print(umap.__file__) early on. It's likely that you're getting a different copy of umap.

You answered your own question, yet it is not a coincidence that you have to change the name of your script.
When Python executes your custom script, it reads:

import umap

Python will first search the local sources (the directory of the script) , and afterwards the installed sources.
In this case Python identifies the local source to be umap.py, and imports the (currently executed) script, while foregoing the installed sources.
Since you are calling functions, not defined in your current script, Python will not be able to fully execute your script.
For more information on importing see Module search path at docs.pyhton.org

It is pretty weird but the only thing that I needed to do is to rename my python script from umap.py to umap_gen.py and right away it started to work. So, I guess there can be such name collisions which we should avoid. Although I would prefer the file to be named as umap.py still...


text:
Seurat installation fails
When I am trying to install Seurat in R markdown file:
install.packages("Seurat")

it fails with the following errors:
Installing package into ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3’
(as ‘lib’ is unspecified)
Warning in install.packages :
  dependency ‘proxy’ is not available
also installing the dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’

trying URL 'https://cran.rstudio.com/src/contrib/igraph_1.2.1.tar.gz'

Content type 'unknown' length 2703426 bytes (2.6 MB)

downloaded 2.6 MB

trying URL 'https://cran.rstudio.com/src/contrib/dtw_1.18-1.tar.gz'

Content type 'unknown' length 849256 bytes (829 KB) downloaded 829 KB trying URL 'https://cran.rstudio.com/src/contrib/diffusionMap_1.1-0.tar.gz'

Content type 'unknown' length 39130 bytes (38 KB)
downloaded 38 KB

trying URL 'https://cran.rstudio.com/src/contrib/Seurat_2.3.0.tar.gz'

Content type 'unknown' length 846065 bytes (826 KB)

downloaded 826 KB
* installing *source* package ‘igraph’ ...
** package ‘igraph’ successfully unpacked and MD5 sums checked
installing to /sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs
** R
** demo
** inst
** preparing package for lazy loading
** help
*** installing help indices
** building package indices
** testing if installed package can be loaded

Error in dyn.load(file, DLLpath = DLLpath, ...) :    unable to load shared object '/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph/libs/igraph.so': libicui18n.so.58: cannot open shared object file: No such file or directory

Error: loading failed Execution halted ERROR: loading failed
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/igraph’ Warning in install.packages :   installation of package ‘igraph’ had non-zero exit status

ERROR: dependency ‘proxy’ is not available for package ‘dtw’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/dtw’ Warning in install.packages :   installation of package ‘dtw’ had non-zero exit status

ERROR: dependency ‘igraph’ is not available for package ‘diffusionMap’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/diffusionMap’ Warning in install.packages :   installation of package ‘diffusionMap’ had non-zero exit status

ERROR: dependencies ‘igraph’, ‘dtw’, ‘diffusionMap’ are not available for package ‘Seurat’
* removing ‘/sfs/qumulo/qhome/nv4e/R/x86_64-pc-linux-gnu-library/3.3/Seurat’ Warning in install.packages :   installation of package ‘Seurat’ had non-zero exit status

How could it be fixed? What to try?
Operating system is CentOS Linux release 7.3.1611 (Core)

Update

I also tried development version of Seurat but it failed too:
http://satijalab.org/seurat/devinstall.html

Update

I tried installing Seurat using conda install -c bioconda r-seurat and it is not working. It is giving the error:

-bash-4.2$ conda install -c bioconda r-seurat
  Fetching package metadata .............
  Solving package specifications:
PackageNotFoundError: Packages missing in current channels:

r-seurat -> r-cowplot
r-seurat -> r-diffusionmap
r-seurat -> r-dtw
r-seurat -> r-ggridges
r-seurat -> r-ica
r-seurat -> r-pbapply
r-seurat -> r-ranger
r-seurat -> r-rcppprogress
r-seurat -> r-sdmtools
r-seurat -> r-tclust

We have searched for the packages in the following channels:

https://conda.anaconda.org/bioconda/linux-64
https://conda.anaconda.org/bioconda/noarch
https://repo.continuum.io/pkgs/main/linux-64
https://repo.continuum.io/pkgs/main/noarch
https://repo.continuum.io/pkgs/free/linux-64
https://repo.continuum.io/pkgs/free/noarch
https://repo.continuum.io/pkgs/r/linux-64
https://repo.continuum.io/pkgs/r/noarch
https://repo.continuum.io/pkgs/pro/linux-64
https://repo.continuum.io/pkgs/pro/noarch

You seem to be missing the libicui18n.so library. Please [edit] your question and tell us what operating system you are using. If Linux, tell us what distribution.

I just updated. It is CentOS. I am working with rstudio on the cluster.

Ah, thanks. So what is the output of yum search libicui18n (if that returns nothing, try yum search libicui or even yum search icui until you get a result). Then, install the package. If that works, let me know and I can post answer (or just post it yourself, up to you).

The output in all the three cases is No matches found. Also, I do not have root privileges obviously, so maybe I should download this libicui18n.so library and install it locally somehow?

Hmm, I haven't used an rpm-based system in ages. You might need to search for libicui* or something. I did check CentOS's repository, and there is such a package. I don't remember how to install rpms locally though and I'm on mobile, so I can't help much. You could try asking how to find and install the relevant package on [unix.se].

As @terdon mentioned in the comments, you need to sudo yum install libicu on the relevant nodes (or have your admin do that for you). Alternatively, seurat can be installed via conda, which means you don't need root access.

Actually I tried installing Seurat using conda install -c bioconda r-seurat  and it is not working. It is giving PackageNotFoundError.

You need -c conda-forge too. You might then need to fix readline with conda upgrade -c conda-forge readline.

If you don't wish to use conda as your r package manager (because, for a start, it will install a new copy of R with its own library), you can install just ICU with conda conda install -c conda-forge icu, and then install seurat in R


text:
How can I apply proportional (p) distances (Nucleotide) using bioPython
I am working on phylogenetic tree generation. I used BioPython and ClustalW 2.x.x for this purpose. I have generate the tree using BioPython but when I try to generate tree using "MegaSoftware GUI" tool my tree does not matched with the output of MegaSoftware.
When I use "No. of difference" model then results are matched with biopython tree Like below image shown.

But whenever I choose "p-distance" model my results are not matched.

Output MegaSoftware using p-distance

Output BioPython

Here is my code for generating phylogenetic tree:
def conversion_fasta_dnd(clustalw, filepath, type, clustering, output, pwmatrix):
    """
        ClustalW tool to conversion the file format.
        :param clustalw: path of the executable programm.
        :param filepath: path of `fasta` extension file.
        :param type: dna, protein
        :param clustering: upgma, nj
        :param output: clustal, phylip, fasta
        :param pwmatrix: blosum, pam, id
        :return: dnd, aln (files)
        """
    # esblish the full path to the clustalW program
    clustalw_exe = clustalw

    # check to make sure the program is there
    assert os.path.isfile(clustalw_exe), "Clustal W executable missing"

    # create the appropriate command line
    clustalw_cline = ClustalwCommandline(clustalw_exe,
                                         infile=filepath,
                                         type=type,
                                         clustering=clustering,
                                         output=output,
                                         pwmatrix=pwmatrix)

    # print the command line
    print(clustalw_cline)

    # execute the command
    stdout, stderr = clustalw_cline()

# This function draw the tree on image file
def simple_tree(filepath):
    """
    Generate the simple tree image.
    :param filepath: path of dnd file format.
    """
    def get_label(leaf):
        return leaf.name

    file = filepath
    tree = Phylo.read(file, 'newick')
    Phylo.draw(tree, label_func=get_label, do_show=False, branch_labels=lambda c: c.branch_length)

    # axis management
    pylab.axis('on')
    pylab.savefig(OUTPUT_IMAGE, format='png', bbox_inches='tight', dpi=100, bbox_extra_artists=[])

In above code I have used two things first BioPhython lib and the other is Clustal W for align the sequences.
Could you please clarify what p-distances are?
I only have these reference for p-distance:

p-distance (Nucleotide).
Proportional (p) distances

I need to add p-distance model in my code, someone please help me out How can I do that.

Guide me on how to write this P-distance code in a matrix form.
I want to add P-distance matrix.

I think you should make a new question with a reproducible example.

Based on the sites you mention, it seems that p-distances are pairwise distances between aligned sequences obtained by counting the number of differences and dividing by the length of the alignment. This is the proportion of differences between two sequences.
I'm not sure exactly what you mean by "add p-distance model in my code", but here is how you could compute those distances from an existing alignment (here I use as an example an alignment I have in fasta format):
from Bio import AlignIO

def p_dist(seq1, seq2):
    ali_len = len(seq1)
    assert ali_len == len(seq2)
    return sum(0 if nuc1 == nuc2 else 1 for (nuc1, nuc2) in zip(seq1, seq2)) / ali_len

ali = AlignIO.read("most_targeted_histones/transcript_alignments/H4_top4_loss_spliced_NT.fa", "fasta")
nb_seq = len(ali)
distances = {}
for i in range(nb_seq - 1):
    rec1 = ali[i]
    j = i + 1
    while j < nb_seq:
        rec2 = ali[j]
        distances[(rec1.id, rec2.id)] = p_dist(rec1.seq, rec2.seq)
        j += 1

print(*distances.items(), sep="\n")

Output:
(('ZK131.4_spliced', 'ZK131.8_spliced'), 0.0)
(('ZK131.4_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.4_spliced', 'ZK131.1_spliced'), 0.0)
(('ZK131.8_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.8_spliced', 'ZK131.1_spliced'), 0.0)
(('T10C6.14_spliced', 'ZK131.1_spliced'), 0.10897435897435898)

If you want to build a tree from those distances, I'm not sure you can do this with BioPython. It seems however that there are some possibilities with DendroPy.
You would need to write the distances in a suitable format, and to choose the method between UPGMA and neighbour-joining.


text:
Low fraction reads in cell metric in snRNA-seq data
For some of our snRNA-seq samples we are finding low fraction of reads detected in single-nuclei rna-seq samples from cellranger, while the other metrics are perfect.
While I understand this could be caused by library preparation and abundance of ambient RNA. Just wondering if this behavior is expected for single-nuclei datasets. My confusion is specifically coming from the fact that 10X website also has such a sample that has low number of fraction reads in cells.
https://cf.10xgenomics.com/samples/cell-exp/7.0.0/5k_mouse_brain_CNIK_3pv3/5k_mouse_brain_CNIK_3pv3_web_summary.html
Would love to know if this behavior is expected or seen often who are more experienced in snRNA-seq. The metric documentation page is here
Here is snapshot from 10X report

What do you mean by low fraction? Fraction of what (i.e. what's the denominator)? What were you expecting? What did you see instead?

Great questions, sorry for being vague, uploaded the screen-shot for cell-ranger. low fraction - 39%, expected according to 10x 70% (if you follow the link). I am not sure how they are computing the metric though, I added the man page from 10X

Yes, "low fraction reads in cells" is more common with snRNA-seq from frozen samples, I wouldn't worry about it by itself (you're just throwing away some sequencing data). You should still check though that the other QC metrics are good. In particular, sometimes such samples have an unclear seperation of cells from empty droplets on the UMI rank plot and the consequent estimated cell count is off, requiring alternative approaches to distinguish cellular droplets from ambient RNA etc.

Thanks a lot Chris. I am thinking about more sophisticated pipeline like CellBender etc. I just wanted opinions for professional more experienced than me because I have not dealt with non-standard snRNA-seq data a lot.


text:
What data set was used Gene Expression Data of Postmortem Tissues?
I want to run the experiments "Sample Application to Genotype-Tissue Expression (GTEx) Project Gene Expression Data of Postmortem Tissues" mentioned in "Enter the Matrix: Factorization Uncovers Knowledge from Omics" by Stein-O'Brien et.al link
The paper metioned

dataset from postmortem samples in the GTEx project [42] (Figure 3).
Specifically, we select a subset of GTEx data containing 12 brain
tissues in the GTEx data for seven individuals for which we had
previously performed only NMF analysis [41] (codes are provided in the
supplemental information online).

I check with https://gtexportal.org/home/datasets and assumed that the paper used the GTEx Analysis V6.
But which files do I need to download? I am computer scientist I am unsure which ones contain the brain  tissues
Further I found
Brain_Anterior_cingulate_cortex_BA24_Analysis_cis-eQTLs.txt.gz  1.5G
Brain_Caudate_basal_ganglia_Analysis_cis-eQTLs.txt.gz   1.9G
Brain_Cerebellar_Hemisphere_Analysis_cis-eQTLs.txt.gz   1.7G
Brain_Cerebellum_Analysis_cis-eQTLs.txt.gz  1.9G
Brain_Cortex_Analysis_cis-eQTLs.txt.gz  1.9G

But I am unsure which one was used in the paper


text:
LD analysis in PLINK based on reference and a SNP list
I'm new to PLINK and genetics, and getting confused with two PLINK commands for LD analysis:

plink --bfile hapmap --r2 --ld-window-r2 --ld-snp-list --ld-window
plink --noweb --bfile --clump hapmap --clump inputfile --clump-field P --clump-p1 --clump-p2 --clump-r2 --clump-kb

It seems to me that for #1, the input files are a reference data and a SNP list file, while for #2, it is a reference data and a .assoc file.
In the output for #2, p-values are included with other results (according to PLINK manual) which is not the case with the output from #1. If I'm interested in finding proxy SNPs for the SNPS on my list in hapmap data, will it be sufficient for me to use the output from #1?
After invoking command #1 using PLINK 1.9, its output gave me 6 columns: "Chromosome number", "chromosomal position", "SNP name" for query/target SNP and 1 column of "R2". I cannot carry out command #1 because I do not have the .assoc file and only have a list of SNPs.
My main question is - what is the difference between these two commands and which of the two should give me the best result for finding proxy SNPs for my list of SNPs.

It looks like the first command is calculating pairwise linkage disequilibrium r2 between your SNPs.
The second command is clumping your dataset based on a p-value over a file of association test results, which I guess is a pre-cursor to computing polygenic risk scores.
Clumping is the process of thinning your SNP list to only keep uncorrelated SNPs so that each representative region has only one marker. In this case you are asking PLINK to keep the marker with the highest p-value.
I'm not sure what your definition of proxy SNPs is but if you are looking for tag SNPs to a SNP of your interest, you should go with command 1 and choose a proper LD r2 threshold.

Many thanks @Vivek for your explanation. I really appreciate it!


text:
How to get all PDB homologs from Uniprot (mapping + BLAST)?
I'd like to create a dataset consisting of all sequences which are either present in the PDB, or whose homolog is present in the PDB. In other words, any sequence in the PDB or any sequence related to it. The similiarity margin is to be very wide, so anything above, for instance, 25% sequence identity is accepted (perhaps there are better criteria, but you get the point).
No further limitations, so any organism, with or without experimental data, not limited to the manually curated Swissprot set.
I have tried the SIFTS database, but that contains, as far as I can tell, only mappings for sequence identity >85%. I would like to go far below that.

Thanks for the edit, that's much clearer! So, next question: how do you define "homologs"? Are you looking for functional homology alone? Is sequence homology enough? I am guessing you will be assuming that homologous sequences will share a structure but that is most certainly not always the case. Is there any reason you don't do it the other way around? Get all sequences from the PDB and then map those to UniProt accessions?

@terdon Homologs are defined here as >25% sequence identity, so very very broadly. I would, in fact, be willing to run a BLAST search for each PDB entry vs Uniprot, but I assume this will take quite some time, so I'm checking to avoid reinventing the wheel.

For the correspondence between PDB and Uniprot entries you can use SIFTS -- a semi-automated mapping between PDB and UniProt maintained by PDBe.
The pipeline that creates the mappings uses BLAST and a few other criteria to decide which UniProt entry should be assigned to each PDB entry. The SIFTS website has all the data in CSV files.

Thanks for the advice! I took a look. In the referenced publication, I only found a mention that they are looking for a 90% sequence identity to accept a cross-reference; this is not enough for me - I need all UniProt entries with, say, sequence identity above 25%.
As I understand it, this is what is contained in the files on the ftp server ftp://ftp.ebi.ac.uk/pub/databases/msd/sifts/ , notably the uniprot_segments_observed.

You can download a BLAST db containing all sequences of proteins in Uniprot and in the PDB. The way I would go about this is first download the databases for uniprot and PDB, then query the PDB database for each sequence from Uniprot. If you get a BLAST hit above a certain threshold (whatever you define as homolog) then add that sequence to a file, and voila, you have a full set of sequences that fulfill your criteria. The only tools you would need for this are the BLAST executables (found here) and the databases (info on how to download can be found here).
I should note that you don't really need any Python to get this to work, only some shell scripting. Although you can run these commands through a Python wrapper if you really want to.
Good luck!

So that's like 60M BLAST queries? How long would that take?


text:
Installing RnBeads via bioconductor - .onLoad failed in loadNamespace(), call: NULL
I am trying to install RnBeads from bioconductor but the installation of its dependancy TxDb.Hsapiens.UCSC.hg19.knownGene fails. The error is:
biocLite("RnBeads")
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) ‘RnBeads’
also installing the dependencies ‘FDb.InfiniumMethylation.hg19’, ‘methylumi’

trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/FDb.InfiniumMethylation.hg19_2.2.0.tar.gz'
Content type 'application/x-gzip' length 38996322 bytes (37.2 MB)
==================================================
downloaded 37.2 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/methylumi_2.24.1.tar.gz'
Content type 'application/x-gzip' length 6178909 bytes (5.9 MB)
==================================================
downloaded 5.9 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/RnBeads_1.10.2.tar.gz'
Content type 'application/x-gzip' length 6437974 bytes (6.1 MB)
==================================================
downloaded 6.1 MB

* installing *source* package ‘FDb.InfiniumMethylation.hg19’ ...
** R
** data
** demo
** inst
** preparing package for lazy loading
Error: package or namespace load failed for ‘TxDb.Hsapiens.UCSC.hg19.knownGene’:
 .onLoad failed in loadNamespace() for 'TxDb.Hsapiens.UCSC.hg19.knownGene', details:
  call: NULL
  error: evaluation nested too deeply: infinite recursion / options(expressions=)?
Error : package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ could not be loaded
ERROR: lazy loading failed for package ‘FDb.InfiniumMethylation.hg19’

The downloaded source packages are in
    ‘/tmp/RtmpY8Yf4l/downloaded_packages’
installation path not writeable, unable to update packages: doMC, doSNOW,
  foreach, iterators, snow, boot, class, cluster, codetools, foreign,
  KernSmooth, lattice, MASS, Matrix, mgcv, nlme, nnet, rpart, spatial, survival

As suggested by the error message, I tried setting options(expressions=500000), which unfortunately did not change the error message. An important note may be that I am under network mount and don't have administrative access rights.
Edit
I don't have a .Rprofile file.
Note about installation path not writeable:
Typing .libPaths() reveals that my home R folder is at the first position (and it is writable). The other paths ("/usr/local/lib/R/site-library"                          , "/usr/lib/R/site-library" and "/usr/lib/R/library") are not writable without sudo.
I also tried to delete my home R folder and to reinstall it. The reinstallation of all packages worked for all packages except for the three listed above (same error message). I also tried biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene", lib.loc=.libPaths()[1], lib=.libPaths()[1]) to specifically address my home R directory, but the same error is showing up (but without the installation path not writeable part).
Here is my sessionInfo():
R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

Matrix products: default
BLAS: /usr/lib/openblas-base/libopenblas.so.0
LAPACK: /usr/lib/lapack/liblapack.so.3.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] BiocInstaller_1.28.0

loaded via a namespace (and not attached):
[1] compiler_3.4.3 tools_3.4.3

Thank you very much!

It  might be better then to raise the error to support.bioconductor.org to let the maintainers know about it.

Ok, I will do. If the error is fixable from user side, I will report the answer here.

I have reported the error on support.bioconductor.org but we were not able to solve the problem. It seems that some dependencies are making problems, but as reinstalling my home R folder did not resolve the issue, I have come up with the following workaround:
I downloaded the latest R version from my favorite CRAN mirror, opened an R session with pathToNewR/bin/R from my terminal and installed RnBeads there. This works fine.
Note: I had to make sure that pathToNewR/library is the only entry in .libPaths(), as per default, my home R folder was also contained in .libPaths().
As a consequence, I am independent from the (maybe broken) R version and packages on my server and can finally work with RnBeads :)

If I come up with another solution that works fine with my server R version, I will post an update here.


text:
How is a principal component analysis for population structure made from a bam or other alignment file?
I'm curious how a PCA for population structure is built from an alignment/BAM or SAM file. How do programs like angsd and ngsCovar actually build a PCA? What exactly makes things come out as different and how can this be determined? If I have 3 alignments that look like this:

Locus1:
>Indiv1
ACGT
>Indiv2
ACGC
>Indiv3
ACGA

Locus2:
>Indiv1
GACT
>Indiv2
GACA
>Indiv3
GAGT

Locus3:
>Indiv1
GGTC
>Indiv2
GGCT
>Indiv3
GCTA

(These example alignments leave out quality scores for simplicity, assume that each has high quality across the board.)
What algorithm or system would be used to graph them on a PCA? How many dimensions would that PCA have? Is it possible to model this by hand so that I can understand what's going on/how it is built?

Can you update your post to somewhere mention that this is related to population structure? There are other uses of PCA using BAM files that use completely unrelated methods (e.g., deepTools).

In general, the process works as such:

Variants are found at a number of loci in each sample (missing loci are often imputed).
Each sample is scored at each locus, often according to the number of copies of the reference allele that it has (so, 0-2).
The resulting matrix is used for PCA.

@terdon "...according to the number of copies of the reference allele that it has (so, 0-2)."    So the PCA is related to heterozygosity?

You said each locus was scored often based on number of copies- what else could it be scored on?

So you build a matrix (often) based on allele copy (0, 1, or 2) then make a covariance matrix and do associated eigendecomposition stuff to generate a PCA?

One could consider making a numeric score for the genotype, such that each non-reference genotype got a different value (instead of them being all lumped into 0/1). I've not seen someone try that, so I don't know whether it's useful. Usually PCA is done via SVD these days, but really however you want to do it should lead to similarish results.


text:
How can I retreive specific qualifiers from genbank files in R?
I'm trying to extract the host field from the features table using the biofiles package. From what I understand, the function filter (from the same package) could be of help but I can't extract the information of interest. This is the simplified series of steps that I'm following:
library(biofiles)
x_search <- reutils::esearch(db = "nuccore", #retmax = 9999,
             term = "((bacteria[ORGN]) AND host[ALL]) AND Loripes[ALL]")

x_fetch <- reutils::efetch(x_search[1:3], "nuccore",rettype = "gbwithparts",
            retmode = "text")

x_rec <- biofiles::gbRecord(x_fetch)

x_ft <- biofiles::getFeatures(x_rec)

x_ft[1]

The last output should show the list of features for that entry, including the source with the qualifier "host", which I would like to retrieve for all the ids.
EDIT: I'm open to solutions outside the biofiles package.

Easy with biopython if Ur open to answers outside of R

@LiamMcIntyre Sure! I would like to know your suggestion.


text:
What is the index fastq file (sample_I*.fastq.gz) generated when demultiplexing Illumina paired-end runs?
What is the index fastq file that comes with some Illumina sequencing datasets? (The samplename_I*.fastq.gz file.)
For example, I recently received some 10X Chromium reads for two libraries sequenced on the same lane. This was a 2x150 sequencing run, so there should be two fastq files. The sequencing center demultiplexed the libraries and generated two separate directories - one for each library.

In each directory there are three fastq files:

Mysample_I1_001.fastq.gz
Mysample_R1_001.fastq.gz
Mysample_R2_001.fastq.gz

I know that the R1 fastq file contains the forward reads and the R2 fastq file contains the reverse reads, but what is the index file? The first few sequences do not match the 10X indices for this library (ACATTACT, CAGCCCAC, GGCAATGG, TTTGGGTA).
>bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | head
NTTGGGTA
NGCAATGG
NAGCCCAC
NAGCCCAC
NCATTACT
NCATTACT
NGCAATGG
NAGCCCAC
NTTGGGTA
NAGCCCAC

I don't know 10x Chromium software, but I know their single cell pipeline, and many versions ago, the software was expecting the read indices to be in a separate file along with the reads themselves.  Cellranger has since changed their pipeline so that it doesn't want these anymore, but whoever gave you your files might have seen that it was a 10xGenomics protocol, and decided to give you the indices like that just in case you needed them,

For what it's worth, 10X doesn't require the index files as input (you can just give it demultiplexed samples) and Illumina's software won't produce those index files by default. Also, you can have up to 2 index files for Illumina data (though you typically won't if you're making standard 10X libraries).

tldr - The I*.fastq.gz file contains the read index sequences.
long explanation
Illumina uses a program called bcl2fastq to demultiplex sequencing runs.
This software takes a list of samples and their associated indices and uses those sequences to make one or more fastq files per sample, binned by one or two index sequences on either end of the sequencing molecule (i5 and i7 indices, see page 6 for HiSeq).
Illumina sequencing isn't perfect though, and there are sometimes errors in reading the index sequence. For example this index, CAGCCCAC can easily have read errors in the A sandwiched between many Cs: CAGCCCAC -> CAGCCCCC.
Instead of throwing out all of the reads with indexing sequencing errors, the bcl2fastq program includes reads that are discernibly derived from a true index for that sample as long as there is no overlap with another sample.
In the example provided above, the four indices for mysample were: ACATTACT, CAGCCCAC, GGCAATGG, TTTGGGTA. If we look at all of the index sequences in the Mysample_I1_001.fastq.gz file, we will see that sequences with the four correct indices are indeed the most prevalent, however there are also reads that have indices with sequencing errors derived from the correct indices.
bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | \
    sort | uniq -c | sort -k1 -nr | head

count      index       source
41362311   CAGCCCAC    True index
37209190   GGCAATGG    True index
36863213   ACATTACT    True index
33674467   TTTGGGTA    True index
 1140358   NAGCCCAC    CAGCCCAC
 1026099   NGCAATGG    GGCAATGG
 1016754   NCATTACT    ACATTACT
  933342   NTTGGGTA    TTTGGGTA
  119626   TTTGGGGA    TTTGGGTA
   98657   GTTGGGTA    TTTGGGTA
   96625   GGCAATGA    GGCAATGG

It took me a while to figure out that the "index" is the same thing as the "barcode" that says which sample each sequence is from on a multiplexed run.
If your data is not demultiplexed (single R1.fastq and R2.fastq files contain the information for multiple samples), then this I1.fastq file is what you use to assign each sequence to a sample (ie to "demultiplex"). The spreadsheet listing which barcode/index is which sample is often called a "mapping file."
If your data is already demultiplexed (you have separate R1.fastq and R2.fastq files for each sample), then you don't really need it... however some people use the index sequences in their quality control workflows. For example, this 2016 paper in BMC Genomics by Wright and Vetsigien "Quality filtering of Illumina index reads mitigates sample cross-talk": https://doi.org/10.1186/s12864-016-3217-x

The term index is almost always used to refer to the 8-bp sequence that the Illumina machine reads to identify the library-of-origin of the individual read, while the term barcode could also mean some sequence within R1 or R2 itself that the user added to label provenance - such as being from a specific cell or specific amplicon.

Yes, in the context of a 10xGenomics single cell run, the "barcode" means the cell barcode, which is not at all the same thing as the sample index.

I received 3 files that exactly match the file names listed by the original poster, except that my data was clearly not demultiplexed yet. My index file referred to what my mapping file listed as the "SampleBarcode".  It's certainly possible that formats are different for single-cell data. For perspective, my data is from a multiplexed 16S rRNA gene amplicon miSeq run. My barcodes identified which sample each sequence was from, and they were a 12 bp sequence attached to the forward primer (so it was not dual-indexed). This is the protocol used by the Earth Microbiome Project.

I looked at 10x data a couple of years ago, and at that time, their protocol would use four separate indexes, with balanced index sequences, as a way of introducing diversity into the library.  These four indexes were separate from the individual oligos used to mark single cells.  For an individual sample, you would need to combine the reads from all four separate indexes, which is how it looks like in your data.  There is a brief explanation here: https://support.10xgenomics.com/single-cell-gene-expression/index/doc/specifications-sample-index-sets-for-single-cell-3


text:
How to generate two fasta files from samtools phase output?
samtools phase is an easy-to-use tool to phase variants from a bam file and a reference.
For a sequence that has two haplotypes, like this:
   position: 123456789
haplotype 0: GACCATATT
       SNPs:  **    *
haplotype 1: GTACATACT

The output looks like this:
CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
CC
CC
PS      frame-1   2    9
M1      frame-1  60    2      A       T       1       3       0       23      0
M1      frame-1  60    3      C       A       2       35      1       24      0
M1      frame-1  60    8      T       C       3       46      0       30      0

How can I take the reference fasta file and the output of samtools phase to generate one fasta file for each of the haplotypes?

Answer January 20th, 2020. If you run samtools phase with the option -b my_prefix it outputs a bam file for both haplotypes, named my_prefix.0.bam and my_prefix.1.bam. If the sequence has been correctly phased, use these bam files to correct the reference using pilon. This is preferable to the method described in the January 14th, 2020 answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that samtools phase found. samtools phase does not find indels.
Answer January 14th, 2020. This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.
Run the program like this:
samtools phase reads.bam | phase2fasta.py ref.fa > haplotypes.fasta.
Here is the code of phase2fasta.py:
#!/usr/bin/env python

import os
import sys

if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
    raise Exception("The first and only argument is the reference fasta.")
allele0 = []
seqname = ""
with open(sys.argv[1]) as f:
    for line in f:
        if line[0] == ">":
            seqname = line[1:].split(" ")[0].strip()
        else:
            for c in line.strip().upper():
                allele0.append(c)
allele1 = [x for x in allele0]
for line in sys.stdin:
    if line[0:2] == "M1":
        fields = line.strip().split()
        if fields[1] == seqname:
            pos = int(fields[3]) - 1
            allele0[pos] = fields[4]
            allele1[pos] = fields[5]
print(">{}_allele0".format(seqname))
print("".join(allele0))
print(">{}_allele1".format(seqname))
print("".join(allele1))

Oh hey, thanks for this! We need this sort of self-answered post!


text:
What is the fastest way to get the reverse complement of a DNA sequence in python?
I am writing a python script that requires a reverse complement function to be called on DNA strings of length 1 through around length 30. Line profiling programs indicate that my functions spend a lot of time getting the reverse complements, so I am looking to optimize.
What is the fastest way to get the reverse complement of a sequence in python? I am posting my skeleton program to test different implementations below with DNA string size 17 as an example.
#!/usr/bin/env python
import random
import timeit

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

By the way, I get output like this. It varies by the call, of course!
naive implementation
       1.83880s total,
       271916.7 strings per second
       -0.7% increase over baseline
global dict implementation
       1.74645s total,
       286294.3 strings per second
       4.3% increase over baseline

Edit: Great answers, everyone! When I get a chance in a day or two I will add all of these to a test file for the final run. When I asked the question, I had not considered whether I would allow for cython or c extensions when selecting the final answer. What do you all think?
Edit 2: Here are the results of the final simulation with everyone's implementations. I am going to accept the highest scoring pure python code with no Cython/C. For my own sake I ended up using user172818's c implementation. If you feel like contributing to this in the future, check out the github page I made for this question.
the runtime of reverse complement implementations.
10000 strings and 250 repetitions
╔══════════════════════════════════════════════════════╗
║ name                       %inc   s total  str per s ║
╠══════════════════════════════════════════════════════╣
║ user172818 seqpy.c        93.7%  0.002344  4266961.4 ║
║ alexreynolds Cython (v2)  93.4%  0.002468  4051583.1 ║
║ alexreynolds Cython (v1)  90.4%  0.003596  2780512.1 ║
║ devonryan string          86.1%  0.005204  1921515.6 ║
║ jackaidley bytes          84.7%  0.005716  1749622.2 ║
║ jackaidley bytesstring    83.0%  0.006352  1574240.6 ║
║ global dict                5.4%  0.035330   283046.7 ║
║ revcomp_translateSO       45.9%  0.020202   494999.4 ║
║ string_replace            37.5%  0.023345   428364.9 ║
║ revcom from SO            28.0%  0.026904   371694.5 ║
║ naive (baseline)           1.5%  0.036804   271711.5 ║
║ lambda from SO           -39.9%  0.052246   191401.3 ║
║ biopython seq then rc    -32.0%  0.049293   202869.7 ║
╚══════════════════════════════════════════════════════╝

You might also get some good advice from CodeReview.SE.

Your implementation of my approach is not doing what I suggested. I suggested working with bytes instead of strings throughout. This means you need your DNAStrings to be in bytes instead of a string and so it would need a separate generation function.

Ah, you meant use them for the entire program. I didn't think to do that. If I were to test that then I would need to convert the entire list of strings to bytestrings before testing, correct?

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

def reverse_complement_table(seq):
    return seq.translate(tab)[::-1]

The thing with hashing is that it adds a good bit of overhead for a replacement set this small.
For what it's worth, I added that to your code as "with a translation table" and here is what I got on my workstation:
global dict implementation
       1.37599s total,
       363374.8 strings per second
       3.3% increase over baseline
naive implementation
       1.44126s total,
       346919.4 strings per second
       -1.3% increase over baseline
with a translation table
       0.16780s total,
       2979755.6 strings per second
       88.2% increase over baseline

If you need python 3 rather than python 2, then substitute tab = str.maketrans("ACTG", "TGAC") for tab = string.maketrans("ACTG", "TGAC"), since maketrans is now a static method on the str type.
For those wondering, using biopython is slower for this (~50% slower than the naive implementation), presumably due to the overhead of converting the strings to Seq objects. If one were already reading sequences in using biopython, though, I wouldn't be surprised if the performance was much different.

@Chris_Rands True, updated with the change needed for python3 (thankfully, it's only a single line difference).

Yeh, it''s str.makestrans though, Biopython's implementation is pretty much the same but with some more alphabet and sanity checks https://github.com/biopython/biopython/blob/master/Bio/Seq.py#L860

Here's a Cython approach that might suggest a generic approach to speeding up Python work.
If you're manipulating (ASCII) character strings and performance is a design consideration, then C or Perl are probably preferred options to Python.
In any case, this Cython test uses Python 3.6.3:
$ python --version
Python 3.6.3 :: Anaconda custom (64-bit)

To install Cython, e.g.:
$ conda install -c anaconda cython

The Cython code below seems to offer about the same speed bump as the translation table — perhaps similar code is run under the hood of that.
Two files are needed, starting with setup.py:
from distutils.core import setup
from Cython.Build import cythonize

setup(
  name = 'Reverse complement C test',
  ext_modules = cythonize("revcomp_c.pyx"),
)

And then a second file called revcomp_c.pyx:
from libc.stdlib cimport malloc

cdef int seq_len = 17
cdef char *seq_dest = <char *>malloc(seq_len + 1)
seq_dest[seq_len] = '\0'

def reverse_complement_c_v1(str seq):
    cdef bytes py_bytes = seq.encode('ascii')
    cdef char *seq_src = py_bytes
    cdef int i = 0
    cdef int d = 0
    for i in range(seq_len):
        d = seq_len - i - 1
        if seq_src[i] == 'A':
            seq_dest[d] = 'T'
        elif seq_src[i] == 'G':
            seq_dest[d] = 'C'
        elif seq_src[i] == 'C':
            seq_dest[d] = 'G'
        elif seq_src[i] == 'T':
            seq_dest[d] = 'A'
    return seq_dest[:seq_len].decode('ascii')

This can be compiled into a Python module like so:
$ python setup.py build_ext --inplace

Then we can modify the test bench to include this Cython module and the relevant test method:
#!/usr/bin/env python
import random
import timeit
from revcomp_c import reverse_complement_c_v1

global complement
complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}

DNAlength = 17

#randomly generate 100k bases
int_to_basemap = {1: 'A', 2: 'C', 3: 'G', 4: 'T'}
num_strings = 500000
random.seed(90210)
DNAstrings = ["".join([int_to_basemap[random.randint(1,4)] for i in range(DNAlength)])
              for j in range(num_strings)]
#get an idea of what the DNAstrings look like
print(DNAstrings[0:5])

def reverse_complement_naive(seq):
    this_complement = {'A': 'T', 'C': 'G', 'G': 'C', 'T': 'A'}
    return "".join(this_complement.get(base, base) for base in reversed(seq))

def reverse_complement(seq):
    return "".join(complement.get(base, base) for base in reversed(seq))

tic=timeit.default_timer()
rcs = [reverse_complement_naive(seq) for seq in DNAstrings]
toc=timeit.default_timer()
baseline = toc - tic
namefunc = {"naive implementation": reverse_complement_naive,
            "global dict implementation": reverse_complement,
            "Cython implementation (v1)": reverse_complement_c_v1}

for function_name in namefunc:
    func = namefunc[function_name]
    tic=timeit.default_timer()
    rcs = [func(seq) for seq in DNAstrings]
    toc=timeit.default_timer()
    walltime = toc-tic
    print("""{}
       {:.5f}s total,
       {:.1f} strings per second
       {:.1f}% increase over baseline""".format(
           function_name,
           walltime,
           num_strings/walltime,
           100- ((walltime/baseline)*100)   ))

A sample run:
$ python ./revcomp_test.py
['ACTGCAATAGACCACGC', 'CAGCTTGAGCCATTAAC', 'GGCCCAAGAGTTCGAAC', 'CGACTGTCTCGAATTGT', 'ATCCGCTATGCGCCGTC']
naive implementation
       1.73295s total,
       288525.0 strings per second
       -2.5% increase over baseline
global dict implementation
       1.70896s total,
       292575.5 strings per second
       -1.0% increase over baseline
Cython implementation (v1)
       0.21458s total,
       2330158.7 strings per second
       87.3% increase over baseline

One easy way to speed this up is to use a static const unsigned char array as an ASCII lookup table, which maps a residue directly to its complement. This would replace the nest of if statements and probably give a nice little boost (and it appears it does, making it among the best performers so far!).
Here is my fast implementation of a reverse complement function in C:
https://gist.github.com/alexpreynolds/4f75cab4350e9d937f4a
You might be able to use this directly in Python via the subprocess library. Outsourcing the reverse complement step to a utility written in C will almost always beat the best that Python can do, and you can do nice and important things like bounds checking etc. without losing much speed.
It's unclear how "pure" the answer needs to be, but making a system call from Python seems fair if you're processing strings and your goal is performance.
Another direction to take may be to look at multithreading, if you don't need ordered output. If you have many thousands of sequences stored in memory, you could split an array of sequences up into smaller arrays by use of offsets or array indices. Each thread would work on "rc"-ing sequences in its own piece of the array.

Nice, +1 for Cython

Nice! Can you compare this to the .translate() version in my code? It'd be interesting to see how much Cython is able to speed things up in comparison given the general python overhead.

@Devon_Ryan: With this test bench, the "Cython implementation (v2)" on my Python 3 setup gave a 91.1% increase over baseline and "table" (translate) gave a 84.6% increase. I'm not sure how a Python 2 Cython setup compares.

The most reliable and simplest way is probably using Biopython:
from Bio.Seq import Seq

def my_reverse_complement(seq):
    return Seq(seq).reverse_complement()

print(my_reverse_complement('ATGCGTA'))
# TACGCAT

As Devon has already said here using Biopython isn't as fast as the naive Python solution, and I also tested that shown here with ipython. However, this is because Biopython's implementation, although similar to the naive approach, includes other features; it can reverse complement RNA as well as DNA and it will tell you if you're mixing DNA and RNA
Note that if you really want a fast way you could look at Cython or another python extension. As I edit this now, there are several nice answers taking this approach from user172818 and Alex Reynolds.

@JackAidley I mentioned in my own reply that biopython is ~50% slower than the naive code in the original post. It's good that this one actually included the code for that, though.

Another python extension but without cython. The source code is available at the bottom of this answer or from this gist. On Mac with Python3:
string
   0.38166s total,
   1310067.4 strings per second
   84.4% increase over baseline
seqpy
   0.20524s total,
   2436182.7 strings per second
   91.6% increase over baseline

On Linux with Python2 (seqpy is the first):
seqpy
   0.10960s total,
   4561842.5 strings per second
   93.3% increase over baseline
string
   0.17811s total,
   2807330.9 strings per second
   89.1% increase over baseline

This is file setup.py:
from setuptools import setup, Extension

seqpy_module = Extension('seqpy', sources = ["seqpy.c"])

setup(
    name = "seqpy",
    description = "seqpy",
    ext_modules = [seqpy_module])

This is file seqpy.c:
#include <stdlib.h>
#include <stdint.h>
#include <Python.h>

static unsigned char seq_comp_table[256] = {
      0,   1,   2,   3,   4,   5,   6,   7,   8,   9,  10,  11,  12,  13,  14,  15,
     16,  17,  18,  19,  20,  21,  22,  23,  24,  25,  26,  27,  28,  29,  30,  31,
     32,  33,  34,  35,  36,  37,  38,  39,  40,  41,  42,  43,  44,  45,  46,  47,
     48,  49,  50,  51,  52,  53,  54,  55,  56,  57,  58,  59,  60,  61,  62,  63,
     64, 'T', 'V', 'G', 'H', 'E', 'F', 'C', 'D', 'I', 'J', 'M', 'L', 'K', 'N', 'O',
    'P', 'Q', 'Y', 'S', 'A', 'A', 'B', 'W', 'X', 'R', 'Z',  91,  92,  93,  94,  95,
     64, 't', 'v', 'g', 'h', 'e', 'f', 'c', 'd', 'i', 'j', 'm', 'l', 'k', 'n', 'o',
    'p', 'q', 'y', 's', 'a', 'a', 'b', 'w', 'x', 'r', 'z', 123, 124, 125, 126, 127,
    128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
    144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159,
    160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175,
    176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191,
    192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207,
    208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223,
    224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239,
    240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255
};

static PyObject *seqpy_revcomp(PyObject *self, PyObject *args)
{
    PyObject *r;
    char *seq, *rev;
    int i, len;
    PyArg_ParseTuple(args, "s#", &seq, &len);
    rev = (char*)malloc(len);
    for (i = 0; i < len; ++i)
        rev[len - i - 1] = seq_comp_table[(uint8_t)seq[i]];
    r = Py_BuildValue("s#", rev, len);
    free(rev);
    return r;
}

static PyMethodDef seqpy_methods[] = {
    {"revcomp", seqpy_revcomp, METH_VARARGS, "Reverse complement a DNA sequence"},
    {NULL, NULL, 0, NULL}
};

#if PY_MAJOR_VERSION >= 3
static struct PyModuleDef seqpy_module = { PyModuleDef_HEAD_INIT, "seqpy", NULL, -1, seqpy_methods };
PyMODINIT_FUNC PyInit_seqpy(void) { return PyModule_Create(&seqpy_module); }
#else
PyMODINIT_FUNC initseqpy(void) { Py_InitModule3("seqpy", seqpy_methods, NULL); }
#endif

To compile and test:
python setup.py build_ext -i
python -c 'import seqpy; print(seqpy.revcomp("GGGTT"))'

Here is a revision of my original Cython answer which incorporates my suggestion to use a char lookup array:
from libc.stdlib cimport malloc

cdef int seq_len = 17
cdef char *seq_dest = <char *>malloc(seq_len + 1)
seq_dest[seq_len] = '\0'

cdef char *basemap = [ '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',   '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  'T', '\0',  'G', '\0', '\0', '\0',  'C', '\0', '\0', '\0', '\0', '\0', '\0',  'N', '\0', '\0', '\0', '\0', '\0',  'A',  'A', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  't', '\0',  'g', '\0', '\0', '\0',  'c', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0', '\0',  'a',  'a' ]

def reverse_complement_c_v2(str seq):
    cdef bytes py_bytes = seq.encode('UTF-8')
    cdef char *seq_src = py_bytes
    cdef int i = 0
    for i in range(seq_len):
        seq_dest[seq_len - i - 1] = basemap[<int>seq_src[i]]
    return seq_dest[:seq_len].decode('UTF-8')

In testing under OS X, Python 3.6.3:
$ python ./revcomp.py
['ACTGCAATAGACCACGC', 'CAGCTTGAGCCATTAAC', 'GGCCCAAGAGTTCGAAC', 'CGACTGTCTCGAATTGT', 'ATCCGCTATGCGCCGTC']
...
Cython implementation (v1)
       0.20664s total,
       2419626.4 strings per second
       88.3% increase over baseline
Cython implementation (v2)
       0.14170s total,
       3528668.5 strings per second
       92.0% increase over baseline

Using my lookup array approach ("v2") adds a very decent performance bump over using if blocks ("v1"), and you can keep everything as a Python string.

Use a bytearray instead of a string and then employ maketrans to translate
You do not need the more advanced string encoding capabilities of string to store a string of bases, but you're still paying for it in performance. Devon Ryan's suggestion of maketrans is the huge improvement, 10x faster than your naive implementation. Using the same approach, but swapping everything out for bytes allows a further 40% speed improvement, however:
# DNA sequence generation also needs to be swapped to bytes
tab_b = bytes.maketrans(b"ACTG", b"TGAC")

def reverse_complement_bytes(seq):
    return seq.translate(tab_b)[::-1]

In my comparison:
Table implementation
       0.24339s total,
       2054315.7 strings per second
       82.5% increase over baseline
Bytes table implementation
       0.14985s total,
       3336755.9 strings per second
       89.2% increase over table

If one needs to convert back to string to interface with the rest of the code, what is the impact on speed ?

@bli It is still about 10% faster if you can work with bytes all the way through and then transfer to a string at the end. If you need to go string->bytes->string then it is about 25-30% slower than staying with strings

Since at least version 1.71 of biopython you can use Bio.Seq.reverse_complement, which also works on plain strings natively (no conversion to Seq objects).  On my mac I get 800k strings converted with that implementation ("biopython just rc") when using the benchmark.
10000 strings and 250 repetitions
                         percent increase over baseline  seconds total strings per second
name
alexreynolds Cython (v2)                          94.0%       0.002491          4015188.7
user172818 seqpy.c                                93.2%       0.002833          3529944.4
alexreynolds Cython (v1)                          90.7%       0.003885          2573863.7
devonryan string                                  86.0%       0.005804          1722983.5
jackaidley bytes                                  85.7%       0.005959          1678085.3
jackaidley bytesstring                            82.1%       0.007459          1340693.2
biopython just rc                                 70.5%       0.012275           814638.6
revcomp_translateSO                               49.9%       0.020822           480254.6
revcom from SO                                    41.7%       0.024253           412320.2
string_replace                                    35.2%       0.026946           371107.3
lambda from SO                                    10.7%       0.037111           269464.0
global dict                                        0.0%       0.041556           240641.7
biopython seq then rc                            -30.9%       0.054430           183721.9
naive (baseline)                                  -0.9%       0.041962           238313.6

def revcomplement(s):
    s = s.upper()
    c = { "A":"T","C":"G","T":"A","G":"C"}
    co = str()
    for i in s:
        co = c[i] + co
    return co

I don't think this piece of code actually "reverts" the sequence but just changes the bases with their complementary bases.

No it reverts the sequence and gives reverse complement.

My bad, co = c[i] + co does the trick. Could be useful if you run a benchmark along with the other offered solutions.

Thanks for joining the community Amardeep. I don't doubt that your code works, but I am a bit sceptical if it answers the original question (seeking for fastest solution). As a matter of fact, your solution is sort of included in the question already (reverse_complement_naive). There are plenty of questions that need to be answered, just make sure that you are really addressing what the person have asked :-)

Sorry.  I implement what you said. Thanks

authors: Костя Дряхлов


text:
How to calculate the memory usage of storing kmers in RAM
I want to write a program in C++ that stores kmers in a hash or in a trie. How can I calculate how much RAM I would need for each type of data structure? For this application my kmers are strand-specific, so I cannot reduce the size complexity by considering canonical kmers.
For a hash, for example, I know that there are $4^k$ possible kmers of length $k$. I also am assuming that DNA character with an alphabet $\in \{G,A,T,C\}$ can be stored in two bits ($2^2 = 4$). Using this logic the kmer hash should occupy $(4^k) \cdot 2 \cdot k$ bits. Plugging this into Wolfram Alpha ((4^13) * 2 * 13 bits to gigabytes), I calculated the amount of RAM that some small kmer sizes would require to store the entire hash in memory in GB.
k = 13 requires    0.2181 GB
k = 15 requires    4.027  GB
k = 17 requires   73.01  GB
k = 19 requires 1306   GB

This escalated rather quickly, and I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?
The kmers that I want to store are sometimes sparse compared to the entire space (for k=13, only 29 million or so out of 67 million possibilities).
Constraints

I do not need to store the number of occurences of a kmer, just its existence
I cannot risk false positives. Bloom filters won't work in this case. :^/

To follow up on Devon Ryan's answer, I thought it would be a little fun to write a Python script that demonstrates using a bit array to maintain a presence/absence table.
Note: I wrote a C++ port that includes a custom bitset implementation that can be sized at runtime. This and the Python script are available on Github: https://github.com/alexpreynolds/kmer-boolean
#!/usr/bin/env python

import sys
import os
import bitarray

# read FASTA
def read_sequences():
    global seqs
    seqs = []
    seq = ""
    for line in sys.stdin:
        if line.startswith('>'):
            if len(seq) > 0:
                seqs.append(seq)
                seq = ""
        else:
            seq += line.strip()
    seqs.append(seq)

# build and initialize bit array
def initialize_bitarray():
    global ba
    ba = bitarray.bitarray(4**k)
    ba.setall(False)
    sys.stderr.write("Memory usage of bitarray.bitarray instance is [%ld] bytes\n" % (ba.buffer_info()[1]))

# process sequences
def process_sequences():
    global observed_kmers
    observed_kmers = {}
    for seq in seqs:
        for i in range(0, len(seq)):
            kmer = seq[i:i+k]
            if len(kmer) == k:
                observed_kmers[kmer] = None
                idx = 0
                for j in range(k-1, -1, -1):
                    idx += 4**(k-j-1) * bm[kmer[j]]
                ba[idx] = True

def test_bitarray():
    test_idx = 0
    for j in range(k-1, -1, -1):
        test_idx += 4**(k-j-1) * bm[test_kmer[j]]
    test_result = ba[test_idx]
    if test_result:
        sys.stdout.write("%s found\n" % (test_kmer))
        sys.exit(os.EX_OK)
    else:
        sys.stdout.write("%s not found\n" % (test_kmer))
        sys.exit(os.EX_DATAERR)

def main():
    global k
    k = int(sys.argv[1])

    global bm
    bm = { 'A' : 0, 'C' : 1, 'T' : 2, 'G' : 3 }

    read_sequences()
    initialize_bitarray()
    process_sequences()

    try:
        global test_kmer
        test_kmer = sys.argv[2]
        if len(test_kmer) == k:
            test_bitarray()
        else:
            raise ValueError("test kmer (%s) should be of length k (%d)" % (test_kmer, k))
    except IndexError as err:
        keys = list(observed_kmers.keys())
        for i in range(0, len(keys)):
            sys.stdout.write("%s found\n" % (keys[i]))

    sys.exit(os.EX_OK)

if __name__== "__main__":
    main()

Note that this doesn't look at canonical kmers, e.g., AG is considered a distinct 2mer from its reverse complement CT.
To use this script, you pipe in your FASTA, specify the k, and an optional kmer that you want to test for presence/absence, e.g.:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 2 CG
Memory usage of bitarray.bitarray instance is [2] bytes
CG found

Or:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 3 AAA
Memory usage of bitarray.bitarray instance is [8] bytes
AAA not found

Or if the optional test kmer is left out:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 5
Memory usage of bitarray.bitarray instance is [128] bytes
CATTC found
ATTCT found
TTCTC found
TCTCG found
CTCGG found
TCGGG found
CGGGA found
GGGAC found
TTATA found
TATAT found
TTTAT found
TTATT found
TATTA found
ATTAG found
TTAGA found
TAGAC found
AGACC found
GACCT found
ACCTC found
CCTCT found

Or for the ~67M kmers in a 13-mer set, for which a roughly 8.4MB bit array is reserved:
$ echo -e ">foo\nCATTCTC\nGGGAC\n>bar\nTTATAT\n>baz\nTTTATTAG\nACCTCT" | ./kmer-bool.py 13
Memory usage of bitarray.bitarray instance is [8388608] bytes
TTTATTAGACCTC found
TTATTAGACCTCT found

I think that my calculations must be wrong. Otherwise, how could programs count kmers in RAM?

Hash table based k-mer counters only keep k-mers seen in data. For $16<k\le32$, you need 64-bit integers to store a k-mer. Given $n$ distinct k-mers, a naive implementation with open addressing hash table would roughly take $2n\cdot 64/8=16n$ bytes. We assume the table is half full – this is where factor 2 comes from. There are various tricks to reduce the memory footprint. You'd better read papers. On a side note, don't use std::unordered_map. It wastes extra memory.
A bitmap described by other answers takes $4^k/8=2^{2k-3}$ bytes, which is a better choice sometimes, but rarely used by k-mer counters as it does not work with long k, either. Trie or full-text indices are overkilling for fixed k-mers and are usually slower.
By the way, this question reminds me of your question on efficient reverse complement. Note that for k-mer counting, you traverse both strands at the same time. You should not reverse complement the whole k-mer.

If you simply need a set of present/absent k-mers and your kmer length is short, then convert the sequence length to an integer (a 13-mer would require 26 bits or 4 bytes) and use that as an index into a bit array. Initialize the array to 0 and use 1 for "k-mer present". If you have ~67 million k-mers then your array is 67 million bits, which is only about 8.5 megs.
This will scale "OK", but eventually balloon to a huge size if you use larger k-mers. At that point you simply have to use other data structures. The reality is that most applications are fairly tolerant to small amounts of error.


text:
Is nanopolish worth it since faster polishing software is available?
For Oxford Nanopore contigs produced by any long-read assembler has anyone performed any benchmarks to compare the polishing tools racon, nanopolish, or something like wtpoa-cns?
The compute time of nanopolish is quite long and I would rather use another tool if it will suffice. After long-read polishing I will polish more with two rounds of pilon using 200x coverage Illumina shotgun data from the same individual.

Nanopolish is necessary if you want to get high-quality consensus. racon etc don't use signal data. They can't achieve high quality. pilon at times doesn't work well. Even if it worked perfectly, it wouldn't help repetitive regions inaccessible by short reads.

You're quite right. I polished twice with racon then twice with pilon. Mapped short read pairs and long reads. All of the gaps in the short read mappings were covered with long reads. Going to use nanopolish now!


text:
What is the most compact data structure for canonical k-mers with the fastest lookup time?
edit: Results are current as of Dec 4, 2018 13:00 PST.
Background
K-mers have many uses in bioinformatics, and for this reason it would be useful to know the most RAM-efficient and fastest way to work with them programmatically. There have been questions covering what canonical k-mers are, how much RAM k-mer storage theoretically takes, but we have not yet looked at the best data structure to store and access k-mers and associated values with.
Question
What data structure in C++ simultaneously allows the most compact k-mer storage, a property, and the fastest lookup time? For this question I choose C++ for speed, ease-of-implementation, and access to lower-level language features if desired. Answers in other languages are acceptable, too.
Setup

For benchmarking:

I propose to use a standard fasta file for everyone to use. This program, generate-fasta.cpp, generates two million sequences ranging in length between 29 and 300, with a peak of sequences around length 60.
Let's use k=29 for the analysis, but ignore implementations that require knowledge of the k-mer size before implementation. Doing so will make the resulting data structure more amenable to downstream users who may need other sizes k.
Let's just store the most recent read that the k-mer appeared in as the property to retrieve during k-mer lookup. In most applications it is important to attach some value to each k-mer such as a taxon, its count in a dataset, et cetera.
If possible, use the string parser in the code below for consistency between answers.
The algorithm should use canonical k-mers. That is, a k-mer and its reverse complement are considered to be the same k-mer.

Here is generate-fasta.cpp. I used the command g++ generate_fasta.cpp -o generate_fasta to compile and the command ./generate_fasta > my.fasta to run it:
return 0;
//generate a fasta file to count k-mers
#include <iostream>
#include <random>

char gen_base(int q){
  if (q <= 30){
    return 'A';
  } else if ((q > 30) && (q <=60) ){
    return 'T';
  } else if ((q > 60) && (q <=80) ){
    return 'C';
  } else if (q > 80){
    return 'G';
  }
  return 'N';
}

int main() {
  unsigned seed = 1;
  std::default_random_engine generator (seed);
  std::poisson_distribution<int> poisson (59);
  std::geometric_distribution<int> geo (0.05);
  std::uniform_int_distribution<int> uniform (1,100);
  int printval;
  int i=0;
  while(i<2000000){
    if (i % 2 == 0){
      printval = poisson(generator);
    } else {
      printval = geo(generator) + 29;
    }
    if (printval >= 29){
      std::cout << '>' << i << '\n';
      //std::cout << printval << '\n';
      for (int j = 0; j < printval; j++){
        std::cout << gen_base(uniform(generator));
      }
      std::cout << '\n';
      i++;
    }
  }
  return 0;
}

Example
One naive implementation is to add both the observed k-mer and its reverse complement as separate k-mers. This is obviously not space efficient but should have fast lookup. This file is called make_struct_lookup.cpp. I used the following command to compile on my Apple laptop (OS X): clang++ -std=c++11 -stdlib=libc++ -Wno-c++98-compat make_struct_lookup.cpp -o msl.
#include <fstream>
#include <string>
#include <map>
#include <iostream>
#include <chrono>
//build the structure. measure how much RAM it consumes.
//then measure how long it takes to lookup in the data structure

#define k 29

std::string rc(std::string seq){
  std::string rc;
  for (int i = seq.length()-1; i>=0; i--){
    if (seq[i] == 'A'){
      rc.push_back('T');
    } else if (seq[i] == 'C'){
      rc.push_back('G');
    } else if (seq[i] == 'G'){
      rc.push_back('C');
    } else if (seq[i] == 'T'){
      rc.push_back('A');
    }
  }
  return rc;
}

int main(int argc, char* argv[]){
  using namespace std::chrono;
  //initialize the data structure
  std::string thisline;
  std::map<std::string, int> kmer_map;
  std::string header;
  std::string seq;
  //open the fasta file
  std::ifstream inFile;
  inFile.open(argv[1]);

  //construct the kmer-lookup structure
  int i = 0;
  high_resolution_clock::time_point t1 = high_resolution_clock::now();
  while (getline(inFile,thisline)){
    if (thisline[0] == '>'){
      header = thisline.substr(1,thisline.size());
      //std::cout << header << '\n';
    } else {
      seq = thisline;
      //now add the kmers
      for (int j=0; j< thisline.size() - k + 1; j++){
        kmer_map[seq.substr(j,j+k)] = stoi(header);
        kmer_map[rc(seq.substr(j,j+k))] = stoi(header);
      }
      i++;
    }
  }
  std::cout << "  -finished " << i << " seqs.\n";
  inFile.close();
  high_resolution_clock::time_point t2 = high_resolution_clock::now();
  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << time_span.count() << " seconds to load the array." << '\n';

  //now lookup the kmers
  inFile.open(argv[1]);
  t1 = high_resolution_clock::now();
  int lookup;
  while (getline(inFile,thisline)){
    if (thisline[0] != '>'){
      seq = thisline;
      //now lookup the kmers
      for (int j=0; j< thisline.size() - k + 1; j++){
        lookup = kmer_map[seq.substr(j,j+k)];
      }
    }
  }
  std::cout << "  - looked at " << i << " seqs.\n";
  inFile.close();
  t2 = high_resolution_clock::now();
  time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << time_span.count() << " seconds to lookup the kmers." << '\n';

}

Example output
I ran the above program with the following command to log peak RAM usage. The amount of time the lookup of all k-mers in two million sequences is reported by the program. /usr/bin/time -l ./msl my.fasta.
The output was:
 -finished 2000000 seqs.
562.864 seconds to load the array.
  - looked at 2000000 seqs.
368.734 seconds to lookup the k-mers.
     1046.94 real       942.38 user        78.96 sys
11680514048  maximum resident set size

So, the program used 11680514048 bytes = 11.68GB of RAM and it took 368.734 seconds to lookup the k-mers in two million fasta files.
Results
Below is a plot of the results from each user's answers.

try std::unordered_map instead of std::map. This should already give you quite a boost.

@PeterMenzel, thanks for your comment. What do you think about submitting this comment as an answer? I can implement it and add the results.

You can try std::unordered_map instead of std::map as a simple improvement (hopefully)

Can you be specific about which you care about more: compactness or speed? If both are equally important, then there will likely be no answer to this question.

Your generate-fasta.cpp does not compile. I also don’t understand its logic; why the multiple random distributions? What are you modelling? The gen_base function confuses me.

YMMV depending on your compiler. It worked for me with g++ Apple LLVM version 10.0.0 (clang-1000.10.44.2). The gen_base function uses a random int between 1-100 as input to select a base. The ranges in the function give it some AT bias, The Poisson distribution makes a bunch of sequences around 59 bp long, but the geometric distribution gives it a long tail of sequence lengths. In this case I was modeling a distribution of some sequences I need to look up kmers for.

You might mark the other answer as the correct one.

The question and the accepted answer are not about k-mer data structure at all, which I will explain in detail below. I will first answer the actual question OP intends to ask.
The simplest way to keep k-mers is to use an ordinary hash table. The performance is mostly determined by the hash table library. std::unordered_map in gcc/clang is one of the worst choices because for integer keys, it is very slow. Google dense, ska::bytell_hash_map and ska::flat_hash_map, tsl::robin_map and absl::flat_hash_map are much faster. There are a few libraries that focus on smaller footprint, such as google sparse and sparsepp, but those can be a few times slower.
In addition to the choice of hash table, how to construct the key is critical. For k<=32, the right choice is to encode a k-mer with a 64-bit integer, which will be vastly better than std::string. Memory alignment is also important. In C/C++, as long as there is one 8-byte member in a struct, the struct will be 8-byte aligned on x86_64 by default. Most C++ hash table libraries pack key and value in std::pair. If you use 64-bit keys and 32-bit values, std::pair will be 8-byte aligned and use 16 bytes, even though only 12 bytes are actually used – 25% of memory is wasted. In C, we can explicitly define a packed struct with __attribute__ ((__packed__)). In C++, probably you need to define special key types. A better way to get around memory alignment is to go down to the bit level. For read mapping, for example, we only use 15–23bp seeds. Then we have 18 (=64-23*2) bits left unused. We can use these 18 bits to count k-mers. Such bit-level management is quite common.
The above is just basic techniques. There are a few other tricks. For example, 1) instead of using one hash table, we can use 4096 (=2^12) hash tables. Then we can store 12 bits of k-mer information into the 4096 part. This gives us invaluable 12 bits in each bucket to store extra information. This strategy also simplifies parallel k-mer insertions as with a good hash function, it is rare to insert into two tables at the same time. 2) when most k-mers are unique, the faster way to count k-mers is to put k-mers in an array and then sort it. Sorting is more cache friendly and is faster than hash table lookups. The downside is that sort counting can be memory demanding when most k-mers are highly repetitive.

The other answer is spending considerable (probably the majority of) time on k-mer iteration, not on hash table operations. The program loops through each position on the sequence and then each k-mer position. For an $L$-long sequence, this is an $O(kL)$ algorithm. It has worse theoretical time complexity than hash table operations, which is $O(L)$. Although hash table operations are slow due to cache misses, a factor of k=29 is quite significant. Another issue is that all programs in the question and in the other answer are compiled without -O3. Adding this option brings the bytell_hash_map lookup time from 314s to 34s on my machine.
The C program at the end of my post shows the proper way to iterate k-mers. It is an $O(L)$ algorithm with a tiny constant. The program keeps track of both forward and reverse k-mers at the same time and update them with a few bit operations at each sequence position. This echoes my previous comment "You should not reverse complement the whole k-mer". On the same machine, the program looks up k-mers in 5.5s using 792MB RAM at the peak. This 6-fold (=34/5.5) speedup mostly comes from k-mer iteration, given that the hash table library in use is known to have comparable performance to bytell_hash_map.

#include <stdio.h>
#include <stdint.h>
#include "khash.h"

static inline uint64_t hash_64(uint64_t key)
{ // more sophisticated hash function to reduce collisions
  key = (~key + (key << 21)); // key = (key << 21) - key - 1;
  key = key ^ key >> 24;
  key = ((key + (key << 3)) + (key << 8)); // key * 265
  key = key ^ key >> 14;
  key = ((key + (key << 2)) + (key << 4)); // key * 21
  key = key ^ key >> 28;
  key = (key + (key << 31));
  return key;
}

KHASH_INIT(64, khint64_t, int, 1, hash_64, kh_int64_hash_equal)

unsigned char seq_nt4_table[128] = { // Table to change "ACGTN" to 01234
  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, 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
};

static uint64_t process_seq(khash_t(64) *h, int k, int len, char *seq, int is_ins)
{
  int i, l;
  uint64_t x[2], mask = (1ULL<<k*2) - 1, shift = (k - 1) * 2, tot = 0;
  for (i = l = 0, x[0] = x[1] = 0; i < len; ++i) {
    int c = (uint8_t)seq[i] < 128? seq_nt4_table[(uint8_t)seq[i]] : 4;
    if (c < 4) { // not an "N" base
      x[0] = (x[0] << 2 | c) & mask;                  // forward strand
      x[1] = x[1] >> 2 | (uint64_t)(3 - c) << shift;  // reverse strand
      if (++l >= k) { // we find a k-mer
        uint64_t y = x[0] < x[1]? x[0] : x[1];
        khint_t itr;
        if (is_ins) { // insert
          int absent;
          itr = kh_put(64, h, y, &absent);
          if (absent) kh_val(h, itr) = 0;
          tot += ++kh_val(h, itr);
        } else { // look up
          itr = kh_get(64, h, y);
          tot += itr == kh_end(h)? 0 : kh_val(h, k);
        }
      }
    } else l = 0, x[0] = x[1] = 0; // if there is an "N", restart
  }
  return tot;
}

#include <zlib.h>
#include <time.h>
#include <unistd.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
  khash_t(64) *h;
  int i, k = 29;
  while ((i = getopt(argc, argv, "k:")) >= 0)
    if (i == 'k') k = atoi(optarg);
  h = kh_init(64);
  for (i = 1; i >= 0; --i) {
    uint64_t tot = 0;
    kseq_t *ks;
    gzFile fp;
    clock_t t;
    fp = gzopen(argv[optind], "r");
    ks = kseq_init(fp);
    t = clock();
    while (kseq_read(ks) >= 0)
      tot += process_seq(h, k, ks->seq.l, ks->seq.s, i);
    fprintf(stderr, "[%d] %.3f\n", i, (double)(clock() - t) / CLOCKS_PER_SEC);
    kseq_destroy(ks);
    gzclose(fp);
  }
  kh_destroy(64, h);
  return 0;
}

Thanks for updating your comments. Doesn't 'reduce collisions' in the hash function imply that collisions are possible? Collisions aren't OK for some applications. How often should this happen and under what circumstances?

@conchoecia you are misunderstanding what collisions mean in hash tables. Ordinary hash functions all have collisions. Perfect hash function doesn't, but constructing it is a non-trivial task and requires dedicated libraries. Perfect hash function is not always practical depending on applications.

What I think you mean by collisions -> "A collision means that two unique elements produce the same hash value. In this case for kmers, two unique DNA k-mer strings would produce the same hash value. So when counting 3-mers in the dataset AAAT, if AAT and AAA have the same hash value, our data structure would show that the k-mers AAT occurs twice and AAA occurs twice."

@conchoecia Hash table keeps the actual keys, not only their hashes. AAT and AAA are distinct keys and will only be counted once. Collision is something that every standard library has to resolve. Every implementation here has collisions. Please read the wiki page to understand how hash table works. PS: collisions have nothing to do with the results. They only affect performance.

The more I go over your code the more I am learning. This is amazing. Thank you so much.


text:
Multithread fastq processing with kseq.h in C++11?
Background
I am using the wonderful kseq.h as my fastq/fasta parser for a C++ project. My actual implementation uses something like this, where I am counting kmers of each sequence and performing a bunch of other processing on every read pair.
Thanks to @user1782818's code my program is very fast, but I want it to be faster. The only way I know how to do this is with multithreading. However, I have never worked with multithreading in C++.
For each read pair I am accessing a vector or map of structs. Within each struct I am incrementing counters or pushing values onto vectors. Essentially I am doing something like this:
//compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
// execute with ./this_prog my_file.fastq.gz
#include <iostream>//cout
#include <map>
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char **argv)
{
  //initialize prefix
  std::map<std::string, int> prefix_counter;
  std::string prefix;
  prefix.resize(2);

  gzFile fp;
  kseq_t *seq;
  int n = 0, slen = 0, qlen = 0;
  fp = gzopen(argv[1], "r");
  seq = kseq_init(fp);
  while (kseq_read(seq) >= 0){
    ++n;
    //COUNT PREFIXES
    std::copy(seq->seq.s, seq->seq.s + 2, begin(prefix));
    if (prefix_counter.count(prefix) == 0){
      prefix_counter[prefix] = 1;
    } else prefix_counter[prefix]++;
  }
  std::cout << "Total n. of sequences: " << n << '\n';
  kseq_destroy(seq);
  gzclose(fp);

  //print counts
  for( const auto& sm_pair : prefix_counter ){
    std::cout << sm_pair.first << ": " << sm_pair.second << '\n';
  }
  return 0;
}

Example output from the above program is something like this for a 100,000-read test fastq file.
Total n. of sequences: 100000
AA: 8561
AC: 5178
AG: 7634
AT: 8523
CA: 7526

Question
Q1: Is it possible to perform a task like the one shown above using multithreading and kseq.h to read in a fastq file and process the reads in chunks? I only want to use standard libraries. C++11 is OK - I know that multithreading support was added in C++11.
Q2: Can you please provide a C++ example of how to convert the above or similar program to use multithreading using kseq.h?

Note: I'm not posting this on SO to avoid "read the manual"-type responses and I believe that others in the bioinformatics SE community will benefit from an answer to this question. Thanks!

You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads.  The library contains an example program to demonstrate how to use it with multiple threads.  It's worth noting that kseq is already very fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).
I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on kseq.h, so omitting third-party libraries must not be a hard constraint.

Thanks for posting this - I will take a look. Yes, the work on each pair is nontrivial so I've been wondering how to best store reads in chunks/process them or have several threads of kseq.h read different portions of a single file.