title: Jitter violin coloured by sample id

text:
# Question
Title: Jitter violin coloured by sample id
I want to have a jitter violin coloured by sample ids but I get this weird plot

```
> str(long)
Classes ‘data.table’ and 'data.frame':  156 obs. of  3 variables:
 $ gene    : chr  "gene2" "gene3" "gene4" "gene1" ...
 $ variable: chr  "1016" "1016" "1016" "1016" ...
 $ value   : num  5.82 8.41 11.22 7.21 5.61 ...
 - attr(*, ".internal.selfref")=<externalptr>
>

> p <- ggplot(long, aes(x=gene, y=value,colour=variable))
> p + geom_violin(width = 1.3) + geom_quasirandom(alpha = 0.2, width = 0.2)+ theme_classic()
There were 21 warnings (use warnings() to see them)

```

I can not figure out where I am doing something wrong that I get this instead of getting a proper violin plot by dots coloured with sample IDs

# Answer
What happens if you group by gene?

```
> long |>
  ggplot() +
  aes(x=gene, y=value, colour=variable, group=gene) +
  geom_violin(width = 1.3) +
  geom_quasirandom(alpha = 0.2, width = 0.2) +
  theme_classic()

```

@llrs raised an insightful point in their comment. There's a good chance that the "weird plot" that you see is a side effect of trying to display individual "grouped-by-variable" boxplots that have too little data. Without knowing what you would like the boxplot to look like, it's hard to create something to improve the presentation.
> 1 votes

---
Tags: r, ggplot2
---


# Question
Title: How to make a violin plot around quasirandom dots
I want to have a plot violin plot around quasirandom dots - like below

This is my data

```
> str(long1)
Classes ‘data.table’ and 'data.frame':  116 obs. of  4 variables:
 $ gene    : chr  "WT1" "KRT7" "MUC16" "EPCAM" ...
 $ variable: chr  "1016" "1016" "1016" "1016" ...
 $ value   : num  5.82 8.41 11.22 7.21 5.4 ...
 $ group   : chr  "Tumour" "Tumour" "Tumour" "Tumour" ...
 - attr(*, ".internal.selfref")=<externalptr>
>

```

By this code I get plot which does not have any violin

```
library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
p <- ggplot(long1, aes(x=group, y=value,fill=group))
p +geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,aes(colour = variable))+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right") +
  theme(legend.title=element_blank())+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 12,face="italic"))

```

I also tried adding geom\_violin like

```
p <- ggplot(long1, aes(x=group, y=value,fill=group))
p +geom_violin()+geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,aes(colour = variable))+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right") +
  theme(legend.title=element_blank())+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 12,face="italic"))

```

But gives a weird plot not similar to what I want

Any help please

***EDITED***

Thanks to @haci

By this code I get this plot but as you see the point color inside the violin is not really similar to what I want

Any suggestion to get a similar color scheme for dots inside the violins

```
p <- ggplot(long1, aes(x=group, y=value,fill=group))+guides(fill = "none")
p +geom_violin(alpha = 0.1)+scale_fill_manual(values=c("#fcf4e0", "#fcf4e0"))+geom_quasirandom(method = 'quasirandom',alpha = 0.2, width = 0.2,size=2,stroke=0.5,aes(colour = variable))+scale_color_hue(l = 5, c = 150)+ labs(title="Supplementary Figure 6B", y = "Expression(logcounts)")+theme_classic()+
  theme(legend.position="right",legend.text=element_text(size=12)) +
  theme(legend.title=element_blank(),legend.text.align = 0,text = element_text(family = "Arial"),
        axis.title.x = element_text(size = 0.00001),
        axis.title.y = element_text(size = 17))+theme(axis.text = element_text(
          color="black",
          size=12,
          face=1)
        )+facet_wrap(~gene, scales = "free")+theme(strip.background = element_blank(),strip.text.x = element_text(size = 17,face="italic"))+ theme(axis.line = element_line(colour = 'black', size = 0.8))

```

# Answer
> 3 votes
You are almost there, if you set the `alpha` parameter within your `geom_violin()` call (to sth like 0.1) you should be able to get a plot similar to the example plot.

The fact that you only have two data points for the `Stroma` group and these two points having similar values, you will end up getting a horizontal-line like violin for this group. Moreover, you do not have that many data points in general anyway so you might consider not using violins at all, especially if the outliers in `Tumour` group is coming from the same (maybe somewhat problematic?) sample, that might need to be dropped from the analysis. The reader should be able to grasp the essence of the data with our without violins (violins could help making the outliers obvious though).

---
Tags: r, ggplot2
---


# Question
Title: R ggtree apply colors to branches with multiple colors
*This question was also asked on StackOverflow*

I have a phylogenetic tree made with a Newick format

```
((a:1,b:1):2, (c:1, d:1):3):1;

```

The output will be

I have drawn this with http://etetoolkit.org/treeview/, but output with R will be same.

What I want to do is apply colors to a certain branch.

That would be

For example,

Since the Newick input doesn't include any information other than the branch (c,d) has length 3, is there any way I can include multiple color information on the branch?

I usually use ggtree for making trees, but other suggestions can be helpful.

The main script for this kind of tree building is like this, but I can't figure out how this works.

(Full script Github link: https://github.com/seongyeol-park/Human\_Lineage\_Tracing/blob/main/Rscripts/HSC\_commitment.R)

```
ggtree(tree) %<+% m_lng_dt + theme_tree2()+
  coord_cartesian(xlim = c(-1,35))+
  geom_text2(aes(subset = is.na(Source_class2) == F, x=x-n_pointmt+3, y = y+0.25, label = Source_class2), size=2.5, color="red")+
  geom_segment2(aes(subset = is.na(VAF1) == F, x=x-n_pointmt, xend = x -n_pointmt +1, y = y, yend = y, color = VAF1))+
  geom_text2(aes(subset = is.na(VAF1) == F, x=x-n_pointmt, y = y, label = VAF1), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF2) == F, x=x-n_pointmt+1, xend = x -n_pointmt +2, y = y, yend = y, color = VAF2))+
  geom_text2(aes(subset = is.na(VAF2) == F, x=x-n_pointmt+1, y = y, label = VAF2), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF3) == F, x=x-n_pointmt+2, xend = x -n_pointmt +3, y = y, yend = y, color = VAF3))+
  geom_text2(aes(subset = is.na(VAF3) == F, x=x-n_pointmt+2, y = y, label = VAF3), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF4) == F, x=x-n_pointmt+3, xend = x -n_pointmt +4, y = y, yend = y, color = VAF4))+
  geom_text2(aes(subset = is.na(VAF4) == F, x=x-n_pointmt+3, y = y, label = VAF4), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF5) == F, x=x-n_pointmt+4, xend = x -n_pointmt +5, y = y, yend = y, color = VAF5))+
  geom_text2(aes(subset = is.na(VAF5) == F, x=x-n_pointmt+4, y = y, label = VAF5), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF6) == F, x=x-n_pointmt+5, xend = x -n_pointmt +6, y = y, yend = y, color = VAF6))+
  geom_text2(aes(subset = is.na(VAF6) == F, x=x-n_pointmt+5, y = y, label = VAF6), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF7) == F, x=x-n_pointmt+6, xend = x -n_pointmt +7, y = y, yend = y, color = VAF7))+
  geom_text2(aes(subset = is.na(VAF7) == F, x=x-n_pointmt+6, y = y, label = VAF7), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF8) == F, x=x-n_pointmt+7, xend = x -n_pointmt +8, y = y, yend = y, color = VAF8))+
  geom_text2(aes(subset = is.na(VAF8) == F, x=x-n_pointmt+7, y = y, label = VAF8), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF9) == F, x=x-n_pointmt+8, xend = x -n_pointmt +9, y = y, yend = y, color = VAF9))+
  geom_text2(aes(subset = is.na(VAF9) == F, x=x-n_pointmt+8, y = y, label = VAF9), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF10) == F, x=x-n_pointmt+9, xend = x -n_pointmt +10, y = y, yend = y, color = VAF10))+
  geom_text2(aes(subset = is.na(VAF10) == F, x=x-n_pointmt+9, y = y, label = VAF10), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF11) == F, x=x-n_pointmt+10, xend = x -n_pointmt +11, y = y, yend = y, color = VAF11))+
  geom_text2(aes(subset = is.na(VAF11) == F, x=x-n_pointmt+10, y = y, label = VAF11), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF12) == F, x=x-n_pointmt+11, xend = x -n_pointmt +12, y = y, yend = y, color = VAF12))+
  geom_text2(aes(subset = is.na(VAF12) == F, x=x-n_pointmt+11, y = y, label = VAF12), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF13) == F, x=x-n_pointmt+12, xend = x -n_pointmt +13, y = y, yend = y, color = VAF13))+
  geom_text2(aes(subset = is.na(VAF13) == F, x=x-n_pointmt+12, y = y, label = VAF13), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF14) == F, x=x-n_pointmt+13, xend = x -n_pointmt +14, y = y, yend = y, color = VAF14))+
  geom_text2(aes(subset = is.na(VAF14) == F, x=x-n_pointmt+13, y = y, label = VAF14), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF15) == F, x=x-n_pointmt+14, xend = x -n_pointmt +15, y = y, yend = y, color = VAF15))+
  geom_text2(aes(subset = is.na(VAF15) == F, x=x-n_pointmt+14, y = y, label = VAF15), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF16) == F, x=x-n_pointmt+15, xend = x -n_pointmt +16, y = y, yend = y, color = VAF16))+
  geom_text2(aes(subset = is.na(VAF16) == F, x=x-n_pointmt+15, y = y, label = VAF16), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF17) == F, x=x-n_pointmt+16, xend = x -n_pointmt +17, y = y, yend = y, color = VAF17))+
  geom_text2(aes(subset = is.na(VAF17) == F, x=x-n_pointmt+16, y = y, label = VAF17), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF18) == F, x=x-n_pointmt+17, xend = x -n_pointmt +18, y = y, yend = y, color = VAF18))+
  geom_text2(aes(subset = is.na(VAF18) == F, x=x-n_pointmt+17, y = y, label = VAF18), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF19) == F, x=x-n_pointmt+18, xend = x -n_pointmt +19, y = y, yend = y, color = VAF19))+
  geom_text2(aes(subset = is.na(VAF19) == F, x=x-n_pointmt+18, y = y, label = VAF19), size=2.5)+
  geom_segment2(aes(subset = is.na(VAF20) == F, x=x-n_pointmt+19, xend = x -n_pointmt +20, y = y, yend = y, color = VAF20))+
  geom_text2(aes(subset = is.na(VAF20) == F, x=x-n_pointmt+19, y = y, label = VAF20), size=2.5)+
  scale_color_gradientn(colors = my_palette, breaks = col_breaks)+
  theme(axis.text = element_text(size=15, angle=45, hjust=1), plot.title = element_text(size=20), plot.margin= unit(c(0.2,0,0.2,0.2),"cm"))

```

# Answer
> 1 votes
The simple solution is Figtree a GUI that easily colours trees on a 'point and click' basis. In linux this can be driven via CLI, which is useful.

The easy coding solution is to use extended Newick format (NHX) to colour the tree.

```
(((a:1,b:1[&NHX:C=2.2.2]):2[&NHX:C=2.2.2], (c:1[&NHX:C=100.100.100], d:1[&NHX:C=100.100.100]):3[&NHX:C=10.10.10]):1);

```

The next thing is a viewer that accepts NHX format and here I used 'Icytree' here

The alternative format a tree viewer might take is:

```
(((a:1,b:1[&&NHX:C=2.2.2]):2[&&NHX:C=2.2.2], (c:1[&&NHX:C=100.100.100], d:1[&&NHX:C=100.100.100]):3[&&NHX:C=10.10.10]):1);

```

I noted you'd used ETE3 toolkit website for your tree. Just to mention ETE3 has two levels,

* the underlying raw Python library
* the website GUI

The underlying raw Python is cool, except it was written to v3.6 because you combine it with your own code, the website is alright, but I think there's better.

Using NHX allows the branch width to be set as well:

```
(((a:1,b:1[&NHX:C=2.2.2:W=3]):2[&NHX:C=2.2.2:W=3], (c:1[&NHX:C=100.100.100:W=3], d:1[&NHX:C=100.100.100:W=3]):3[&NHX:C=10.10.10:W=3]):1);

```

Sadly you can't set the font (which you obviously can in `ggtree`

NHX normally used for auxiliary information, such as species name, or a particular phenotype.

In summary, I do prefer NHX if you are coding over `R`.

**Note** For tree colouring you have to set the use NHX format BTW in the options, a tree browser will not do it automatically.

---

The basic `ggtree` first step is to make the tree ...

```
library("ggtree")
tree <- read.tree(text = "(((a:1,b:1):2, (c:1, d:1):3):1);") # my bracketing might be wrong
p<-ggtree(tree)
p+geom_tiplab()

```

Then use `ggplot2` style functions to overlay the presentation, I'm assuming you do understand `ggplot2`. What the code you showed was doing was individually changing every font and the position of every branch in the tree. This is definitely not cool code and I wouldn't consider it further, i.e. its not even worth bothering with, that ain't how graphs nor trees are supposed to be represented.

In any case, overall `ggplot2` style code is always fiddly and to be honest I gave up using it, but you increment/or layer the changes you want as the code you supplied describes. Honestly I don't want to get into the fiddly incremental layered changes that `ggplot2` graphs require.

`ETE3` is good and does do coloured trees, which are described here, http://etetoolkit.org/docs/latest/tutorial/tutorial\_drawing.html. The only hassle it was written for Python 3.6 - that long ago, but you can access it via `conda` and set the Python version.

`ETE4` will be very cool but its badly delayed in its release date, they did get funding for this so the delay is not cool, maybe never who knows. Getting it up to Python 3.10 would be useful at the very least.

My personal take is ETE3 is useful for automating the tree output so you can manually assess the iterations, its quite critical for e.g. testing and debugging. Whether you are turning that into publication quality - which is universally how graphs are handled - I personally find excessive. The difference is that complex graphs represent loads of information which is not always immediately obvious to a casual observer, therefore gamma (transparency), plot type loads of stuff are relevant. A tree is a tree and 'dressing up a tree' makes not a jot of difference for an experienced observer. If someone wasn't quite clear about a tree, they'd just ask you for the tree file and file it through 'Figtree'.

---

# Answer
> 1 votes
See change colors or line types of arbitrarily selected branches in `ggtree` documentation for details. Here's my R code for applying a color to a specific branch in `ggtree`:

```
suppressPackageStartupMessages(library(ggtree))
# read tree data
tree <- read.newick('tree.nwk')
# create a dataframe of the default color (black) for branches
branchColor <- data.frame(node = 1:Nnode2(tree), colour = 'black')
# set a color for the specifc branch
branchColor[7, 2] <- 'blue'
branchColor[3, 2] <- 'orange'
branchColor[6, 2] <- 'purple'
branchColor[c(1, 2), 2] <- 'red'
# plot tree
ggtree(tree, size = 1) %<+%
  branchColor +
  aes(colour = I(colour)) +
  geom_tippoint(size = 2, color = 'black') +
  geom_tiplab(size = 5, color = 'black') +
  geom_treescale(x = 0, y = 0, width = 0.5, fontsize = 5, linesize = 1)

```

The tree looks like this:

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


# Question
Title: Is there a command line tool to split a SAM/BAM file by CB (cell barcode) tag?
I have a BAM file from a single cell sequencing experiment. Each read has had the cell barcode annotated in the `CB` tag. Some reads do not have a `CB` tag. For each of an arbitrary set of barcodes, I would like to split out the reads for each of those barcodes into a separate file. Is there a tool that does this?

# Answer
I wrote a package called sinto (https://github.com/timoast/sinto) that contains a function to split a BAM file into different files based on the cell barcodes.

The installation can be done via anaconda:

```
conda create -p conda/envs/sinto python=3.8 sinto

```

And a usage example from the documentation is:

```
sinto filterbarcodes -b my_bam.bam -c my_barcodes.txt --barcodetag "CB"

```

Where the `my_barcodes.txt` file contains cluster names in a tab separated list format:

```
TGGCAATGTTGAAGCG-1  A
GACCAATCACCATTCC-1  A
CAGGATTCAGAACTTC-1  B
GAACCTAAGAGAGGTA-1  B,A
..

```
> 1 votes

# Answer
I wrote a python script based on pysam that is free for anyone to use:

<pre><code>"""
MIT License
Copyright (c) 2020 Warren W. Kretzschmar
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

from pathlib import Path

import pysam

class BamWriter:
    def __init__(self, alignment, barcodes, prefix):
        self.alignment = alignment
        self.prefix = prefix
        self.barcodes = set(barcodes)
        self._out_files = {}

    def write_record_to_barcode(self, rec, barcode):
        if barcode not in self.barcodes:
            return
        if barcode not in self._out_files:
            self._open_file_for_barcode(barcode)
        self._out_files[barcode].write(rec)

    def _open_file_for_barcode(self, barcode):
        self._out_files[barcode] = pysam.AlignmentFile(
            f"{self.prefix}_{barcode}.bam", "wb", template=self.alignment
        )

def main(input_bam, barcodes_file, output_prefix, contigs):
    """Split a 10x barcoded sequencing file into barcode-specific BAMs

    input:
    barcodes_file: a file containing barcodes or a single barcode
    contigs: '.' for all contigs, 'chr1' for the contig 'chr1',
    or '1-5' for chromosomes 1, 2, 3, 4, and 5
    """
    alignment = pysam.AlignmentFile(input_bam)
    if Path(barcodes_file).is_file():
        with open(barcodes_file, "r") as fh:
            barcodes = [l.rstrip() for l in fh.readlines()]
    else:
        barcodes = [barcodes_file]
        print(f"Extracting single barcode: {barcodes}")
    writer = BamWriter(alignment=alignment, barcodes=barcodes, prefix=output_prefix)
    if contigs == ".":
        print("Extracting reads from all contigs")
        recs = alignment.fetch()
    else:
        if "-" in contigs:
            start, end = contigs.split("-")
            print(f"Extracting reads from contigs {start} to {end}")
            recs = (alignment.fetch(str(contig)) for contig in range(start, end + 1))
        elif "," in contigs:
            contigs = contigs.split(",")
            print(f"Extracting reads from contigs {contigs}")
            recs = (alignment.fetch(str(contig)) for contig in contigs)
        else:
            print("Extracting reads for one contig: {contigs}")
            recs = (alignment.fetch(c) for c in [contigs])

    for rec in recs:
        try:
            barcode = rec.get_tag("CB")
            writer.write_record_to_barcode(rec=rec, barcode=barcode)
        except KeyError:
            pass

if __name__ == "__main__":
    import sys

    main(
        input_bam=sys.argv[1],
        barcodes_file=sys.argv[2],
        output_prefix=sys.argv[3],
        contigs=sys.argv[4],
    )
```
</code></pre>
> 4 votes

# Answer
You can do it with samtools:

```
# Save the header lines
samtools view -H $BAM_FILE > SAM_header

# Filter alignments using filter.txt. Use LC_ALL=C to set C locale instead of UTF-8
samtools view $BAM_FILE | LC_ALL=C grep -F -f filter.txt > filtered_SAM_body

# Combine header and body
cat SAM_header filtered_SAM_body > filtered.sam

# Convert filtered.sam to BAM format
samtools view -b filtered.sam > filtered.bam

```

And 10x Genomics has developed subset-bam a tool that can split a BAM file by barcode.
> 0 votes

# Answer
If you have access to a lot of cores:

10x has subset-bam a pretty fast way to pull out all the cells with a specific CB/BX if you only need a small number.
> 0 votes

---
Tags: sam, single-cell, data-preprocessing
---


title: Problem with sequence alignment

text:
# Question
Title: Problem with sequence alignment
I'm trying to find if the peptides in column A are present within the longer peptides of column B.

Pairwise alignment tools take the data from both columns as one string (i.e., TIIDYTNNHYTNNHLEII... and PFYMNITGYTDAYNTDIIDYIIDYTNNHLEIITLFN...) which leads to a number of peptides from column A left unaligned even though the peptide is present within one of the longer peptides in column B.

I would like to find a method that works for detecting such peptides.

# Answer
> 3 votes
I don't think that you want an *alignment* so much as you want to *test for the presence of an exact substring* (column A) *within larger strings* (column B).

There are simple excel utilities that do this, for example the `SEARCH()` function. There is no need to use alignments.

---
Tags: proteins, sequence-alignment, sequence
---


# Question
Title: Any way to filter out highly correlated genes with limma linear model?
I have Affymetrix gene level expression matrix (genes in the rows and sample ID on the columns), and I have annotation data of microarray experiment observation where sample ID in the rows and description identification on the columns.

**example data**

here is Affymetrix gene level expression matix:

```
> dim(exprs_mat)
[1] 112830   735
> exprs_mat[1:4, 1:3]
        Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at           6.062215        6.125023        5.875502
10_at          3.796484        3.805305        3.450245
100_at         5.849338        6.191562        6.550525
1000_at        3.567779        3.452524        3.316134

```

here is annotation data which contain experiment observation:

```
> dim(ano)
[1] 735   6
> head(ano)
                       SampleID   GA Batch     Set Train Platform
Tarca_001_P1A01 Tarca_001_P1A01 11.0     1 PRB_HTA     1    HTA20
Tarca_013_P1B01 Tarca_013_P1B01 15.3     1 PRB_HTA     1    HTA20
Tarca_025_P1C01 Tarca_025_P1C01 21.7     1 PRB_HTA     1    HTA20
Tarca_037_P1D01 Tarca_037_P1D01 26.7     1 PRB_HTA     1    HTA20
Tarca_049_P1E01 Tarca_049_P1E01 31.3     1 PRB_HTA     1    HTA20
Tarca_061_P1F01 Tarca_061_P1F01 32.1     1 PRB_HTA     1    HTA20

```

I intend to see how the genes in each sample are correlated with GA value of corresponding samples in the annotation data. How can I get my expected correlation matrix and filter out the genes by correlation value? any idea to make this happen correctly?

**my attempt**

```
library(limma)
fit <- limma::lmFit(exprs_mat, design = model.matrix( ~ 0 + t(ano$GA))
fit <- eBayes(fit)
topTable(fit, coef=2)

```

but this attempt gave me dimension error and I still didn't get highly correlated genes from `exprs_mat` with each `ano$GA`in `ano`data. I just removed `0+` and ran the code, then I got this error:

```
row dimension of design doesn't match column dimension of data object

```

I am pretty sure I did `df[!is.na(df$GA), ]`, removed `NA` in both `exprs_mat` and `ano` metadata already, but the same error still persists. I also tried `df[complete.cases(df),]` to make sure I removed all `NA` values in my data.

How can I get those highly correlated genes correctly? any possible way to make this happen in R? any thought?

How can I cut down `exprs_mat` with the shortlist of genes? any idea?

**desired output**

I intend to filter out genes in the gene expression matrix `expr_mat` by using correlation value (only keep highly correlated genes), expected sub expression matrix should have same data structure as `expr_mat`. any way to get this done?

# Answer
The problem is that you are transposing the vector of GA values with `t(ano$GA)`. Why would you do that? It produces a row matrix that is inappropriate for input to `model.matrix`.

All you need is

```
design <- model.matrix(~ano$GA)

```

and everything will work fine.
> 2 votes

---
Tags: r, microarray, gene-expression, limma
---


# Question
Title: Comparing differential expression across samples - is batch effect correction needed?
I have a bulk RNA-seq dataset made up of control and treatment conditions for a range of cell lines. This dataset was generated in two batches, such that the cell lines are split between batches but both the treatment and control for each cell line are within the same batch. As the cell lines are very different, I'm looking to do DE analysis for each cell line and then compare these DE responses between cell lines.

I'm wondering whether batch effect correction is useful here since I'm only ever comparing DE values and not the raw counts themselves. Also, the baseline expression values for the cells lines are so different that they already constitute a strong biological batch effect that is perfectly correlated with the technical batch effect, which batch correction wouldn't be able to account for.

# Answer
> 2 votes
There is no need for batch correction because all your comparisons are within batch. As long as you are making comparisons for each line, and treatment and control for the same line are in the same batch, then any batch effects are already accounted. for.

# Answer
> 0 votes
If your treatment and control are prepared in different batches, it's impossible to distinguish between treatment-specific effects, and batch-specific effects. DESeq2 will complain when you attempt to carry out such correction in the statistical model for this reason.

Baseline differences in expression are a large reason for carrying out batch correction in the first place; it makes differential expression analyses *more* powerful, not less. I'd recommend consulting with a statistician (prior to carrying out an experiment / sequencing) regarding the best way to split samples into different batches so that batch effects can be discovered and accounted for. As you mention, there are multiple effects that should be considered (e.g. cell lines, sequencing run).

---
Tags: rna-seq, statistics, differential-expression, batch-effects
---


title: Basic RNA Differential Expression in R

text:
# Question
Title: Basic RNA Differential Expression in R
I have two matrices, one for individuals before treatment and one for the same individuals after treatment. Both matrices are raw read counts of RNA expression.

```
        Treated_1  Treated_2  Treated_3  Treated_4  Treated_5
RNA_1   105        283        64         155        51
RNA_2   359        27         47         348        84
RNA_3   99         10         89         345        77
RNA_4   48         100        77         74         83

        Untreated_1  Untreated_2  Untreated_3  Untreated_4  Untreated_5
RNA_1   100          200          50           130          33
RNA_2   200          10           30           300          75
RNA_3   90           10           60           320          60
RNA_4   66           50           33           60           55

```

How can I tell which RNA's are significantly differentially expressed in R using `edgeR` or `DESeq2`? I'm thinking first normalize using TMM then compare the mean of the normalized values for treated vs. untreated?

# Answer
See here for a Differential Expression guide which discusses carrying out differential expression using DESeq2:

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

It starts with a standard "quick start" workflow, and goes into quite a lot of detail explaining all the different things that can be done:

```
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_trt_vs_untrt")

```
> 3 votes

# Answer
> I'm thinking first normalize using TMM then compare the mean of the normalized values for treated vs. untreated?

That is a bit too simple. You need to analyse your experiment as a paired comparison. You need to combine the treated and untreated count matrices into one and form factors for individual and treatment:

```
Counts <- cbind(Counts.treated, Counts.untreated)
Individual <- gl(5,1,length=10)
Treatment <- gl(2,5,length=10,labels=c("Treated","Untreated"))
Treatment <- relevel(Treatment, ref="Untreated")
design <- model.matrix(~Individual+Treatment)

```

Then follow a standard workflow: https://bioconductor.org/packages/release/workflows/vignettes/RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html

```
library(edgeR)
y <- DGEList(counts=Counts)
y <- calcNormFactors(y)
y <- estimateDisp(y,design)
fit <- glmQLFit(y,design)
test <- glmQLFTest(fit)
topTags(test)

```
> 1 votes

---
Tags: r, deseq2, differential-expression, rna, edger
---


# Question
Title: Aligning scRNA-seq fastq to .bam without cell barcodes
I would like to convert the fastq files from this dataset: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE98679 to .bam files. The group did not provide the scripts they used to to their alignment, I know they already posted their reads but I need the .bam files for my pipeline. Is it possible to do this without knowing the cell barcodes? If so with what tools?

# Answer
The linked dataset is based on the Smart-seq2 technology. This is a plate-based assay in which individual cells are first sorted into individual wells of a microwell plate via FACS. The cells are then individually lysed and libraries are prepared. Since this happens per well one assigns a single sequencing barcode to each well/cell, hence no further cellular barcoding as in droplet-based methods such as 10x Chromatium is required. Long story short, each fastq file that one obtains from the sequencing represents one individual cell. Hence, one can use any standard RNA-seq alignment or quantification technique such as `STAR` or `salmon` to obtain alignments/quantifications. One would then need to parse the individual quantifications into a single count matrix. If one used aligners such as `STAR` this could be done with high-performance quantifiers such as `featureCounts`. If using pseudo- or selective aligners such as `salmon` or `kallisto` I would use `tximport` to aggregate the per-transcript quantifications to the gene level while importing into R to get a gene-level count matrix. As each file is a single cell there is no need for more complicated solutions such as `CellRanger` which do not apply here.
> 1 votes

---
Tags: sequence-alignment, scrnaseq, sam, fastq
---


# Question
Title: Existing tool for converting gff3 to genbank (gbk)
I want to convert my gff3 annotation files to genbank format for use in Mauve. I found the seqret tool here https://www.ebi.ac.uk/Tools/sfc/emboss\_seqret/ which can perform this task, but my files (bacterial genomes) are too big.

Can anyone suggest a simple method to make this conversion?

# Answer
I believe `seqret` is already the simplest approach. If you install the emboss command line utility, genome (file) sizes will not be a problem. You can even install it with conda, see https://anaconda.org/bioconda/emboss.

After that, you can use it as the example:

```
seqret -sequence {genome file} -feature -fformat gff -fopenfile {gff file} -osformat genbank -osname_outseq {output prefix} -ofdirectory_outseq gbk_file -auto

```

Hope it helps
> 3 votes

# Answer
You can use BioPerl to do this. Read in the GFF and make gene objects and write out those gene objects in Genbank format.

Perhaps this tool (doing the opposite) will get you started: https://metacpan.org/pod/bp\_genbank2gff3.pl

Here is a script doing what you ask but it doesn't use any libraries (so is probably buggy): https://gist.github.com/avrilcoghlan/5386810

Hope that helps, Dan.
> 1 votes

# Answer
Here's a python script based on the bcbio-gff (reading gff3) and biopython (writing genbank): https://github.com/chapmanb/bcbb/blob/master/gff/Scripts/gff/gff\_to\_genbank.py
> 1 votes

---
Tags: format-conversion, gff3, genbank
---


# Question
Title: How can I remove FASTA sequences with an asterisk inside the sequence?
I have several protein sequences with an asterisk in the middle of the sequence, and I'd like to remove those sequences.

In the example below, I have three sequences:

From this example, I would like to obtain only the following output:

# Answer
Try `seqkit grep`

```
seqkit grep --by-seq --invert-match --pattern '*' seq.fa

```

`--by-seq` search by seq

It will output sequences without an asterisk.
> 4 votes

# Answer
It’s good you brought up the question.

### `*` are stop codons

Don't simply parse out the `*` symbol. If done incorrectly, it could constitute **'tampering'** and that’s a *red line* in biology.

The `*` means a stop codon in the middle of a protein sequence. There are quite a number of reasons why there is an `*`, i.e., a stop codon in the sequence, for example, in no particular order:

* the open reading frame is out of frame which is part of the natural molecular biology, i.e., there can be multiple reading frames in a gene
* the open reading frame is out of frame due to a sequencing error and poor annotation (which is common);
* the wrong reading frame was used (very, very common);
* the locus isn't translated, e.g., it’s an intron or UTR (untranslated region)

The list goes on, e.g., it could be expression regulation ...

---

1. If it is not an intron/UTR the way you might 'remove' `*` is by changing the reading frame, e.g., try using the second codon position (then the third) as the opening reading frame: if the `*` 'disappear' it was very likely just a reading frame issue. How to do that is a separate question. The resulting amino acid sequence will of course be very different.
2. There's a second possibility if the reading frame is correct, which is to treat each locus between the `*` as a separate protein. That way the `*` is 'removed' because their a multiple proteins within the locus. In that case, a `\n` would replace the `*` (Mac and Linux), i.e., carriage return.
> 4 votes

# Answer
You could use Biopython for this. It makes FASTA parsing really easy.

```
'''
Python version: 3.10.0
Biopython version: 1.79
'''
from Bio import SeqIO
for seq_record in SeqIO.parse("input.fasta", "fasta"):
    if "*" not in (seq_record.seq):
        print(">",seq_record.id)
        print(seq_record.seq)

```
> 1 votes

---
Tags: fasta, phylogenetics, proteins, phylogeny, sequence
---


# Question
Title: Example of a small but moderately complex PDB file?
I am developing a bioinformatics framework only for learning purposes. I want to calculate --

1. dihedral angles
2. detect H, E, and C components
3. detect Hydrogen bonds
4. etc.

Therefore, I want a PDB file that represents a protein with a small chain but is complex enough for me to experiment with.

Can you suggest such a PDB file?

# Answer
Although Insulin is a great, well studied example of a small structure it does have the slight downside of not having any sheet structure in the monomer. There is a small cross-chain sheet between residues 24-26 in both the B and D chains (from https://www.ebi.ac.uk/pdbe/entry/pdb/1zni/ ).

Another possibility might be crambin (http://www.cathdb.info/version/latest/domain/1ejgA00) which does have a small antiparallel sheet, although it is tiny at only 2 residues per strand.

A more general approach might be to use something like the top2018 structure set from the Richardson lab (https://zenodo.org/record/5777651#.YtQB59LMJH4) . I used the top500 back in the day, and they are a high-quality set of structures.
> 3 votes

# Answer
Insulin, perhaps?

https://www.rcsb.org/structure/1trz

From here:

> Insulin is a perfect molecule for exploring protein structure. It is small enough that you can display all of the atoms and still have a picture that is not too confusing. Human insulin is pictured here, using entry 1trz. This structure, like many other insulin structures, includes a hexameric complex of the hormone, and only one monomer (composed of one A-chain and one B-chain) is shown here. In the structure, you can see many of the key features that stabilize protein structure. Carbon-rich amino acids, like leucine and isoleucine, cluster in the middle of the molecule, forming a hydrophobic core, and the surface is covered with charged amino acids like arginine and glutamate that interact favorably with the surrounding water. Also notice the three disulfide bridges between cysteine amino acids, which stabilize this tiny protein.
> 2 votes

---
Tags: pdb
---


# Question
Title: Statistical significance of pairwise protein sequence alignment
I have identified a conserved region in a protein I'm studying, by doing multiple sequence alignments against a dozen of its orthologs. I've then found that this region also aligns relatively well with a functional region (binding motif) in a completely unrelated protein, and I'm trying to estimate whether it's reasonable to believe the two regions might be bound by the same protein. blastp doesn't identify this match, would I then assume that this is insignificant?

Let's say I have a 15 aa long peptide sequence **X**, which I align against a longer 200 aa protein sequence **Y**. I find some best match (gaps allowed).

What statistical test could I now perform to answer the question:

*Is the match I found for sequence **X** in sequence **Y** significantly better than what I would expect to find between **X** and a random 200 aa sequence **Z**?*

# Answer
> 1 votes
Just formalise the answer, its a straight use of the Blast E-value and I wouldn't personally advise against anything beyond that (explained below) within molecular evolution.

So the Google gives the definition

> The BLAST E-value is the number of expected hits of similar quality (score) that could be found just by chance. E-value of 10 means that up to 10 hits can be expected to be found just by chance, given the same size of a random database.

Thus the smaller the E-value the better the hit.

Using a local blast server will provide a specific query verse reference E-value.

```
blastp -query gene.fa -subject genome.fa

```

The broader question concerns the whether the Blast hit is a genuine homology and that is very complicated. This is simply because amino acid loci could have arisen by convergent evolution rather than direct descendancy which are homologues and this can become controversial.

When trying to assign a function to a amino acid locus this starts becoming relevant because e.g. the selection pressures of convergence could very difference.

Identifying homologues would require a phylogenetic tree which is compatible with the established phylogenetic tree, e.g. for house-keeping genes as a starting point. It gets complex quickly beyond that. Functional domains within the Blast variant (Blast-Psi) could be of use, its field specific however.

Thus the simple answer is to use Blast E-value (bit-score) and it is a universally recognised metric, the reader can then make up their own mind what the value represents. You could simply simple state the definition of the E-value in the report.

---
Tags: sequence-alignment, statistics
---


# Question
Title: Common tools for generating results and visualizations for bioinformatics and genomics publications
I am a computational scientist with background in computer science and engineering however very new to the field of bioinformatics and computational genomics. Looking at papers in nature, PLOS, Cell, PNAS etc. within the topic of computational genomics, I see a very similar set of very sophisticated graphical results, ranging from clustered/sorted heat map representation of gene expressions to log2 fold changes, volcano plots, circular graph network visualizations etc. Coming from an engineering background I would often generate customized tools for visualizing results myself as part of each paper, albeit not at the sophistication level of results I see in bioinformatics. But here, it seems there is a set of tools standardized where many are using them to display their results. Is there a summary of these result generation/visualization tools along with the references to their implementation/usage (opensource I assume)? Something like that could be very helpful to newcomers like myself. Thanks in advance.

Below is an example /elaboration based to guide the question further into specifics:

Consider the example paper:

https://doi.org/10.1016/j.cell.2018.04.018

This paper contains a principle set of results/visualizations I see almost universally in many nature/cell computational genomics papers. I am wondering if there exists an **integrated library** either in python, R (or Matlab for example when they create toolboxes and integrate a commonly used subset of tools in an area) that everyone taps into for generating these type of results. The "principal set" that I am referring to include:

1. Genome browser track snippets: Figure 1A
2. log2Fold x-y scatter plots : Figure 1B
3. Hierarchically clustered/sorted gene exp heatmaps: Figure 1C,1F
4. Venn diagrams : Figure 4F
5. Hi-C type square ligation heat-maps: Figure 6A

# Answer
> 2 votes
Figure 1A: These genome browser snippets are either plain screenshots from applications like IGV or created from within R with software like gviz.

Figure 1B: This looks like base R plotting routine. Like `plot(x,y)` and then some coloring and legending. No particular "tool" here other than really just base R or ggplot2.

Figure 1C/F: Heatmap packages include ComplexHeatmap and pheatmap. I personally use the former as it is versatile with extensive documentation.

Figure 4F: Venn diagrams, just google "Venn" and "R" for several packages that do it in R. Same for other languages. There are plenty of packages. No preference.

Figure 6A: These plots are often returned by specialized HiC software such as HiCexplorer which is powerful with lots of documentation.

---
Tags: bioconductor, visualization
---


# Question
Title: plotting gene expression after EdgeR DE analysis using RUVg (RUVseq) covariates
I have used the empirical RUVg method (from RUVseq) to estimate the unwanted variation of my dataset (consisting of several public datasets analysed together, with controls and case samples in different datasets.

Afterwards, I used EdgeR to perform my DE analysis with 2 RUV covariates:

```
designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))

y <- DGEList(counts=counts(set2), group=pheno$condition4)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, designRUV)
y <- estimateGLMTagwiseDisp(y, designRUV)
fit <- glmFit(y, designRUV)

```

I got my results after making all the contrast that I needed but I would like to plot the expression of some genes of interest. One plot per gene with all the conditions in `condition4` has 9 categories (expression level or norm counts on the Y axis).

My doubt is which counts to use. I have the normalised ones from RUVg, but it states they should only use them for exploratory analysis (maybe this is the case??). These counts as I understand are the residuals of the model, so not really sure what to do.

Another option is to use the normalised counts produced by EdgeR, but these ones doesn't account for the unwanted variation (W\_1, W\_2 covariates) (right???)

Any advice or nice solution?? :=) Thanks in advance!!

# Answer
> 2 votes
You could correct your logCPMs (`edgeR::cpm(y, log=TRUE)`) with the covariates from RUVseq using limma and then use these for plotting:

```
set2 <- RUVg(...)

designRUV <- model.matrix(~0 + condition4 + W_1 + W_2, data = pData(set2))

y <- DGEList(counts=counts(set2), group=pheno$condition4)
y <- calcNormFactors(y, method="upperquartile")
y <- estimateGLMCommonDisp(y, designRUV)
y <- estimateGLMTagwiseDisp(y, designRUV)
fit <- glmFit(y, designRUV)

# correction with limma
logcpm <- cpm(y, log=TRUE)
design <- model.matrix(~0+condition4, pData(set2))
corrected <- removeBatchEffect(logcpm, design=design, covariates=as.matrix(pData(set2)[,c(W_1,"W_2"])

```

---
Tags: rna-seq, normalization, edger
---


# Question
Title: Converting Chip-seq fastq format ({$x_k$} list of seq segments) to genome browser {sequence $x_i$, count-number $y_i$} format
Very new to bioinformatics. I am trying to figure out the steps required to convert a chip-seq fastq file into a genome browser track, with the end goal of being able to access the sequence (x) vs. activity level/count (y) data associated.

Any genome browser activity/count track data has the following functional format  $y=F(x)$

Here x is the sequence axis, and y is the activity/count axis. Typically the x(i) is the nucleotide associated with the precise genomic coordinate i (an integer number) depending on the genome used to align (human, mouse, etc.). What I need is a format allowing direct access to the nucleotide representation x(i), associated with the index i, as well as the count levels at y(x(i)). What is the closest data file type that I need? And how do I get to it?

In other words, what are the steps for a fastq file associated with a chip-seq data to be transformed to this x(i),y(i) format, or best yet if there is an example Python or R code out there that could illustrate the steps of going from file1.fastq to something like to two files **file1\_x\_sequences** and **file2\_y\_counts**? I only know the first step is an alignment to a genome, but stuck on the next steps to complete the picture. Thank you.

Ultimately I need a list of sequences (tags) lets say 250bp long, and a count number that represents the chip-seq enrichment for the part of the genome (the 250bp region). Further reading, indicates I need to detect/call "peaks" with MACS2/homer to generate a BAM file and ultimately a Fasta file for the sequence lists, but I don't know yet how the 'counts' come in (in what file format and where from). Any hints there would complete the picture I think. Ultimately this is training data for x=250bp sequence to y=count/chip-seq enrichment activity model.

# Answer
I have no idea about what you are actually trying to achieve, but getting these browser tracks is simple. First you align your fastq files to a reference genome, e.g. with `bowtie2`, this returns a SAM file. Sort the SAM file, convert to BAM and sort + index it with `samtools`. Then use something like `bamCoverage` from the deeptools suite to get your browser tracks.

A word of caution. Since you seem to be new to OMICS and bioinformatics in general: NGS data, especially applications that are as noisy as ChIP-seq have lots of pitfalls. Be sure to go through the literature to find references for what you are trying to develop, at least the basic idea and discuss with experienced folks whether it makes sense. It sounds like you want counts of the peak summits, is that correct? In any case, be sure not to reinvent the wheel and to do a meaningful analysis. Discussion with experienced peers is key, that may save you from pitfalls and rookie mistakes. This entire choice of words "sequence axis", "activity level" is very uncommon, I have never heard of this before despite working on NGS data for years, for what it's worth.
> 5 votes

---
Tags: r, python, bioconductor, chip-seq
---


# Question
Title: How to extract DNA sequence from browser track input files (BigWig, bed etc.) files
A BigWig file for something like Chip-seq, or ATAC data has the following format for each line i:

chrom\_i - start\_i - end\_i - value\_i

I need to convert for each line the (start\_i -\> end\_i) reference to the actual DNA sequence 'ATATCCTAG...' so that at the end I have a list of DNA sequences like 'ATATCCTAG...' and the Chip-seq activity value associated with that sequence segment. What are the ways to do this? For example, what tool(s) can I use to generate the (sequence, value) pairs from browser tracks for the Chip-seq or ATAC etc. And in fact, going from BigWig browser track files for Chip-seq the most direct way to produce this kind of (seq,value) data? Thank you.

# Answer
\[Here is the link.\]\[1\]. Instructions are very clear at the website
\[1\]: https://bedtools.readthedocs.io/en/latest/content/tools/getfasta.html
> 1 votes

# Answer
Bedtools can take a list of coordinates and a fasta and return a multifasta of the specified regions.
> 1 votes

---
Tags: r, python, chip-seq, atac-seq, macs2
---


# Question
Title: Determination of transmembrane domain limits using multiple prediction tools
Integral membrane proteins are proteins that have at least one transmembrane (TM) segment. I performed predictions of TM domains using three tools: ProtScale, ΔG predictor 1.0, TOPCONS.

For example, the results in terms of aminoacid positions for the TM1 domain were as follows:

* **ProtScale**: 37-55
* **ΔG predictor**: 37-55
* **TOPCONS**: 36-56

The limits of the predicted TM domains can be used for drawing the topology of this protein.

For this purpose, we need to have representative limits from these results.

**Statistically speaking:**

* which values should be chosen?
* should we calculate a mean : (37+37+36)/3=36.67 knowing that these values correspond to aminoacid positions and 36.67 doesn't mean anything?

How can I use these results (37, 37, 36) to determine the limit corresponding to the 1st aminoacid of TM1 segment?

# Answer
The variation you are seeing are minor differences in estimates in the transmembrane domain, these are usually modelled by hydrophobicity and are accurate predictions. If you look at the secondary/tertiary structure you might see alpha helices which will 'burrow' into the membrane. The region in this particular question is a small domain however, but the position is 'about right'.

Thus question is whether 1 amino acid at 5' and 3' site is inside or outside the membrane - I think this is overthinking. Edge-effects are common in any estimate and 1 amino isn't bad variation at all. To double check have a look at the Kyte-Doolittle index on the Expasy website (which is one of the first), and TMHMM. Finally, you need to check the crystal structure if available.

What does puzzle me is the apparent absence of a 3' GPI anchor, but I dunno the biological system so ain't for me to comment.
> 0 votes

---
Tags: membrane-protein, tools
---


# Question
Title: How can I generate ASV file from nanopore sequencing data?
I am converting the fast5 file to fastq by using guppy basecaller after that by using kraken2 classified the sequence.

Now I am trying to generate ASV file.

Is this possible to generate ASV file from kraken output?

Or is this any other way to generate ASV file from nanopore sequencing?

I am new to nanopore sequencing data. I need help or suggestions to generate ASV file.

Thank you

# Answer
I am not sure that kraken/kraken2 does anything similar to this.

Usually to get OTUs/ESVs I think that you want to use a tool that clusters your sequences into groupings that represent the same sequence, allowing for technical error.

Examples might include dada2 or dadasnake, or qiime. This is after the technical steps for adapter/trimming/denoising steps.

kraken2 is more for classifying into known taxonomic categories based on a reference sequence, which is a related but distinct task.

I don't think that I understand what the ESV file format is supposed to look like. This resource seems to have some info on processing data into ESVs, but it is not for nanopore. Possibly it is still useful for the final step of creating that file format. A linked repository has a ton of one-off scripts that may be helpful but are not very well documented. For example, here is a script that they use to generate TSV files which *may* be the "ASV"/"ESV" format:

```
#!/bin/bash
source activate qiime2-2018.8

for item in `ls */*table.qza`
    do
    name=`echo $item | cut -d\/ -f1`
    qiime tools export --input-path $item --output-path EXPORTS-GTDB-tax/$name
    mv EXPORTS-GTDB-tax/$name/feature-table.biom EXPORTS-GTDB-tax/$name.biom
    rmdir EXPORTS-GTDB-tax/$name
done

qiime tools export --input-path 11_ANT-28-5-SPA_seqs_classified_against_GTDB_r86.1_807F-1050R_primer_sliced/classification.qza --output-path EXPORTS-GTDB-tax/

sed -i '1c#OTUID    taxonomy    confidence' EXPORTS-GTDB-tax/taxonomy.tsv

#Tutorial here: https://forum.qiime2.org/t/exporting-and-modifying-biom-tables-e-g-adding-taxonomy-annotations/3630
for item in `ls EXPORTS-GTDB-tax/*biom`
    do
    filestem=`basename $item .biom`
    biom add-metadata -i $item -o EXPORTS-GTDB-tax/$filestem.with-tax.biom --observation-metadata-fp EXPORTS-GTDB-tax/taxonomy.tsv --sc-separated taxonomy
done

for item in `ls EXPORTS-GTDB-tax/*.with-tax.biom`
    do
    filestem=`basename $item .biom`
    biom convert -i $item -o EXPORTS-GTDB-tax/$filestem.tsv --to-tsv --header-key taxonomy
done

```
> 1 votes

---
Tags: ngs, nanopore
---


# Question
Title: Snakemake scatter-gather with wildcard AmbiguousRuleException
My problem is when using Snakemake scatter-gather feature the documentation is basic and i modified my code according to mentioned in this link

```
rule fastq_fasta:
    input:rules.trimmomatic.output.out_file
    output:"data/trimmed/{sample}.fasta"
    shell:"sed -n '1~4s/^@/>/p;2~4p' {input} > {output}"

rule split:
    input:
        "data/trimmed/{sample}.fasta"
    params:
        scatter_count=config["scatter_count"],
        scatter_item = lambda wildcards: wildcards.scatteritem
    output:
        temp(scatter.split("data/trimmed/{{sample}}_{scatteritem}.fasta"))
    script:
        "scripts/split_files.py"

rule process:
    input:"data/trimmed/{sample}_{scatteritem}.fasta"
    output:"data/processed/{sample}_{scatteritem}.csv"
    script:
        "scripts/process.py"

rule gather:
    input:
        gather.split("data/processed/{{sample}}_{scatteritem}.csv")
    output:
        "data/processed/{sample}.csv"
    shell:
        "cat {input} > {output}"

```

I added wildcard option but, i got `AmbiguousRuleException: Rules fastq_to_fasta(which is previous rule) and split are ambiguous for the file data/trimmed/Ornek_411-of-81-of-81-of-81-of-81-of-81-of-81-of-81-of-81-of-8.fasta`

I tried lots of things but either rules are not calling or take AmbiguousRuleException. What am i missing, can someone help?

# Answer
It's not a problem with the scatter-gather feature per-se; Snakemake is just confusing itself with what text should match what wildcards with these filenames. You could fix this by defining a wildcard constraint for your samples, and/or changing your filename patterns to clarify what should match what. The docs on scatter-gather show an example with different directory names and file suffixes that prevent this from happening. Probably the simplest fix would be to follow that idea and use different directory names for the split parts, like `trimmed-split` and `processed-split`.

A hint of what's happening is the repeating `...-of-81-...` thing, and the mention of ambiguity between rules. (See also the `--debug-dag` option.) If you have a sample named Ornek\_411 and you want to make data/processed/Ornek\_411.csv, Snakemake will look for data/processed/Ornek\_411\_1-of-N.csv and so on, and try to figure out how to make those. But a filename like that could be the output of either the process rule (`sample=Ornek_411` and `scatteritem=1-of-N`) or the gather rule itself again (`sample=Ornek_411_1-of-N`). Obviously you want the former and not the latter, but it doesn't know that automatically.
> 1 votes

---
Tags: python, snakemake
---


title: Multifasta protein alignment

text:
# Question
Title: Multifasta protein alignment
I have around 4000 orthologous proteins for each of 5 species, I'd like to align them and to concatenate alignments. Which tools can I use? And how?

Marco

# Answer
4000 isn't much. You'll need to start with an aligner. For aligners the common options are in no particular order:

The current one that 'outperforms' the others strictly on a combination of two published benchmarks comprising 'accuracy' and 'speed' is MAFFT. FAMSA isn't relevant here unless you've a really small computer, like an old laptop, it excels at speed.

However, the situation could have changed recently with the BioRXIV release of Muscle v5. Please note, I have not used it. Here's the Abstract:

> On Balifam, ensembles generated by MUSCLE are shown to align an average of 59% of columns correctly, 13% better than Clustal-omega (52% correct) and 26% better than MAFFT (47% correct). The ensemble bootstrap is applied to a previously published tree of RNA viruses ...

On this benchmark they are suggest Clustal Omega is better than MAFFT. Its not been the run of play that is certain MAFFT rose to dominance in recent history. Des Higgins Clustal W and X was historically the aligner of choice by a long shot, then Muscle, then MAFFT ... now its context specific.

*Muscle v5* Its difficult to really get to grips with a BioRXIV ms as an algorithm unless you're a specialist in aligner construction. Please note, the 'RNA viruses' is complicated. As an algorithm this author has kudos and has extended that kudos beyond Muscle. The BioRXIV would have been subject to peer-review, but I am not aware Muscle v5 has been accepted for publication.

Historically, Muscle made early advances in accuracy which were duplicated by Clustal Omega at a much later date (2017?).

Essentially, here your choice is between MAFFT and Muscle v5 for this data set (context specific).

`concatenate` this is less clear but it appears you are performing 4000 alignments with 5 taxa in each? You'll need a data pipeline.

Trimming and concatenation are easily performed in BioPython using `BioAlignIO` but you would need coding because that can't be done manually. Computationally, `BioAlignIO` objects are annoyingly expensive - they are `SeqIO` objectives inside a `MultipleAlignment` object. However in your case it matters not a jot and would the recommended route all things considered (unless your coding expertise is Java).
> 1 votes

---
Tags: proteins, sequence-alignment, multi-fasta
---


# Question
Title: Is there a way to extract spike sequence from a Sars-CoV-2 complete genome (preferably in Python)?
I have a complete genome sequence of a Sars-CoV-2 variant, but I am only interested in the Spike sequence?

I should note that I am not a biologist (nor bioinformatics).

The complete genome looks like the following:

```
AAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAAAATCTGTGTGGCTG
TCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGGACAC
GAGTAACTCGTCT...

```

Its length is 30230.

My question is, is there a way to know where exactly the spike sequence lies among this complete genome sequence, and if possible, are there already existing tools to extract it programmatically (in Python). Thanks in advance.

# Answer
Spike Wuhan is here

Using Python then use the BioPython `Bio.pairwise2` program here

Yours sequence is the reference, the shorter Wuhan query is the query. The genomic positions `position1` and `position2` will then be defined via the output of `pairwise2`, then cut it by using the genomic position -1 as follows:

```
spike = completegenome[int(position1)-1:int(position2)-1]

```

`-1` is needed because lists start at 0 rather than 1.
> 0 votes

---
Tags: genome, covid-19, sars-cov-2, python
---


# Question
Title: Pymol: select low confidence regions from AlphaFold pdb file
I have downloaded a predicted structure from AlphaFold as a pdb file (https://alphafold.com/entry/O75376) and loaded it into Pymol (2.3.0). There is quite a large portion of the structure that was modelled with very low confidence (pLDDT \< 50). Is there a way to select these regions in the structure so that I can assign a different representation to them? Thank you.

# Answer
> 2 votes
*Answer from @matteo-ferla, converted from comment:*

They are stored as b-factors. This is well described in the PyMol wiki.

*\[please edit to improve this answer, if possible\]*

# Answer
> -1 votes
use: Spectrum b, Red\_white\_blue

---
Tags: protein-structure, pymol, structural-biology
---


# Question
Title: How to identify genomic regions / peaks associated with enhancers (TF binding sites)? Is there a tool or a formal recipe?
I am attempting to identify DNA sequences/regions associated with enhancers in my NGS sequencing data. How can I identify (from peak files/counts for chip-seq, atac etc.) genomics regions/peaks associated with enhancers? I know proximity to genes and possibly histone (i.e. h3k27ac) activity near known genes is possibly an indicator, but is this the only way and how is the computational process formalized/captured? What are the specific set of tools involved? i.e. Peak annotations etc.? Thank you.

As I understand it in my limited bio knowledge, enhancers are at varying distances from genes and sometimes they are multiples (super-enhancers). Are all enhancers already "known/described" (akin to genes)? What is/was the process of finding them (human manual procedure)? It's that more manual/heuristic process (without ML) I am looking for so we can consider formalizing and automating it.

Based on a very insightful answer: I think some input/clarity on the following bit would make the answer(s) even more valuable:

Possible list of steps towards this include following tools/functions:

1. Bowtie (to align fastq to a genome, mouse, human etc.)
2. Homer/makeTagDirectory (Sep files, associate chr/coords with tags/sequence reads)
3. Homer/findPeaks (identify and list genomic regions/sequences with count values above a threshold)
4. Homer/pos2bed (convert peaks to bed file, genome browser compatible form)
5. Homer/annotatePeaks.pl (Runs on \*bed file - Don't know much about this, but intuitively associate a label with peaks found in 4?)

Question: Is there a step 6? Where do we go after annotatedPeaks for example to get to enhancers?

# Answer
> 2 votes
If you have reads from chip-seq or atac you can use tools like homer or macs to identify peaks in your data. Once you have the peaks, these programs allow you to do motif discovery to annotate your peaks with e.g. what transcription factor binding sites they contain.

---
Tags: r, rna-seq, chip-seq, machine-learning, atac-seq
---


# Question
Title: PCA analysis of samples in a phylogenetic tree
I have a phylogenetic tree.

Each branch ends have samples (s1, s2, ~ s16).

What I want to do is, I want to make PCA analysis plot for each sample.

First I thought each sample has each lineage (expressed in image as (1,1,1,1)~(2,2,2,2)).

So with that coordinates, I can do PCA analysis.

But on the second thought, the difference between (1,1,1,1) and (1,1,1,2) is much less than that of (1,1,1,1) and (2,1,1,1), but in the analysis, the difference between them will be same.

I am not sure if this kind of approach is suitable for making a PCA plot with a phylogenetic tree.

How can I make this work?

# Answer
> 1 votes
There appear two basic issues here:

1. The PCA would represent arbitrary lineage designation not a genetic quantification.
2. The PCA will weight the first, second ... fourth position equally, i.e. it is not hierarchical. This is very important for the analysis.

*Summary* The hierarchical -\> unweighted PCA will cause analytical problems and should not be used, given it is ultimately arbitrary. If there was a tree specialist reviewing your work it would not be accepted.

The reason is that position one in the hierarchical cluster defines the entire data set - it cuts it in half, whilst position 4 defines a small a very small subtree. PCA will not see the hierarchy it will see position 1 as of equal worth to position 4. Therefore PCA will start associating arbitrary designated position 4 e.g. branch 2, with each other. Thus 1,1,1,2 is a long way from 2,1,1,2 but PCA will miss that it will think the shared '2' means their closer than they really are and find it of equal distance to e.g. 1,2,1,2 - which or course it isn't.

*Fixing it* You ask how to fix it

1. You can go back to the original data and use that to create a PCA based on SNPs ... this has been done for other genetic data types, but cannot done reinterpreting a tree. Using this approach there will be a problem with back mutation, but its a parsimony analysis and you would simply state your assumption - ultimately the influence of back mutation concerns the genetic distance between all taxa in the analysis and if thats not big, the error will be small.
2. Alternatively if you had a phylogram you could construct a genetic matrix and use genetic distance instead. Thats a separate question because its quite complicated.

tSNE is often used to refine a multivariate analysis and I think thats what you getting at, PCA followed by tSNE is valid. A tree isn't any regular multivariate analysis, its a tightly defined approximation within clear mathematical boundaries measuring/reconstructing evolutionary divergence. A tree specialist will see this, whilst a general multivariate stats specialist will miss it.

---
Tags: phylogenetics, pca, tsne
---


# Question
Title: how to retain reads with low mapping quality (MAPQ) scores when using samtools view -q
I have been using the `-q` option of `samtools view` to filter out reads whose mapping quality (MAPQ) scores are below a given threshold when mapping reads to a reference assembly with either `bwa mem` or `minimap2`. However, using `samtools idxstats` to count total mapped reads and unmapped reads indicates that these reads with lower MAPQ scores are actually being discarded altogether, since the number of unmapped reads reported (last row) is always zero when the BAM file is created with a `-q` argument, even though fewer reads are being mapped to the assembly.

For example:

```
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b | samtools sort --reference assembly.fasta > nanopore_mapped_no_q_thresh.bam
samtools idxstats nanopore_mapped_no_q_thresh.bam

```

gives

```
contig_1      4751250 9897    0
contig_2      37697   528     0
contig_3      9131    2366    0
contig_4      98997   701     0
contig_5      4539    365     0
contig_6      3933    3164    0
*       0       0       142

```

while

```
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b -q 2 | samtools sort --reference assembly.fasta > nanopore_mapped_q2.bam
samtools idxstats nanopore_mapped_q2.bam

```

gives

```
contig_1      4751250 9716    0
contig_2      37697   528     0
contig_3      9131    1415    0
contig_4      98997   700     0
contig_5      4539    365     0
contig_6      3933    3161    0
*       0       0       0

```

Is there an alternative way to filter mapped reads by their MAPQ scores so that the filtered reads are still retained in the BAM file?

# Answer
> 2 votes
Looking at the samtools view docs, I believe that you need to also set the `-p` flag if you want to retain "filtered" reads in the file and simply have them be unmapped:

> The -p, option sets the UNMAP flag on filtered alignments then writes them to the output file.

---
Tags: sam, samtools, genomics, read-mapping
---


# Question
Title: PWM (4xM matrices) to Sequence Logo Visualizations in MATLAB or Python
Is there a MATLAB, or else Python (not R) tool for visualizing sequence Logos from probability weight matrices (PWMs)?

This seems like a basic and simple tool to implement, yet just above the threshold in the space of tedious to write my own from scratch - I was surprised I couldn't find a basic function like a

```
visualizePWM(pwm_matrix)

```

that would take a pwm\_matrix that is a (4 by m) matrix and spit out the visualization of an m-length sequence motif/logo like:

I see all kinds of over-complicated bloated tools that I thought I ask in case I am missing something.

# Answer
> 1 votes
*answer from @M\__, converted from comments:*

This has been answered on stackoverflow using matplotlib for seq logo plot:

https://stackoverflow.com/questions/42615527/sequence-logos-in-matplotlib-aligning-xticks

IMO its only useful for presentation, pretty pictures, thus there's no need for someone to write a library (which is what you're asking). If I'm using a mutation matrix (which I do), I want the raw numbers and significance values I don't want pretty graphs getting in the way.

Conceptually you would move from 2D to 1D, without subplots (which is doable), i.e. its a matrix. My view on code availability is not universal. I've had requests for separate CLIs (useful for someone clearly) ... which just never happened.

# Answer
> 1 votes
After much searching and trial and error, located possibly the optimal tool for my question:

**seqlogo**

https://github.com/betteridiot/seqlogo

Surprisingly not only seqlogo is an elegant and minimal solution in Python (and R) I have discovered it has a relative implementation in MATLAB! It installed and worked the first try - beautiful output with all the output options I needed - rock solid.

**Python** one liner (not counting the example pwm input generation)

random\_ppm = np.random.dirichlet(np.ones(4), size=16)

ppm = seqlogo.Ppm(random\_ppm)

**seqlogo.seqlogo(ppm ,format='png',filename='out.png')**

Produces the following ('out.png')

And here is **MATLAB** Port:

https://www.mathworks.com/help/bioinfo/ref/seqlogo.html

---
Tags: rna-seq, python, visualization, motifs, matlab
---


# Question
Title: Retrieve chromosome number with just RSID?
I have ~500,000 rsids - but want which chromosome each SNP resides on. How can I do this?

```
MarkerName A1 A2
rs2326918 a g
rs7929618 c g
rs66941928 t c

```

I want the chromosome number associated with each MarkerName. So the output would be like this:

```
MarkerName A1 A2 Chr
rs2326918 a g     6
rs7929618 c g     6
rs66941928 t c    6

```

How can i do this?

# Answer
> 2 votes
Download dbsnp VCF file from https://ftp.ncbi.nih.gov/snp/organisms/human\_9606/VCF/. This contains rs IDs and chromosome info (among other info). Can match your rs IDs to the rs IDs in dbsnp VCF, and extract chromosome from there.

See https://www.ncbi.nlm.nih.gov/variation/docs/human\_variation\_vcf/ for more information on the dbsnp files.

---
Tags: bioconductor
---


# Question
Title: snakemake multiple parameters for multiple input and single output in snakemake. ConbineGVCFs gatk problem
I have written a rule for CombineGVCFs in gatk4. The rule is as follow

```
all_gvcf = get_all_gvcf_list()

rule cohort:
  input:
    all_gvcf_list = all_gvcf,
    ref="/data/refgenome/hg38.fa",
    interval_list = prefix+"/bedfiles/hg38.interval_list",
  params:
    extra = "--variant",
  output:
    prefix+"/vcf/cohort.g.vcf",
  shell:
    "gatk CombineGVCFs -R {input.ref} {params.extra} {input.all_gvcf_list} -O {output} --tmp-dir=/data/tmp -L {input.interval_list}"

```

all\_gvcf is the dataset for all gvcf files which will be combined. But the problem is I need to add --variant parameter before every input. Command I'm getting right now is as follows

```
gatk CombineGVCFs -R /data/refgenome/hg38.fa --variant /data/prjna644607/vcf/SRR12165216_HC.g.vcf /data/prjna644607/vcf/SRR12165217_HC.g.vcf /data/prjna644607/vcf/SRR12165218_HC.g.vcf /data/prjna644607/vcf/SRR12165219_HC.g.vcf -O /data/prjna644607/vcf/cohort.g.vcf --tmp-dir=/data/tmp -L /data/prjna644607/bedfiles/hg38.interval_list

```

The command I want to achieve is as follows

```
gatk CombineGVCFs -R /data/refgenome/hg38.fa --variant /data/prjna644607/vcf/SRR12165216_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165217_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165218_HC.g.vcf --variant /data/prjna644607/vcf/SRR12165219_HC.g.vcf -O /data/prjna644607/vcf/cohort.g.vcf --tmp-dir=/data/tmp -L /data/prjna644607/bedfiles/hg38.interval_list

```

How can I add this extra "--variant" tag before every input? I have added it in the get\_all\_gvcf\_list() function. But then snakmake gives me inputfiles not found problem.

# Answer
> 5 votes
Found out the problem. Turns out I can write a lambda function as follows

```
params:
        extra=lambda wildcards, input: ' -V '.join(input.all_gvcf_list)

```

and add '-V' before {params.extra} in the "shell" sction. That solves the problem

# Answer
> -1 votes
well I came up with the solution. I think it is a little unconventional.

```
    rule gatk_CombineGVCFs:

        input:
            vcf_dummy = expand("gvcf/{sample}.g.vcf.gz", sample = SAMPLES), # a dummy vcf to connect this rule to gatk_HaplotypeCaller
            fa = "data/genome.fa",

        output:
            gvcf = "gvcf_combined/all.g.vcf.gz",
            #index = directory("gvcf_index")                                 #index generated without seperate output
        params:
            list = "gvcfs.list"
            log: "logs/CombineGVCFs.log"
        shell:
            "ls gvcf/*.gz > gvcfs.list && gatk --java-options -Xmx7g CombineGVCFs -R {input.fa} --variant {params.list} -O {output.gvcf} "
             #Two bash comands running sequentially to make a gvcf.list as during my testing the combineGVCFs tool takes lists as input.

```

---
Tags: python, snakemake, gatk
---


title: My model doesn't work in test dataset

text:
# Question
Title: My model doesn't work in test dataset
1. I got the data of stomach cancer from TCGA, I found some different expression genes between cancer and normal sample.
2. I put aside my normal sample, and divide my cancer sample into train(80%) and test(20%) groups.
3. I did Unicox (with survival package), lasso (with glmnet package) and multicox ( with survival package) on train group, to find out which genes are more related to survival.
4. I use coefficients of each genes (that I got from multicox), to built my model and calculate the risk score. Like that: (expression of gene1×coefficient gene1 + expression of gene2×coefficient gene2 +...).
5. I divided my sample in to high risk and low risk groups, to compare survival between them, using Kaplan-Meier.

In training group, Kaplan-Meier shows significant difference between high and low risk group. However when I use my model on my test group, it is not significant. Why this problem occurs? I assumed my model is overfit of my train group, but how can I fix it? Why does this happen even when I use lasso?

# Answer
There is no overtraining here, it isn't deep learning and there appears no parameterisation. If parameterisation was occurring on the training set that can and does cause over-tightening.

Essentially, you haven't included a data split for parameterisation, i.e. there is no parameterisation. The theory of a parameterisation is it must be distinct from training (and obviously) testing. It allows parameter tuned without overtightening. Lasso regression will require parameterisation.

I don't know about the survival regressions for machine learning, but lasso/ridge regression are standard approaches.

---

So just to explain, what you do is have a split for parameters. Thus instead of 80:20 split it could be 70:20:10 split where the 10% is used for parameter optimisation. Bad code risks leakage from the parameter split and that can result in over-tightening.
> 0 votes

---
Tags: r
---


# Question
Title: MetaQuast for assembling samples from complex communities
I'm working with whole genome metagenomic samples from human skin, and I'm using MEGAHIT for assembly and MetaQuast for evaluation. However, MetaQuast requires a list of reference genomes for the analysis (if I understood correctly). When I perform the evaluation without the references, it uses BLASTN for aligning contigs to SILVA 16S rRNA database. Since I'm working with whole genome sequences, this is not an optimal choice. As expected, MetaQuast is identifying only two species in my contigs file, although I have manually checked through BLASTN that there are contigs that match other species' genomes. And this evaluation makes my assembly looks terrible. Since I'm working with samples from human skin, it wouldn't be viable to download a database of reference genomes and perform MetaQuast using a local alignment. So, what would be the best approach to perform the evaluation correctly in my case?

Looking forward to discussing it!

# Answer
> 1 votes
I would suggest first using a binning approach such as MetaBAT2 (or an equivalent tool, there are many) to separate your metagenomic assembly into bins roughly corresponding to genomes. I would then run CheckM to assess the bins' quality/completeness/purity.

If you still want to use quast, you could then run metaquast with a reference genome appropriate to each bin independently. I would suggest GTDB-tk to perform a classification and identify such an appropriate reference genome.

---
Tags: assembly, metagenome, wgs
---


# Question
Title: Python/Biopython - Replace amino acid residue on MSA with "z" from a list of unaligned positions
I'm trying to programmatically replace a set of amino acid residues on an MSA with a "Z" from a list of unaligned positions. Any ideas on how I could do this?

Input: a list of unaligned positions for each sequence + an MSA
Output: an MSA with those unaligned positions on the aligned MSA replaced with a Z for each sequence.

I have a multiple sequence alignment i'm working with and a list of sites on that alignment that are glycosylated. I was trying to run a script that would take the input of the positions as a list of characters along with the MSA to replace the Nth character with a "Z" to easily visualize the sites. So essentially if I feed it the following sequence

```
M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K

```

and positions

```
2, 5, 10, 11, 14

```

it should spit back

```
M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

```

# Answer
> 1 votes
The code is

```
import re

pos = [2, 5, 10, 11, 14]
newseq = []
prot = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        newseq.append('-')
    else:
        if (n + 1 - count) in pos:
            s = re.sub(r'[A-Za-z]', 'Z', s)
        newseq.append(s)
print(''.join(newseq))

```

output

> M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

Precisely as requested, apologies for the delay.

Happy to explain why it works - its just a couple of Python tricks and a method. It must be run under Python 3 (unlikely to work in Python 2) and probably \>3.6.

BTW I normally avoid `re` for complicated reasons but its easy, hence its here.

---

As @gaspanic has pointed out `re` is not only unnecessary it will slow the code (`re` is very inefficient). It is in fact one of the 'legendary inefficiencies' of Python. Thus removing it ...

```
# variable declarations omitted
newseq = list()
count = 0
for n,s in enumerate(prot):
    if str(s) == '-':
        count += 1
        newseq.append('-')
    else:
        if (n + 1 - count) in pos:
            s = 'Z'
        newseq.append(s)
print(''.join(newseq))

```

Same output.

# Answer
> 1 votes
Another alternative, without the need for `re`:

```
seq = 'M----KFLAFLCLLGFANAQ-------------D-GKCG----TLSNKSPS--------------------K'
pos = [2, 5, 10, 11, 14]

new_seq = []
count = 0
for s in seq:
    if s == "-":
        new_seq.append(s)
    else:
        count += 1
        if count in pos:
            new_seq.append('Z')
        else:
            new_seq.append(s)

print(''.join(new_seq))

```

Result:

```
M----ZFLZFLCLZZFAZAQ-------------D-GKCG----TLSNKSPS--------------------K

```

---
Tags: fasta, python, proteins
---


# Question
Title: How to identify which isoforms of a gene are actually expressed in my data?
I am new to bioinformatics and have been assigned the task of discovering which isoforms of a certain protein-coding gene are present in the RNA-seq data that I have. There are several isoforms of this gene listed on Entrez & Ensembl and I want to know which of these isoforms are expressed in the RNA-seq data that I have been given. I am not sure how to begin to tackle this issue so If anyone has any ideas or advice, it would be greatly appreciated.

# Answer
Tools like kallisto or salmon can perform a process call pseudo-alignment that matches reads from RNA-seq data to transcript isoforms present in a reference genome. You can then look at the output of these tools to see which isoforms for a given gene appear to be present in your RNA-seq data. These tools answer the question "how abundant is each transcript in my data?".

But if you want to answer a question like "are there differences in which transcripts are present between experimental conditions?", I'd look for tutorials on differential expression analysis, differential splicing analysis, or differential transcript usage.
> 4 votes

# Answer
If you do as @james-hawley rightly suggests, keep in mind that pseudo-alignment could detect a lot of transcripts with low TPM that are actually false positives, as reported here (it's a paper focused on scRNA-seq, but many statements are also valid for bulk). So, you may want to try to remove false positive transcripts, that will likely result with very low read count and TPM. Salmon's decoy function could help here.

As an alternative to pseudo-alignment, you could run StringTie on your BAM file to detect transcript isoforms, using a reasonable threshold of coverage over splice junctions.
> 3 votes

---
Tags: rna-seq, ensembl, entrez, isoform
---


# Question
Title: What do "PnBS" and "PnMS" mean in FCS flowCore parameters?
I'm using the flowCore packages to analyze cytometry data, and I am trying to understand all the parameters stored in the FCS file. I've figured out most of the elements, but despite a thorough search I haven't been able to find documentation on some of the parameters.

For each channel, there are 9 associated parameters. Here's P1 and its associated values from one of my files :

* `$P1N` FSC-A
* `$P1R` 262144
* `$P1B` 32
* `$P1E` 0,0
* `$P1V` 225
* `$P1G` 1
* `$P1DISPLAY` LIN
* `$P1BS` 0
* `$P1MS` 0

I've figured out the following from various sources (notably the FCS data standards, given here) or from the fields themselves:

* `$PnN` Short name for parameter n.
* `$PnR` Range for parameter number n.
* `$PnB` Number of bits reserved for parameter number n.
* `$PnE` Amplification type for parameter n.
* `$PnV` Detector voltage for parameter n.
* `$PnG` Amplifier gain used for acquisition of parameter n.
* `$PnDISPLAY` the transformation to apply when displaying the data (e.g. LIN, LOG)

However, these 2 don't appear in the file standard description, and flowCore does not have documentation on the parameters either.

* `$PnBS`
* `$PnMS`

I've seen `$PnBS` set at -1, 0, 48, 57, 54, 52 or 25 for various channels in my data, while `$PnMS` always seems to be at 0.

Does someone know what it is, or where to look for a definition? If that helps, my FCS files were generated through BD FACSDiva Software Version 9.0.1.

# Answer
> 0 votes
I found my answer in an old mailing list. Credit to Josef Spidlen who authored the mail.

*Short answer*: both parameters are proprietary keywords from the BD FACSDiva software for storing scaling values.

*Long answer* from the mailing list:

> **PnBS** Indicates the automatically set Biexponential scale value for Instrument Parameter (-1 indicates an invalid scale value). When importing an FCS file, if the FCS file has valid Biexponential scales, they should be applied to the data. If the FCS file doesn't have valid biexponential scales then the scales are calculated from the entire data. If biexponential scales exist, then if AUTOBS is set, the PnBS values are used and applied to the data. If AUTOBS is not set, the PnMS values are used and applied to the data. The information represented by this keyword becomes part of a Tube's attributes, and does not appear in the Tube's list of keywords. It can be viewed in the Biexponential inspector' when the auto scaling option is selected.
>
> **PnMS** Indicates the manually set Biexponential scale value for Instrument Parameter (-1 indicates an invalid scale value). When importing an FCS file, if the FCS file has valid Biexponential scales, they should be applied to the data. If the FCS file doesn't have valid biexponential scales then the scales are calculated from the entire data. If biexponential scales exist, then if AUTOBS is set, the PnBS values are used and applied to the data. If AUTOBS is not set, the PnMS values are used and applied to the data. The information represented by this keyword becomes part of a Tube's attributes, and does not appear in the Tube's list of keywords. It can be viewed in the Biexponential inspector' when the manual scaling option is selected.

---
Tags: r, flow-cytometry
---


# Question
Title: Gene/protein expression specific to a group in omics
I am wondering what is the significance of finding a particular protein specific to a disease or control group? when we detect 1000s of proteins in a proteomics experiment, how can one be sure that the specific proteins found in say, a disease group is really absent in control group than by random chance?

# Answer
You can never be 100% sure that a difference is real and not a statistical fluke. This is the entire concept underlying Type I errors - these are errors where the null hypothesis (usually, that there is no difference) is incorrectly rejected (a difference is called significant, when in fact none truly exists). Especially when performing many hypothesis tests, it's possible that a variable with truly identical distributions between two groups appears to have very different distributions simply due to sampling variability.

You can reduce the chance of making such an error by reducing the level of alpha used to determine significance, and by performing appropriate multiple hypothesis correction. With stringent significance thresholds and conservative correction, you can make it extremely unlikely that a "significant" difference could be due to random chance, but it's not usually possible to conduct an experiment with a Type I error rate of 0%.
> 2 votes

---
Tags: rna-seq, genomics, proteomics, mass-spectrometry
---


title: Reading lines from a fasta format file

text:
# Question
Title: Reading lines from a fasta format file
```
file_name = "consen1.txt"
    with open(file_name) as infile:
        new_dict = {}
        for line in infile:
            if line.startswith(">"):
                key = line.strip('\n')
            else:
                line = ''.join(x.strip('\n') for x in infile)
                new_dict[key] = line
    print(new_dict)
    for key,values in new_dict.items():
        gc = values.count('G') + values.count('C')
        content = (gc/len(values)) *100
        new_dict[key] = content
    print(new_dict)

```

I am reading from a text file, and everything works fine but for some reason it is skipping the line of sequence from the text file. Here's how the text file looks like,

```
>Rosalind_4478
GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGAC
CGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGC
TTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGAT
CGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCAC
TCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTC
CTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAA
TGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGC
TGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGT
ATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCC
GTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCA
AGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCT
CGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTA
CTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCT
TTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTG
TGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAG
TA

```

Expected output:

```
 {'>Rosalind_4478': 'GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGACCGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}

```

Real Output:

```
{'>Rosalind_4478': 'CGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}

```

I do realize that I can do this kind of computation easily with modules like biopython, But I wanted to try doing this from scratch. Thanks in advance for the help.

# Answer
Your approach is a bit risky code ... Anyway not tested but I think this will work:

```
with open(file_name) as infile:
    new_dict = {}
    for l in infile:
        line = l.rstrip('\n')
        if line.startswith(">"):
            k = line.lstrip('>')
            new_dict[k] = [] # must initiate
            # new_dict = {k:[]} is possible alternative initiation
        else:
            new_dict[key].append(line)

```

I not sure without testing which dictionary initiation will work. People do use `defaultdict`

```
from collections import defaultdict
new = defaultdict(def_value)

```

This circumvents the initiation issue to a large extent (if the dictionary key ain't correctly initiated it will throw a value error, or some sort of bug). However it isn't actually needed and I find it unnecessary.

What I am certain about is the following should be avoided because your code depends on the exact line by line order between the fasta header and the code you are associating it with. You likely break that association with the code below,

```
line = ''.join(x.strip('\n') for x in infile) # you are calling the entire infile not the line

```

BTW

```
from collections import Count ()
new = Count()

```

Is a very good approach for the next part of your script, i.e. counting %GC
> 1 votes

# Answer
needed to edit accepted code to make it works:

```
file_name = "consen1.txt"

with open(file_name) as infile:
    new_dict = {}
    for l in infile:
        line = l.rstrip('\n')
        if line.startswith(">"):
            k = line.lstrip('>')
            new_dict[k] = str() # must initiate like string !!!!!!!!!!!!!
            # new_dict = {k:[]} #is possible alternative initiation
        else:
            new_dict[k] += line

print(new_dict)

for key,values in new_dict.items():
    # print(key, values)
    gc = values.count('G') + values.count('C')
    content = (gc/len(values)) *100
    new_dict[key] = content

print(new_dict)

```

output :

```
{'Rosalind_4478': 'GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGACCGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}
{'Rosalind_4478': 49.55654101995565}

```

works using this snippet too. but feel like is not orthodox python:

```
file_name = "consen1.txt"

with open(file_name) as infile:
        new_dict = {}
        for line in infile:
            if line.startswith(">"):
                key = line.strip('\n')

                print(key)

                line_2 = None

            else:

                if not line_2:

                    line_2 = ''.join(line.strip('\n'))

                else:

                    line_2 += ''.join(line.strip('\n'))

                # print(line_2)

                new_dict[key] = line_2

print(new_dict)

for key,values in new_dict.items():
    gc = values.count('G') + values.count('C')
    content = (gc/len(values)) *100
    new_dict[key] = content

print(new_dict)

```

output:

```
>Rosalind_4478
{'>Rosalind_4478': 'GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGACCGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}
{'>Rosalind_4478': 49.55654101995565}

```

be sure to check if they are both working with your input
> 2 votes

---
Tags: python
---


# Question
Title: How To display bar graph with error bars?
I am using the below code to get the bar plot. I want to get the bar plot with error bars displayed.

```
import seaborn as sns

def plot_coefficients(target, coeffs_df, n=10):
    if coeffs_df.shape[0]==0:
        return None # nothing to plot here!
    data = coeffs_df.copy()
    data = data.sort_values(('Coefficient','mean'), ascending=False).head(n)
    g = sns.barplot(
        data=data,
        y='Gene',
        x=('Coefficient','mean'),
        hue='Feature Type',
        hue_order=[
            h for h in [
                'Genomic Variant',
                'Gene Expression',
                'Proteomics',
                'Methylation',
            ]
        ],
        dodge=False
    )
    g.set(
        xlim=(0,1),
        xlabel='Importance [0,1]',
        ylabel='',
        title=f"Top {target} predictor features"
    )
    return g

```

Coefficient here is a score generated from predictive modelling that ranges from 0 to 1. 2. data here is entire data and I created a mean column. Ideally I wanted to calculate it in seaborn plot rather than creating a mean column separately before.

# Answer
> 1 votes
I did put forward a comment, it was removed I think the idea is its a full post or zero. The problem proving a negative is tricky.

I don't think its doable within the exact system specified, i.e. a Seaborn bar chart. Again I could be wrong.

If the OP went to a Seaborn dot line formation - that would work `sns.barplot` is changed to `sns.pointplot` but then all the nice colouring disappears. Standard deviation is one error range if all the values are present: but I think the OP will need to request it (`.sd()` ?? can't remember).

In `matplotlib` there is an error bar method where the coder specifies the error bars and that is `plt.errorbar`. An example from the official docs is here. Thus call a bar chart and then add the error bars onto it.

**Summary** The reasons I don't think its doable in Seaborn:

1. I'm pretty certain Seaborn has no independent errorbar method unlike matplotlib
2. Seaborn can do error bars given the full data set, but those are created automatically for a given method. Maybe there is a "sns.barplot-errorbar" method but I don't know it.

The Seaborn approach does carry the advantage of not having to work out the error bar independently.

---
Tags: python, visualization
---


# Question
Title: Bacterial genomes and snpeff warnings 'WARNING_TRANSCRIPT_NO_START_CODON'
I annotated a bunch of genomes with snpeff. And although the original reads' qualities and the matches with the reference genomes were quite high, I see the above mentioned warning in almost every variant in my vcf file. Thus, I wonder how this could happen. Since I'm using bacterial genomes with a manually built database (from these reference genomes) and there are quite some genes (especially in gene clusters) that have the alternative start codon GTG, I suspect this causing the warning. Any idea how to test if this is really the reason for the warning?

Thanks a lot in advance!

Edit: As requested, here is the html snpEff report. And I've also uploaded an example vcf file after annotation here

# Answer
> 2 votes
Bacteria have a different translation table from eukaryote nuclear DNA; it may be that SnpEff is expecting human genes in the genome. Have you followed the build steps outlined in the documentation?

https://pcingola.github.io/SnpEff/se\_build\_db/

> If your genome, or a chromosome, uses non-standard codon tables you must update snpEff.config accordingly. A typical case is when you use mitochondrial DNA. Then you specify that chromosome 'MT' uses codon. Invertebrate\_Mitochondrial codon table. Another common case is when you are adding a bacterial genome, then you specify that the codon table is Bacterial\_and\_Plant\_Plastid.

# Answer
> 1 votes
According to your report, you are using human HG18 as the reference genome version. Therefore, the report will be wrong along with the wrong codons warnings. It is unfruitful to try with a human database with bacteria.

---
Tags: genome, bacteria, codon, snpeff
---