# Question
Title: Using pysam with cython: htslib/kstring.h not found
I'm trying to learn to use cython to speed up some code based on pysam. My issue is not strictly speaking about bioinformatics, but rather about building tools using a bioinformatics library. I hope this is still relevant in this site (I hesitated posting it as a github issue to pysam, alternatively, but I was not sure that this was a problem with pysam).

In my project's folder, the `bam25prime/libbam25prime.pyx` file has the following content:

```
# cython: language_level=3
from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

```

In my `setup.py`, the `setup` function is called with this argument: `ext_modules = cythonize("bam25prime/libbam25prime.pyx")`.

When I try to do the build (running `python3.6 setup.py build_ext`), cython seems to miss some include paths, because I get the following error:

```
[1/1] Cythonizing bam25prime/libbam25prime.pyx
running build_ext
building 'bam25prime.libbam25prime' extension
gcc -pthread -Wno-unused-result -Wsign-compare -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -O3 -march=native -fomit-frame-pointer -fPIC -I/home/bli/.local/lib/python3.6/site-packages/pysam -I${HOME}/include -I/home/bli/include/python3.6m -I/home/bli/.local/include -c bam25prime/libbam25prime.c -o build/temp.linux-x86_64-3.6/bam25prime/libbam25prime.o
bam25prime/libbam25prime.c:497:28: fatal error: htslib/kstring.h: No such file or directory

```

If I manually re-run the gcc command, adding `-I/home/bli/.local/lib/python3.6/site-packages/pysam/include/htslib`, the header file gets found.

**How do I get cython to correctly determine the path to the necessary header files?**

I guess another way to formulate the question would be the following: **How is cython supposed to determine what include flags to use?**

Or maybe there is a problem with my pysam installation?

If I remember well, I installed it using `pip3.6 --user install pysam`. I have a `~/.pydistutils.cfg` file containing the following settings

```
[install]
optimize=1
[build_ext]
include_dirs=${HOME}/include
library_dirs=${HOME}/lib
rpath=${HOME}/lib
user=1

```

Any advice appreciated.

# Answer
Based on a suggestion by Devon Ryan, I searched how to set `C_INCLUDE_PATH for` cython, and found this issue. It refers to a (currently not working) feature present in cython's documentation, that should let libraries inform cython of the correct localisation for the header files in the module's `setup`:

```
ext_modules = cythonize(
    "bam25prime/libbam25prime.pyx",
    include_path=pysam_get_include())

```

As noted in the issue discussion, this seems to have no effect, but a workaround exists, consisting in using the `Extension` class to store the header localization information.

So the general answer is: **Some python libraries (including pysam) provide a function that returns the paths to their include directories. This can be used to inform cython, but has to be done "manually".**

My `setup.py` now looks as follows:

```
from setuptools import setup, find_packages
# If you have .pyx things to cythonize
from distutils.extension import Extension
from Cython.Build import cythonize
from pysam import get_include as pysam_get_include

name = "bam25prime"
__version__ = "0.1"

# https://github.com/cython/cython/blob/master/docs/src/reference/compilation.rst#configuring-the-c-build
extensions = [
    Extension(
        "bam25prime.libbam25prime", ["bam25prime/libbam25prime.pyx"],
        include_dirs=pysam_get_include()),
    ]

setup(
    name=name,
    version=__version__,
    # [description, author, email, licence, ...]
    packages=find_packages(),
    scripts=["bin/%s" % name],
    ext_modules = cythonize(extensions),
    # ext_modules = cythonize(
    #     "bam25prime/libbam25prime.pyx",
    #     include_path=pysam_get_include()),
    install_requires=["pysam", "pybedtools"])

```

Then, `python3.6 setup.py build_ext` succeeds.

And for whoever might be interested, here are the other files of the project:

`bin/bam25prime` is a link to `bam25prime/bam25prime.py`, which looks like this:

```
#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
"""
When executed, this module collapses reads in a bam file to their 5-prime end
(reducing them to length 1) and writes the result in bed format.
"""

import argparse
from pybedtools import BedTool, Interval
from pysam import AlignmentFile
from .libbam25prime import collapse_ali

# pybedtools.Interval | pysam.AlignedSegment
# chrom               | reference_name
# start               | reference_start
# end                 | reference_end
# bool(int(strand))   | is_reverse
#
# Interval(reference_name, reference_start, reference_end, strand=strand)

def collapse_and_sort(alis):
    """
    Collapse elements from *ali* to their 5 prime end.

    The resulting :class:`pybedtool.BedTool` is sorted.

    :param alis: An iterable of :class:`pysam.AlignedSegment`
    :rtype: :class:`pybedtool.BedTool`

    :class:`pybedtool.BedTool` is an iterable:
    It can be iterated over in a for loop, or produce an iterator
    using :func:`iter` on it
    (see https://stackoverflow.com/a/28353158/1878788).
    """
    return BedTool(
        "\n".join(map(collapse_ali, alis)),
        from_string=True).sort(stream=True)

def main():
    """Main function of the program."""
    parser = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.ArgumentDefaultsHelpFormatter)
    parser.add_argument(
        "-b", "--bamfile",
        required=True,
        help="Sorted and indexed bam file containing the aligned reads.")
    parser.add_argument(
        "-o", "--out_bed",
        required=True,
        help="Bed file in which the results should be written.")
    args = parser.parse_args()
    with AlignmentFile(args.bamfile, "rb") as bamfile:
        collapse_and_sort(bamfile.fetch()).saveas(args.out_bed)
    return 0

if __name__ == "__main__":
    exit(main())

```

`bam25prime/libbam25prime.pyx`:

```
#!/usr/bin/env python3
# vim: set fileencoding=<utf-8> :
# cython: language_level=3
"""
This library contains a cythonized version of a function to collapse aligned
segments to their 5' end.
"""

from pysam.libchtslib cimport BAM_FREVERSE
from pysam.libcalignedsegment cimport AlignedSegment

cdef object ccollapse_ali(AlignedSegment ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    if (ali.flag & BAM_FREVERSE) != 0:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_end - 1),
            str(ali.reference_end),
            ".",
            ".",
            "-"])
    else:
        return "\t".join([
            ali.reference_name,
            # five prime position
            str(ali.reference_start),
            str(ali.reference_start + 1),
            ".",
            ".",
            "+"])

def collapse_ali(ali):
    """
    Return the bed information corresponding to the 5' collapse
    of AlignedSegment *ali*.
    """
    return ccollapse_ali(ali)

def main():
    """Example use case."""
    print("No examples to show here.")
    return 0

if __name__ == "__main__":
    import sys
    sys.exit(main())

```

(Any advice or suggestions welcome: I'd like to improve my python package-making skills.)

---

For those interested, to use cython with `pybedtools` I had to define the following `Extension` (where `glob` is `glob.glob` and `/home/bli/src/bedtools2-2.27.1/` is the location of the source code of bedtools).

```
    Extension(
        "bam25prime.libcollapsebed", ["bam25prime/libcollapsebed.pyx"],
        include_dirs=glob("/home/bli/src/bedtools2-2.27.1/src/utils/*"),
        language="c++"),

```
> 6 votes

---
Tags: python, pysam
---


# Question
Title: How do I filter a GFF file by gene type?
I have a file with following information:

```
chr1    HAVANA  gene    135141  135895  .       -       .       gene_id "ENSG00000268903.1"; gene_type "processed_pseudogene"; gene_name "AL627309.7"; level 2; havana_gene "OTTHUMG00000182518.2";
chr1    HAVANA  gene    137682  137965  .       -       .       gene_id "ENSG00000269981.1"; gene_type "processed_pseudogene"; gene_name "AL627309.8"; level 2; havana_gene "OTTHUMG00000182738.2";
chr1    HAVANA  gene    139790  140339  .       -       .       gene_id "ENSG00000239906.1"; gene_type "antisense_RNA"; gene_name "AL627309.2"; level 2; havana_gene "OTTHUMG00000002481.1";
chr1    HAVANA  gene    141474  173862  .       -       .       gene_id "ENSG00000241860.6"; gene_type "processed_transcript"; gene_name "AL627309.5"; level 2; tag "ncRNA_host"; havana_gene "OTTHUMG00000002480.3";
chr1    ENSEMBL gene    157784  157887  .       -       .       gene_id "ENSG00000222623.1"; gene_type "snRNA"; gene_name "RNU6-1100P"; level 3;
chr1    HAVANA  gene    160446  161525  .       +       .       gene_id "ENSG00000241599.1"; gene_type "lincRNA"; gene_name "AL627309.4"; level 2; havana_gene "OTTHUMG00000002525.1";
chr1    HAVANA  gene    182696  184174  .       +       .       gene_id "ENSG00000279928.2"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.2"; level 2; havana_gene "OTTHUMG00000191962.1";
chr1    HAVANA  gene    185217  195411  .       -       .       gene_id "ENSG00000279457.4"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.1"; level 2; havana_gene "OTTHUMG00000191963.1";

```

In the above file there is a `gene_type` with different names like processed\_pseudogene, lincRNA, snRNA, unprocessed\_pseudogene, antisense\_RNA, processed\_transcript.

I want to create a new file based on `gene_type`. Among all the gene\_types I want the new file to be with certain names of `gene_type` like only processed\_pseudogene, lincRNA, unprocessed\_pseudogene, antisense\_RNA.

The new file should look like below:

```
chr1    HAVANA  gene    135141  135895  .       -       .       gene_id "ENSG00000268903.1"; gene_type "processed_pseudogene"; gene_name "AL627309.7"; level 2; havana_gene "OTTHUMG00000182518.2";
chr1    HAVANA  gene    137682  137965  .       -       .       gene_id "ENSG00000269981.1"; gene_type "processed_pseudogene"; gene_name "AL627309.8"; level 2; havana_gene "OTTHUMG00000182738.2";
chr1    HAVANA  gene    139790  140339  .       -       .       gene_id "ENSG00000239906.1"; gene_type "antisense_RNA"; gene_name "AL627309.2"; level 2; havana_gene "OTTHUMG00000002481.1";
chr1    HAVANA  gene    160446  161525  .       +       .       gene_id "ENSG00000241599.1"; gene_type "lincRNA"; gene_name "AL627309.4"; level 2; havana_gene "OTTHUMG00000002525.1";
chr1    HAVANA  gene    182696  184174  .       +       .       gene_id "ENSG00000279928.2"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.2"; level 2; havana_gene "OTTHUMG00000191962.1";
chr1    HAVANA  gene    185217  195411  .       -       .       gene_id "ENSG00000279457.4"; gene_type "unprocessed_pseudogene"; gene_name "FO538757.1"; level 2; havana_gene "OTTHUMG00000191963.1";

```

# Answer
> 4 votes
```
grep -E 'gene_type "(processed_pseudogene|lincRNA|unprocessed_pseudogene|antisense_RNA)";' input_file.txt

```

# Explanation

* `grep -E`: use extended regular expression so you don't need to escape things.
* This kind of expression `(foo|bar|baz)` allows grep to include either of these strings. Hence, the total expression will match `gene_type "processed_pseudogene";'`, `gene_type "lincRNA";'`, etc.

---
Tags: linux, shell
---


title: Installing a R package from git

text:
# Question
Title: Installing a R package from git
I am trying to install STEMNET R package but a permanent error stops me, what can I do please?

```
> > library(devtools)
> > devtools::install_git("https://git.embl.de/velten/STEMNET/", build_vignettes=TRUE) Downloading git repo
> https://git.embl.de/velten/STEMNET/ Installation failed: Error in
> 'git2r_clone': SSL error: syscall failure: Connection reset by peer
>
> >
> > devtools::install_git("https://git.embl.de/velten/STEMNET/", build_vignettes=FALSE) Downloading git repo
> https://git.embl.de/velten/STEMNET/ Installation failed: Error in
> 'git2r_clone': SSL error: syscall failure: Connection reset by peer
>
> >
> > library(git2r)
>
> Attaching package: ‘git2r’
>
> The following objects are masked from ‘package:Biobase’:
>
>     content, notes
>
> The following object is masked from ‘package:igraph’:
>
>     tree
>
> >
> > libgit2_features function ()  {
>     .Call(git2r_libgit2_features) } <bytecode: 0x1d276d918> <environment: namespace:git2r>
> > devtools::install_git("https://git.embl.de/velten/STEMNET/", build_vignettes=FALSE) Downloading git repo
> https://git.embl.de/velten/STEMNET/ Installation failed: Error in
> 'git2r_clone': SSL error: syscall failure: Connection reset by peer
>
> >
> > install_git("c97sr/idd") Downloading git repo c97sr/idd Installation failed: Error in 'git2r_clone': unsupported URL protocol
>
> >
> > install_git() Error in lapply(url, git_remote, subdir = subdir, branch = branch, credentials = credentials) :    argument "url" is
> missing, with no default
> >

```

# Answer
> 3 votes
> UPDATE 2018/04/22: Some users had problem with this command, a way to install that works for everyone is to download the tar.gz archive, unpack to a directory of your choice (e.g. /path/to) and inside R run devtools::install\_local("/path/to/stemnet\_dir", build\_vignettes=T)

https://git.embl.de/velten/STEMNET

---
Tags: r, scrnaseq
---


# Question
Title: How to assess the quality of assembled .fasta genome files?
I have assembled 3 `.fasta` files from contigs in`fastq` format of 3 different Homo sapiens.

I would like to see if the assembled genomes which are in `.fasta` format are in good quality or not. In addition, I would be happy to see heritages of genes from the maternal or the parental genomes.

More specifically, I have 3 `.fasta` files of a mother, father and daughter.

1. How can I see they are indeed the biological parents?
2. How can I study what genes were inherited from whom of the parents?
3. How can I assess the quality of the assembled genomes? Is there a standard way to do that?
4. Is there a simple frequency analysis of bases n-grams I can use?

# Answer
> 2 votes
To expand on the comment by Devon Ryan, what you really need are (A.) bam alignment files for each individual and (B.) VCF files for each individual showing all of their genomic variants (with respect to the reference human genome). You mentioned that you had fastq files and that you converted them into fasta files. Do you know if the bam files and VCF files are also available for download anywhere?

If the bam files are not available you will have generate them using sequence aligner software (for example: bwa-mem). With the bam files you can calculate the average depth of coverage (and many other things).

If the VCF files are not available then you will have to generate them from the bam files using variant calling software (for example, the GATK package from the Broad Institute).

In terms of assessing the quality of the sequence runs, a commonly applied tool is name FastQC, but once again, it uses fastq formatted files as input, not fasta files.

The main use for fasta formatted files these days seems to be for compact storage of various reference genomes from different genome sequencing projects (using gzip or similar type of compression).

With the VCF files for each individual in hand you should be able to begin addressing your other questions regarding the frequency of shared variants, etc.

---
Tags: genome, python, assembly, phylogenetics, phylogeny
---


# Question
Title: Where to get '--fldMean' and '--fldSD' for single-end Salmon run
I need to run `single-end` `bulk-RNAseq` `salmon` \- based alignment and I am not sure what I should supply as `--fldMean` and `--fldSD` which seem to be crucial for the `single-end` experiment alignment. Can they be somehow derived from the `FASTQ` files that I have, or should I be provided with them by someone who did the experiment?

# Answer
> 1 votes
From looking at the salmon documentation here, it is correct that both of these are essential for single end experiments:

> Since the empirical fragment length distribution cannot be estimated from the mappings of single-end reads

From the discussion here it would appear that you can't estimate these parameters at all via the FASTQ files alone.

The fragment length refers to the physical molecule of D/RNA captured in the library prep stage which you have sequenced (in part) on your sequencer. Thus ideally this sort of data should come from a bioanalyzer trace from whoever performed the library prep or similar QC step, this should ideally show the mean fragment length and standard deviation.

---
Tags: rna-seq, salmon
---


title: Are gene names same across species?

text:
# Question
Title: Are gene names same across species?
I have a bunch of gene names of *Apis mellifera* (specifically 194). I used these gene names as an input on STRING database to create a network for *Drosophila melanogaster*. 54 of those genes were also present in *D. melanogaster* and I received a second network. Is this a credible method of creating homologous gene network. Or should I use the sequences of the 194 genes and BLAST it to the geneset of the *D*. *melanogaster* and do it that way?

# Answer
You could try using BioMart in Ensembl Metazoa. Filter by your list of gene names in honeybee, then get the fly homologues as attributes. The homologues are calculated by sequence comparison and clustering.
> 4 votes

# Answer
Gene names are not consistent across species. Most of the time genes are automatically annotated and the algorithms try to give names that are similar to their matches in other organisms, but not always. It's also more likely that genes will be automatically named/annotated based on their sequence, and manually named/annotated based on their function.

The further away you get from *Homo sapiens* & *Mus musculus* (and maybe *Caenorhabditis elegans*, *Arabidopsis thaliana* & *Danio rerio*), the less likely it is to find consistently named genes.
> 10 votes

---
Tags: gene, networks, sequence-homology
---


# Question
Title: Changing the colour of each cell in tSNE plot
I have plotted a tSNE plot of my 1643 cells from 9 time points by seurat like below as 9 clusters. But, you know I should not expected each cluster of cells contains only cells from one distinct time point. For instance, cluster 2 includes cells from time point 16, 14 and even few cells from time point 12. My question is, how I can change the colour of cells here so cells from each time point has a distinct colour does not matter some cells are being clustered somewhere else because this is unavoidable. I have named my cells already so I know cells from each time point but I want to colour cells to nine colours based on the time point they have been sequenced not the cluster they are being grouped by seurat.

I first created a file in one column is name of my cells and in second colour the corresponding time point

Like so

```
    time
s0.1    h0
s0.2    h0
s0.3    h0
s0.4    h0
s0.5    h0
s0.6    h0
s0.7    h0
s0.8    h0
s0.9    h0
s0.10   h0
s0.11   h0
s0.12   h0
s0.13   h0
s0.14   h0
s0.15   h0
s0.16   h0
s0.17   h0
s0.18   h0
s0.19   h0
s0.20   h0
s0.21   h0
s0.22   h0
s0.23   h0
s0.24   h0
s0.25   h0

```

I read my file

```
time <- read.table("time.txt",header = TRUE)
seurat@meta.data$time <- time
TSNEPlot(object = seurat, group.by = "time")

```

Sorry in my code colour of tSNE changed also the position of cells changed that I don't want that. I want to change the color of cells but they are still in the same cluster

says that

```
> TSNEPlot(object = al)
Error in seq.default(h[1], h[2], length.out = n) :
  'to' must be a finite number

```

when I changed `names(identity)=sampleInfo$cellName` to `names(identity)=src$cellName`

that worked but again the position of cells changed that I don't want that :(

# Answer
Seurat puts the label in the tSNE plot according to the **@ident** slot of the Seurat object.

By default, if you do the tSNE without computing the clusters and you have the correct metadata in the object, the labels should be pointing to your timepoints not to the clusters.

If it is not the case, what you have to do its to change this slot to reflect to which timepoint belong each cell you are plotting.

For instance if your Seurat object is called src:

```
identity=as.factor(sampleInfo$timepoints)
names(identity)=sampleInfo$cellName
src@ident=identity

```

And then plot the tSNE. Something like this should work after adapting it to your metadata.
> 2 votes

---
Tags: r, scrnaseq, seurat
---


title: Genebank Submission

text:
# Question
Title: Genebank Submission
I am not sure this kind of question is allowed here. I am asking because some of you might have already experienced this.

I submitted some virus genomes in the Genebank. The genomes are on hold and are not public yet.My PI thinks thank we should rename the isolates. I have sent two emails to the admin asking that I want to rename the isolates. I sent an Excel file containing accession number, current isolate name and new name that I want to assign. It has been almost a month and I haven't heard back. So, I don't know if this is the right approach.

What should I do now ? Should I just resubmit the genomes ? Or is there any other way to do that ?

Thanks

# Answer
> 1 votes
I usually hear back from the GenBank curators within a single day.

That said, I’d still suggest you wait: it’s unlikely that they missed both of your emails (but maybe double-check the address and whether they’re likely to be caught by a spam filter): they’re probably just swamped with other submissions.

Furthermore, it’s holiday season so they might operate at drastically decreased capacity.

---
Tags: phylogenetics
---


title: RCurl issue when installing SCnorm

text:
# Question
Title: 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.

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

---
Tags: r, scrnaseq, rcurl
---


# Question
Title: Interpretation of the violin plots from sc-RNA-seq
I'm confused about the meaning of the black dots and the red shape in the violin plots from the seurat tutorial:

# Answer
The black dots represent the values for individual cells. The "nGene" plot (the first one) shows the number of detected genes for every cell.

The red shape shows the distribution of the data. For the "nGene" plot, you can see that the average number of genes per cell is about 900 and most of the cells have roughly around 700-1100 genes.
> 1 votes

# Answer
The plot includes the data points that were used to generate it, with jitter on the x axis so that you can see them better. In red you see the actual violin plot, a vertical (symmetrical) plot of the distribution/density of the black data points.
> 0 votes

---
Tags: scrnaseq
---


# Question
Title: Find a map/correspondence between two versions of a genome
I am working with two versions of the *C. elegans* genome. I am finding interesting regions (specifically, tRNA genes) in version 1 and then I would like to know if version 2 also has a tRNA gene in that region. For this, I need to know the correspondence between coordinates in genome 1 and coordinates in genome 2. I think a whole-genome alignment programme might help me, e. g. nucmer or mauve. I have run nucmer on the two genomes and got a delta file but I don't know how to interpret it. I have run `show-coords` on the delta file and obtained a .coords file. It seems more human-friendly. The first few lines look like this:

```
/Users/user/Documents/project/output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS245.genomic.fa /Users/user/Documents/project/output/genome/celegans/strain/seq/bristol/scaffolds.fa
NUCMER

    [S1]     [E1]  |     [S2]     [E2]  |  [LEN 1]  [LEN 2]  |  [% IDY]  |  [LEN R]  [LEN Q]  |  [COV R]  [COV Q]  | [TAGS]
===============================================================================================================================
       1      432  |      498      929  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I
       1      432  |     1002     1433  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I

```

Does this mean that chromosome I starting at 1 and ending at 432 of version 1 corresponds to chromosome I starting at 498 and ending at 929 and that there is 100% identity between these two 432-long sequences? In this case, it's easy to find a correspondence between bases in version 1 and bases in version 2. However, some lines do not align sequences of the same length:

```
   1   237139  |     4022   241181  |   237139   237160  |    99.99  | 15072434 15371751  |     1.57     1.54  | I  I

```

Also, in this sequence position 1 of chromosome I is aligned again? So is it aligning to more than one position in version 2 of the genome? Also one sequence is 237139 bases long and the other 237160 bases long. How do I know the internal structure?

Can anybody guide me, please?

I've also aligned the two genomes using Mauve:

```
progressiveMauve c_elegans.PRJNA13758.WS245.genomic.fa scaffolds.fa --output scaffolds.mauve > scaffolds.txt

```

I've tried to import the .mauve file into Python using Biopython but I got an error:

```
alignment = SeqIO.read('scaffolds.mauve', "mauve")

Traceback (most recent call last):
  File "code/mauve.py", line 4, in <module>
    alignment = SeqIO.read('scaffolds.mauve', "mauve")
  File "/Users/user/anaconda/lib/python3.6/site-packages/Bio/SeqIO/__init__.py", line 712, in read
    raise ValueError("More than one record found in handle")
ValueError: More than one record found in handle

```

Note that the alignment opened in the Mauve GUI looks fine. Unfortunately, the picture won't load here.

I am working on Mac OSX.

I would appreciate your help in understanding the delta file, resolving the Python error or any other support in finding a function between coordinates in version 1 of the genome and coordinates in version 2 of the genome.

# Answer
A relatively convenient method would be to create a .chain file, as used by tools like `liftOver` and `crossmap` (and probably remap from NCBI). For your use-case, you can simply use the DoSameSpeciesLiftOver.pl command as described in the link. The resulting chain files are pretty self-explanatory and can be used to map BED files of regions of interest.
> 1 votes

---
Tags: sequence-alignment, python, mauve
---


# Question
Title: Filtering the expression values of a gene from GEO
I'm filtering the gene expression value of a particular gene using Bioconductor in R. Following the instruction given here, I could filter the gene expression value for a gene from the dataset GDS3782.

The output is as follows for the gene of my interest,

```
GSM524151  "23.82"
GSM524152  "27.49"
GSM524153  "24.26"
GSM524154  "17.98"
GSM524155  "25.87"
GSM524156  "26.16"
GSM524157  "26.02"
GSM524158  "25"
GSM524159  "29.85"
GSM524160  "18.25"
GSM524161  "28.86"
GSM524162  "27.33"
GSM524163  "27.42"
GSM524164  "28.23"
GSM524165  "22.09"
GSM524166  "24.85"
GSM524167  "27.11"
GSM524168  "32.01"
GSM524169  "28.97"
GSM524170  "30.98"

```

The above data corresponds to disease and control which I get using

```
Meta(gds3782)$description

```

I wish to see the description of each value. For example, from the webpage, I have seen that the first ten values correspond to "disease-state" and the last ten correspond to "control". How do I view this classification while processing the data using Bioconductor in R?

# Answer
You're looking for the `Columns()` function:

```
> Columns(gds3782)
     sample        disease.state           description
1  GSM524151 non-diabetic control  Value for GSM524151: beta-cells_non-diabetic condition_donor1; src: human beta-cells_non-diabetic
2  GSM524152 non-diabetic control  Value for GSM524152: beta-cells_non-diabetic condition_donor2; src: human beta-cells_non-diabetic

```

Note that it's probably easier to work with a more common structure like and expressionSet:

```
eset = GDS2eSet(gds3782)

```

You'd then use the `pData()` function to return the same thing.
> 3 votes

---
Tags: bioconductor
---


title: SC-RNA seq percent.mito

text:
# Question
Title: SC-RNA seq percent.mito
Is that possible that the percent.mito level is zero? I have done vlnplot and it shows that its percent.mito level is zero. Does that mean my mitochondrial gene didn't detect?

It is violin plot. And the command is:

# Answer
I've just come across a similar issue. It seems the mouse and human mitochondrial genes use mixed capitalizations. See the following PR for details:

https://github.com/Hoohm/dropSeqPipe/pull/43

Bottom line is, if you ignore the case in the grep pattern it might work (assuming you use human/mouse):

```
mito.genes <- grep("^mt-", rownames(shox2crepos_mtmg@data), value=TRUE, ignore.case=TRUE)

```
> 1 votes

# Answer
The simplest answer is that the bit of code you copied from the Seurat vignette is not working on the gene names you have. Look at the code. If none of your gene names begin with "MT-", none will be identified as mitochondrial.
> 1 votes

---
Tags: scrnaseq
---


title: Comparing the gene expression data

text:
# Question
Title: Comparing the gene expression data
I'm trying to compare the gene expression data reported in the studies documented in GEO for a specific gene expressed in a tissue of my interest in Homo sapiens. I compared the values reported in 4 studies, in two studies the expression value of gene A lies between 15-30, in the other two studies the range lies between 1-4.

I would like to ask for advice on whether it is meaningful to compare the data from different studies. What are the factors that one should keep in mind while making such comparisons?

Any help would be highly appreciated

# Answer
There's going to be a batch effect between studies, so unless you have some sort of constant group that you can use for between-study normalization you're going to completely muck up your results from using data from multiple studies.

More specifically, this is what packages like `sva` and the `combat()` command in it are used to correct. Note that you should go back to the raw .CEL or .fastq files if you want to process multiple studies, in order to at least ensure that everything is processed in the same way.
> 1 votes

# Answer
It doesn't make sense to compare gene expression data from different studies unless you have a negative spiked-in control such as sequins.

The most obvious thing you will need to do is normlizate the data for various bias such as batch effects, sequencing platforms, sequencing depth, bioinformatic workflows and more ... Your task will be quite **complicated** as you must justify how your normalization would work for data from different studies.

You should check how the raw data was generated (methods section), convince yourself, and then apply gene quantification yourself from the raw files.
> 1 votes

---
Tags: gene, gene-expression, geoquery
---


# Question
Title: RNA-Seq: clustering/treatment of genes with low expression
I have some RNA-Seq data from leukaemia patients. I want to do unsupervised clustering on them with some other published leukaemia RNA-Seq data and see how they cluster. There are a few problems I encountered while doing this.

1. I read mix messages of whether using log2(TPM) or rlog(counts) (e.g. Deseq2 rlog or limma voom transformation) for clustering. Which one should I choose?
2. If I am filtering out genes with low counts, should I do it prior normalisation to library size or after?
3. I tried using filterByExpr from edgeR for filtering but it removed many of the genes that are only expressed by subtypes with small sample number. If I am to filter the counts "manually", is there a recommendation of how the threshold should be set?
4. If I should use log2(TPM) for unsupervised clustering, can I treat low TPM (e.g. \<2) as 0? As this paper suggested TPM \< 2 basically means no expression, is it alright to treat them as 0?

# Answer
1. `rlog(normalized counts)` is going to be more robust than `log2(TPM)`, so use it instead.
2. Do it afterward, keeping the low counts can be helpful for library-size normalization.
3. This is very much more of an art than a science, though have a look at the `kOverA()` function from the `genefilter` package. That will give you a bit more fine-grained control on filtering.
4. You'd be better off using something like zFPKMs if you want to filter by "expressed or not", since there's no really great definition of that.
> 3 votes

---
Tags: rna-seq, clustering
---


title: how to filter out cell after doing QC

text:
# Question
Title: how to filter out cell after doing QC
I was wondering that the threshold that we set after doing QC in Seurat. there is an example

```
# GenePlot is typically used to visualize gene-gene relationships, but can be used for anything
# calculated by the object, i.e. columns in object@meta.data, PC scores etc.
# Since there is a rare subset of cells with an outlier level of high mitochondrial percentage
# and also low UMI content, we filter these as well
par(mfrow = c(1, 2))
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "percent.mito")
GenePlot(object = pbmc, gene1 = "nUMI", gene2 = "nGene")`

```

Then its filter cell command is under:

```
# We filter out cells that have unique gene counts over 2,500 or less than 200
# Note that low.thresholds and high.thresholds are used to define a 'gate'
# -Inf and Inf should be used if you don't want a lower or upper threshold.
pbmc <- FilterCells(object = pbmc, subset.names = c("nGene", "percent.mito"), low.thresholds = c(200, -Inf), high.thresholds = c(2500, 0.05))`

```

So I want to know whether the threshold setting is depending on the picture I show?

And the below is the Geneplot I have got base on my data. How can I set my thresholds?

# Answer
Yes, you end up placing filters where things become inconsistent or similar to background. In the example, which appears to be from the Seurat tutorial, there's a downward plunge in `nGenes vs. nUMI` at around 200 genes and few believable cells with more than 2000-2500 genes, so the tutorials sets filters around there. The 0.05 mitochondrial contamination is a more arbitrary round number and is less supported by the plots (one would use more like 0.06 or 0.07 if one based this solely on the plots).
> 1 votes

---
Tags: scrnaseq
---


# Question
Title: Minimap2 -ax map-pb doesn't output tlen field
I have used minimap2 to map some pacbio reads to a reference genome. I would like to know the "insert size" (true length of sequence) relative to the reference. More specifically, I want to know the first base position in the reference where the read maps and the last base position in the reference where the mapping ends.

As I understand it I need the TLEN field from bam files (sam format field 9). However, all the alignment files from minimap2 have `0` in the TLEN field.

Does anyone know how to get pacbio alignments with the TLEN field? I am also happy using `bwa` or another aligner.

The command I used was `minimap2 -ax map-pb {fasta} {reads} | samtools view -h -b - | samtools sort - > myfile.bam`.

# Answer
The TLEN for single-end reads is defined as 0 in the SAM specification, which is why this happens. You can still get the end position of the alignment by parsing the CIGAR string. An example using the pysam package in python is below:

```
import pysam

bam = pysam.AlignmentFile("myfile.bam")
for b in bam.fetch():
    print("This alignment spans {}:{}-{}".format(b.reference_name, b.reference_start+1, b.reference_end))

```

The `reference_end` value is automatically calculated from the CIGAR string.
> 2 votes

---
Tags: sequence-alignment, phylogenetics, bwa, minimap2
---


title: Visualize PLINK pedigree files

text:
# Question
Title: Visualize PLINK pedigree files
Is there a pedigree viewer that works with PLINK pedigree files? I saw an old post by Pierre, but I don't see the script on his blog.

# Answer
We use `kinship2` R package to visualize pedigree files (.fam plink files). For some reasons, some characters from plink ped aren't supported (like "-9" for unknown genre or affected status).
Here is the R script we use:

```
args<-commandArgs(TRUE)
pedPath<-args[1]
pdfFile<-args[2]

library(kinship2)
pedfile <- read.table(pedPath,sep="\t")
pedfile$V5 <- apply(pedfile,1,function(data){if(data[5] == 1 || data[5] ==2){data[5]}else{3}})
pedfile$affected <- apply(pedfile,1,function(data){if(as.numeric(data[6]) == 1){res <- data[6]}else{if(as.numeric(data[6]) == 2){res <- data[6]}else{res <- NA}}})

familyIDs <- as.character(unique(pedfile$V1))

pdf(file=pdfFile)
for(famID in familyIDs){
    ped <- pedigree(id=as.character(pedfile$V2[pedfile$V1==famID]),dadid=as.character(pedfile$V3[pedfile$V1==famID]),momid=as.character(pedfile$V4[pedfile$V1==famID]),sex=as.numeric(pedfile$V5[pedfile$V1==famID]),affected=as.numeric(pedfile$affected[pedfile$V1==famID]))
    plot.pedigree(ped,cex=0.7,align=TRUE)
    title(famID)
}

```

And we call it this way:

```
Rscript <kinshipScript.R> <plink ped file> <pdf output>

```
> 3 votes

---
Tags: plink, pedigree
---


# Question
Title: DiffBind: Model must include count data for contrasts.
I am trying to use `DiffBind` for my `ATACseq` experiment going through the `vignette` for `DiffBind`. I was able to import the data and make a heatmap with `plot(data)`. But somehow calling `data <- dba.contrast(data, categories=DBA_TISSUE)` function is not working for me giving a warning message:

> Warning message: Model must include count data for contrasts.

I tried running beforehand: `data <- dba.count(data)` but it is giving an error:

> Error: Error processing one or more read files. Check warnings(). In addition: Warning messages: 1:
> 2:
> 3:
> 4:
> 5:
> 6:
> 7:

My `sample_sheet.csv`looks like that:

I just need to compare `Dura` vs `Microglia` which is `Tissue` column. How could this error be fixed? Any suggestions would be greatly appreciated.

# Answer
The problem was that I was using wrong files. I ran a pipeline from here and it generates several folders with many files. The `.bam` files that should be used are located in the `aligned_mm10` folder and are with the extension: `.pe.q10.sort.rmdup.bam`.
> 0 votes

---
Tags: atac-seq
---


# Question
Title: Do people today with a disease caused by single gene mutation that is inheritable - all share a common ancestor?
First question here so sorry if it is about both genetics and evolution.

So, if I get a disease which can be clearly identified on the gene level and all my descendants inherit it, after 500 years do all people with the disease share the common ancestor (given that it's only inheritable)?

I mean, my clear guess first would be "yes" but I really want to check it since I haven't read it anywhere.

Feel free to post materials on this too. Thanks!

# Answer
That can happen, but it's more likely that the disease will be present in more than one family. Diseases with a single, small genetic cause will be more likely to flare up elsewhere (assuming the disease isn't *too* bad, and it wouldn't be if it's been in the family for 500 years).

However, what you will find is that in *your* family that disease is almost always present together with some other inherited genetic feature that is nearby on the genome (i.e. an inherited haplotype block). This association is less likely to be shared with other families.
> 2 votes

---
Tags: phylogenetics
---


# Question
Title: How to map reads shorter than 32bp with minimap2?
I mapped fasta sequences to a reference. The fasta sequences range in length from 17bp up to several hundred bp. I used `minimap2` with the following command: `minimap2 -ax sr ref.fa seqs.fa`, et cetera. Only the sequences greater than 32bp len mapped though.

```
samtools view -F 4 filename.sorted.bam | \
  cut -f10 | awk '{print length($1)}' | \
  sort | uniq -c | sort -k2 -n

   2304 32
  34931 33
   5443 34
   5343 35
   8811 36
   5485 37
   5417 38

```

Is it possible to make minimap2 map sequences shorter than 32bp? Do I need to use `bwa aln` or some other mapper?

# Answer
> 2 votes
Minimap2 is not the most optimal aligner in this case, it was designed for long reads, of which 32bp reads are clearly not. You could try lowering the mnimizer k-mer length and with `-k` the default being 15, also try lowering the minimzer window size with `-w`. However I suggest using `bwa aln` and `bwa samse` for all reads under 32bp and then merging the alignment later. Although I struggle to see what the utility of these short reads will be against a backdrop of longer reads?

# Answer
> 3 votes
Minimap2 doesn't work well with 32bp reads by design. Older mappers like bowtie1 and bwa-aln will be better.

---
Tags: fasta, read-mapping, phylogenetics, bwa, minimap2
---


# Question
Title: Efficiently aligning a lot of reads on the same small reference sequence
The context: I have a DNA-sequence coding for a protein, about 1500 bp in length. Using NGS, a lot of reads of (mutants of) this same sequence were acquired. All of these reads need to be aligned to the reference. We're talking about a lot of reads (100,000 - 1,000,000) of not too short length (120 - 300 bp). All of these reads can belong to different mutants of the template sequence, so the alignment is necessary to determine the exact sequence of every single mutant.

Currently, I'm just using Smith-Waterman-based local alignment to align every single read to the reference one by one. Yet I can't help but feel there might be a more computationally/time-efficient solution to this specific problem.

Maybe there exists an algorithm that's not very efficient for most alignment problems, but that becomes very worthwile if it has to map a ton of reads to the same place over and over again. For example, it might do some time-consuming operations on the short template that make it fast to align reads to it, but these operations aren't worth it if only aligning a couple of reads. That's just an idea, I don't know all the different techniques that are out there.

**So, to recap**: I want to align a lot (100,000 - 1,000,000) of 120-300 bp reads to the same short 1500 bp reference sequence. If anyone has any suggestions about an algorithm, or just a specific workflow, that is particularly suited to do this, it would be appreciated. I work in R so I can implement some things myself, it doesn't have to be a ready-to-use software package or anything like that. Thanks in advance!

# Answer
> 4 votes
I would suggest using a whole genome aligner like BWA or BOWTIE2 and swapping out the whole genome reference with your specific reference.

These alignment algorithms are designed to align large numbers of NGS reads against a single reference. They can be used with any reference including short ones like your example or small viruses.

To ensure that local alignment is correct, keep in mind that `bowtie2` needs the `--local` while `bwa` should be run with on the `BWA-MEM` modality

---
Tags: sequence-alignment, ngs, algorithms
---


# Question
Title: Calculating read average length in a Fastq file with bioawk/awk
I found here this awk script:

```
BEGIN {
  headertype="";
  }
{
  if($0 ~ "^@") {
    countread++;
    headertype="@";
  }
  else if($0 ~ "^+") {
    headertype="+";
  }
  else if(headertype="@") { # This is a nuc sequence
    len=length($0);
    if (len>4) {
      readlength[len]++;
    }
  }
}
END {
  for (i in readlength){
    countstored+=readlength[i];
    lensum+=readlength[i]*i;
    print i, readlength[i];
  }
  print "reads read = "countread > "/dev/stderr";
  print "reads stored = "countstored > "/dev/stderr";
  print "read average length = "lensum/countstored > "/dev/stderr";
}

```

and I just wonder if it is possible to shorten it with bioawk?

# Answer
This can also be done with regular awk.

```
awk '{if(NR%4==2) {count++; bases += length} } END{print bases/count}' <fastq_file>

```

The `NR%4==2` count the second line out of every block of 4. `length` is a built-in that defaults to the length of the line, same as `length($0)`. In this case you can inject you custom printing to the `END{}` block but `countread` and `countstored` will always be the same since the averaging is done on the fly.
> 8 votes

# Answer
This script is wrong because a quality string may start with `@`. With bioawk, it can be simplified to:

```
bioawk -c fastx '{ readlength[length($seq)]++; countread++ } END{...}'

```

The `END{}` block is the same as your original version.

EDIT: forgot `-c fastx` is the original answer. Thank @MattBashton for pointing this out.
> 3 votes

---
Tags: fastq, awk, bioawk
---


title: fuse fastq files with multiple records

text:
# Question
Title: fuse fastq files with multiple records
I'd like to fuse ~50 pairs of corresponding fastq files *on per sequence bases*. I'd be happy for a solution as a Linux script, R script (if feasible with huge files), or any free software that runs preferably under Debian Linux (otherwise Windows). Aim is preprocessing for a Linux pipeline for NGS data. Ideally both the sequence and quality data are fused.

The **aim** of this request is the following: there is an existing pipeline for downstream mapping of immunoreceptor sequences (migec/mixcr/vdjtools etc.). Previously we got Illumina reads (read 1, read 2 = reverse), with the index and UMI sequences attached to one of the reads =\> can be alligned and processed to habe 1 biological useful sequence with UMI and index at one end. NOW we get read 1, read 2 and seperately "read 3" which is just the UMI+Index reads, and which was attached by the core facility to read 1 using software we do not have access to. Since the previous pipeline works well, we prefere to use it in the future, but are open for suggestions how to use the index reads with migec/mixcr/vdjtools otherwise.

Each of the 2 files contains the same number (millions) of short (~300 bp) and corresponding DNA reads:

File\_A.fastq:

* @record1 / aaa... / +... / ...
* @record2 / ...
* @record3 / ...
* ... (millions)

File\_B.fastq:

* @record1 / ttt... / +... / ...
* @record2 / ...
* @record3 / ...
* ... (millions)

The target file should look like:

* @record1 / aaa...ttt... / +... / ...
* @record2 / ...
* @record3 / ...
* ... (millions)

(of course the aaa...ttt... is just an example)

I was looking for a solution on Stackexchange and Google, but the solutions, I found do not work for me because

* manual processing is not an option given the number of sequences
* I need to fuse the files on a **per sequence** basis, i.e. a simple Linux command like cat file1 file2 is not what I need (as far as I understand, this would give me (record1\_A record2\_A.... record1B record2\_B).
* R, package shortread, readFastq + writeFastq would be perfect, but cannot handle such large numbers of input lines

Maybe I did not use the right search terms, but in any case I'd be happy for a solution

# Answer
Based on the comments of the question I thought to put together this answer using GNU coreutils which should exist on pretty much all systems.

The idea is to paste the 2 files together with `paste` (assuming they are in the same order and each line correspond to the same line in the other file).

Since only the sequence and quality of both files should be kept, the lines starting with `@` and `+` need to be made empty for one of the files. This can be achived with sed, e.g. `sed 's/^@.*//'` empties all lines starting with `@` but keeps the empty line.

To modify the second file on the fly, one can using process substitution `<()`. Everything strung together should be something like this:

```
file1=File_A.fastq
file2=File_B.fastq
paste -d "" $file1 <(cat $file2 | sed 's/^@.*//' | sed 's/^+.*//') > target.fastq

```

I haven't tested it thoroughly, but could be a good starting point.

If you have a lot of files to paste, you might however want to process them first (not on the fly) using:

```
for file in $(ls *.fastq); do
  cat $file | sed 's/^@.*//' | sed 's/^+.*//' > ${file%.fastq}_cleaned.fastq
done

```

This creates `File_A_mod.fastq` from `File_A.fastq` etc. Best to move them in a separate directory (make sure the first one is the original instead of the \_cleaned one since it serves as the base). You should have the following in the folder: `File_A.fastq File_B_cleaned.fastq File_C_cleaned.fastq ...`. Those can be than pasted with `paste -d "" *.fastq > target.fastq`.
> 1 votes

# Answer
What you actually want to do is overlap the reads, for which you can use tools like FLASH or BBMerge (I've had better luck with FLASH).
> 0 votes

---
Tags: fastq
---


# Question
Title: Why is cell count stated relative to animal size?
Why do some single-cell RNA-seq papers give the number of cells sampled *as a fraction of the size of the organ, animal, or embryo*, as opposed to just saying how many cells they sampled?

For instance, Wagner et al.'s paper "Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo" says "For different developmental stages, we sampled 0.17x to 0.97x of the total cells per embryo, sufficient to detect cell states as rare as 0.1-0.5% of all cells." A pollster would never write "We sampled 4 out of every 10,000 registered voters in the USA"; they'd just give the number. Statistical results on cluster detection (e.g. Disentangling Gaussians by Adam Tauman Kalai, Ankur Moitra, and Gregory Valiant) are also phrased in terms of number of samples, not this ratio of sample size to organism size.

It's worth mentioning that these studies use low-yield techniques. In order to sequence 0.97 times the number of cells in a zebrafish embryo, Wagner et al. had to sequence at least 50-100 embryos. So, it's more like sampling 97 of 100 cells with replacement and less like sampling 97 of 100 cells without replacement.

# Answer
The typical goal of scRNA-seq is to assess the possible cell-types and states that exist in an organism/tissue/etc. Given that, it can be useful to know what portion of the possible cells have been assayed, since that sets bounds on the number and size of cell populations **not** assayed. If the goal of a study is to assess the variability in a single cell-type, then you're correct that mentioning "0.8% of the population was sequenced" isn't actually informative. Note that "in a single cell-type" is important there, since if one is assaying multiple cell-types then it's important that enough of each are assayed given their prevalence in the relevant population.

Practically speaking, numbers like this are only ever written in papers if they look impressive, so no one is going to write about sequencing 0.000001% of the cells in an organism, since the reader/reviewer will then start to question the relevance of the results.
> 1 votes

---
Tags: single-cell
---


# Question
Title: At what stage of a transcriptome assembly is it better to perform read contaminant filter?
I'm trying to assemble a bivalve transcriptome. Since bivalves are filter feeders, their transcriptomes tend to be highly contaminated by bacteria, algae and whatnot. Since I pooled several transcriptomes, I have a high amount of reads (\>2B reads). I also have a not-so-high quality reference genome.

Now, I'm not sure when to perform a BLAST-based filter on these reads:

* Do I filter reads with BLAST before assembling them?
* Do I assemble the transcriptome and then filter out contaminant contigs?
* Since I have a reference genome, do I simply not care about contaminant reads because they will not align to it?

Having this amount of reads, I'd be inclined to first filter out those reads that don't have a match against Metazoa, and then to perform a second filter after assembly in order to retain only those contigs that match against Lophotrochozoa.

If you know resources on how to manage read contamination, feel free to share them.

# Answer
> 3 votes
First, to answer your question about mapping to a low-quality reference:

### 1. Mapping

For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome assembly does not contain any scaffolds from contaminating species like bacteria, algae, et cetera.

### 2. Reasons to not BLAST your raw reads

There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

* BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
* Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
* Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in **#4**, then reassembling, may work better.

### 3. Using multiple individuals of highly heterozygous species in assembly.

Most importantly for non-model organisms and highly heterozygous species - **do not assemble a transcriptome using data from more than one individual**. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.

### 4. A strategy to filter transcripts based on mapping and BLAST

Given **#1**, **#2**, and **#3** above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.

I would try this:

* Assemble the transcriptome of each individual separately.
    + This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.
* Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.
    + Use somthing like `cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta`.
* Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
* Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.
    + Use `samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam` to only keep the transcripts that **did not** mapped to the reference. `-f 4` is the magic flag here that keeps things that didn't map.
* Make a list of all the transcripts that **did not** map to the reference...
    + `samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt`
* ...and get only the transcripts from your fasta file that are potential contaminants.
    + `seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta`
* BLAST all of the `unmapped.fasta` and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
* Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this `known_contaminants.fasta`
* Map **all reads** from **all individuals** to the `known_contaminants.fasta`.
* Get all the reads that map to known contaminating transcripts, `known_contaminants.fasta` using:
    + `samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt`
* Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
* Reassemble the filtered fastq files.

This strategy has the following benefits:

* Does not require weeks/months of compute time to BLAST 2 billion reads.
* Ensures that all reads from the reference genome are left untouched.
* Ensures that all reads from known contaminants are removed.
* Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
* Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
* Leverages contaminants from different individuals to filter reads from all individuals.

---
Tags: transcriptome, reference-genome, read-correction
---


# Question
Title: Why DoHeatmap Does not show all genes in genes.use?
I am heatmaping a list of genes by DoHeatmap function in Seurat R package. I am sure I have 212 genes but heat map shows only a few of my genes

```
> DoHeatmap(
+        object = seurat,
+        genes.use = genes,
+        slim.col.label = TRUE,
+        remove.key = TRUE)
> dim(as.matrix(seurat@data))
[1] 12293   209

> length(genes)
[1] 212
> class(genes)
[1] "character"

> class(seurat)
[1] "seurat"
attr(,"package")
[1] "Seurat"
>

```

Please look at heat map, I was trying to plot 212 genes but you seeing only few genes are among my 212 genes

Please somebody save me from this confusion

# Answer
https://github.com/satijalab/seurat/issues/287

In this post there is a answer for my question, however when I set scale to FALSE again the number of genes did not change but the shape of plot dramatically changed. So, maybe the issue is because of scaling that did not solve my problem
> 1 votes

---
Tags: r, scrnaseq, seurat
---


# Question
Title: Get results of keyword search on Pfam via python script
I'm interested in all proteins that are in any way associated with *Danio rerio*. I decided to look them up at Pfam data base and when I just make a keyword search, I get a a nice list which looks like this http://pfam.xfam.org/search/keyword?query=danio+rerio&submit=Submit. Important thing to note is that most of the associations of these protein families can be found in **seq\_info** column.

My question is: how to extract that list from a python script? Moreover, I want to make a script that when run, will fetch the latest release, filter out the proteins that cannot be found in *D. rerio* and do some additional processing.

I've been searching for ways to directly download files from Pfam ftp server and mine for the proteins that are associated with *D. rerio*, however, I didn't find a single way to do it. I simply cannot see where do they store the information about species in which the protein can be found.

I would like to retrieve everything about protein families of interest, because the next step I will need to do is to find the DNA binding proteins and for finding a way to filter out the ones of no interest, I would liketo have as much data as possible.

# Answer
If I understand your question correctly, then you can use the Pfam proteom search to get the data you are looking for.

Go to the

> Pfam website -\> Browse -\> select the letter D from proteoms category -\> D. rerio

There you can find a download link under the section `Domain composition`:

> You can download a list of all regions for this proteome from our FTP site.

The download link seems to be bound to the species ID and if you want to automate your data retrieval you may be able to simply download the file from there in your python script.
> 1 votes

---
Tags: python, public-databases
---


# Question
Title: Are these standard species abbreviations and how to look up others?
Reactome uses a three letter code for each species which are listed below.

Are the letter codes used by Reactome a standard?

If so is there a web page or download file that can validate the code and translate it to a common name?

```
 HSA  Homo sapiens
 SSC  Sus scrofa
 GGA  Gallus gallus
 DME  Drosophila melanogaster
 RNO  Rattus norvegicus
 DDI  Dictyostelium discoideum
 MMU  Mus musculus
 ATH  Arabidopsis thaliana
 TGU  Taeniopygia guttata
 DRE  Danio rerio
 OSA  Oryza sativa
 BTA  Bos taurus
 XTR  Xenopus tropicalis
 CFA  Canis familiaris
 SCE  Saccharomyces cerevisiae
 PFA  Plasmodium falciparum
 CEL  Caenorhabditis elegans
 SPO  Schizosaccharomyces pombe
 ROT  Rotavirus
 FLU  Influenza A virus
 HCV  Hepatitis C Virus
 HBV  Hepatitis B virus
 HER  Human herpesvirus 1
 NME  Neisseria meningitidis serogroup B
 NME  Neisseria meningitidis
 CTR  Chlamydia trachomatis
 STY  Salmonella typhimurium
 STY  Salmonella enterica
 ECO  Escherichia coli
 OCU  Oryctolagus cuniculus
 XLA  Xenopus laevis
 CAE  Cercopithecus aethiops
 CGR  Cricetulus griseus
 HER  Human herpesvirus 8
 MCV  Molluscum contagiosum virus subtype 1
 MCV  Molluscum contagiosum virus
 HCY  Human cytomegalovirus
 HIV  Human immunodeficiency virus 1
 PCH  Penicillium chrysogenum
 TAE  Triticum aestivum
 VRA  Vigna radiata var. radiata
 VRA  Vigna radiata
 CFS  Crithidia fasciculata
 FCA  Felis catus
 AMA  Arenicola marina
 MGA  Meleagris gallopavo
 NGO  Neisseria gonorrhoeae
 SAU  Staphylococcus aureus
 MVI  Measles virus
 HAL  Human alphaherpesvirus 2

```

---

EDIT

Not a complete answer as this only covers a subset of all organisms but at present better than nothing.

Using KEGG Organisms

https://www.genome.jp/kegg-bin/get\_htext

Enter letter code, e.g. GGA

---

EDIT - After reading answer by conchoecia

The graph database schema for Reactome only has two labels that use the 3 letter code which is named abbreviation in the schema:

Species and ReferenceTherapeutic

Species has taxid (Taxon Identifier)

```
abbreviation     taxId       displayName
AMA              6344        Arenicola marina
ATH              3702        Arabidopsis thaliana
BAN              1392        Bacillus anthracis
BTA              9913        Bos taurus
CAE              9534        Cercopithecus aethiops
CAL              5476        Candida albicans
CBO              1491        Clostridium botulinum
CEL              6239        Caenorhabditis elegans
CFA              9615        Canis familiaris
CFS              5656        Crithidia fasciculata
CGR              10029       Cricetulus griseus
CPH              10703       Corynephage beta
CPO              10141       Cavia porcellus
CSB              60711       Chlorocebus sabaeus
CTE              1513        Clostridium tetani
CTR              813         Chlamydia trachomatis
CVI              10243       Cowpox virus
DDI              44689       Dictyostelium discoideum
DME              7227        Drosophila melanogaster
DRE              7955        Danio rerio
ECO              562         Escherichia coli
FCA              9685        Felis catus
FLU              11320       Influenza A virus
GGA              9031        Gallus gallus
HAL              10310       Human alphaherpesvirus 2
HAM              6706        Homarus americanus
HBV              10407       Hepatitis B virus
HCV              11103       Hepatitis C Virus
HCV              31649       Hepatitis C virus genotype 2a
HCV              31646       Hepatitis C virus subtype 1a
HCY              10359       Human cytomegalovirus
HER              10298       Human herpesvirus 1
HER              37296       Human herpesvirus 8
HGA              10376       Human gammaherpesvirus 4
HIV              11676       Human immunodeficiency virus 1
HPV              333760      Human papillomavirus type 16
HPV              337041      Alphapapillomavirus 9
HSA              9606        Homo sapiens
IBR              11120       Infectious bronchitis virus
LMO              1639        Listeria monocytogenes
LMO              1906951     Listeria monocytogenes serovar 1/2a
LPN              446         Legionella pneumophila
MCV              10280       Molluscum contagiosum virus subtype 1
MCV              10279       Molluscum contagiosum virus
MGA              9013        Meleagris gallopavo
MML              9544        Macaca mulatta
MMU              10090       Mus musculus
MTU              1773        Mycobacterium tuberculosis
MTU              83332       Mycobacterium tuberculosis H37Rv
MVI              11234       Measles virus
NGO              485         Neisseria gonorrhoeae
NME              491         Neisseria meningitidis serogroup B
NME              487         Neisseria meningitidis
OAR              9940        Ovis aries
OCU              9986        Oryctolagus cuniculus
OSA              4530        Oryza sativa
PCH              5076        Penicillium chrysogenum
PFA              5833        Plasmodium falciparum
RNO              10116       Rattus norvegicus
ROT              10912       Rotavirus
SAU              1280        Staphylococcus aureus
SCE              4932        Saccharomyces cerevisiae
SEV              11191       Sendai virus
SPO              4896        Schizosaccharomyces pombe
SSC              9823        Sus scrofa
STY              90371       Salmonella typhimurium
STY              28901       Salmonella enterica
TAE              4565        Triticum aestivum
TGO              5811        Toxoplasma gondii
TGU              59729       Taeniopygia guttata
VRA              3916        Vigna radiata var. radiata
VRA              157791      Vigna radiata
VVI              10245       Vaccinia virus
XLA              8355        Xenopus laevis
XTR              8364        Xenopus tropicalis

```

While ReferenceTherapeutic does not.

The species can be identified using the NCBI taxon identifier, e.g.

https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=6344

---

EDIT - After reading answer by user172818

Searching the file https://www.uniprot.org/docs/speclist.txt for 6344

```
AREMA E    6344: N=Arenicola marina
                 C=Lugworm
                 S=Lumbricus marinus

```

or using the URL: https://www.uniprot.org/taxonomy/6344

and with reference to

> I think databases should stop inventing new species abbreviations.

The Wikipedia article on Lugworm shows the many Taxon identifiers in multiple databases.

# Answer
No, there are no standardized three-letter abbreviations for species names. Three-letter abbreviations with roman letters is 26^3 combinations, which is only 17,576 possible combinations. As a result there are not enough three-letter abbreviations for all living species.

All three-letter abbreviations are either specific to a database or perhaps are shared between a database or two. There is currently no substitute for a binomial name (although NCBI taxid is a convenient parsable substitute for the species included in NCBI to date).
> 7 votes

# Answer
These are not. The closest to "standard" is the 5-character abbreviation by Swissprot. It names 25,886 species and has been used for decades. It is also easy to remember for some common species such as HUMAN, RAT, MOUSE and HORSE. I think databases should stop inventing new species abbreviations.
> 6 votes

---
Tags: database, data-mining, taxonomy
---


# Question
Title: How many reads do I need for hybrid assembly
I have Illumina and PacBio reads and I would like to use dbg2olc for hybrid assembly.

Part of dbg2olc is `SelectLongestReads` which select reads that sum to total\_length of bases.

How many bases do I need for a 1.8G allotetraploid to remain?

# Answer
> 1 votes
I have never used dbg2olc. However other OLC assemblers managed to do a solid assembly from 20x of corrected long reads. But obviously, the more the better.

In you case of allotetraploid, you need to figure out if the heterozygosity is sufficient to separate haplotypes of not and adjust parameters of the assembler. For instance in bdelloid rotifers the tetraploid genome has an AA'BB' structure, all the assemblies were able to separate A and B, but not A and A' and B and B'. I guess an allotetraploid could have a similar structure.

It would be good to know what you are dealing with and what is your goal in the assembly. Do you want a single haplotype?

-- edit --

You can check the genome structure using short reads by smudgeplot, a tool I made (yes, it's a shameless plug).

---
Tags: assembly, illumina, pacbio, coverage
---


# Question
Title: A strange error while clustering of single cell RNA-seq using URD package
I am using URD R package on my single cell RNA-seq data hoping to be able to trace gene expression changes during the development. After dimension reduction step by PCA and diffusion map whatever I am trying to cluster some cells, but whatever I do I am getting this error

```
> object <- graphClustering(object, num.nn = c(5,8,10,15), method="Louvain", do.jaccard = TRUE)
Error in graphClustering(object, num.nn = c(5, 8, 10, 15), method = "Louvain",  :
  unused arguments (num.nn = c(5, 8, 10, 15), method = "Louvain", do.jaccard = TRUE)

> print(object.6s)
URD object: 12451 genes x 1556 cells.

object.6s.mnn <- graphClustering(object.6s)
Show Traceback Rerun with Debug Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘graphClustering’ for signature ‘"URD"’

```

I tried to change Louvain, jaccard, ... but I am still getting this error. How can I get rid of it?

# Answer
This solved the error

```
object <- URD::graphClustering(object)

```

***By this function I am telling R to use the graphClustering() function from the URD namespace. I probably had loaded another package with a different graphClustering() function - or defined one myself - which is taking precedence over URD's function.***
> -1 votes

---
Tags: r, scrnaseq
---


# Question
Title: Counting letters in phylip alignment columns with Biopython
I have been using python 3.6 and biopython 1.72 to work with protein data files. I am using a protein sequence file (phylip format), for example:

```
    14    678
Zebrafish  LSSCGVVSGD LISVILPASS LEETQTSSAA AHQTHTDQQA GGSHVSSSSS
Fugu       LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS-
Chicken    LASCGIVSGD LIRVALIRAA LMAADAPDRA DDGGGHSEQA VAMATNQVS-
Zebra      LASYGIVPGD LICLLLEETD AKSSLPPPSS SSPPSLQNGR ETSTLTPKNQ

```

I want to manipulate the file in a way so that;

1. I may get the number of columns which have all the same letters (amino acids)
2. The number of columns which have at least 2 different amino acids (with no dashes).
3. The number of columns which have any single dash.

I have iterated over the columns of the file with:

```
alignment = AlignIO.read(open("example.phy"), "phylip")
def column(alignment):
    n = len(alignment[0].seq)

    for i in range(n):
        c = []

        for record in alignment:
            c.append(record.seq[i])

        yield c

```

This returns the columns of the file, but I don't understand how I can manipulate them further.

Can I get some helpful suggestions?

# Answer
I suggest you use `zip()` to transform your data into columns and then apply the following counts similar to heathobrien suggests:

```
from Bio import AlignIO

alignment = AlignIO.read(open('example.phy'), 'phylip')

var1 = var2 = var3 = 0
for col in zip(*alignment):
    num_unique = len(set(col))
    if num_unique == 1 and col[0] != '-':
        var1 += 1
    elif num_unique > 1 and '-' not in col:
        var2 += 1
    elif '-' in col:  # assumes 1 or more dashes
        var3 += 1

print('Number of columns with the same amino acid: {}\n'
      'Number of columns with at least 2 amino acids (no gaps): {}\n'
      'Number of columns with at least one gap: {}'
       .format(var1, var2, var3))

```

For your example data, this outputs:

```
Number of columns with the same amino acid: 9
Number of columns with at least 2 amino acids (no gaps): 40
Number of columns with at least one gap: 1

```
> 3 votes

# Answer
using something like `for col in column(alignment):`

1. I may get the number of columns which have all the same letters (amino acids)

```
len(set(col)) == 1

```
2. The number of columns which have at least 2 different amino acids (with no dashes):

```
len(set(col)) > 1 and '-' not in col

```
3. The number of columns which have any single dash:

```
'-' in col

```

or if you really want columns with a single dash:

```
    col.count('-') == 1

```
> 4 votes

---
Tags: sequence-alignment, python, phylogeny, phylogenetics
---


# Question
Title: Get a list of PDB ids from a list of sequence protein number?
I am working on a project about phosphorylation, and I am creating a database of some proteins. I have a list of protein accession ids for which I need to get the corresponding PDB file for each one to map it to PSM matrix. I want to get the PDB file using Python or BioPython tools. Any suggestions will be appreciated.

# Answer
A typical Uniprot URL would look like the example here:

https://www.uniprot.org/uniprot/?query=ID:P26378&format=tab&columns=id%2Cdatabase(PDB)

The response will be:

> Entry Cross-reference (PDB)
>
> P26378 1FXL;1G2E;

We want our query to include the protein ID, the output format should be tab separated and the only contain the columns `id` and `database (PDB)`.

To get the PDB IDs you would need to split the response, take the last line, split this line after each tab and split it again after each semi-colon.

Using Python's Requests library the whole code would be.

```
import requests
import time

uniprot_ids = ['P26378', 'O35433', 'Q02910']
url = 'https://www.uniprot.org/uniprot/'

protein_to_pdb = {}
for protein in uniprot_ids:
    params = {
        'format': 'tab',
        'query': 'ID:{}'.format(protein),
        'columns': 'id,database(PDB)'
    }
    contact = ""  # Please set your email address here.
    headers = {'User-Agent': 'Python {}'.format(contact)}
    r = requests.get(url, params=params, headers=headers)

    protein_to_pdb[protein] = str(r.text).splitlines()[-1].split('\t')[-1].split(';')
    protein_to_pdb[protein].pop(-1)
    time.sleep(1)  # be respectful and don't overwhelm the server with requests

print(protein_to_pdb)

```

---

* To accelerate things, you could store your `protein_to_pdb` dict in a file and only retrieve the new ones
* Checking if the returned ID matches your input ID is always a good idea
> 2 votes

# Answer
You could use mapping between UniProt and PDB entries from the SIFTS project. The `uniprot_pdb.csv` file there has UniProt accession names followed by a semicolon-separated list of PDB four letter codes.
> 0 votes

---
Tags: python, pdb
---


# Question
Title: Convert TPM-normalized matrices back to UMI in python
I want to process a plenty of `scRNAseq` datasets in `python`, and I want to run `TMM` normalization on all of them. Unfortunately some of them are already normalized with `TPM` method. Is there a way of either converting it back to `UMI`matrix, so that `TMM` normalization could be applied to it or could we just apply `TMM` directly to `TPM`-normalized matrices? Would that work? Any suggestions would be greatly appreciated.

# Answer
> 1 votes
It depends what your downstream requirement is. If you want to do count based statisitcs, then you definately shouldn't just do TMM normalisation on the TPM data as count statistics only work on real count data.

If you are intersted in doing something like clustering or dimensionality reduction, then I wouldn't use TMM, I'd probably use rLog, which includes normalisation.

You can have a crack at converting the TPMs back to something that looks like counts, but I not sure how well its going to work. To convert them back to REAL counts, you are going to need to know: the total number UMIs in the original dataset and the precise gene length used for each gene. Its posible the data set provides this ... Otherwise you are going to have problems.

$\frac{TPM\_i}{1,000,000}$ is the fraction of all transcripts in a sample that belong to a given gene.

If $$fL\_i = \frac{length\_i}{\Sigma\_{j\in G} length\_j} $$ is the fraction of the transcribed genome in a gene i. Thus $fR\_i = fL\_i \times \frac{TPM\_i}{1,000,000}$ should be the fraction of all counts that come from a gene $i$. Multiply this by the total number of reads to get the original counts.

The total number of reads may or may not be critcal, I think this might come out in the wash of the normalization, but it might mess with the disperision estimates, I'm not sure. But it is absolutely critical that

1) The TPM dataset was counted across exactly the same annotation as the count datasets

2) The gene lengths used are identical.

This could be untrue for several reasons:

* Different annotation or annotation version (ENSEMBL, UCSC, REFSEQ, CCDS etc).
* Different subset of transcripts used, or different way of combinng transcripts.
* Different ways of calculating "effective length", used to correct for things like which parts of the gene reads can or can't map to (reads too close to ends, unmappable regions of the gene etc).

Much safer to convert the count sets to TPM, but this would preclude doing count based statisitcs (e.g. DESeq2, edgeR, voom).

However I'd also worry that if these datasets were processed at this point differently, they are going to have been processed in different ways, and that the data is basically uncomparable.

---
Tags: python, scrnaseq, normalization, tpm, tmm
---


# Question
Title: How I can test my hypothesises computationally
I have single cell RNA-seq data on about 2000 cells in 9 time point. I have clustered my cells in each time point by Seurat. I am seeing in some time points I have 3 clusters while in another time points I have 2 clusters of cells (please look at the plot).

One hypothesis would be that, a new cell type is being emerged at last time point. But maybe if we had another further time point this third cluster being disappeared because that has been arose just because of undifferentiated cells from the later time point (second hypothesis).
Another hypothesis would be this new third cluster is a new cells type that only temporarily expresses a set of genes but if we had cells for another further time point, we would not see the expression of these genes anymore and we had again 2 clusters. I am really struggling with this how I could test my hypothesises.
The final hypothesis would be that the third cluster is really a new cell type born from old cell types during the differentiation.

Like UAD algorithm addresses my question. I just tried FateID algorithm on my data but got stopped by

```
> pr <- prcurve(y,fb,dr,k=2,m="tsne",trthr=0.4,start=2)
Error in svd(xstar) : a dimension is zero

```

Any suggestion would be highly appreciated.

# Answer
The answer is behind these papers:
The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution
Single-cell mapping of gene expression landscapes and lineage in the zebrafish embryo

Also this.
> 1 votes

# Answer
This is what packages like FateID and Monocle are for, namely taking single-cell RNAseq data and inferring differentiation trajectories from it. Don't try to reinvent the wheel on this, there are a number of packages out there to do this sort of thing and they're going to have pretty complicated methods that you're really not going to want to reinvent unless you absolutely have to.
> 3 votes

---
Tags: r, scrnaseq, seurat
---


# Question
Title: Finding the location and unit length of repetitive sequences within a long read
After discovering a few difficulties with genome assembly, I've taken an interest in finding and categorising repetitive DNA sequences, such as this one from *Nippostrongylus brasiliensis* \[each base is colour-coded as A: green; C: blue; G: yellow; T: red\]:

\[FASTA file associated with this sequence can be found here\]

These sequences with large repeat unit sizes are only detectable (and assembleable) using long reads (e.g. PacBio, nanopore) because any subsequence smaller than the unit length will not be able to distinguish between sequencing error and hitting a different location within the repeat structure. I have been tracking these sequences down in a bulk fashion by two methods:

1. Running an all-vs-all mapping, and looking for sequences that map to themselves lots of times
2. Carrying out a compression of the sequence (e.g. bzip2), and finding sequences that have a compression rate that is substantially higher than normal

After I've found a suspicious sequence, I then want to be able to categorise the repeat (e.g. major repeat length, number of tandem repeats, repetitive sequence). This is where I'm getting stuck.

For doing a "look, shiny" demonstration, I currently have a very manual process of getting these sequences into a format that I can visualise. My process is as follows:

1. Use LAST to produce a dot plot of self-mapping for the mapping
2. Visually identify the repetitive region, and extract out the region from the sequence
3. Use a combination of `fold -w <width>` and `less -S` to visually inspect the sequence with various potential repeat unit widths to find the most likely repeat unit size
4. Display the sequence in a rectangular and circular fashion using my own script, wrapping at the repeat unit length

But that process is by no means feasible when I've got thousands of potential repetitive sequences to fish through.

Is there any better way to do this? Given an arbitrary DNA sequence of length \>10kb, how can I (in an automated fashion) find both the location of the repeat region, and also the unit length (bearing in mind that there might be multiple repeat structures, with unit lengths from 30bp to 10kb)?

An example sequence can be found here, which has a ~21kb repeat region with ~171bp repeat units about 1/3 of the way into the sequence.

## A Kmer-based Analysis

I've now seen human sequences with repetitive regions in excess of 10kb (i.e. out of the range of most linked-read applications). My current idea is centred around creating hash tables of short sequences (currently 13-mers) and tracking their location:

1. Process the sequence, storing the location of each kmer
2. For each kmer, find out how many times it appears in the sequence
3. For repeated kmers, find out how much of a gap there is between the next time that kmer occurs
4. Report the median and modal gap length of repeated kmers, and statistics associated with their frequency in the sequence

Some local repetitive regions may be lost in the statistics with this approach, it's hard to tell if there are multiple repetitive regions within a single sequence, and if the repeat units are themselves slightly repetitive (enough that a kmer is duplicated within a repeat unit), then the algorithm will under-report the repetitiveness (see step 3).

# Answer
> 3 votes
I've sorted the visualisation out. Here are three alternative representations of repetitive structures for the same sequence:

These were generated using the same R script, callable from the command line:

```
 $ repaver.r -style dotplot -prefix dotplot MAOA_chrX_43645715-43756016.fa
 $ repaver.r -style profile -prefix profile MAOA_chrX_43645715-43756016.fa
 $ repaver.r -style semicircular -prefix semi MAOA_chrX_43645715-43756016.fa

```

More details about this are in the presentation I gave at Queenstown Research Week, 2018. I also wrote a chapter in a peer-reviewed eBook here.

This is fast enough that I can run it on the nanopore C. elegans genome in about half an hour, producing these plots for each contig. I don't quite have a method to iterate through this plot and pick out the dominant repeat length at each location, but that's a \[relatively\] simple extension of what I've already done.

With a lot of optimisations for speed and memory consumption, I've now been able to run it on the full human genome. It takes a couple of days on my desktop computer (64GB RAM + SSD swap) to categorise the 100bp repeat structure of the assembled T2T/CHM13-v1.0 chromosomes.

Code here.

# Answer
> 2 votes
It could be an idea to fragment the long reads into small sequences, like simulating Illumina reads of 150 bp, and then map these small sequences against the original long reads and extract regions with a high coverage?

---
Tags: nanopore, long-reads, repeat-elements
---


title: bioawk removed part of FASTQ header

text:
# Question
Title: bioawk removed part of FASTQ header
I used

```
bioawk -cfastx 'length($seq) > 1 {print "@"$name"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

```

in order to keep a particular read length, but this command shortened the header from

```
@A00199:161:HF3JLDMXX:1:1101:5882:1063 2:N:0:CTGTAGCC+GAACAGGC

```

to

```
@A00199:161:HF3JLDMXX:1:1101:5882:1063

```

How is it possible to fix the `bioawk` command to remain the full command?

# Answer
```
bioawk -cfastx 'length($seq) > 1 {print "@"$name $comment"\n"$seq"\n+\n"$qual}' in.fq.qz | gzip > out.fq.qz

```

Note the addition of `$comment`, which is what follows the name on the first line of each entry.
> 5 votes

---
Tags: bioawk
---


# Question
Title: Counting the characters in a row using python
I am using python 3 and biopython 1.72. I have been using a protein sequence file in clustal format.

```
CLUSTAL 2.1 multiple sequence alignment

Human           ------------------------------------------------------------
Chimpanzee      ------------------------------------------------------------
Dog             ------------------------------------------------------------
Mouse           ------------------------------------------------------------
Xenopus         ------------------------------------------------------------
Amphioxus       MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA

Human           ------------------------------------------------------------
Chimpanzee      ------------------------------------------------------------
Dog             ------------------------------------------------------------
Mouse           ------------------------------------------------------------
Xenopus         ------------------------------------------------------------
Amphioxus       HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK

Human           ------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS
Chimpanzee      ------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS
Dog             ------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS
Mouse           ------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS
Xenopus         ------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST
Amphioxus       INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD
                                  *:***::  :.  :::     *** **: *        * .

Human           NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI
Chimpanzee      NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI
Dog             DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL
Mouse           DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL
Xenopus         EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ
Amphioxus       DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH
                :  * *:** :*.* **::.* . **:***** ::* .   .               *

```

I need to iterate over the rows of my file so that I could count the occurrence of characters in each row. I tried using the code:

```
file_var = input("File: ")
count = 0
with open (file_var , 'r') as fo:
    for i in fo:
        count += 1
        print (i)
    print (count)

```

But this results just counting the total number of rows in a file`(count variable)`. If I iterate over the file using biopython, it does'nt work also.

```
for seq_record in SeqIO.parse("sample.aln","clustal"):
    print (seq_record.id)
    print (repr(seq_record.seq))

```

**seq\_record.id** results the names of the sequences used in the file. **repr(seq\_record.seq)** results the sequences of each specie alongwith description.

While I need to

1. get the rows only with amino acids (exclusing the symbols `'*,.,:'` line) and nothing else from the row and for each row I need to count the number of characters (exclusing `'-'` dashes)
2. the number of characters of that row be placed at the end of each of the row.

for example: for above example file, the counter should work like;

```
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42

```

# Answer
Here's one way with pure Python:

```
#!/usr/bin/env python3
import sys
import re

with open (sys.argv[1] , 'r') as fo:
    for line in fo:
        fields = line.split()
        if len(fields)==2:
            print(fields[1],len(fields[1].replace('-','')))

```

Running this on your file gives:

```
$ foo.py file.aln
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS 42
------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS 42
------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS 42
------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST 42
INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD 60
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL 45
EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ 45
DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH 60

```

---

Some other options:

1. `awk`

```
$ awk 'NF==2{l=$2;gsub(/-/,"",$2); print l,length($2)}' file.aln
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS 42
------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS 42
------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS 42
------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST 42
INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD 60
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL 45
EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ 45
DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH 60

```
2. Perl

```
$ perl -lane '$l=$F[1];$l=~s/-//g; $#F==1 && print "$F[1] ", length($l)' file.aln
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
MQWTGFRVSMTTLMMIMGVVAVLIALLPAKAQQPHDKSLRTTSTLTDTGASADEADMGSA 60
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
------------------------------------------------------------ 0
HVELLDGDDDVGNGSDQMMVTLHLQSIFQCIRRPCEKVDRAIDPVTQRWRTANTRNDYQK 60
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSHLRQSLLCTWGYSS 42
------------------MRLRVRLLKRTWPLEVPETEPTLGHLRSRLRQSLLCTWGYSS 42
------------------MKLRVRLQKRTWPLDLPDAEPTLGQLRAHLSQALLPSWGFGS 42
------------------MKLRVRLQKRTQPLEVPESEPTLGQLRAHLSQVLLPTLGFSS 42
------------------MKLRVRVRKQTNRLELEAESPTLGDLRSKLSSVTLPALGYST 42
INVCVVPAYDVSLSTGVRMKLRVKISGQKTRVDVGQDCHTLGTLRTLLAPVLGEQYGLGD 60
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
NTRFTITLNYKDPLTGDEETLASYGIVSGDLICLILQDDIPA---------------PNI 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLILEDAIPA---------------PNL 45
DTRFAITLNNKDALTGDEETLASYGIVSGDLICLVLEDDMPA---------------PNL 45
EANFTITLNGKDALTGDQNTLESAGIISGDLIVVVLPDSQLP---------------PPQ 45
DMPFEISLNGRDALLGDDKPLSDLGIVSGDLIHILLASVDQPSTNHNTQQQGQHPSSPEH 60

```
> 5 votes

---
Tags: python
---


title: How is prior.count used by edgeR's cpm

text:
# Question
Title: How is prior.count used by edgeR's cpm
edgeR's cpm function has an argument called `prior.count`. Based on my understanding of the documentation, it is supposed to be adding a fixed number per sample which is proportional to the library size of said sample. Average of all numbers added to all samples would be equal to `prior.count`.

However, looking at actual data this does not seem to be the case. Given an imaginary data frame of

```
df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)

```

We can try to calculate log cpms by doing

```
logCPM = cpm(df,log = TRUE,prior.count = .5)

```

We can also calculate regular cpms by doing

```
CPM = cpm(df)

```

If we wanted to see, what is the number that is added to the CPMs before getting logged, we could do

```
difference = 2^logCPM - CPM

```

To make it look pretty lets use tibble

```
tibble::as.tibble(difference)

# A tibble: 4 x 3
                     a                    c                    b
                 <dbl>                <dbl>                <dbl>
1  25641                25641                25641
2  12821                12821                12821
3 -    0.0000000000582 -    0.0000000000582 -    0.0000000000582
4  38462                38462                38462

```

Here we see that the number that is added has no relation to the library size of the sample. What I want to learn is, how is this number that is added to each cpm is calculated based on prior.count.

I have been digging through the code but that part goes into C++ territory which hinders my understanding of what is going on

# Answer
> 7 votes
The prior count ends up getting scaled by the ratio of a library size to the average library size and then multiplied by 2 before getting added to each library size (I'm sure there's a good reason for that, but I don't know what it is). Using the example `df` data frame in your post, let's walk step-by-step through what `cpm()` is doing:

```
df = data.frame(a = c(1,2,3,0),c = c(1,2,3,0)*3, b = c(1,2,3,0)*2)
prior.count = 0.5
# First, we need to calculate a library size
lib.size = colSums(df)
# Calculate the average library size and the adjusted priors
ave.lib = mean(lib.size)
adjusted.prior = prior.count * lib.size / ave.lib
# Update the library sizes
adjusted.lib.size = lib.size + 2*adjusted.prior
# Now we can compute the CPM
customCPM = t((log(t(df) + adjusted.prior) - log(adjusted.lib.size) + log(1000000))/log(2))

```

The matrix transposition looks really ugly in the last step, but it's needed so R adds things like `adjusted.prior` to each row of `df`.

Note that the prior count is only used if you compute log(CPM), it's unused otherwise.

---
Tags: r, edger
---