# Question
Title: How do you generate read-length vs read-quality plot for long-read sequencing data (e.g., MinION)?
How do you generate read-length vs read-quality plot (heat map with histograms in the margin) for long-read sequencing data from the Oxford Nanopore Technologies (ONT) MinION? The MinKNOW software from ONT provides a plot like this during base calling.

This would also be very helpful for PacBio data.

# Answer
I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using:

* A fastq file (optionally compressed)
* A bam file
* The sequencing\_summary.txt file generated by albacore

Using optional flags you can:

* Log transform the read lengths
* Use aligned reads rather than sequenced reads
* Downsample the reads
* Set a maximum read length

I've added an example below, plotting log transformed read length vs average read quality (using a kernel density estimate). More examples can be found in the gallery on my blog.

I welcome all feedback and suggestions!
> 6 votes

# Answer
It's important to always consider read length and quality jointly with high-error read data, and current long-read technologies (e.g., MinION and PacBio) have high error rates. ***Considering read length and quality jointly will help you determine how successful the run was, how many reads were 'high quality', whether the longer reads are 'real' (or just pore noise), etc.***

I've had a recent spike of interest in similar plots and came across a project called pauvre (french for 'poor', play on 'pore') through the Oxford Nanopore Technologies (ONT) community that I think is even better than MinKNOW's base calling plot. Plus, you can generate these plots from a fastq file when ever you want to, unlike with MinKNOW.

***\[Note: I am not the original author, but I am now contributing because I liked (and needed) it.\]***

Pauvre will also report useful statistics:

```
fastq stats for fastq_runid_bb8b8ddedb22bdd6802b2bfa2b4e424c92c30d28_0.fastq
numReads: 2164829
numBasepairs: 4970615217
meanLen: 2296.077527139557
medianLen: 1495.0
minLen: 5
maxLen: 392031
N50: 3450
L50: 402786

                               Basepairs >= bin by mean PHRED and length
minLen          Q0          Q5         Q10         Q15       Q17.5         Q20      Q21.5  Q25  Q25.5  Q30
     0  4970615217  4970611559  4835461771  3889995868  2900103275  1087779109  165637656  429      0    0
 50000    11531044    11531044      270324      160128       50729       50729          0    0      0    0
100000     6260554     6260554           0           0           0           0          0    0      0    0
150000     3504240     3504240           0           0           0           0          0    0      0    0
200000     2501101     2501101           0           0           0           0          0    0      0    0
250000     1609592     1609592           0           0           0           0          0    0      0    0
300000     1033423     1033423           0           0           0           0          0    0      0    0
350000      392031      392031           0           0           0           0          0    0      0    0

                    Number of reads >= bin by mean Phred+Len
minLen       Q0       Q5      Q10      Q15    Q17.5     Q20  Q21.5  Q25  Q25.5  Q30
     0  2164829  2164605  2083436  1626706  1183812  435687  77341    1      0    0
 50000      109      109        5        3        1       1      0    0      0    0
100000       36       36        0        0        0       0      0    0      0    0
150000       15       15        0        0        0       0      0    0      0    0
200000        9        9        0        0        0       0      0    0      0    0
250000        5        5        0        0        0       0      0    0      0    0
300000        3        3        0        0        0       0      0    0      0    0
350000        1        1        0        0        0       0      0    0      0    0

```

***These plots and stats would be equally useful with PacBio, but that's not super easy (though it is possible) with current raw output from the Sequel sequencer: Which quality score encoding does PacBio use?***

Pauvre currently uses Biopython to parse the fastq and `matplotlib` for the actual plot, and will let you choose the output image format (e.g., .png, .pdf, etc.). You can also choose whether the background is transparent or white (for .png output).

The parser is currently super slow because it's using `SeqIO.parse`, but we're changing parsers to speed that up. We're also adding some extra features (e.g., choose whether to include y-axes in margin histograms, print some stats directly to the plot for documentation, etc.)

Purple is currently the only color choice (which I personally love), but adding options to change that will be super easy.
> 4 votes

---
Tags: python, ngs, fastq, nanopore, pacbio
---


# Question
Title: Running differential expression analyses on count matrices with many zeroes
Note: I also posted this issue (with less context) in the bioconductor support site: https://support.bioconductor.org/p/97424/

I'm working on a snakemake workflow that identifies various small RNA species in *C. elegans* small RNA-seq libraries. Some are endogenous siRNAs supposedly generated from RNA templates (through RNA-dependant RNA polymerases (RdRP)) that can be variously classified (protein-coding genes, transposons, and other repeat types).

I count such small RNA reads in a set of libraries and try to identify differentially producing sources. For this, I use DESeq2 (that I run using rpy2 from within snakemake).

I'm not sure DESeq2 is always appropriate for this kind of data, but so far, the analyses would at least complete. However, I recently added new potential types of small RNA source (simple repeats and satellites), and these happen to have low counts. I'm not 100% sure, but I suspect these low counts are the reason for failures during DESeq2 analyses (for debugging purposes, I ran this manually in R):

```
> dds <- DESeq(dds, betaPrior=T)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  :
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)
> warnings()
Warning messages:
1: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
2: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
3: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
4: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
5: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
6: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
7: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
8: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
9: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
10: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
11: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight
12: In lfproc(x, y, weights = weights, cens = cens, base = base,  ... :
  procv: no points with non-zero weight

```

The count matrix indeed contains a high proportion of zeroes, the rest are mostly ones, but there are some sources that seem a little more significantly producing small RNAs:

```
> mean(counts_data == 0)
[1] 0.7488889
> mean(counts_data == 1)
[1] 0.1507407
> max(counts_data)
[1] 34

```

**How would you recommend handling such kind of data?**

Should I for instance discard rows without enough counts? If so, what would be a reasonable threshold?

### Edits: trying to get a dispersion plot:

The same error occurs when trying to estimate dispersion:

```
> dds <- estimateSizeFactors(dds)
> dds <- estimateDispersions(dds)
gene-wise dispersion estimates
mean-dispersion relationship
-- note: fitType='parametric', but the dispersion trend was not well captured by the
   function: y = a/x + b, and a local regression fit was automatically substituted.
   specify fitType='local' or 'mean' to avoid this message next time.
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  :
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)

```

This seems to prevent the generation of dispersion plots:

```
> plotDispEsts(dds)
Error in min(py[py > 0], na.rm = TRUE) :
  invalid 'type' (list) of argument
In addition: Warning message:
In structure(x, class = unique(c("AsIs", oldClass(x)))) :
  Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

```

Using `fitType="local"` also fails:

```
> dds <- DESeq(dds, betaPrior=T, fitType="local")
using pre-existing size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
Error in lfproc(x, y, weights = weights, cens = cens, base = base, geth = geth,  :
  newsplit: out of vertex space
In addition: There were 12 warnings (use warnings() to see them)

```

But `fitType="mean"` works:

```
> dds <- DESeq(dds, betaPrior=T, fitType="mean")
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

```

The dispersion plot then looks as follows:

**What do these `fitType` options mean, and are there good reasons why they would have to be changed in order to accomodate for low counts?**

# Answer
> 6 votes
I don't think that the issue is the low counts, but rather the number of features without any real variance (the black dots at the bottom).

# So what the heck is the dispersion plot and why does one need to fit it anyway?

In a typical RNAseq experiment, one measures many thousands of genes with only a few replicates per biological group. This then leads to issues when performing testing, since your variance measurement for any given gene will be larger or smaller than reality, simply due to the lack of a large number of replicates. Consequently, at least since the `limma` package came out, most RNAseq packages have tried to fit the mean-dispersion with some sort of line and then used that as a background distribution toward which the variance assigned to each gene should be "shrunken" (i.e., moved toward). Getting this at least approximately correct is important, since if you have few replicates and the final "shrunken variance" is significantly wrong, then your p-values will also be wrong.

# What are the various fitType options and what are they doing?

The three options that one can use to fit this relationship are `parametric`, `local`, and `mean`. `parametric` is the default and basically results in fitting to a bent line. This usually works quite well, but you can see that you have two clouds of data, one shaped like a bent line (at the top) and another (the black dots at the bottom) with no variance. That's causing this fit type to fail. Likewise, `local` is doing a local fit and failing for the same reason. `mean` is just taking the mean (average) of all of the estimates. This is fool-proof, but you're then regressing to something that really doesn't fit your data. `mean` is really a last choice if nothing else works.

For your particular dataset, my main concerns would be figuring out (1) why the maximum mean counts are around 20 (that's really low) and (2) why so many of these features have incredibly low variance.

# Answer
> 5 votes
Remove low-count features in advance. This is standard for most tools including DESeq2 and edgeR (see section 2.6).

This will keep you from testing a lot of features that cannot be differentially expressed (using NB model you need at least four reads - below that the uncertainty is too high to ever be DE) and you will make the false discovery correction unnecessarily hard.

---
Tags: r, bioconductor, differential-expression, deseq2
---


# Question
Title: Can you tell me where these RNAi’s cleave on the Drosophila melanogaster CG2316 mRNA?
How can I work out at which sites on the Drosophila melanogaster CG2316 mRNA the RNAi\`s:

* 12170/FBti0089992,
* 12168/FBtp0030589,
* 107343/FBtp0042163 and
* 41984/FBst0034349

(stock numbers/flybase ID) cleave?

I attempted to work this out myself some time ago but with limited success. I found the Fasta sequence of CG2316 and searched for the two complementary sites of the RNAi hairpin using the word doc find function. I feel this rudimentary method was flawed on a few levels, not least because I only got a single hit. As such I have posted the problem here.

# Answer
Firstly, I would download some sort of editor that is built to handle nucleotide sequences. There are a handful out there, but APE is a fairly popular software designed for plasmid editing. It is much more efficient at searching for nucleotide strings as it can search for reverse complements and to a certain extent, mismatches. It's also free!

http://biologylabs.utah.edu/jorgensen/wayned/ape/

For your question, I'm not exactly sure what you are asking, but I would use the BLAST alignment tool to try to see similarity between your RNAi sequences and the target mRNA. You input fasta files and the program searches for any possible overlap between the two inputs. The higher the "max-score" the more likely the two strands (or subsections of strands) have high similarity. This program also considers the complement

https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE\_TYPE=BlastSearch&LINK\_LOC=blasthome

You can use the results from BLAST to annotate in APE using different colors which is extremely helpful. You can also use BLAST in the user interface of APE.

Also keep in mind RNAi's are not always effective, even though they are predicted to bind at a given location, that may not be a biological reality

Good luck!
> 1 votes

---
Tags: fasta, transcriptome-regulation
---


# Question
Title: How can I dock a protein to a nucleic acid?
I have a protein of interest and I would like to now how it interacts with RNA. I have structures of both molecules.

What tool can I use?

# Answer
I'm not sure what you meant but you can take a look at NPDock (disclaimer we wrote that tool). If you have a structure of your protein of interest, you can dock it to the structure of your DNA/RNA of interest. Mind that this is a rigid body docking which means that the structure will not change upon binding.
> 3 votes

---
Tags: proteins, software-recommendation, rna, docking
---


# Question
Title: What are the labels in my rstats distance matrix?
I have a distance matrix generated by hierfstat, thusly, with a link to the fasta file here:

```
library(adegenet)
library(hierfstat)
snps <- fasta2DNAbin('test.fa', chunkSize = 50)
gi <- DNAbin2genind(snps)
# manually define the populations for passing to pairwise.fst
p <- c('Botswana', 'Botswana', 'Botswana', 'Botswana', 'France', 'France', 'Vietnam', 'Vietnam', 'Uganda', 'Uganda', 'Uganda', 'Uganda', 'Vietnam', 'Vietnam', 'Laos', 'Laos', 'Laos', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam', 'Vietnam')
f <- pairwise.fst(gi, p, res.type = c('dist', 'matrix'))
as.matrix(f)

```

which produces

```
    1   2   3   4   5
1   0   0.2189008   0.225567409 0.1821518   0.259409722
2   0.2189008   0   0.130736953 0.1648034   0.191050772
3   0.2255674   0.130737    0   0.1669077   0.006396789
4   0.1821518   0.1648034   0.166907683 0   0.203931457
5   0.2594097   0.1910508   0.006396789 0.2039315   0

```

Is it safe to assume that 1 = Botswana, 2 = France, 3 = Vietnam, 4 = Uganda and 5 = Laos? i.e. that the distance matrix rows/columns follow the order that each population first appeared in `p`?

Is there a way to determine for certain, which row/column in the distance matrix corresponds to which population?

# Answer
Disclaimer: I don't have adagenet and hierfstat installed, I just looked at the source code. It seems like pairwise.fst (the deprecated code from the adagenet package, which you seem to be using) already returns a matrix. Set res.type="matrix" and try calling `rownames(f)` and `colnames(f)` afterwards.
> 2 votes

---
Tags: r
---


# Question
Title: How to convert a phylogeny to a dendrogram in R
I want to convert my phylogeny into a dendrogram so I can use it with dendextend in R to produce a tanglegram. I have made some progress but I keep encountering errors, see below:

```
library(ape)
library(dendextend)

Tree <- rtree(10, rooted=F)

Tree <- read.tree(file="clipboard", text=NULL)
# test <- as.dendrogram(Tree)
## Error in ape::as.hclust.phylo(object) : the tree is not ultrametric

Tree.ultra <- chronos(Tree)
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not rooted

Tree.ultra$root.edge <- 0
# test <- as.dendrogram(Tree.ultra)
## Error in ape::as.hclust.phylo(object) : the tree is not binary

```

This thread seems to be the closest to answering my issues but the workflow is different and Im not sure how to address this binary issue. I have tried using as.hclust, as.hclust.phylo and hclust.

# Answer
After searching high and low I have found an answer from this thread

Workflow goes like:

```
library(DECIPHER)
dend1 <- ReadDendrogram(file="clipboard")

# Alternatively to load from existing environment
Tree$node.label <- NULL # Need to remove node labels
dend1 <- ReadDendrogram(textConnection(write.tree(Tree)))

```

This object can then be used in the dendextend package and tanglegram function.
> 1 votes

# Answer
For your example to work you should replace the `rtree` function with the `rcoal` function, see example below.

```
library(ape)
library(dendextend)
Tree <- rtree(10, rooted=F)
is.ultrametric(Tree)
[1] FALSE
is.binary.tree(Tree)
[1] TRUE

Tree_rcoal <- rcoal(10, rooted=F)
is.ultrametric(Tree_rcoal)
[1] TRUE
is.binary.tree(Tree_rcoal)
[1] TRUE

Tree_rcoal2 <- rcoal(10, rooted=F)

tanglegram(Tree_rcoal, Tree_rcoal2)

```

To make use of a file in newick format, you can try to convert with the `chronos` function.

```
s <- "owls(((Strix_aluco:4.2,Asio_otus:4.2):3.1,Athene_noctua:7.3):6.3,Tyto_alba:13.5);"
cat(s, file = "ex.tre", sep = "\n")
tree.owls <- read.tree("ex.tre")
is.ultrametric(tree.owls)
[1] FALSE
dendrogram <- chronos(tree.owls)

Setting initial dates...
Fitting in progress... get a first set of estimates
         Penalised log-lik = -24.04992
Optimising rates... dates... -24.04992
Optimising rates... dates... -24.04908

Done.
is.ultrametric(dendrogram)
[1] TRUE

```
> 2 votes

---
Tags: r, phylogenetics, phylogeny
---


# Question
Title: How do I generate a color-coded tanglegram?
I want to compare two phylogenies and colour the association lines based on some metadata I have. I have been using ape cophyloplot but I have not had any success in getting the lines to colour accurately according to my data (see previous question).

Note that in my actual work flow I define the colour scheme using a palette to control the colour outcome.

I want a means to make a tangle using phylogenies which I can format. Preferably in R. I like to get an output like this:

# Answer
> 3 votes
Thanks everyone for the suggestions @b.nota answer is useful to colour according to clade groups but does not address my actual question on using metadata to colour the lines. For this answer see below:

dendextend states for `color_lines` a vector of colors for the lines connected the labels. If the colors are shorter than the number of labels, they are recycled (and a warning is issued). The colors in the vector are applied on the lines from the bottom up.

```
# Extract labels from dendrogram on the left
labels <- dendA %>% set("labels_to_char") %>% labels

#Using a metadata table with colours create a vector of colours
labels <- as.data.frame(labels)
labels2 <- merge(labels, metadata, by.x="labels", by.y="Sample.name", sort=F)
cols <- as.character(labels2$Colours)

# Make tanglegram
tanglegram(dendA, dendC, color_lines = cols)

```

Side note if anyone with a phylogeny needs to convert it to a dendrogram to try this on their data see my other question.

# Answer
> 11 votes
I think you can try dendextend, in this manual there is an example of coloring the branches. I don't think it is exactly like your coloring, but with a little tweaking you might get your colorscheme in there.

The manual mentions an argument called `color_lines` for the function `tanglegram()`:

```
# The `which` parameter allows us to pick the elements in the list to compare
iris_dendlist %>% dendlist(which = c(1,4)) %>% ladderize %>%
   # untangle(method = "step1side", k_seq = 3:20) %>%
   set("rank_branches") %>%
   tanglegram(common_subtrees_color_branches = TRUE)

```

According to the manual, the code above produces an image like this:

Try to see if you can fit your color palette in there.

# Answer
> 3 votes
To get your colors in the original order, you can make a dataframe with the labels in the right order with your metadata. You can use `merge` for this.

So from your example code it would be something like this:

```
site <- structure(list(name = structure(c(1L, 3L, 4L, 5L, 6L, 7L, 8L,9L, 10L, 2L), .Label = c("t1", "t10", "t2", "t3", "t4", "t5","t6", "t7", "t8", "t9"), class = "factor"), site = c(1L, 1L,1L, 2L, 2L, 3L, 1L, 3L, 2L, 2L)), .Names = c("name", "site"), row.names = c(NA,10L), class = "data.frame")

library(ape)
library(dendextend)

t1 <- rcoal(10)
t2 <- rcoal(10)

str(site)
# The name are factors now, so make chr first
site$name <- as.character(site$name)

# make a data.frame of your labels
labels_df <- data.frame(t1$tip.label)

#merge the 2 data.frames in the right order (hence sort=F)
colors <- merge(labels_df,site,by.x="t1.tip.label", by.y="name",all.x=T, all.y=F,sort=F)

tanglegram(t1, t2, color_lines=colors$site)

```

# Answer
> 2 votes
As an alternative (though I realise this doesn't really answer the OPs question directly), Dendroscope and Treemap can also produce these types of plots. I don't know if they'll do coloured connections off the top of my head, but I think they can differentially colour the tree branches.

Just in case that is of use!

---
Tags: r, phylogeny, software-recommendation, phylogenetics
---


# Question
Title: Trouble using biomaRt to retrieve hgnc symbols from Ensembl transcript ids
I have a matrix of gene counts which I'm going to use as input for DESeq. Right now, each gene is labeled by its Ensemble transcript ID, but I'd like to convert these to their HGNC symbols before I input them into DESeq for analysis. I'm attempting to do this conversion using biomaRt, however, this is my first time using the program, and I'm running into issues. Here is my current code:

```
library("biomaRt")
mart <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", dataset = "mmusculus_gene_ensembl")
transcript.ids <- rownames(txi.kallisto$counts)
hgnc_symbols <- getBM(attributes = "hgnc_symbol", filters = "ensembl_transcript_id", values = transcript.ids, mart = mart)

```

Printing hgnc\_symbols shows that it isn't what I want:

```
>ghnc_symbols
[1] hgnc_symbol
<0 rows> (or 0-length row.names)

```

Lastly, here is what the input transcript IDs look like:

```
> head(rownames(txi.kallisto$counts))
[1] "ENSMUST00000178862.1" "ENSMUST00000178537.1" "ENSMUST00000196221.1"
[4] "ENSMUST00000179664.1" "ENSMUST00000177564.1" "ENSMUST00000179520.1"

```

# Answer
You need to specify the number without the version. Instead of "ENSMUST00000178862.1" just "ENSMUST00000178862":

You can do this with one more line:

```
g <- gsub("\\..*", "", rownames(txi.kallisto$counts))
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = g, mart = mart))
##  hgnc_symbol chromosome_name ensembl_transcript_id
##1          NA              14    ENSMUST00000177564
##2          NA               6    ENSMUST00000178537
##3          NA               6    ENSMUST00000178862
##4          NA              12    ENSMUST00000179520
##5          NA              14    ENSMUST00000179664
##6          NA              14    ENSMUST00000196221

```

With these transcripts it doesn't seem that there is an equivalent hgnc\_symbol. But is recognized, otherwise they are not recognized:

```
(hgnc_symbols <- getBM(attributes = c("hgnc_symbol", "chromosome_name", "ensembl_transcript_id"), filters = "ensembl_transcript_id", values = rownames(txi.kallisto$counts), mart = mart))
## [1] hgnc_symbol           chromosome_name       ensembl_transcript_id
##<0 rows> (or 0-length row.names)

```

However you might be interested in doing something like this:

```
human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
getLDS(attributes = "ensembl_transcript_id",
       filters = "ensembl_transcript_id", values = g,mart = mouse,
      attributesL = "hgnc_symbol", martL = human)
## The query to the BioMart webservice returned an invalid result: the number of
## columns in the result table does not equal the number of
## attributes in the query. Please report this to the mailing list.

```
> 8 votes

---
Tags: r, bioconductor, deseq2, biomart, ensembl
---


# Question
Title: Is there any way of using biopython to write Swissprot files?
I have the following mwe for filtering a Swissprot file based on a certain feature, in this case, transmembrane regions.

```
from Bio import SeqIO
records=[]
for record in SeqIO.parse("Input.txt", "swiss"):
    transmembrane_protein=False
    print record.id
    for i, feature in enumerate(record.features):
        if feature.type == "TRANSMEM":
            transmembrane_protein=True
    if transmembrane_protein==True:
        records.append(record)
SeqIO.write(records, "Output.txt", "swiss")

```

The script works when `SeqIO.write(records, "Output.txt", "swiss")` becomes `SeqIO.write(records, "Output.txt", "fasta")`

However, such a method is not yet supported.

```
ValueError: Reading format 'swiss' is supported, but not writing

```

From the docs, I see that writing is not supported for swiss:

> Note that while Bio.SeqIO can read all the above file formats, it cannot write to all of them.

Is there any unofficial way of using biopython/python to write swissprot files from parsed swissprot files?

# Answer
Using `SeqIO.index` rather than `SeqIO.parse` lets you read all the records into a `dict`, from which you can then extract the raw entry:

```
from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
for key in record_dict:
    print(record_dict.get_raw(key).decode())

```

Now you should be able to apply your test for a transmembrane protein to each entry, and write out only the ones you want to keep.

```
from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
out = open('Output.txt', 'wb')
for key in record_dict:
    record = record_dict[key]
    for i, feature in enumerate(record.features):
    if feature.type == "TRANSMEM":
        out.write(record_dict.get_raw(key).decode())

```
> 6 votes

# Answer
I don't know how to do this with biopython, but it's simple enough in Perl using Swisskinfe:

```
#!/usr/bin/env perl
use strict;
use warnings;
use SWISS::Entry;

my $flat=$ARGV[0]||die("Need a flat file\n");

## Change the line termination string so we read an
## entire entry at a time
local $/ = "\n//\n";

my $fh;
## Deal with compressed or uncompressed files
if($flat=~/\.gz$/){
    open($fh,"zcat $flat |")|| die("cannot open flat file $flat : $! \n");
}
else{
    open($fh,"$flat")|| die("cannot open flat file $flat : $!: $@\n");
}
## Parse the file.
entry: while (my $current_entry = <$fh>) {
    ## Read the current entry as an "entry" object
    my $entry = SWISS::Entry->fromText($current_entry);
    ## Iterate over the entry's features
    foreach my $feature ($entry->FTs->elements()) {
        ## If any of the features are 'TRANSMEM', print the entry
        if ($feature->[0] eq 'TRANSMEM') {
            print "$current_entry\n";
            ## Skip to the next entry
            next entry;
        }
    }
}

```

You need to install the Swissknife module first:

1. Download Swissknife\_1.73.tar.gz from here (SourceForge link).
2. Install

```
tar xvzf Swissknife.tar.gz
cd Swissknife_1.73/
perl Makefile.PL
sudo make install

```

Now, make the script above executable and run with:

```
foo.pl file.flat

```
> 2 votes

---
Tags: database, python
---


# Question
Title: How to find hairpin loop structure in a large set of long non coding RNA transcripts
I have done de novo assembly of pair end raw read sequences and resulted transcripts sequences were separated based on coding potential into two categories: long non-coding RNA transcripts and coding RNA transcripts.

I want to find how many of long non coding can make hairpin loop structure and how can I annotate these miRNA which are generated from long non coding RNA in plants.

# Answer
For the first question: "how many of long non coding can make hairpin loop structure?"

The easiest thing to would be to run your sequences through an RNA secondary structure prediction tool. There's a few tools for doing this but the most commonly used are RNAfold from the ViennaRNA package and mfold.

This will give you an output that looks like this, where the parentheses indicate which bases are paired with each other.

`GGGCUAUUAGCUCAGUUGGUUAGAGCGCACCCCUGAUAAGGGUGAGGUCGCUGAUUCGAAUUCAGCAUAGCCCA
(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).`

To see if there's hairpins, you'll need to annotate the elements within it. Take a look at this question on biostars for an example of how to do this. In short, you can use the forgi library to produce an annotation like this:

`(((((((..((((.........)))).(((((.......))))).....(((((.......)))))))))))).
sssssssmmsssshhhhhhhhhssssmssssshhhhhhhsssssmmmmmssssshhhhhhhsssssssssssse`

If there's an `h` in the output annotation, you have a hairpin. The remainder of your question about annotating miRNAs will be left to somebody else.
> 4 votes

---
Tags: bioconductor, rna-structure
---


# Question
Title: Obtaining data table headers from GEO using GEOquery
For a study in GEO, I would like to obtain the data table header descriptions, specifically the "VALUE" column for all samples in the study.

If you go here: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE99511) and then scroll down and click one of the samples: let's choose "GSM2644971". Then scroll down and you should see "Data table header descriptions" and below that you should see "**VALUE** Normalized (provided the normalization method) Average Beta". That information is what I want.

I tried using assayData() from the Biobase package but I didn't know whether the method takes a sample, a matrix of samples, or something else as a parameter. Apparently it takes a S4 and then returns some random text but I don't know how to use that random text.

EDITED:

For example, if I do:

```
library(Biobase)
library(GEOquery)

getgeo<-getGEO("GSE99511")
assayData(getgeo$GSE99511_series_matrix.txt.gz)

```

When I use assayData() I would want it to return "Normalized (provided the normalization method) Average Beta".

Instead it returns:

```
<environment: 0x110674178>

```

Please let me know if the question doesn't make sense or if there's another method that I should be using.

# Answer
The short answer is that if you are seeing information in a GEO sample, then it is the sample that you need to download and access using the appropriate method. In this case, that method is `Columns()`:

```
library(GEOquery)
gsm <- getGEO("GSM2644971")
Columns(gsm)

          Column      Description
1         ID_REF
2         VALUE       Normalized (provided the normalization method) Average Beta
3         Detection   Pval

```

The longer answer: you will really benefit from sitting down with the GEOquery documentation and trying to understand the objects that it uses and their methods. It is not something you can work out by trial and error for every case-by-case basis. The tutorial is a really good starting point (you will find `Columns()` in there, for example).
> 3 votes

---
Tags: r, database, bioconductor, text
---


# Question
Title: How are PDO and PDX used in computational and predicative models for tumour biology?
There are wet methods: Patient-derived models: Patient Derived Xenograft (PDX) and Patient Derived Organoids (PDO) to reflect tumor biology.

Are there any databases or computational tools that use the outcomes from PDO/PDX experiments to create a predictive computational model for cancer or other diseases?

# Answer
> 1 votes
**Potential pitfall!**

I'm not sure about predictive models, but you need to be aware of a potential pitfall in blindly aligning PDX or PDO based sequencing data without first removing contaminating host organism reads, as otherwise these will lead to a lot of false positive variants caused by miss alignment. In my experience even a small mount of host material can lead to a ten fold increase in called variants due to miss aligned host reads looking like true variants. I'd recommend using Xenome, source.

Note Xenome was designed for DNA sequencing it's not splice-aware, so for RNA-Seq there is no optimal solution other than removing reads which align to the host organisms genome. Although the issue here is that in conjunction with a splice-aware read aligner synteny between the host and grafted genomes might create some interesting problems. However for RNA-Seq provided FACS or similar confirms low-levels of host contamination I expect levels of expression will not be badly affected. Although this really needs investigating.

Finally and rather annoyingly Xenome produces none-standard FASTQ so you'll need to fix it's output with:

`awk '{if (NR % 4 == 1) print \"@\"$0; else if (NR % 4 == 3) print \"+\"$0; else print $0 }'` as reported here on seqanswers.

---
Tags: cancer, drugs, modelling
---


# Question
Title: Run kallisto iteratively across many samples
I am on a Mac using UNIX. I am trying to use the kallisto quant command on all files in a directory (instead of manually entering them). Because I'm running the analysis against the same index file, I first enter the following:

```
./kallisto index -i --index --make-unique index.fa

```

This successfully creates an index file. Then, I tried this:

```
for file in *.fasta; do kallisto quant --single -l 200 -s 0.1 -o $file-aligned
"$file"; done

```

While the $file-aligned folders are created, they're empty. I get an error that states that the index file is missing. So I assumed I would need to specify the index file in the command line. I then tried this:

```
for file in *.fasta; do kallisto quant -i index.fa --single -l 200 -s 0.1 -o
$file-aligned  "$file"; done

```

I get an error that says "Error: incompatible indices. Found version 0, expected version 10.

I have now also tried the following:

```
for file in *.fasta; do ./kallisto quant -i index --single - 200 -s 0.01 -o
"${file}-aligned" "${file}";

```

but I get an error saying the index file is not found.

---

I am now running

```
or file in *.fasta; do ./kallisto quant -i index --single -l 200 -s 0.1 -o $file-aligned

```

"$file"; done

This sort of works. The index file error messages are gone, but it runs quant on the same file 10 times instead of each file once. If I run echo $file it spits back the file that my command is running repeatedly.

If I close terminal, start a new session and run the same command, I get

```
Error: file not found ""

```

# Answer
**Running the jobs in parallel**

Iteratively is one solution, but since you're not using bootstrapping, which can be run multithread (and is not needed if you're taking the popular tximport route) you can do better using GNU parallel, this will enable you to run as many jobs as you have execution threads on your system simultaneously. Note that Kalliso only runs multithreaded with `-t` if you're using bootstrapping.

1. You'll need to make the index first:

```
kallisto index --make-unique -I index.fa.idx index.fa

```
2. Next lets define the input FASTQ files and the sample IDs you want to use using a tab-delimited input file like so:

```
sample_1    sample_1_R1.fastq.gz
sample_2    sample_2_R1.fastq.gz
sample_3    sample_3_R1.fastq.gz

```

Here the first column is the sample ID which will be used for output and the second separated by a tab character is the input file.

You can produce this using a text editor but as a short-cut you can feed the list of fastq.gz you have into a list via:

```
ls -1 *.fastq.gz > editme

```
3. Next we create a list of jobs to run on our system using a utility script `Make_job_list_kallisto.sh` which consumes our tab delimited input file list above. I've used your `-l` and `-s` parameters below, note that with paired end data normally these are estimated for you.

```
#!/bin/bash -eu

[ $# -ne 2 ] && { echo -en "\n*** This script generates jobs for GNU parallel. *** \n\n Error Nothing to do, usage: < input tab delimited list > < output run list file >\n\n" ; exit 1; }
set -o pipefail

# Get command-line args
INPUT_LIST=$1
OUTPUT=$2

# Set counter
COUNT=1
END=$(wc -l $INPUT_LIST | awk '{print $1}')

echo " "
echo " * Input file is: $INPUT_LIST"
echo " * Number of runs: $END"
echo " * Output job list for GNU parallel saved to: $OUTPUT"
echo " "

# Main bit of command-line for job
CMD="kallisto quant -i index.fa.idx --single -l 200 -s 0.1"

# Main Loop
[ -e $OUTPUT ] && rm $OUTPUT
while [ $COUNT -le $END ];
do
    LINE=( $(awk "NR==$COUNT" $INPUT_LIST) )
    # Make file list
    echo "Working on $COUNT of $END Sample ID: ${LINE[0]}, Files ${LINE[@]:1}"
    echo "$CMD -o ${LINE[0]} ${LINE[@]:1}" >> $OUTPUT
    ((COUNT++))
done

```

We use the script like so:

```
 Make_job_list_kallisto.sh input_file_list.txt job_list.txt

```
4. Now we have a job list for GNU parallel we can run our jobs in parallel like so:

```
parallel --progress --jobs 4 --joblog kallisto_joblog.txt < job_list.txt

```

This will run 4 of the kallisto jobs simultaneously (Increase this if you have more threads - don't mind consuming all the IO/CPU threads on your system) producing a nice progress report of those completed, additionally a job log will be written to the file `kallisto_joblog.txt`.

For paired end data I already have these scripts + similar instructions on GitHub which people might find useful.
> 7 votes

# Answer
For generating the index, you should be using one of `-i` or `--index`, but not both:

```
./kallisto index --make-unique -i index.fa.idx index.fa

```

You need to point kallisto to the generated index file (usually has a `.idx` extension), rather than the fasta file:

```
for file in *.fasta;
  do kallisto quant -i index.fa.idx --single -l 200 -s 0.1 -o "${file}-aligned" "${file}";
done

```
> 5 votes

---
Tags: rna-seq
---


title: Workflow of metabarcoding analysis

text:
# Question
Title: Workflow of metabarcoding analysis
For a diet analysis of an insect-eating animal, all species in a sample shall be identified. For this a sequencing of the metagenomic sample was done, where the COI/COX region was used as a barcode and amplified. This resulted in 20 million reads (after quality control).

While looking for an apropriate workflow, I am faced with these issues:

* For many tools 20 million reads are to much to handle (BLAST for instance)
* For the majority of metagenomic tools the focus is on micro-organisms, not insects
* The database to search against must contain COI/COX data, mainly from insects, however, the animal could also eat snails or frog eggs...
* Since the reads are amplicons some analysis are not recommended

Leaving me with these questions:

1. What are the typical analysis steps in a workflow like this?
2. Which tools could I use for my purpose?
3. Where could I get a database of COI/COX data?

# Answer
> 0 votes
I'm working on something similar as a side project. I'll detail the steps that I'm using, but I'm curious if others can suggest a better way.

* You don't mention if your data are paired-end, but in my case I have paired-end sequences that overlap in the middle. I used FLASh to merge the reads and seqtk to convert the reads to fasta format.
* I used the cluster\_fast algorithm in vsearch to cluster sequences at 100% identity, then filtered out rare sequence types (\<10 sequences) to account for sequencing errors.
* This resulted in a manageable number of unique sequences, which I searched against a custom database using vsearch, but which could also have been searched against the nr database using blast.

# Answer
> 0 votes
Regarding your third question ("Where could I get a database of COI/COX data?"), you can apparently search and download fasta sequences from BOLD.

However, I don't see a "select all" button to select all records found by one search in a single fasta file.

There seems to exist an R API: https://cran.r-project.org/web/packages/bold/bold.pdf which refers to this: http://www.boldsystems.org/index.php/resources/api#sequenceParameters\>.

From this I seem to understand that one can directly query using urls:

```
wget -O Gastropoda.fa http://www.boldsystems.org/index.php/API_Public/sequence?taxon=Gastropoda&marker=COI

```

And I get a fasta file containing COI sequences for Gastropoda. There may be ways to refine the query geographically, using the `geo` parameter.

---
Tags: metagenome, barcode
---


# Question
Title: How to append numbers only on duplicates sequence names?
I have a reference database with contains 100s of sequences in fasta format. Some of these sequences have duplicate names like so:

```
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName
atgc
>1_duplicateName
atgc

```

Is is possible to run through a large file like this and change the names of only the duplicates?

```
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName_1
atgc
>1_duplicateName_2
atgc

```

# Answer
Sure, this little Perl snippet should do it:

```
$ perl -pe 's/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/; ' file.fa
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc_2
>1_duplicateName
atgc_3
>1_duplicateName_2
atgc_4

```

Or, to make the changes in the original file, use `-i`:

```
perl -i.bak -pe 's/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/; ' file.fa

```

Note that the first occurrence of a duplicate name isn't changed, the second will become `_2`, the third `_3` etc.

### Explanation

* `perl -pe` : **p**rint each input line after applying the script given by `-e` to it.
* `++$seen{$_}>1` : increment the current value stored in the hash `%seen` for this line (`$_`) by 1 and compare it to `1`.
* `s/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/` : if the current line starts with a `>` and the value stored in the hash `%seen` for this line is greater than 1 (if this isn't the first time we see this line), replace the end of the line (`$`) with a `_` and the current value in the hash

---

Alternatively, here's the same idea in `awk`:

```
$ awk '(/^>/ && s[$0]++){$0=$0"_"s[$0]}1;' file.fa
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName
atgc
>1_duplicateName_2
atgc

```

To make the changes in the original file (assuming you are using GNU `awk` which is the default on most Linux versions), use `-i inplace`:

```
awk -iinplace '(/^>/ && s[$0]++){$0=$0"_"s[$0]}1;' file.fa

```

### Explanation

In awk, the special variable `$0` is the current line.

* `(/^>/ && s[$0]++)` : if this line starts with a `>` and incrementing the value stored in the array `s` for this line by 1 evaluates to true (is greater than 0).
* `$0=$0"_"s[$0]` : make the current line be itself with a `_` and the value from `s` appended.
* `1;` : this is just shorthand for "print this line". If an expression evaluates to true, awk will print the current line. Since `1` is always true, this will print every line.

---

If you want all of the duplicates to be marked, you need to read the file twice. Once to collect the names and a second to mark them:

```
$ awk '{
    if (NR==FNR){
        if(/^>/){
            s[$0]++
        }
        next;
    }
    if(/^>/){
        k[$0]++;
        if(s[$0]>1){
            $0=$0"_"k[$0]
        }
    }
    print
}' file.fa file.fa
>1_uniqueGeneName
atgc
>1_anotherUniqueGeneName
atgc
>1_duplicateName_1
atgc
>1_duplicateName_2
atgc

```

---

**IMPORTANT**: note that all of these approaches assume you don't already have sequence names ending with `_N` where `N` is a number. If your input file has 2 sequences called `foo` and one called `foo_2`, then you will end up with two `foo_2`:

```
$ cat test.fa
>foo_2
actg
>foo
actg
>foo
actg
$ perl -pe 's/$/_$seen{$_}/ if ++$seen{$_}>1 and /^>/; ' test.fa
>foo_2
actg
>foo
actg
>foo_2
actg

```

If this can be an issue for you, use one of the more sophisticated approaches suggested by the other answers.
> 5 votes

# Answer
And here’s a solution using R (with the Bioconductor):

```
fa = ShortRead::readFasta(infile)
ids = as.character(ShortRead::id(fa))
fa@id = Biostrings::BStringSet(make.unique(ids, sep = '_'))
ShortRead::writeFasta(fa, outfile)

```
> 7 votes

# Answer
Sure, using biopython:

```
#!/usr/bin/env python
from Bio import SeqIO

records = set()
of = open("output.fa", "w")
for record in SeqIO.parse("foo.fa", "fasta"):
    ID = record.id
    num = 1
    while ID in records:
        ID = "{}_{}".format(record.id, num)
        num += 1
    records.add(ID)
    record.id = ID
    record.name = ID
    record.description = ID
    SeqIO.write(record, of, "fasta")
of.close()

```

Change `output.fa` and `foo.fa`. One doesn't need to explicitly change the `.name` and `.description`, but that's handy to prevent the original ID from not appearing still (after a space).

# Explanation

* `records = set()` This will create a lookup table of all IDs **written**
* `for record in SeqIO.parse("foo.fa", "fasta"):` Open `foo.fa` as a fasta file and iterate over the **entries** in it. These are objects with an ID, name, description, and sequence attribute (the name and description are the same as the ID if not present).
* `ID = record.id` Memoize the entry ID (e.g., `1_duplicateName`)
* `while ID in records:` As long as the ID has already been seen, keep looping.
* `ID = "{}_{}".format(record.id, num)` Start adding increasing numbers after the ID, such as `1_duplicateName_1` and `1_duplicateName_2`. This will continue until the ID has not been seen.
* `records.add(ID)` Add the unseen ID to the set.
* `record.id = ID` Update the ID, the `.name` and `.description` are the same. If you don't do that, then you get output like `>1_duplicateName_1 1_duplicateName`. That's not really a problem, but it's excessive.
* `SeqIO.write(record, of, "fasta")` Write the record to the output file in fasta format.

A benefit of biopython is the easy ability to change formats, so instead of fasta one could have substituted Genbank or another format if needed.
> 6 votes

---
Tags: fasta, text
---


# Question
Title: Handling sample identity in aggregated 10x libraries?
`cellranger aggr` can combine multiple libraries (samples), and appends each barcode with an integer (e.g. AGACCATTGAGACTTA-1). The sample identity is not recorded in the combined `matrix.mtx` file.

https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/aggregate#gem\_groups

How can we keep and assign sample information to each cell after loading the data into `R`? (e.g. `Seurat::Read10X()`)

# Answer
> 1 votes
### An example solution for Seurat:

Retrieve sample IDs from the .csv used with `cellranger`:

```
samples <- read.csv(file.path("/path/to/csv", "nameof.csv"), stringsAsFactors=F)

```

Load the 10x dataset and initialize the Seurat object:

```
cells.data <- Read10X("path/to/filtered_gene_bc_matrices")

cells <- new("seurat", raw.data=cells.data)

```

Get barcodes and suffix:

```
cellcodes <- as.data.frame(cells@raw.data@Dimnames[[2]])
colnames(cellcodes) <- "barcodes"
rownames(cellcodes) <- cellcodes$barcodes

cellcodes$libcodes <- as.factor(gsub(pattern=".+-", replacement="", cellcodes$barcodes))
cellcodes$samples <- as.vector(samples$library_id[cellcodes$libcodes])

```

Create dataframe for `meta.data` argument and set up object:

```
sampleidentity <- cellcodes["samples"]

cells <- Setup(cells,
               meta.data=sampleidentity,
               min.cells=3,
               min.genes=200,
               do.logNormalize=T, total.expr=1e4, project="projectname")

```

Use the `group.by` argument of various functions.

# Answer
> 1 votes
You have to keep track of the file order you used for `aggr`. The suffix number represents which represents the aggregated sample. You could store this information in a text file and load it into R independently of the `Read10X` function and combine it with the cell names to get their sample names. From the page you link:

> Gem Groups
>
> This is an integer that is appended to each barcode in the gene-barcode matrix. For example, *AGACCATTGAGACTTA-1* and *AGACCATTGAGACTTA-2* are distinct cell barcodes from different libraries, despite having the same nucleotide sequence.
>
> The numbering of the GEM groups will reflect the order that the libraries were provided in the Aggregation CSV.

---
Tags: rna-seq, r, scrnaseq, 10x-genomics
---


# Question
Title: How to obtain .bed file with coordinates of all genes
I want to get a .bed file with the genes' names and canonical coordinates, also I would like to have coordinates of exons, too. I can get the list from UCSC, however, if I choose UCSC Genes - knownCanonical, I can not extract coordinates of exons. If I use other options - I am getting coordinates of as many transcriptional isoforms as were detected while I need only one canonical form.

How can I get such BED file?

# Answer
> 16 votes
Via Gencode and BEDOPS `convert2bed`:

```
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff - \
    > genes.bed

```

You can modify the `awk` statement to get exons, by replacing `gene` with `exon`.

To get HGNC symbol names in the ID field, you can add the `--attribute-key="gene_name"` option to v2.4.40 or later of `convert2bed`. This slight modification extracts the `gene_name` attribute from the annotation record and puts it in the fourth (ID) column:

```
$ wget -qO- ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "gene"' - \
    | convert2bed -i gff --attribute-key="gene_name" - \
    > genes.bed

```

BEDOPS: https://github.com/bedops/bedops

This is based off an answer I wrote on Biostars, which includes a Perl script for generating a BED file of introns from gene and exon annotations: https://www.biostars.org/p/124515/#124522

# Answer
> 4 votes
You can indeed get this information from the UCSC Table Browser. Select knownGene as your primary table, make a filter, add knownCanonical to as a linked table to filter on, then in the free-form query section add "1" without the quotes. Then click submit and select Bed output, where you can choose "exons plus" as an output option. This will lead to output like the following:

```
chr1    17368   17436   uc031tla.1_exon_0_0_chr1_17369_r    0   -
chr1    29553   30039   uc057aty.1_exon_0_0_chr1_29554_f    0   +
chr1    30563   30667   uc057aty.1_exon_1_0_chr1_30564_f    0   +
chr1    30975   31097   uc057aty.1_exon_2_0_chr1_30976_f    0   +
chr1    30365   30503   uc031tlb.1_exon_0_0_chr1_30366_f    0   +
chr1    34553   35174   uc001aak.4_exon_0_0_chr1_34554_r    0   -
chr1    35276   35481   uc001aak.4_exon_1_0_chr1_35277_r    0   -
chr1    35720   36081   uc001aak.4_exon_2_0_chr1_35721_r    0   -

```

For more step-by-step information, please check the answer from the mailing list archives here (in particular the Exon Method 1 section): https://groups.google.com/a/soe.ucsc.edu/d/msg/genome/BJ-6DlaZNCY/grgGIpuJAwAJ

Also note that UCSC provides several mailing lists for support, if you have futher questions please send them there:

* General questions: genome@soe.ucsc.edu
* Questions involving private data: genome-www@soe.ucsc.edu
* Questions involving mirror sites: genome-mirror@ose.ucsc.edu

# Answer
> 2 votes
You can download a list of transcript annotations as a flat file from UCSC:

http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz

It is not exactly a BED file, but it does contain information about the known transcripts for this assembly (hg19, in this case):

```
585 NR_046018   chr1    +   11873   14409   14409   14409   3   11873,12612,13220,  12227,12721,14409,  0   DDX11L1 unk unk -1,-1,-1,
585 NR_024540   chr1    -   14361   29370   29370   29370   11  14361,14969,15795,16606,16857,17232,17605,17914,18267,24737,29320,  14829,15038,15947,16765,17055,17368,17742,18061,18366,24891,29370,  0   WASH7P  unk unk -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,

```

The gene name can be found in column 13. The exon start and end positions are in columns 10 and 11, respectively. More information can be found in the other columns:

1. bin
2. name
3. chrom
4. strand
5. txStart
6. txEnd
7. cdsStart
8. cdsEnd
9. exonCount
10. exonStarts
11. exonEnds
12. score
13. name2
14. cdsStartStat
15. cdsEndStat
16. exonFrames

(The column information was obtained from the sql file here)

---
Tags: bed, public-databases
---


# Question
Title: Genome databases of microbial eukaryotes
I am interested in collecting the reference genomes of marine protists and would like to know which databases can be used for such tasks, and is there a single, centralized database of protistological genomic data?

# Answer
> 2 votes
You can find many protist genomes on ensembl

http://protists.ensembl.org/index.html

---
Tags: database, genome
---


# Question
Title: How to isolate genes from whole genomes for phylogenetic tree analysis?
I have 446 whole Klebsiella Pneumoniae genomes I want to build a phylogenetic tree from. After reading about constructing phylogenetic trees it seems the only option for large numbers of genomes is to isolate a gene with low variability from generation to generation and use this gene to build a tree. For example Lars Jensen recommends using "16S rRNA \[or\] all ribosomal-protein-coding genes" https://www.biostars.org/p/1930/. What program isolates these genes of interest from the whole genome fasta files and can put them into a multiple alignment file? Or outputs them in a formate ready for a multiple alignment program such as Muave? The reason I say a multiple alignment file is because this it the type of file most phylogenetic tree programs take (I.E. clonalframe).

# Answer
There are a lot of ways to do this. I suggest using Prokka/Roary to produce a core genome alignment. There's a useful tutorial on the Roary website:

```
for file in *
do
    prokka --kingdom Bacteria --outdir "${file%%.*}"  --genus Listeria --locustag "${file%%.*}" "$file"
    mv "${file%%.*}"/PROKKA_07052017.gff GFF/"${file%%.*}".gff # use current date
done
roary -f Alignment -e -n -v GFF/*.gff

```

Alignment/core\_gene\_alignment.aln can be used as input for phylogenetic analyses
> 3 votes

# Answer
Extract desired gene sequences using standalone blast

Simply provide a reference database with your desired output. Set up your command and away you go. You can set the search up with a for loop for a batch of sequences. Command may look like

```
for f in *.fasta; do
   f=$(basename $f .fasta)
   blastn \
   -outfmt "6 sseqid qseq %" \
   -query $f.fasta \
   -subject reference.fna \
   > out/$f.fas
done

```

Watch the output as blast will spit out the gene as detected in + or - sense. If you want to only gather positive sense use `-strand`option. The default output I have here is tab output which requires a few sed commands to make into fasta.

```
sed -i \
-e 's/\s*$//g' \
-e 's/^/>/g'  \
-e 's/\s\+/\n/g' \
*.fas

```

Online alignment servers are an easy way to align small datasets e.g. EBI
> 3 votes

---
Tags: gene, phylogenetics, phylogeny
---


# Question
Title: PASA pipeline: compare experimental transcripts to the reference annotation
I would like to ask if anyone has experience in running a subset of the PASA pipeline, in particular for the reconciliation of some experimental 'transcripts' with the reference annotation.

In more detail, I am working with RNA-seq data from *D. melanogaster*. I have reconstructed the 'transcripts' using Trinity. I have aligned these 'transcripts' to the reference genome using GMAP. Now I would like to match these 'experimental transcripts' with the reference annotation, to see how they compare.

I was wondering if I can run just a subset of the whole PASA pipeline, so basically skipping the step of alignment to the reference and using the alignment file that I generated externally.

Also, I would be glad to receive more general redirecting to any released software allowing the comparison of experimental transcripts to the reference annotation.

Cross-posted on the PASA Google Group pasapipeline-users

Cross-posted on biostars

# Answer
Answer from the authors:

> PASA ends up using gff3 format instead of sam format for uploading the alignments. If somebody has the gmap gff3 output, then can be done to just upload that directly using the custom alignment importer, but not with sam. They suggest just having PASA rerun GMAP as part of its regular routine.
> 2 votes

---
Tags: rna-seq, assembly, transcriptome, isoform
---


# Question
Title: using python to write bioinformatics pipelines tutorial
I was wondering if there is a tutorial or a small code snippet to understand how to write bioinformatics pipeline using python, for example

1. use a aligner (say hisat)
2. get the output and process it using samtools

I was able to use subprocess from python2.7 for this purpose using samtools but i am not able to link both the processes.i.e given path(which I can use argparse) for directory with fastq files output would be processed bam.

sample code for samtools sam to bam :

```
import subprocess
subprocess.run(['samtools','view', '-bS',
                '../some.sam',
                '>',
                '../some.bam'])

```

# Answer
> 4 votes
I agree that using a specialized tool is probably a good idea.

Nevertheless if you want to stick with Python, I suggest using `plumbum` instead of `subprocess`. It has a very nice syntax for this kind of problems.

```
from plumbum import local

# load a command
samtools = local['samtools']

# you can then easily redirect
(samtools['view', '-bS', '../some.sam'] > '../some.bam')()

# it is also easy to pipe
hisat2 = local['hisat2']
chain = hisat2['-x', 'genome.index', '-1', 'reads_R1.fq.gz', '-2', 'reads_R2.fq.gz'] | samtools['sort'] > 'reads_vs_genome.bam'

chain()

```

# Answer
> 13 votes
Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and runs it robustly.

* **Snakemake** is my favourite. It's (nearly) pure Python and can generate reports.
* **Nextflow** is growing in popularity and is pretty straightforward
* **Ruffus** used to be reasonably popular, seemed fine when I used it
* **Bpipe** is for bioinformatics
* **Airflow** is a more industrial "big" solution but in wide use
* See this big list of pipeline systems

I'm sure you can find something there.

# Answer
> 5 votes
BioPython has some good tools for processing reads and alignments. http://biopython.org/DIST/docs/tutorial/Tutorial.html

There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls https://pysam.readthedocs.io/en/latest/

I would use subprocess to call the aligner and specify the output to a bam file that you have named and then read that bam file with pysam to do the analysis that you are interested in.

If you provide more specifics about the analysis or aligners I can help you more but you are in the right track just need to connect the two calls together.

example:

```
import subprocess
import pysam
import os

for root, dirs, files in os.walk(rootPath):
    for filename in fnmatch.filter(files, ".fastq"):

        subprocess.run([`hisat`, `-f`, filename `-o`, `some.bam`])

        #get depth
        depth = pysam.depth('some.bam')

```

# Answer
> 4 votes
If you're just doing alignment and conversion to a sorted BAM file, there's no need to run it through python. A simple pipe on the unix command line works just as well (and probably runs faster):

```
hisat2 -x genome.index -1 reads_R1.fq.gz -2 reads_R2.fq.gz |
  samtools sort > reads_vs_genome.bam

```

# Answer
> 0 votes
Have a look at luigi for writing pipelines in python. You can found BioLuigi if you like :-P

---
Tags: rna-seq, samtools, python, genomics
---


# Question
Title: scanBam from Rsamtools is not importing one of my reads into R
I have this read in my BAM file. It maps on chromosome 1.

I open this BAM file in IGV, and I can see the alignment on chromosome 1.

But when I open this file in R with Rsamtools:

```
bamContigsCel <- Rsamtools::scanBam('output/alignment/pacbio/bwa/ref     /bristolAssemblySorted.bam', param = Rsamtools::ScanBamParam(what = Rsamtools::scanBamWhat(), flag = Rsamtools::scanBamFlag(isMinusStrand = FALSE), tag = bamTags))[[1]]

```

I then check if the read maps to chromosome 1 in my R object but I cannot find it.

```
bamContigsCel$rname[bamContigsCel$qname == '000000F|arrow']
[1] II II IV
Levels: I II III IV MtDNA V X

```

But if I look at the BAM file, it's there. Why isn't Rsamtools importing my read into R?

```
$ samtools view bristolAssemblySorted.bam | grep -n '000000F' | head  -c80
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   0 0   GAAGCTGTCTAAACTTTGGC

```

Cross-posted on biostars

# Answer
Note the flag; that read is mapped in a reverse-complemented manner, so `isMinusStrand = FALSE` is filtering it out. I tested this by making a BAM file with only that read:

```
@SQ SN:I    LN:1000000
000000F|arrow   2064    I   336331  7   2926310H260M1774267H    *   00  *   *

```

Then in R:

```
> library(Rsamtools)
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = TRUE)))[[1]]
> blah$qname
[1] "000000F|arrow"
> blah = scanBam("foo.bam", param=ScanBamParam(what=scanBamWhat(), flag=scanBamFlag(isMinusStrand = FALSE)))[[1]]
> blah$qname
character(0)

```
> 4 votes

---
Tags: r, sam, samtools
---


# Question
Title: COG Annotation - Dealing with genes assigned to two or more COG categories
I am dealing with a list of genes which have been selected from a gene enrichment analysis. In order to see what kind of genes are overrepresented, I ran eggnog-mapper to do an orthology assignment of each gene against a bacterial database (provided in eggnog).

After running eggnog you get a .csv file containing for each gene which is the closest ortholog, and a list of GO, KEGG and COG terms associated.

From there, I would like to do a Fisher-test (also considering Bonferroni correction) to highlight COG categories overrepresented in our dataset compared to our reference genome.

Finally my question: How do you deal with genes assigned to multiple COG categories (same would apply for GO terms). You count +1 to each of the COG categories assigned? Do you just choose one category per gene?

# Answer
> 3 votes
You add 1 to the count of each COG category. You are looking for over represented *COG categories* so you must count them all. Many genes will be assigned to multiple categories. In fact, a majority of human genes are given multiple GO cellular function annotations. I have worked considerably in this field and one of the surprises I found is how common protein multifunctionality is. But I digress.

The point here is that you want to see whether any COG categories are over represented in your test set. Whether a gene is annotated to one or to a hundred categories is irrelevant. You are counting genes per COG, so you need to count each of the COGs assigned to a gene separately.

---
Tags: sequence-annotation
---


# Question
Title: Y Chromosome Aligned Reads in scATAC-seq data from a female-derived cell line?
I'm working with scATAC-Seq data on the K562 cell line, which is supposed to be derived from a female patient. While following the scATAC-seq data analysis pipeline, after performing bowtie alignment they recommend filtering out all reads aligned to mitochondrial DNA and on the Y chromosome.

Out of curiosity, I decided to count how many reads aligned to chrY and found that it can be quite high - chrY has approximately 10% as many reads as a similarly-sized chromosome (chr19 in hg19).

Is it normal to find reads mapping to chrY when the cell line is female-derived, or did I mess up somehow?

# Answer
> 7 votes
There are homologous regions between X an Y chromosomes: https://en.wikipedia.org/wiki/Pseudoautosomal\_region

It is therefore normal to have some female-derived reads mapping in Y chromosome.

You should probably check what proportion of such reads fall in other parts of the Y chromosome than pseudoautosomal regions.

---
Tags: sequence-alignment
---


# Question
Title: Software to produce a table of post-translational modifications from a peptide list
Does anyone know if there is a program/library/script in R or Python that takes as input a list of proteins/peptides and a list of post-translational modifications (PTMs; like oxidation of methionine and acetylation of cysteine), and as output returns a matrix/list with the imported peptides and all possible locations of the selected modifications?

# Answer
> 2 votes
I don't know of a tool that can do exactly what you want. However, there are various tools that can predict specific classes of post translational modifications. For example:

1. Phosphorylation: GPS

> Computational prediction of phosphorylation sites with their cognate protein kinases (PKs) is greatly helpful for further experimental design. Although ~10 online predictors were developed, the PK classification and control of false positive rate (FPR) were not well addressed. Here we adopted a well-established rule to classify PKs into a hierarchical structure with four levels.
2. Mitochondrial import: Mitoprot

> MitoProt calculates the N-terminal protein region that can support a Mitochondrial Targeting Sequence and the cleavage site. A complete description of the method to make the prediction is available in: M.G. Claros, P. Vincens. Computational method to predict mitochondrially imported proteins and their targeting sequences. Eur. J. Biochem. 241, 779-786 (1996).
3. Acetylation: "PAIL: Prediction of Acetylation on Internal Lysines"

> In this work, we present a novel online predictor for protein acetylation sites prediction of PAIL, Prediction of Acetylation on Internal Lysines. We have manually mined scientific literature to collect 249 experimentally verified acetylation sites of 92 distinct proteins. Then the BDM (Bayesian Discriminant Method) algorithm has been employed.

You can use tools like the above and then combine the results to get the sort of matrix you are asking for.

---
Tags: r, proteins, python
---


# Question
Title: How can I compute gene expression for a set of RNA reads?
I'm trying to compute a gene expression profile for an organism. I have gene nucleotide sequences of the mentioned organism stored in a fasta file and a set of paired reads stored in two separate files with the same fasta format. Now I want to compute the coverage of every gene by the given reads. That is, on average how many reads cover a given segment of the sequence.

Is there any tool for this specific purpose? I'm trying to use the velvet sequencer but I have encountered a few problems.

# Answer
> 8 votes
It's unclear if your paired-end reads are actually in fasta format, I'll presume that they're in fastq instead.

The easiest tool to use is salmon, which nicely deals with things like multimapping. If you're trying to judge the quality of an assembly, then I recommend having a look at transrate, which uses some related methodology for assessing contig quality.

# Answer
> 3 votes
In addition to salmon that was already mentioned, you can also try kallisto or RSEM, which are also fairly popular/respectable and will work with a transcriptome FASTA.

# Answer
> 2 votes
It's not possible to compute absolute expression from RNASeq reads if they are processed in the usual way, where a sequencer produces the same number of reads regardless of the input RNA amount. At best, RNASeq will give you an indication of proportional expression within a single sample. For this reason, relative expression (i.e. that used by differential expression tests) is easier to determine than absolute expression.

The closest approximation to absolute expression is to generate an expression relative to the average expression of a set of housekeeping genes, but I don't think there's a universal set that has been decided on. Gene expression, even for common housekeeping genes, can vary depending on the environmental conditions of the cell. For example, GAPDK is involved in immune cell activation.

However, as long as experimental conditions are similar, and you're not planning on looking for statistical significance, the proportional expression can still give qualitative insights into how cell populations behave in relation to other populations. DESeq2 provides a variance-stabilising transformation function that minimises variation for small-count genes, assuming that each sample has roughly the same total expression. I have found that I get better outcomes / comparisons from this transformation when carrying out a further adjustment to account for gene length (i.e. divide by the length of the longest transcript for each gene). See our Th2 paper, section "Read mapping and differential expression analysis" for more information. The "transcripts-per-million" values produced by Kallisto and Salmon provide measures similar to this.

If, on the other hand, you were able to modify the experimental design, single cell sequencing (or "known cell count" sequencing) can be used for determining absolute expression: use a spiked-in transcript that is added in proportion to the *cell count*, so that results can be compared in proportion to the expression for that transcript.

---
Tags: rna-seq, rna-alignment
---


# Question
Title: Selecting sites from VCF which have an alt AD > 10
I have high-depth variant calling created using the HaplotypeCaller with `--output_mode EMIT_ALL_SITES` I'm interested in finding all sites (regardless of genotype call heterozygous or homozygous) where at least one of the **alternative** alleles have an `AD` value (Allelic Depth) greater than 10, *I*.*e*. are supported by more than 10 reads. Also ideally I want back more than just the first alternative allele. Note that I don't want back lines of VCF were we only see an AD count for the ref allele only.

So in the example VCF snippet below I'm wanting to select lines: 6,7,8,12,13 and 14, which have GT:AD values `1/1:1,988:989` `0/1:116,92` `0/1:220,234` `0/1:62,611` `1/1:0,109` respectively.

```
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  12908_DIAG
3   187446740   .   T   .   Infinity    .   AN=2;DP=1095;MQ=60.00   GT:AD:DP    0/0:1095:1095
3   187446741   .   C   .   Infinity    .   AN=2;DP=1117;MQ=60.00   GT:AD:DP    0/0:1117:1117
3   187446752   .   A   .   Infinity    .   AN=2;DP=1297;MQ=60.00   GT:AD:DP    0/0:1297:1297
3   187446763   .   C   .   Infinity    .   AN=2;DP=1494;MQ=60.00   GT:AD:DP    0/0:1494:1494
3   187451574   .   C   .   Infinity    .   AN=2;DP=1493;MQ=60.00   GT:AD:DP    0/0:1493:1493
3   187451609   rs1880101   A   G   39794.03    .   AC=2;AF=1.00;AN=2;BaseQRankSum=1.859;ClippingRankSum=0.000;DB;DP=995;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;MQRankSum=0.000;QD=24.56;ReadPosRankSum=0.406;SOR=8.234  GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0.000;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0.00;MQ=59.97;MQRankSum=0.000;ReadPosRankSum=-2.523;SOR=6.357   GT:AD:DP:GQ:PL  0/0:211,23:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0.000;DB;DP=215;ExcessHet=3.0103;FS=1.110;MLEAC=1;MLEAF=0.500;MQ=59.97;MQRankSum=0.000;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809  GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803671 .   C   A   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-0.880;ClippingRankSum=0.000;DP=450;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.953;SOR=0.572    GT:AD:DP:GQ:PL  0/0:445,2:447:99:0,1272,15958
4   1803681 .   T   C   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.654;ClippingRankSum=0.000;DP=483;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=-0.422;SOR=0.664    GT:AD:DP:GQ:PL  0/0:479,2:481:99:0,1408,18538
4   1803703 .   A   G   0   .   AC=0;AF=0.00;AN=2;BaseQRankSum=-1.704;ClippingRankSum=0.000;DP=458;ExcessHet=3.0103;FS=0.000;MLEAC=0;MLEAF=0.00;MQ=60.00;MQRankSum=0.000;ReadPosRankSum=0.299;SOR=0.497 GT:AD:DP:GQ:PL  0/0:454,2:456:99:0,1325,18095
4   1803704 rs2234909   T   C   6676.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0.000;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=14.71;ReadPosRankSum=0.324;SOR=0.849   GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.60 .   AC=1;AF=0.500;AN=2;BaseQRankSum=8.083;ClippingRankSum=0.000;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.500;MQ=60.00;MQRankSum=0.000;QD=16.51;ReadPosRankSum=0.180;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 .   AC=2;AF=1.00;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=29.22;SOR=9.401    GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

```

A dropbox link for file

I'd initially considered using GATK's SelectVariants but I'm not sure JEXL has the ability to select out what I want specifically other than a blanket AD \> 10 which will give me both ref and alt alleles with AD \> 10. Perhaps there is a bioawk solution or something more elaborate with coreutils which could successfully return sites with an alt AD count \> 10?

# Answer
> 4 votes
This now works with the development version of Bcftools v1.5 (commit 4f134df). Thanks to Petr Danecek for adding the feature. I expect this feature to make its way into the next release of Bcftools:

```
git clone git://github.com/samtools/htslib.git
git clone git://github.com/samtools/bcftools.git
(cd bcftools; make)

bgzip Test.vcf
./bcftools/bcftools index Test.vcf.gz
./bcftools/bcftools filter -i 'AD[0:1-] > 10' Test.vcf.gz

```

Output without header (I have modified the second line to be tri-allelic to demonstrate the filtering works):

```
3   187451609   rs1880101   A   G   39794   PASS    AC=2;AF=1;AN=2;BaseQRankSum=1.859;ClippingRankSum=0;DB;DP=995;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;MQRankSum=0;QD=24.56;ReadPosRankSum=0.406;SOR=8.234   GT:AD:DP:GQ:PL  1/1:1,988:989:99:39808,2949,0
4   1803279 .   T   G   0   PASS    AC=0;AF=0;AN=2;BaseQRankSum=-6.652;ClippingRankSum=0;DP=245;ExcessHet=3.0103;FS=89.753;MLEAC=0;MLEAF=0;MQ=59.97;MQRankSum=0;ReadPosRankSum=-2.523;SOR=6.357 GT:AD:DP:GQ:PL  0/0/0:211,3,34:234:99:0,364,6739
4   1803307 rs2305183   T   C   2486.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-5.049;ClippingRankSum=0;DB;DP=215;ExcessHet=3.0103;FS=1.11;MLEAC=1;MLEAF=0.5;MQ=59.97;MQRankSum=0;QD=11.95;ReadPosRankSum=-0.045;SOR=0.809   GT:AD:DP:GQ:PL  0/1:116,92:208:99:2494,0,3673
4   1803704 rs2234909   T   C   6676.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-2.605;ClippingRankSum=0;DB;DP=456;ExcessHet=3.0103;FS=1.753;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=14.71;ReadPosRankSum=0.324;SOR=0.849  GT:AD:DP:GQ:PL  0/1:220,234:454:99:6684,0,6366
4   1803824 rs2305184   C   G   2030.6  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=8.083;ClippingRankSum=0;DB;DP=124;ExcessHet=3.0103;FS=6.128;MLEAC=1;MLEAF=0.5;MQ=60;MQRankSum=0;QD=16.51;ReadPosRankSum=0.18;SOR=0.096    GT:AD:DP:GQ:PL  0/1:62,61:123:99:2038,0,1766
4   1805296 rs3135883   G   A   3876.03 PASS    AC=2;AF=1;AN=2;DB;DP=110;ExcessHet=3.0103;FS=0;MLEAC=2;MLEAF=1;MQ=60;QD=29.22;SOR=9.401 GT:AD:DP:GQ:PL  1/1:0,109:109:99:3890,326,0

```

**Edited** (Dec 2020): bcftools filter syntax has been updated, answer now reflects changes found in version v1.11 (not sure at what point these were introduced). The old query `'AD[1-] > 10'` would produce the error `The FORMAT tag AD can have multiple subfields, run as AD[sample:subfield]`. A discussion of these changes can be found here.

# Answer
> 5 votes
using vcfilterjs

and the following script:

```
function accept(vc)
    {
    var i,j;
    for(i=0;i< vc.getNSamples();++i)
        {
        var genotype = vc.getGenotype(i);
        if(!genotype.hasAD()) continue;
        var ad = genotype.getAD();
        /* loop over AD starting from '1' =  first ALT */
        for(j=1  ;j< ad.length ;++j)
            {
            if(ad[j]>10) return true;
            }
        }
    return false;
    }
accept(variant);

```

usage:

```
java -jar dist/vcffilterjs.jar -f script.js Test.vcf

```

# Answer
> 4 votes
You can do this in Hail:

```
from hail import *
hc = HailContext()
(hc.import_vcf('test.vcf')
 .filter_variants_expr('gs.exists(g => g.ad[1:].exists(d => d > 10))')
 .export_vcf('filtered.vcf'))

```

This works with any number of samples and will keep the variants where at least one sample has a genotype with an alternate allele support by more than 10 reads.

To verify we got the expected 6 variants:

```
>>> hc.import_vcf('filtered.vcf').count()
(1L, 6L)

```

`count` returns the number of samples (1) and number of variants (6).

Take a look at the getting started page or tutorials if you want to try it out!

---
Tags: variant-calling, vcf, bioawk
---


# Question
Title: Using data mining of papers in order to derive genomic connections
Are there any examples of researchers data mining articles and papers (e.g., from pubmed or google scholar) in order to derive relations between genes and diseases?

# Answer
> 4 votes
There are many databases that have used publication scraping for oncogenic gene fusions. There are publications for the individual methods they used for their scraping and aggregation.

# Answer
> 3 votes
This existed as a closed silo, at least in 2015. Qiagen has a team of hired students and Post-Docs for collating research papers into their Knowledge Base, an extensive database that is integrated into a few of their commercial products. Qiagen's claim is that by providing a consistently-structured and well-formatted database, the process of research discovery is accelerated. I've got a few notes from a seminar they gave us a couple of years ago:

* The curation is better than a free / public database because free databases typically have wrong or inconsistent data
* Their knowledge base is an army of hired MDs and PhDs, 5.3M "findings" as of Nov 05 2015
* Scientific articles have unstructured data, so expert knowledge is required to process the data into a "finding" for their knowledge base

I was already doubtful about their approach, and this doubt was magnified when they started getting into the details of their process. I started to lose trust and write my own thoughts about it:

* The knowledge base they have created is private
* The approach is to use a function (i.e. MDs/PhDs) to convert research articles into scientifically-validated associations. However, the function is effectively a black box; Qiagen allow the experts to decide for themselves how a particular article should be converted into a finding.
* The findings in the knowledge base have no strength or likelihood.
* The Experts were encouraged to only enter findings that were "generally accepted by the scientific community. but it seemed like only a single expert was required to get a finding into the knowledge base (i.e. findings weren't curated by other people).
* When using their software (e.g. IPA), p-values of 10^-15 and 10^-22 are reported; these numbers are meaningless.
* The team didn't look at methylation, phosphorylation, etc.; plans are in place to hire additional experts to look at this at a later date.

After listening to their seminar, I realised that Qiagen's idea of a general database that attempts to capture everything isn't going to work well. Experts will disagree about whether a particular association is valid, relevance (or correctness) will differ depending on the area of interest, and the established correctness of a piece of information can change over time.

It would be possible to create something that might work by adding a lot more metadata to all the "findings" in such a database, but it would need to be a huge, public, collaborative effort. Something on the scale of SRA/ENA, with encouragement from the scientific community to enter details into the database prior to publications being accepted/indexed. If the people who wrote papers entered their findings in such a database, other researchers wouldn't need to be hired to try to understand papers and guess at the findings of published research.

---
Tags: gene, genome, cancer, data-mining
---


# Question
Title: Is Wikidata a reliable KB for genomic and drugs annotations?
I'm looking for a Linked Open Data approach to annotate a dataset with human genome and drugs information.

According to the Linked Open Data cloud, there are a lot of interconnected `RDF` vocabularies devoted to life sciences but, unfortunately, a lot of them seem out of date.

I recently found Wikidata:

> Wikidata is a free and open knowledge base that can be read and edited by both humans and machines.

It seems exactly what Linked Data should be and it seems continuously updated. At a first sight, it covers very well genomic and drugs knowledge.

Does Anyone know if Wikidata is a reliable source for genomic and drugs annotation?

# Answer
> 1 votes
The question is what you mean with "reliable". For bacterial genome annotations there is one Wikidata based platform called WikiGenomes (publications).

---
Tags: sequence-annotation, public-databases, data-management
---


# Question
Title: alternatives to MEDIPS to analyse MeDIP datasets
MEDIPS is an established tool with functions for the quality control and analysis of data derived from immunoprecipitation (IP)-seq samples, like Methylation IP sequencing datasets.

I would like to know if there are any other tools I should consider nowadays that could do a better job and hence replace MEDIPS for the analysis of MeDIP datasets.

# Answer
> 2 votes
I guess QSEA should be an obvious answer:

> QSEA (quantitative sequencing enrichment analysis) was developed as the successor of the MEDIPS package for analyzing data derived from methylated DNA immunoprecipitation (MeDIP) experiments followed by sequencing (MeDIP-seq).

See http://bioconductor.org/packages/release/bioc/vignettes/qsea/inst/doc/qsea\_tutorial.html

---
Tags: methylation
---


# Question
Title: Where can I find a database that has phenotype information together with associated SNPs?
I have been searching for a long time and the furthest I got is some database with the functional description for genes. Then, I have to parse these descriptions manually to figure out the association. Moreover, for some species, there are not even any functional descriptions.

I need this information for mice, rice, and Arabidopsis Thaliana. Can anyone help me?

# Answer
> 8 votes
International Mouse Phenotyping Consortium is building a database of phenotypes and knock-outs of mouse. I believe that this database will be fairly complete (20000 knock-outs), but these are knock-outs, not SNPs... There are several mouse GWAS studies, but I am not aware of a database that would pull all the results together.

*Arabidopsis* big GWAS project is called 1001 genomes project. A lot more information about *Arabidopsis* from individual experiments is collected and maintained by The Arabidopsis Information Resource, however this resource is subscription based, therefore what you need might not be available for free.

As mentioned, SNPs of rice were identified and linked to different rice groups in The 3,000 rice genomes project. There is also GWAS done on 512 individuals, not sure if it is sufficient for reliable conclusions (even they claim that they found supersigificant links).

Is this what you were searching for?

# Answer
> 4 votes
A few days ago, I was trying to find some GWAS datasets to download. I hope this site for 3000 rice will help: http://snp-seek.irri.org/. You can click Genotype in the index page for "Query for SNPs from the 3000 genome project".

**UPDATE:** With the help of user manual, I can get IRIS ID from Genotype. According to this paper, iris id may associate with some phenotype information. But the website of IRIS is under maintenance, I can not get more information about phenotype.

# Answer
> 1 votes
Have you looked at Ensembl? At least for mice, they have pretty rich data for SNPs e.g. http://www.ensembl.org/Mus\_musculus/Variation/Explore?db=core;r=2:3225281-3226281;v=rs27096498;vdb=variation;vf=802773

Assumes you have known dbSNP ids (or 'rs' numbers) for your SNPs. If not you could use their VEP tool to annotate them: http://www.ensembl.org/Homo\_sapiens/Tools/VEP?db=core

VEP is also available for rice and Arabidopsis: http://plants.ensembl.org/Oryza\_sativa/Tools/VEP?db=core

---
Tags: gene, gwas
---


# Question
Title: Converting gene names from one public database format to another
*This is a question from /u/apivan19 on reddit. The original post can be found here.*

I have some proteomics data that was given to me with the UniProt gene identifiers in column 1. I've been trying to convert these to normal gene symbols using various programs, but it is proving to be difficult.

The Uniprot website does it fairly decently, but it is not able to convert all of them and then adds some unknown genes into my list.

For example, I will give it 5439 genes in UniProt notation, and will say "5420 of 5439 UniProt identifiers have been converted to 5450 gene symbols"... which is ridiculous.

I tried using David to change the symbols, but it returns them to me in some ridiculous, random order and there's no way I can sort... actually there might be but it'll take a second.

What are some of the easiest ways to do this? It's already very time consuming and am looking for simpler solutions

# Answer
> 9 votes
I tend to use Ensembl Biomart for such queries since there are APIs for various programming languages, e.g. biomaRt, and, maybe more interestingly, via a REST API (although it’s a pretty terrible one).

To translate identifiers from different databases, proceed as follows:

1. Choose database “Ensembl genes”
2. Choose dataset *your desired oganism*
3. Go on “Filters” › “Gene:” › “Input external reference ID list”
    1. Select the chosen source database
    2. Provide a list of IDs, delimited by newline
4. Go to “Attributes” › “Gene:” › untick “Transcript stable ID”
    1. If Ensembl IDs are desired, leave “Gene stable ID” ticked …
    2. Otherwise untick it; go to “External:”, tick your desired identifier format
5. Click “Results” at the top left. This gives a preview that can be exported into various formats; alternatively the top centre buttons “XML” and “Perl” provide the query in XML (for SOAP/REST requests) and as a (horrendously formatted) executable Perl script.

# Answer
> 8 votes
If you're comfortable doing a little programming, check out mygene.info (web services for gene annotations of all sorts). ID translation is specifically one of the use cases addressed in the bioconductor client (see the vignette), and there is a python client as well available through pypi. The documentation for mygene can be found here.

# Answer
> 4 votes
You can do the same using AnnotationDbi Package from Bioconductor. Download the organism specific annotation file like org.Mm.eg.db for mouse and map current gene ids to the gene names/gene symbols.

# Answer
> 3 votes
My favourite gene database conversion site is db2db. You provide a list of IDs in one of a large number of different public formats, and can select one or more IDs as translation targets. It will then walk through various known paths to do the translation, picking what it determines to be the most reliable route to get the information you have asked for. Results appear in the browser as a table, but can also be exported as an Excel file, or as a tab-separated text file.

Note that the mapping of genes from one database to another is not a one-to-one mapping. It is likely the case that there will be some genes in the source database that map to multiple genes in the target database (and *vice versa*), and some genes that aren't present in the target database. These phenomena probably account for the "ridiculous" results that have been seen here.

# Answer
> 3 votes
I'm not a huge fan of the Ensembl BioMart system because I find it difficult to use. The Synergizer has a very straightforward interface and works pretty well for most lists. Note: it hasn't been updated in a while.

# Answer
> 2 votes
By far my preferred option for doing this manually is PICR: http://www.ebi.ac.uk/Tools/picr/

BTW it is not "ridiculous" to get different numbers of genes reported for a given set of proteins. For several reasons:

1. Uniprot IDs can disappear, merge or be split
2. not all uniprot and gene IDs have a 1-to-1 relationship
3. depending on species some gene symbols can be ambiguous or synonymous.

---
Tags: database, gene, conversion
---


title: Rule Extraction from nnet results

text:
# Question
Title: Rule Extraction from nnet results
I used a script in R language that uses nnet library to predict promoter bacteria and i would like to know how to extract rules from this neural network results.

As my input of the neural network i have n positive examples of promoter and n false examples. I use as input these examples and the FASTA file with the genome of the bacteria. As a result, i have a value of 0 to 1 for each network tested corresponding to its learning. I want to discover roles that can improve my network in future experiments in the same bacteria.

Which algorithms or softwares could i use?

# Answer
Short answer: you can't.

Neural networks use positive and negative examples to add weights to the neural network architecture that is provided to it. Trying to deconvolve the meaning behind the weights is nigh-on impossible except for the simplest perceptron.

This is common to many machine learning algorithms: they work very much like black boxes. One exception is decision trees. They will report what features are used to classify the positive vs negative datasets.

However, as with all ML methods, you need to be very careful of your training datasets. This is especially true of motif searching and especially for negative datasets. It's quite easy to find a negative dataset which is completely inappropriate for learning and give you misleading accuracies.

Like @user172818, I would try traditional methods as they can work well, again, if given appropriate data. The MEME-suite would be a good start.
> 4 votes

# Answer
R's nnet package only supports fully connected neural networks with one hidden layer. This is the most primitive type of network. I doubt it will work well for promotors finding. In addition, such network won't give you useful interpretations.

If you want to explore neural networks, you should use a one-dimension convolution layer as is described in this paper. This layer effectively represents a positional weight matrix (PWM). You can know which motifs are used. With backtracking, you can also identify the coordinates of motifs. To deploy such a model, you need to be fairly familiar with deep learning and to learn one deep learning framework. For framework, you may start with keras. Once you get used to keras, you can implement your full promotor finder in less than 100 lines of python code. Alternatively, you may try dragonn. It is supposed to simplify deploying models for DNA sequences. I have no experience with it, though.

For promotor finding, it is also worth trying traditional methods. It is interesting that few/no neuralNet-based publications have evaluated traditional methods, probably because many ML people know little about classical motif finding. These methods could work well, too.
> 7 votes

---
Tags: r, fasta, machine-learning
---


# Question
Title: Reordering scaffolds according to a reference without a genetic map
I am trying to reorder scaffolds of a rice species, but no genetic map is available right now. *Oryza sativa Japonica* is a close relative of this rice species. Mummer was used to do a whole genome alignment, and I am trying to reorder scaffolds according the mummer result. But I found it's too complicated to do this by myself, many scaffolds may align to multiple Oryza sativa Japonica chromosomes or multiple loci in the same chromosome.

What software can I use to handle this task?

# Answer
> 2 votes
I've been using ATAC to successfully align different versions of the rice genome assembly. There is a chance it will work for you between two closely related species (potato and tomato work well for example).

It's not massively user friendly, but it's very fast, and it produces long, 'clean' alignments (i.e. almost all 1 to 1 alignments). The output format is a bit cryptic (zero-based start position and length instead of the more usual 1-based start and end positions), but as I said, it should basically output the order for you in about 30 minutes.

# Answer
> 5 votes
Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need.

Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a reference and/or b) rearrange contigs to improve assembly quality.

I use the tutorial as above which uses the GUI for Mauve. I

---
Tags: sequence-alignment, assembly, software-recommendation, scaffold
---


# Question
Title: How to interpret methylation calls from Bismark on opposite strands?
I'm looking at Reduced representation bisulfite sequencing (RRBS) data from ENCODE, and to align the FASTQ files I've used Bismark with Bowtie 1. When I load the resulting BAM file into `R` with `Rsamtools` and `GenomicRanges`, I get something like this:

```
> reads
GRanges object with 2 ranges and 5 metadata columns:
      seqnames         ranges strand |                     seq                                   XM          XR          XG        NM
         <Rle>      <IRanges>  <Rle> |          <DNAStringSet>                          <character> <character> <character> <integer>
  [1]     chrM [16523, 16558]      - | ATAAAACCTA...CCCTTAAATA .....h........h.......z.............          CT          GA         3
  [2]     chrM [16524, 16559]      + | TAAAGTTTAA...TTTTAAATAA .....hh.......hhh.h.z...hhhh........          CT          CT        11
  -------
  seqinfo: 25 sequences from an unspecified genome

```

For the read on the - strand (which is aligned to the "G-\>A"-converted reference), how should I compare methylation calls to those on the + strand, since they don't line up?

For example, the `z` in the methylation string of the - strand is one position ahead of the `z` in the + strand (which makes sense because of symmetric CpG methylation). But how should I determine whether these two methylation calls are "essentially the same" or not?

# Answer
The strand the bismark reports is related to the strand from which the read originated, not necessarily how it's aligned. So, alignments on the + strand shouldn't have calls overlapping those on the - strand, since you can't have a C in the same place on the same strand. One should often see Z/z next to each other on opposite strand, like in your example, since these are CpG (so it should be a Z/z on the + strand and then a Z/z on the following base on the - strand). Thus, in your example you have two reads supporting unmethylation for the CpG as a whole (1 read for each of the Cs).

The most confusing thing about BSseq is that reads (or paired-end reads) are only ever informative for a single strand. This is actually true for all sequencing with Illumina instruments (single-stranded fragments are loaded, after all), but since strands are almost always complementary we can usually ignore this fact.

As an aside, you might find MethylDackel useful (full disclosure, I wrote it). It'll be much faster at extracting methylation calls than bismark and supports nice things like excluding regions of bias methylation and likely variant positions.
> 4 votes

---
Tags: r, samtools, methylation
---