title: Coexpression network analysis

text:
# Question
Title: Coexpression network analysis
I have a question related with the most popular frameworks for network construction: WGCNA for weighted networks and igraph package for unweighted networks.

What are the advantages or disadvantages of applying these approaches and when it is recommended to use one over another.

# Answer
> 3 votes
`igraph` is a serious general network theory framework for any given data science application. It permits full directionality and describes vertices (nodes), edges (branches) and does cool stuff like De Bruijn graphs and Laplacian graphs. Understanding graph theory would be helpful and will be a steep learning curve. It's a serious data science package.

`WGCNA` is a specific ready made RNA expression package which appears to be network theory and is undirectional, i.e. no directionality which is fine for RNA expression. It will have far less overall power but if you've a load of RNA-seq data and want 'plug and play' type stuff WGCNA is the only way. It's a correlation network, based on *r*2 (Pearson's correlation), to assess co-expression. I would assume the weighting is geared towards identifying co-expression. *There is an error in Fig. 2 legend* BTW Langfelder & Horvath (2008): the legend states r^2 what they mean is *r*2. I guess they messed up on their LaTeX.

If this is a one of coexpression analysis `WGCNA`. If this is a big project, graph theory is central to that effort and you're looking towards developing a data science tool set then its `igraph`. There are other graph theory packages BTW.

---

In case my answer wasn't clear `igraph` will do weighted edges via `is_weight`. The information is here, basically it will do anything.

If your question is "what does weighting contribute to coexpression analysis", that is a good question. I know graph theory, I dunno the biology of gene expression versus the maths. What `WGCNA` is doing is statistically is simple if weighting is excluded ... my reasoning is it must be pivotal to its success.

---

**To address the OPs comments.**

*Statistical power* In the pure sense, I was incorrect to use this term because technically 'lack of power' is the inability to reject a hypothesis when it is false. For example, parametric statistics are more powerful than non-parametric statistics for this reason.

What I meant was two things:

1. In graph theory you can do anything within the "network", what is called a "graph", e.g. find the centre of the graph, minimal distance across the "graph", cyclisation, anything you can think of.

1a. `WGCNA` is fixed it only does what is described in Langfelder & Horvath (2008). Nothing can be done outside that without rewriting the underlying code and it would be better to use `igraph` instead if that was the objective in my opinion.

Thus `igraph` is more "powerful" because it will do a lot more stuff.

2. My concern of `WGCNA` is its over-reliance on *r*2. This is a parametric statistic and may not always be applicable but from my reading of the paper (albeit briefly) it was entirely dependent on that. To centre everything on *r*2, personally I would at a minimum use *r* in addition (Spearman's correlation). This is extremely complex data and there is a risk of over-simplification. *r* and *r*2 can certainly give different results.

The answer to the question however it is not what I think but everyone else thinks and Langfelder & Horvath (2008) has almost 10000 citations. I am not aware that the calculation has been overturned, e.g. DSeq2 has been overturned and demonstrated to be inappropriate. So I would assume it must be capturing the overall picture, so the positives outweigh potential negatives.

---
Tags: r, networks, wgcna, coexpression
---


# Question
Title: Struggling to install scvi-tools on Mac M1
I was wondering if somebody here has experience with installing scvi-tools on a Mac with M1 chip. I tried the following command:

```
conda install scvi-tools -c conda-forge

```

But got the following error:

```
Collecting package metadata (current_repodata.json): failed
CondaSSLError: OpenSSL appears to be unavailable on this machine. OpenSSL is required to
download and install packages.
Exception: HTTPSConnectionPool(host='conda.anaconda.org', port=443): Max retries exceeded with url: /conda-forge/osx-64/current_repodata.json (Caused by SSLError("Can't connect to HTTPS URL because the SSL module is not available."))

```

Then I tried to install OpenSSL with following command:

```
conda install -c anaconda openssl

```

and I got the exact same error as the conda install scvi-tools -c conda-forge command.

I do not understand why collection of metadata is failing for openssl. Any inputs would be highly helpful. Thanks!

# Answer
> 2 votes
Did you try:

```
pip install scvi-tools

```

?

---
Tags: software-installation
---


title: Definition of genotype in demuxlet

text:
# Question
Title: Definition of genotype in demuxlet
I am reading the Online methods of the demuxlet paper.

The genotype $g$ is taken from the set $\\{0, 1, 2\\}$, defined as "the true genotype of the sample corresponding to $c$-th droplet at $v$-th variant". What is the mapping between these values and the corresponding genotype?

Each variant is modeled using a binomial with $n=2,p=\frac{g}{2}$, so I would assume that 1 is heterozygous, 0 is homozygous REF and 2 homozygous ALT. Is this assumption correct?

# Answer
The genotypes $\\{0,1,2\\}$ are mappings to the genotypes AA (homref), AB (heterozygous), BB (homalt) (from supplementary resources). This is also mapped to the number of alterations (e.g., 0 -\> 0 alterations -\> homref).
> 1 votes

---
Tags: single-cell, genotyping, demultiplexing
---


# Question
Title: How to extend the x axis in Dimplot Seurat
I have performed a Seurat PCA via Dimplot.

How do I extend the x axis?

As you can see in my figure the double x axes overlap.

# Answer
Looking at the source code for Seurat's `DimPlot()` it's based on ggplot2 graphics (https://github.com/satijalab/seurat/blob/HEAD/R/visualization.R#L794), so you can manipulate the dimplot object and save it using ggplot2 commands. Using `ggsave()` you can specify the figure dimensions; if you increase the width of the final figure it should 'fix' your label overlap issue, i.e.

```
#BiocManager::install("Seurat")
library(Seurat)
#> Attaching SeuratObject
library(ggplot2) # for plotting
library(cowplot) # for the theme used by Seurat

# example data
data("pbmc_small")
p <- DimPlot(object = pbmc_small, split.by = "ident")

# example dimplot
p

```

```
# increase font size to make the labels overlap
p +
  theme_cowplot(font_size = 28) +
  theme(strip.background = element_blank())

```

```
# save the 'overlapping label' figure and specify a relatively large width to 'space out' the three facets
ggsave("dimplot_test.png", width = 40, height = 10, units = "cm")

```

---

You could also decrease the font size and keep the original dimensions to fix the problem, or alter the breaks on the x axis, e.g.

```
# reduce the size of the x axis text labels
p +
  theme_cowplot() +
  theme(strip.background = element_blank(),
        axis.text.x = element_text(size = 8))

```

```
# change the number of breaks
p +
  theme_cowplot() +
  theme(strip.background = element_blank()) +
  scale_x_continuous(breaks = scales::pretty_breaks(n = 2))

```

```
# Specify the breaks you want to use
p +
  theme_cowplot() +
  theme(strip.background = element_blank()) +
  scale_x_continuous(breaks = c(-20, 0, 20, 40))

```

Created on 2022-12-01 with reprex v2.0.2
> 1 votes

# Answer
As a direct answer to your question, restrict the output to producing a single column:

```
Dimplot(..., ncol=1)

```

But this doesn't seem like a useful answer. Please provide more context to your question if you want more help than this.
> 0 votes

---
Tags: seurat, single-cell
---


# Question
Title: Tutorial for the WGCNA: changes in heatmap colours
I am trying to reproduce the results of the R tutorial of the WGCNA package.

In section I number 5, when generating the heatmap it is quite similar to the one provided by the pdf but the colours of the heatmap differ from the document. Is it possible to change it?

I attach the image together with the code used.

\[

Image Obtained

Image of tutorial

```
nSelect = 400
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
sizeGrWindow(9,9)
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

```

# Answer
> 1 votes
You can add a color palette to the plot.

e.g.

```
my_palette <- colorRampPalette(c("red", "orange", "yellow","white"))(n = 100);
TOMplot(plotDiss, selectTree, selectColors, main ="Network heatmap plot, selected genes", col=my_palette)

```

---
Tags: rna-seq, wgcna
---


# Question
Title: Maximum likelihood estimation for tree construction
Could someone direct me to a resource (textbook chapter, lecture notes, online video, etc.) that demonstrates, in a mathematically or computationally rigorous way, how to use maximum likelihood estimation to generate a phylogenetic tree, given a set of aligned sequences?

I've found several resources which describe the procedure for maximum likelihood estimation, such as:

1. giving the formula and the parameters that need to be determined (which I already know);
2. how to use it for aligned sequences (I already know this), I am seeking a worked-out example.

Alternatively, if someone is willing to work out a small example, how would you use maximum likelihood estimation to generate a tree for the following three sequences?

```
organism_1 ATCG
organism_2 ATCC
organism_3 ATTT

```

# Answer
@DataScienceNovice ... here's a rough explanation. Just to note a tree with three taxa doesn't really work as a tree - you need 4 taxa.

1. If we have **tree 1** ((A,B),(C,D)) the probability of the tree is expressed as the summed log likelihood for each position 1, 2, 3 .... in the alignment. The probability is derived from a mutation matrix, thus certain mutations e.g. transitions `A<->G`; `C<->T` are more probable (high likelihood), e.g. transversions `A<->C`; `A<->T`; `G<->C`; `G<->T` are less probable in accordance with the matrix. This hypothetical matrix (most matrices) is scoring transitions as much more probable than transversions.
2. The 'total probability' expressed log likelihood for each column, i.e. each alignment position is added together leveraging an interesting property of logs (see final sentence of this paragraph). That becomes the 'score' or 'likelihood' for that tree. Note in the example you provided you have 4 alignment positions, so 4 probabilities of the log-likelihood will be added. Note, the probabilities are actually multiplied but log-likelihoods added is the same thing.
3. We will now rearrange the tree lets try **tree 2** ((C,D),(B,G)), again the hypothetical matrix scores transition/transversions, thus essentially if this tree has higher transitions than the tree in point 1. that will provide a higher likelihood. The tree that maximises the transitions **as sister groups** will maximise the likelihood. Equally maximising transversions as the most **divergent groups** will maximise the likelihood of that tree.
4. Re-arrange the tree again **tree 3** ... ((C,B),(D,A)) ... which tree has the highest transitions vs transversions? ...
5. In this example Tree 2 has the highest likelihood. Thus that is the maximum likelihood. Simple?

In the example you have the outcome will depend on the matrix. Alignment position 3 comprises a transition, whilst position 4 comprises a transition and a transversion. Generally transversions are 1:2 to 1:4 less likely even though there are twice as many. So the example you gave the 'tree' (you can't really have a tree with 3 taxa) would be ((`organism2`,`organism3`),`organism1`). This is against the total number of mutations because `organism2` and `organism3` has two mutations whilst `organism2` and `organism1` only have one, but the mutation types are different (also see point 3 below). In other words the occurrence of a transversion between `organism1` against `organism2`,`organism3`, mean `organism2` and `organism3` are more likely **sister groups** even though they comprise 2 mutations (transitions) against the single mutation between `organism1` and `organism2`.

A good maximum likelihood program for phylogenetics is **RAxML**.

---

1. The reality is it's not that simple because the mutation matrix is a function of the tree, which is a function of the mutation matrix and in technical tree building the matrix is also calculated alongside the total probability of the tree via maximum likelihood. Just to be clear it can be calculated as a separate maximum likelihood calculation and then fixed for the tree calculation, but it's not a good idea because of non-linear interactions.
2. In addition, mutation rate (heterogeneity) is also calculated via maximum likelihood.
3. The final thing to note is back-mutations are also calculated and that refers to past posts on this site, which I will link if I find them.

**Complexity** If the parameters are calculated via maximum likelihood simultaneous to the tree the parameter space becomes EXTREMELY large and the calculation is searching the tree landscape against the parameter landscape which becomes a multi-dimensional calculation.

**Reality** In reality much of a maximum likelihood calculation becomes about the search strategy because its impossible to traverse all of the parameter space so a heuristic estimate needs to be performed and the calculation can get stuck on lower plateaus of probability and fail to identify the tree maximum likelihood. This is the nuts and bolts of real world tree building.

**Note** However if you fix the matrix and the rate, which a lot of people do, then points 1 to 5 describe the maximum likelihood tree.
> 1 votes

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


# Question
Title: Getting translation table of organisms through Entrez direct (or other CLIs)
I have a long list of NCBI taxonomic ids for organisms and I want to write a script to get their corresponding translation tables. Is there a way to do this through Entrez?

(packages in R, Python etc. would be fine too)

# Answer
> 1 votes
I came up with this idea. So, the information for the translation table of any particular organism can be found on its page on the Taxonomy Browser: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi

If I search a taxid into the website, eg 7227, the page for *Drosophila melanogaster* will come up but the URL will remain unmodified. However, you can get a modified URL specific to a taxid by appending ?mode=Info&id=*your\_taxid* to the above URL. Eg for *D. melanogaster* again, it has this unique URL for its Taxonomy browser page: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=7227

Since you have a list of taxids, you can programmatically access Taxonomy browser pages by modifying the URL in the way I said. To programmatically get the translation tables of *D. melanogaster* with this approach, use the urllib package in Python.

```
import urllib.request
url = "https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=7227"
weburl = urllib.request.urlopen(url)
data = weburl.read().decode("utf-8")  # decode because page returned in binary

data_s = data.split('<')

for d in data_s:
    if "Translation table" in d:
        print(d.split('>')[1])

```

The output is:

```
Translation table 1 (Standard)
Translation table 5 (Invertebrate Mitochondrial)

```

This will work for one taxid. But you have a list of them. However, with this approach in mind, it's elementary to construct a Python function for getting the Translation table of all your taxids.

```
import urllib.request

my_taxids = ['7227','9606','10090']

def get_translation_table(taxid):
    url = "https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?mode=Info&id=" + taxid
    weburl = urllib.request.urlopen(url)
    data = weburl.read().decode("utf-8")

    data_s = data.split('<')

    tr_table = []

    for d in data_s:
        if "Translation table" in d:
            tr_table += [d.split('>')[1]]

    return tr_table

your_tables = list(map(get_translation_table,my_taxids))

print(your_tables)

```

The output is:

```
[['Translation table 1 (Standard)', 'Translation table 5 (Invertebrate Mitochondrial)'], ['Translation table 1 (Standard)', 'Translation table 2 (Vertebrate Mitochondrial)'], ['Translation table 1 (Standard)', 'Translation table 2 (Vertebrate Mitochondrial)']]

```

And if you want the results nicely formatted into a pandas dataframe, just add:

```
import pandas as pd

```

To the top of the script, and

```
t = pd.DataFrame(your_tables,columns=["Nuclear","Mitochondrial"])
print(t)

```

To the bottom. Output of this print (if you have too many taxids use the head method instead so not everything prints):

```
                          Nuclear                                     Mitochondrial
0  Translation table 1 (Standard)  Translation table 5 (Invertebrate Mitochondrial)
1  Translation table 1 (Standard)    Translation table 2 (Vertebrate Mitochondrial)
2  Translation table 1 (Standard)    Translation table 2 (Vertebrate Mitochondrial)

```

# Answer
> 0 votes
You can find them here at least: https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi

You can download the complete list of codes (using FTP) here.

However, to link a specific organism to a specific codon usage table, you'll need to do a bit of extra work. Perhaps the script here can give you a starting point?

---
Tags: phylogenetics, entrez
---


title: How to make bar bigger in UpsetR plot

text:
# Question
Title: How to make bar bigger in UpsetR plot
I have tried to use UpsetR to visualize the input file which can be found here.

```
library("UpSetR")
orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)

#All species
selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"

head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367

ncol(orthogroups_df)
[1] 14

orthogroups_df[orthogroups_df > 0] <- 1
# we only show intersections of interest
intersections=list(list(selected_species),
                   list("Nlab", "NQLD", "Ngla"),
                   list("Nlab", "NQLD", "Nsyl"),
                   list("Nlab", "Ngla", "Nsyl"),
                   list("NQLD", "Nsyl", "Ngla"))
upset(orthogroups_df,
      text.scale = c(1.4, 1.4, 1.4, 1.4, 1.4, 1.4),
      sets=rev(selected_species),
      nsets = ncol(orthogroups_df),
      #keep.order=T,
      #mb.ratio=c(0.5,0.5),
      #order.by='degree',
      #order.by='freq',
      intersections = intersections,
      sets.x.label="Total nbr of orthogroups",
      mainbar.y.label = "Nbr of orthogroups",
      #empty.intersections = "on",
      group.by = "intersections")

```

How is possible to make the bar bigger of the double digits values or is a way to rotate the plot?

# Answer
I'm not quite sure what you mean. You can alter the % of height that the bottom part of the graph takes up, that will make all the vertical bar parts proportionately taller, but I don't think you should because you have so many groups to display in the bottom part.

I think what you'd want to do is to alter the y-axis so that it stops being totally proportional, to compress the size of the tallest bar, so you can better see the difference in bar sizes between the other comparisons?

I'm not sure you can do that. Making the y axis logged with "scale.intersections" might help.

You'd probably be best off making a second graph that omits that "common to all" intersection. Once that is out of the picture, the difference between the others will be more visible.
> 1 votes

---
Tags: r, ggplot2, upsetr
---


# Question
Title: Concatenate multiple fasta files into one file for MLST
I have seven fasta files, one per gene, with more than 400 fasta entries per file, like this: Input:

Gene1.fasta

```
>1721_1 gene name
ATG
>10_14 gene name
GTT

```

Gene2.fasta

```
>1721_1 gene name
TGA
>10_14 gene name
GAT

```

Output:

```
>1721_1
ATGTGA
>10_14
GTTGAT

```

What I want to do is create one fasta (continuously) file containing seven genes per sample. As you can see, the sample name is not the same before the _ in the fasta headers.

Is there any way to do this with awk or other command-line tools? Thank you very much.

# Answer
> 2 votes
The `concat` command from seqkit could concatenate sequences with same ID from multiple files.

```
seqkit concat [file1] [file2] ...

```

# Answer
> 0 votes
I think it can be done using python dictionary. You can split the file using the header and then make the fasta id as your key and the sequence as your value. Then using simple loop and conditions you can concatenate them.

# Answer
> 0 votes
Assuming your files are small enough to allow you to fit one of them into memory, you can do this:

```
$ awk '{
        if(/^>/){ name=$0 }
        else{
            if(NR==FNR){seq[name]=$0 }
            else{
                print name;
                printf "%s%s\n", seq[name],$0
            }
        }
      }' file1.fa file2.fa
>1721_1 gene name
TGA
>10_14 gene name
GTTGAT

```

If you cannot fit the files in memory, you can try this:

```
$ join -o1.1,1.2,2.2 -t $'\t' <(FastaToTbl file1.fa ) <(FastaToTbl file2.fa) | sed 's/\t//2' | TblToFasta
>1721_1 gene name
ATGTGA
>10_14 gene name
GTTGAT

```

Note that this approach requires the `FastaToTbl` and `TblToFasta` scripts I have posted here: https://bioinformatics.stackexchange.com/a/2651/298.

Or, if brevity is more your thing, the same approach but more compact:

---
Tags: fasta, linux, awk
---


# Question
Title: How to make unrooted tree for Likelihood mapping result by using IQ2-tree?
I am a biologist, and I do not fully understand the tree topology of the experimental species. I used four-taxon set (4 sequences) to identify the Four-cluster Likelihood-Mapping by using

```
iqtree2 -s ***.fa -lmap 2000 -n 0

```

I have obtained the svg file - triangle (attached - 1), however there were no data mentioned about T1, T2 and T3 unrooted tree based on my input 4 sequences (Haur\_2, Shad\_6, Equa\_oga2 and apa). I used single orthologous copy genes as the supermatrix for this data.

Therefore, I do not understand which topology of tree is reliable depends on the sequence data. I knew that my data set showed similar result topology from Iq2 tree output. There is no species sequence data available for the species that I am looking for a reliable topology of those species. If anyone knows how to understand the assigned sequence on the unrooted tree based on Four-cluster Likelihood-Mapping by using iqtree2? Could you kindly answer this question? I want to know the orange unrooted tree (image 2) for my analyzed data. Thank you very much.

output \***.iqtree Quartet support of regions a1, a2, a3 (mainly for clustered analysis):

```
     #quartets    a1 (% a1)        a2 (% a2)        a3 (% a3)    name
-----------------------------------------------------------------------------
1   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  Haur_2
2   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  Shad_6
3   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  Equa_oga2
4   2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )  apa
-----------------------------------------------------------------------------
2000     648     (32.40 ) 712     (35.60 ) 640     (32.00 )

```

# Answer
You have provided two different trees. The second was the Arthropoda and labeled "image 2".

The first example the explanation is that it is an exact polytomy. I know this because the values are equally split between the 3 possible alternatives for a 4-taxon tree. Thus there is no bifurcated solution to the data, any bifurcating combination is equally valid therefore it is simply expressed as all 4-taxa sharing the same node (hypothetical common ancestor). In other words there are 3 equally valid solutions.

Thus, the tree is 4-taxa tree with 1-node. Thats the solution. A normal 4-taxon tree will have 2-nodes, this is a star-tree, i.e. one node, technically a polytomy.

The example below concerns Arthropoda and demonstrates the top tree ((Chelicerta, Myriapoda),(Hexapoda, misc.)) is the correct tree. That 90.4 value demonstrates that. What this says is that neither myriapods or Chelicerta share a direct common ancestor with hexapods, moreover "misc" (whatever that is) shares a direct ancestor with the hexapods. "Image 2" is a classical bifurcating tree and represents only one solution.

---

So to specifically answer the question. The tree for arthropoda is:

and the other tree is,

I should have removed the scale bar in this example, because I didn't have the data, but not when you have genetic quantification of the distances.
> 1 votes

---
Tags: rna-seq, phylogenetics, phylogeny, transcriptome, orthologues
---


# Question
Title: Downloading genes from NCBI in fasta format
I'm pretty new in bioinformatics.

I need to download FASTA sequences of several genes. The list of genes I have assembled consists of 140 genes, so I'd rather do this through via code than download each gene manually from the NCBI database.

The genes belong to 2 different organisms (bacteria) I have chosen, most of them are found in both organisms and some are specific to one or the other bacterium.

I've tried running this code which I found here: https://www.ncbi.nlm.nih.gov/books/NBK179288/

```
esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efetch -format fasta

```

The gene of interest being aminoacyl-tRNA hydrolase here. But I don't know how to filter it for the specific organism, Neisseria meningitidis. Tried running this, without luck:

```
esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efilter -organism [Neisseria meningitidis] | efetch -format fasta

```

Appreciate any help :)

# Answer
If I put the species name in quotes, I get expected (by me) output:

```
% esearch -db nucleotide -query "aminoacyl-tRNA\ hydrolase" | efilter -organism "Neisseria meningitidis" | efetch -format fasta
# outputs:
>NZ_JWNU01000010.1 Neisseria meningitidis strain Nm440716 Scaffold10, whole genome shotgun sequence
GACCGCCTATTCGCTTGCCGCCGGCGGCCCCATCATGCAGGCAGGATTACACGCCTTCACACTCGTACCC
# [snip, tons of output]

```

Note that filtering down to only *N. meningitidis* still leaves all of the different sequenced accessions of this species, such that you will still have a lot of different annotated copies of this gene. You will need to filter more specifically to reduce that.
> 2 votes

---
Tags: fasta, phylogenetics, data-download, sequence
---


title: Can I run dotPlotly from Rstudio?

text:
# Question
Title: Can I run dotPlotly from Rstudio?
I'd like to run dotPlotly from RStudio but I can't find any description/help for this. Is it at all possible? Another side to this question: does anyone know any tool for dot plots for genome comparison apart from pafr (doesn't provide colour legend) and d-genies (we have problems with installation)? Is there a way to make dot plots for a minimap2 output with ggplot2?

To specify:

1. I need help in running dotPlotly (preferably from Rstudio)
2. I'd like to build dot plots for genome comparisons showing % of identity

Thanks

# Answer
> 1 votes
For the issue of just making a dotplot in ggplot2 from a PAF file, it is pretty trivial to do a simple one, though it is necessarily not the most beautiful.

I grabbed 2 random E. coli genomes and mapped them:

```
% minimap2 ~/Downloads/GCF_000157115.2_Escherichia_sp_3_2_53FAA_V2_genomic.fna.gz ~/Downloads/GCF_000308975.1_ASM30897v2_genomic.fna.gz > 2_ecolis.paf

```

Then I opened Rstudio and ran the following (which is only slightly modified from what I was first shown by an old colleague, Zev Kronenberg):

```
ecoli_paf = read.table("~/2_ecolis.paf", sep="\t", fill=T)
ggplot() +
geom_segment(
  data=ecoli_paf, mapping=aes(
    x=ecoli_paf$V3, y=ecoli_paf$V8, xend=ecoli_paf$V4, yend=ecoli_paf$V9, color=ecoli_paf$V6
  )
)

```

I am far from a ggplot2 expert as you can tell, and there are a few issues here (reverse complement is shown opposite to what you'd expect, linear interpolation of mappings from mapped endpoints), but it makes a somewhat informative dotplot:

I am confident that with minimal fiddling you can make this look a lot better. I will try to revisit and beautify.

---
Tags: r, ggplot2
---


title: How to automate NCBI genome download

text:
# Question
Title: How to automate NCBI genome download
I need to download all the **completely assembled cyanobacterial genome**'s GenBank file(.gbff) from NCBI(RefSeq or INSDC ftp data).

For this I think, the steps are:

1. Need to find the completely assembled genomes.
2. find the GenBank file URL based on the taxonomic name.
3. download the GenBank file(.gbff file)

Is their any way to do this by using any python module or any other idea??

# Answer
> 3 votes
Outline of solution:

Edit:
Since I will need something like that soon, I made a Perl script for downloading genomes.
For Cyanobacteria, you would do `./download-refseq-genomes.pl 1117`

# Answer
> 2 votes
As an alternative, you can use the NCBI Datasets CLI, in a two step process to generate a dehydrated dataset for later hydration:

```
# Download metadata for all Cyanobacteria
# And pull out accession for 'Complete' genomes
./datasets assembly-descriptors taxon Cyanobacteria --limit ALL \
     | jq -r '.assemblies[].assembly | "\(.assembly_accession) \(.assembly_level)"' \
     | grep Complete \
     | cut -f 1 -d ' ' > accs.txt

# Download a dehydrated dataset for the provided set of accessions
./datasets download assembly -i accs.txt -f 1117_complete.zip \
     --dehydrated \
     --exclude-protein \
     --exclude-rna \
     --exclude-seq \
     --exclude-gff3 \
     --include-gbff

```

You may be interested in some of the other files that are excluded.

Once you have the dehydrated dataset, you can rehydrate it when you are ready to analyze the data with:

```
unzip 1117_complete.zip -d 1117_complete
./datasets rehydrate -f 1117_complete

```

# Answer
> 0 votes
Use the NCBI Datasets command line tool to download genomes by taxon (NCBI Taxonomy ID, scientific or common name at any tax rank).

```
datasets download genome taxon cyanobacterial \
  --assembly-source refseq --assembly-level complete_genome \
  --exclude-genomic-cds --exclude-gff3 --exclude-protein \
  --exclude-rna --exclude-seq --include-gbff \
  --filename genomes.zip

```

---
Tags: python, phylogenetics
---


# Question
Title: heatmap of specific sequence motifs in aligned fasta files
I have a collection of fasta files, each containing three aligned sequences. I am interested in understanding the distribution of a specific sequence motif PGP, RGP and KGP in all the alignments. I was wondering if there is possible way in R to create a heatmap highlighting the location of these motifs.

I have tried the code in this post to create a heatmap but it is not working as expected. One issue is that ggmsa does not import the fasta file. The fasta file has to be copied into the installation directory for it to work. I have close to 50 aligned files, so not very efficient. The other major problem is it creates a heatmap compared to a reference sequence and not the motifs. Can't figure how to change this.

Any help is much appreciated. Two fasta files containing aligned sequence are as follows,

```
>chain A CO1A1_HUMAN/1-1464 Collagen alpha-1(I) chain
MFSFVDLRLLLLLAATALLTHGQEEGQVEGQDEDIPPITCVQNGLRYHDRDVWKPEPCRICVCDNGKVLCDDVICDETKNCPGAEVPEGECCPVCPDGSESPTDQETTGVEGPKGDTGPRGPRGPAGPPGRDGIPGQPGLPGPPGPPGPPGPPGLGGNFAPQLSYGYDEKSTGGISVPGPMGPSGPRGLPGPPGAPGPQGFQGPPGEPGEPGASGPMGPRGPPGPPGKNGDDGEAGKPGRPGERGPPGPQGARGLPGTAGLPGMKGHRGFSGLDGAKGDAGPAGPKGEPGSPGENGAPGQMGPRGLPGERGRPGAPGPAGARGNDGATGAAGPPGPTGPAGPPGFPGAVGAKGEAGPQGPRGSEGPQGVRGEPGPPGPAGAAGPAGNPGADGQPGAKGANGAPGIAGAPGFPGARGPSGPQGPGGPPGPKGNSGEPGAPGSKGDTGAKGEPGPVGVQGPPGPAGEEGKRGARGEPGPTGLPGPPGERGGPGSRGFPGADGVAGPKGPAGERGSPGPAGPKGSPGEAGRPGEAGLPGAKGLTGSPGSPGPDGKTGPPGPAGQDGRPGPPGPPGARGQAGVMGFPGPKGAAGEPGKAGERGVPGPPGAVGPAGKDGEAGAQGPPGPAGPAGERGEQGPAGSPGFQGLPGPAGPPGEAGKPGEQGVPGDLGAPGPSGARGERGFPGERGVQGPPGPAGPRGANGAPGNDGAKGDAGAPGAPGSQGAPGLQGMPGERGAAGLPGPKGDRGDAGPKGADGSPGKDGVRGLTGPIGPPGPAGAPGDKGESGPSGPAGPTGARGAPGDRGEPGPPGPAGFAGPPGADGQPGAKGEPGDAGAKGDAGPPGPAGPAGPPGPIGNVGAPGAKGARGSAGPPGATGFPGAAGRVGPPGPSGNAGPPGPPGPAGKEGGKGPRGETGPAGRPGEVGPPGPPGPAGEKGSPGADGPAGAPGTPGPQGIAGQRGVVGLPGQRGERGFPGLPGPSGEPGKQGPSGASGERGPPGPMGPPGLAGPPGESGREGAPGAEGSPGRDGSPGAKGDRGETGPAGPPGAPGAPGAPGPVGPAGKSGDRGETGPAGPAGPVGPVGARGPAGPQGPRGDKGETGEQGDRGIKGHRGFSGLQGPPGPPGSPGEQGPSGASGPAGPRGPPGSAGAPGKDGLNGLPGPIGPPGPRGRTGDAGPVGPPGPPGPPGPPGPPSAGFDFSFLPQPPQEKAHDGGRYYRADDANVVRDRDLEVDTTLKSLSQQIENIRSPEGSRKNPARTCRDLKMCHSDWKSGEYWIDPNQGCNLDAIKVFCNMETGETCVYPTQPSVAQKNWYISKNPKDKRHVWFGESMTDGFQFEYGGQGSDPADVAIQLTFLRLMSTEASQNITYHCKNSVAYMDQQTGNLKKALLLQGSNEIEIRAEGNSRFTYSVTVDGCTSHTGAWGKTVIEYKTTKTSRLPIIDVAPLDVGAPDQEFGFDVGPVCFL
>chain B CO1A1_HUMAN/1-1464 Collagen alpha-1(I) chain
MFSFVDLRLLLLLAATALLTHGQEEGQVEGQDEDIPPITCVQNGLRYHDRDVWKPEPCRICVCDNGKVLCDDVICDETKNCPGAEVPEGECCPVCPDGSESPTDQETTGVEGPKGDTGPRGPRGPAGPPGRDGIPGQPGLPGPPGPPGPPGPPGLGGNFAPQLSYGYDEKSTGGISVPGPMGPSGPRGLPGPPGAPGPQGFQGPPGEPGEPGASGPMGPRGPPGPPGKNGDDGEAGKPGRPGERGPPGPQGARGLPGTAGLPGMKGHRGFSGLDGAKGDAGPAGPKGEPGSPGENGAPGQMGPRGLPGERGRPGAPGPAGARGNDGATGAAGPPGPTGPAGPPGFPGAVGAKGEAGPQGPRGSEGPQGVRGEPGPPGPAGAAGPAGNPGADGQPGAKGANGAPGIAGAPGFPGARGPSGPQGPGGPPGPKGNSGEPGAPGSKGDTGAKGEPGPVGVQGPPGPAGEEGKRGARGEPGPTGLPGPPGERGGPGSRGFPGADGVAGPKGPAGERGSPGPAGPKGSPGEAGRPGEAGLPGAKGLTGSPGSPGPDGKTGPPGPAGQDGRPGPPGPPGARGQAGVMGFPGPKGAAGEPGKAGERGVPGPPGAVGPAGKDGEAGAQGPPGPAGPAGERGEQGPAGSPGFQGLPGPAGPPGEAGKPGEQGVPGDLGAPGPSGARGERGFPGERGVQGPPGPAGPRGANGAPGNDGAKGDAGAPGAPGSQGAPGLQGMPGERGAAGLPGPKGDRGDAGPKGADGSPGKDGVRGLTGPIGPPGPAGAPGDKGESGPSGPAGPTGARGAPGDRGEPGPPGPAGFAGPPGADGQPGAKGEPGDAGAKGDAGPPGPAGPAGPPGPIGNVGAPGAKGARGSAGPPGATGFPGAAGRVGPPGPSGNAGPPGPPGPAGKEGGKGPRGETGPAGRPGEVGPPGPPGPAGEKGSPGADGPAGAPGTPGPQGIAGQRGVVGLPGQRGERGFPGLPGPSGEPGKQGPSGASGERGPPGPMGPPGLAGPPGESGREGAPGAEGSPGRDGSPGAKGDRGETGPAGPPGAPGAPGAPGPVGPAGKSGDRGETGPAGPAGPVGPVGARGPAGPQGPRGDKGETGEQGDRGIKGHRGFSGLQGPPGPPGSPGEQGPSGASGPAGPRGPPGSAGAPGKDGLNGLPGPIGPPGPRGRTGDAGPVGPPGPPGPPGPPGPPSAGFDFSFLPQPPQEKAHDGGRYYRADDANVVRDRDLEVDTTLKSLSQQIENIRSPEGSRKNPARTCRDLKMCHSDWKSGEYWIDPNQGCNLDAIKVFCNMETGETCVYPTQPSVAQKNWYISKNPKDKRHVWFGESMTDGFQFEYGGQGSDPADVAIQLTFLRLMSTEASQNITYHCKNSVAYMDQQTGNLKKALLLQGSNEIEIRAEGNSRFTYSVTVDGCTSHTGAWGKTVIEYKTTKTSRLPIIDVAPLDVGAPDQEFGFDVGPVCFL
>chain C CO1A2_HUMAN/1-1366 Collagen alpha-2(I) chain
----------------------------------------------------------------------------------------MLSFVDTRTLLLLAVTLCLATCQSLQEETVRKGPAGDRGPRGERGPPGPPGRDGEDGPTGPPGPPGPPGPPGLGGNFAAQYDGKGVGLGPGPMGLMGPRGPPGAAGAPGPQGFQGPAGEPGEPGQTGPAGARGPAGPPGKAGEDGHPGKPGRPGERGVVGPQGARGFPGTPGLPGFKGIRGHNGLDGLKGQPGAPGVKGEPGAPGENGTPGQTGARGLPGERGRVGAPGPAGARGSDGSVGPVGPAGPIGSAGPPGFPGAPGPKGEIGAVGNAGPAGPAGPRGEVGLPGLSGPVGPPGNPGANGLTGAKGAAGLPGVAGAPGLPGPRGIPGPVGAAGATGARGLVGEPGPAGSKGESGNKGEPGSAGPQGPPGPSGEEGKRGPNGEAGSAGPPGPPGLRGSPGSRGLPGADGRAGVMGPPGSRGASGPAGVRGPNGDAGRPGEPGLMGPRGLPGSPGNIGPAGKEGPVGLPGIDGRPGPIGPAGARGEPGNIGFPGPKGPTGDPGKNGDKGHAGLAGARGAPGPDGNNGAQGPPGPQGVQGGKGEQGPPGPPGFQGLPGPSGPAGEVGKPGERGLHGEFGLPGPAGPRGERGPPGESGAAGPTGPIGSRGPSGPPGPDGNKGEPGVVGAVGTAGPSGPSGLPGERGAAGIPGGKGEKGEPGLRGEIGNPGRDGARGAPGAVGAPGPAGATGDRGEAGAAGPAGPAGPRGSPGERGEVGPAGPNGFAGPAGAAGQPGAKGERGAKGPKGENGVVGPTGPVGAAGPAGPNGPPGPAGSRGDGGPPGMTGFPGAAGRTGPPGPSGISGPPGPPGPAGKEGLRGPRGDQGPVGRTGEVGAVGPPGFAGEKGPSGEAGTAGPPGTPGPQGLLGAPGILGLPGSRGERGLPGVAGAVGEPGPLGIAGPPGARGPPGAVGSPGVNGAPGEAGRDGNPGNDGPPGRDGQPGHKGERGYPGNIGPVGAAGAPGPHGPVGPAGKHGNRGETGPSGPVGPAGAVGPRGPSGPQGIRGDKGEPGEKGPRGLPGLKGHNGLQGLPGIAGHHGDQGAPGSVGPAGPRGPAGPSGPAGKDGRTGHPGTVGPAGIRGPQGHQGPAGPPGPPGPPGPPGVSGGGYDFGYDGDFYRADQPRSAPSLRPKDYEVDATLKSLNNQIETLLTPEGSRKNPARTCRDLRLSHPEWSSGYYWIDPNQGCTMDAIKVYCDFSTGETCIRAQPENIPAKNWYRSSKDKKHVWLGETINAGSQFEYNVEGVTSKEMATQLAFMRLLANYASQNITYHCKNSIAYMDEETGNLKKAVILQGSNDVELVAEGNSRFTYTVLVDGCSKKTNEWGKTIIEYKTNKPSRLPFLDIAPLDIGGADQEFFVDIGPVCFK

```

and

```
>chain A CO4A1_HUMAN/1-1669 Collagen alpha-1(IV) chain
---------MGPRLSVWLLLLPAALLLHEEHSRAAA--KGGCAGSGC-GKCDCHGVKGQKGERGLPGLQGVIGFPGMQGPEGPQGPPGQKGDTGEPGLPGTKGTRGPPGASGYPGNPGLPGIPGQDGPPGPPGIPGCNGTKGERGPLGPPGLPGFAGNPGPPGLPGMKGDPGEILGHVPGMLLKGERGFPGIPGTPGPPGLPGLQGPVGPPGFTGPPGPPGPPGPPGEKGQMGLSFQGPKGDKGDQGVSGPPGVPGQA-------QVQEKGDFATKGEKGQKGEPGFQGMPGVGEKGEPGKPGPRGKPGKDGDKGEKGSPGFPGEPGYPGLIGRQGPQGEKGEAGPPGPPGIVIGTGPLGEKGERGYPGTPGPRGEPGPKGFPGLPGQPGPPGLPVPGQAGAPGFPGERGEKGDRGFPGTS-LP-GPSGRDGLPGPPGSPGPPGQPGYTNGIVECQPGPPGDQGPPGIPGQPGFIGEIGEKGQKGESCLICDIDGYRGPPGPQGPPGEIGFPGQPGAKGDRGLPGRDGVAGVPGPQGTPGLIGQPGAKGEPGEFYFDLRLKGDKGDPGFPGQPGMPGRAGSPGRDGHPGLPGPKGSPGSVGLKGERGPPGGVGFPGSRGDTGPPGPPGY---GPAGPIGDKGQAGFPGGPGSPGLPGPKGEPGKIVP--------------------LPGPPGAEGLPGSPGFPGPQGDRGFPGTPGRPGLPGEKGAVGQPGI-GFPGPPGPKGVDGLPGDMGPPGTPGRPGFNGLPGNPGVQGQKGEP---GVGLPGLKGLPGLPGIPGTPGEKGSIGVPGVPGEHGAIGPPGLQGIRGEPGPPGLPGSVGSPGVPGI-GPPGARGPPGGQGPPGLSGPPGIKGEKGFPGFPGLD-MPGPKGDKGAQGLPGITGQSGLPGLPGQQGAPGIPGFPGSKGEMGVMGTPGQPGSPGPVGAPGLPGEKGDHGFPGSSGPRGDPGLKGDKGDVGLPGKPGSMDKVDMGSMKGQKGDQGEKGQIGPIGEKGSRGDPGTPGVPGKDGQAGQPGQP-GPKGDPGISGTPGAPGLPGPKGSVGGMGLPGTPGEKGVPGIPGPQGSPGLPGDKGAKGEKGQAGPPGIGIPGLRGEKGDQGIAGFPGSPGEKGEKGSIGIPGMPGSPGLKGSPGSVGYPGSPGLPGEKGDKGLPGLDGIPGVKGEAGLPGTPGPTGPAGQKGEPGSDGIPGSAGEKGEPGLPGRGFPGFPGAKGDKGSKGEVGFP-GLAGSPGIPGSKGEQGFMGPPGPQGQPGLPGSPGHA-TEGPKGDRGPQGQPGLPGLPGPMGPPGLPGIDGVKGDKGNPGWPGAPGVPGPKGDPGFQGMPGIGGSPGITGSKGDMGPPGVPGFQGPKGLPGLQGIKGDQGDQGVPGAKGLPGPPGPPGPYDIIKGEPGLPGPEGPPGLKGLQGLPGPKGQQGVTGLVGIPGPPGIPGFDGAPGQKGEMGPAGPTGPRGFPGPPGPDGLPGSMGPPGTPSVDHGFLVTRHSQTIDDPQCPSGTKILYHGYSLLYVQGNERAHGQDLGTAGSCLRKFSTMPFLFCNINNVCNFASRNDYSYWLSTPEPMPMSMAPITGENIRPFISRCAVCEAPAMVMAVHSQTIQIPPCPSGWSSLWIGYSFVMHTSAGAEGSGQALASPGSCLEEFRSAPFIECHG-RGTCNYYANAYSFWLATIERSEMFKKPTPSTLKAGELRTHVSRCQVCMRRT
>chain B CO4A1_HUMAN/1-1669 Collagen alpha-1(IV) chain
---------MGPRLSVWLLLLPAALLLHEEHSRAAA--KGGCAGSGC-GKCDCHGVKGQKGERGLPGLQGVIGFPGMQGPEGPQGPPGQKGDTGEPGLPGTKGTRGPPGASGYPGNPGLPGIPGQDGPPGPPGIPGCNGTKGERGPLGPPGLPGFAGNPGPPGLPGMKGDPGEILGHVPGMLLKGERGFPGIPGTPGPPGLPGLQGPVGPPGFTGPPGPPGPPGPPGEKGQMGLSFQGPKGDKGDQGVSGPPGVPGQA-------QVQEKGDFATKGEKGQKGEPGFQGMPGVGEKGEPGKPGPRGKPGKDGDKGEKGSPGFPGEPGYPGLIGRQGPQGEKGEAGPPGPPGIVIGTGPLGEKGERGYPGTPGPRGEPGPKGFPGLPGQPGPPGLPVPGQAGAPGFPGERGEKGDRGFPGTS-LP-GPSGRDGLPGPPGSPGPPGQPGYTNGIVECQPGPPGDQGPPGIPGQPGFIGEIGEKGQKGESCLICDIDGYRGPPGPQGPPGEIGFPGQPGAKGDRGLPGRDGVAGVPGPQGTPGLIGQPGAKGEPGEFYFDLRLKGDKGDPGFPGQPGMPGRAGSPGRDGHPGLPGPKGSPGSVGLKGERGPPGGVGFPGSRGDTGPPGPPGY---GPAGPIGDKGQAGFPGGPGSPGLPGPKGEPGKIVP--------------------LPGPPGAEGLPGSPGFPGPQGDRGFPGTPGRPGLPGEKGAVGQPGI-GFPGPPGPKGVDGLPGDMGPPGTPGRPGFNGLPGNPGVQGQKGEP---GVGLPGLKGLPGLPGIPGTPGEKGSIGVPGVPGEHGAIGPPGLQGIRGEPGPPGLPGSVGSPGVPGI-GPPGARGPPGGQGPPGLSGPPGIKGEKGFPGFPGLD-MPGPKGDKGAQGLPGITGQSGLPGLPGQQGAPGIPGFPGSKGEMGVMGTPGQPGSPGPVGAPGLPGEKGDHGFPGSSGPRGDPGLKGDKGDVGLPGKPGSMDKVDMGSMKGQKGDQGEKGQIGPIGEKGSRGDPGTPGVPGKDGQAGQPGQP-GPKGDPGISGTPGAPGLPGPKGSVGGMGLPGTPGEKGVPGIPGPQGSPGLPGDKGAKGEKGQAGPPGIGIPGLRGEKGDQGIAGFPGSPGEKGEKGSIGIPGMPGSPGLKGSPGSVGYPGSPGLPGEKGDKGLPGLDGIPGVKGEAGLPGTPGPTGPAGQKGEPGSDGIPGSAGEKGEPGLPGRGFPGFPGAKGDKGSKGEVGFP-GLAGSPGIPGSKGEQGFMGPPGPQGQPGLPGSPGHA-TEGPKGDRGPQGQPGLPGLPGPMGPPGLPGIDGVKGDKGNPGWPGAPGVPGPKGDPGFQGMPGIGGSPGITGSKGDMGPPGVPGFQGPKGLPGLQGIKGDQGDQGVPGAKGLPGPPGPPGPYDIIKGEPGLPGPEGPPGLKGLQGLPGPKGQQGVTGLVGIPGPPGIPGFDGAPGQKGEMGPAGPTGPRGFPGPPGPDGLPGSMGPPGTPSVDHGFLVTRHSQTIDDPQCPSGTKILYHGYSLLYVQGNERAHGQDLGTAGSCLRKFSTMPFLFCNINNVCNFASRNDYSYWLSTPEPMPMSMAPITGENIRPFISRCAVCEAPAMVMAVHSQTIQIPPCPSGWSSLWIGYSFVMHTSAGAEGSGQALASPGSCLEEFRSAPFIECHG-RGTCNYYANAYSFWLATIERSEMFKKPTPSTLKAGELRTHVSRCQVCMRRT
>chain C CO4A2_HUMAN/1-1712 Collagen alpha-2(IV) chain
MGRDQRAVAGPALRRWLLLGTVTVGFLAQSVLAGVKKFDVPCGGRDCSGGCQCYPEKGGRGQPGPVGPQGYNGPPGLQGFPGLQGRKGDKGERGAPGVTGPKGDVGARGVSGFPGADGIPGHPGQGGPRGRPGYDGCNGTQGDSGPQGPPGSEGFTGPPGPQGPKGQKGEP-YALPKEERDRYRGEPGEPGLVGFQGPPGRPGHVGQMGPVGAPGRPGPPGPPGPKGQQGNRGLGFYGVKGEKGDVGQPGPNGIPSDTLHPIIAPTGVTFHPDQYKGEKGSEGEPGIRGISLKGEEGIMGFPGLRGYPGLSGEKGSPGQKGSRGLDGYQGPDGPRGPKGEAGDPGPPGLP--AYSPHPSLAKGARGDPGFPGAQGEPGSQGEPGDPGLPGPPGLSIGDGDQRRGLPGEMGPKGFIGDPGIPALYGGPPGPDGKRGPPGPPGLPGPPGPDGFL-FGLKGAKGRAGFPGLPGSPGARGPKGWKGDAGECRCTEGDEAIKGLPGLPGPKGFAGINGEPGRKGDRGDPGQHGLPGFPGLKGVPGNIGAPGPKGAKGDS-RTITTKGERGQPGVPGVPGMKGDDGSPGRDGLDGFPGLPGPPGD-GIKGPPGDPGYPGIPGTKGTPGEMGPPGLGLPGLKGQRGFPGDAGLPGPPGFLGPPGPAGTPGQIDCDTDVKRAVGGDRQEAIQPGCIGGPKGLPGLPGPPGPTGAKGLRGIPGFAGADGGPGPRGLPGDAGREGFPGPPGFIGPRGSKGAVGLPGPDGSPGPIGLPGPDGPPGERGLPGEVLGAQPGPRGDAGVPGQPGLKGLPGDRGPPGFRGSQGMPGMPGLKGQPGLPGPSGQPGLYGPPGLHGFPGAPGQEGPLGLPGIPGREGLPGDRGDPGDTGAPGPVGMKGLSGDRGDAGFTGEQGHPGSPGFKGIDGMPGTPGLKGDRGSPGMDGFQGMPGLKGRPGFPGSKGEAGFFGIPGLKGLAGEPGFKGSRGDPGPPGPP-PVILPGMKDIKGEKGDEGPMGLKGYLGAKGIQGMPGIPGLSGIPGLPGRPGHIKGVKGDIGVPGIPGLPGFPGVAGPPGITGFPGFIGSRGDKGAPGRAGLYGEIGATGDFGDIGDT-INLPGRPGLKGERGTTGIPGLKGFFGEKGTEGDIGFPGITGVTGVQGPPGLKGQTGFPGLTGPPGSQGELGRIGLPGGKGDDGWPGAPGLPGFPGLRGIRGLHGLPGTKGFPGSPGSDIHGDPGFPGPPGERGDPGEANTLPGPVGVPGQKGDQGAPGERGPPGSPGLQGFPGITPPSNISGAPGDKGAPGIFGLKGYRGPPGPPGSAALPGSKGDTGNPGAPGTPGTKGWAGDSGPQGRPGVFGLPGEKGPRGEQGFMGNTGPTGAVGDRGPKGPKGDPGFPGAPGTVGAPGIAGIPQKIAVQPGTVGPQGRRGPPGAPGEMGPQGPPGEPGFRGAPGKAGPQGRGGVSAVPGFRGDEGPIGHQGPIGQEGAPGRPGSPGLPGMPGR-SVSIGYLLVKHSQTDQEPMCPVGMNKLWSGYSLLYFEGQEKAHNQDLGLAGSCLARFSTMPFLYCNPGDVCYYASRNDKSYWLSTTAPLP--MMPVAEDEIKPYISRCSVCEAPAIAIAVHSQDVSIPHCPAGWRSLWIGYSFLMHTAAGDEGGGQSLVSPGSCLEDFRATPFIECNGGRGTCHYYANKYSFWLTTIPEQSFQGSPSADTLKAGLIRTHISRCQVCMKNL

```

# Answer
> 1 votes
Here's a solution based on Jesse's comment. The steps are:

1. read in MSA with `Biostrings::readAAStringSet()`
2. right-fill sequences so they are all the same length (assumes left alignment)
3. search for motif substrings in each amino acid sequence with `gregexpr()`. If a motif in the motif set is found, record position of center of motif.
4. bind the results from each sequence to a single table

```
find_motifs <- function(msa, motifs) {

  require(tibble)
  require(Biostrings)

  # read in fasta as Biostrings object
  seqs <- Biostrings::readAAStringSet(filepath = msa, format = "fasta")
  maxseqlen <- max(width(seqs))

  # initiate results tibble
  results <- tibble::tibble(position = seq(1:maxseqlen))

  # readAAMultipleAlignment requires equal length of sequences
  # instead, use readAAStringSet and fill "-" on right
  # assumes left-aligned MSA
  for (i in 1:length(seqs)) {

    seqlen <- nchar(seqs[i])

    if (seqlen < maxseqlen) {
      fill <- maxseqlen - seqlen
      seqs[[i]] <- paste0(as.character(seqs[[i]]),
                          paste(rep("-", fill), collapse = ""))
    }

    # initiate vector of zeros to store hit positions
    positions <- rep(0, maxseqlen)

    lapply(motifs, function(motif) {

      # factor to add to gregexpr output to center hit
      # for motifs with an even number of characters, use left center
      shift <- ceiling(nchar(motif)/2)

      hits <- unlist(gregexpr(motif, seqs[[i]])) + shift

      lapply(hits, function(hit) {
        positions[hit] <<- 1
      }) # end apply over hits

    }) # end apply over motifs

    results <- dplyr::bind_cols(results, positions)

  } # end loop over sequences

  names(results) <- c("position", paste0("seq", seq(1:length(seqs))))
  results

} # end function

```

Running this function on the provided fastas, we get a tables where positions in sequences matching motifs are recorded as `1`.

```
setwd("/path/to/files")
ngp <- c("PGP","RGP","KGP")
msa1 <- find_motifs(msa = "msa1.fasta",
                    motifs = ngp)
msa2 <- find_motifs(msa = "msa2.fasta",
                    motifs = ngp)

dplyr::slice(msa1, 120:130)

# A tibble: 11 × 4
   position  seq1  seq2  seq3
      <int> <dbl> <dbl> <dbl>
 1      120     0     0     0
 2      121     0     0     0
 3      122     1     1     1
 4      123     0     0     0
 5      124     0     0     0
 6      125     1     1     0
 7      126     0     0     0
 8      127     0     0     0
 9      128     0     0     1
10      129     0     0     0
11      130     0     0     0

```

And here's a function to make the plots.

```
heatmap_separate <- function(data) {
  require(ggplot2)
  require(tidyr)
  numseqs <- ncol(data) - 1
  numrows <- nrow(data)
  toplot <- data |> tidyr::pivot_longer(cols = -position,
                                        names_to = "sequence",
                                        values_to = "hits")
    ggplot(toplot) +
    geom_tile(mapping = aes(x = position,
                            y = sequence,
                            fill = as.factor(hits),
                            color = as.factor(hits))) +
    scale_x_continuous(limits = c(1,numrows),
                       n.breaks = 10,
                       expand = c(0,0)) +
    scale_y_discrete(expand = c(0,0),
                     limits = rev) +
    scale_fill_manual(values = c("#f2f0f7","#54278f")) +
    scale_color_manual(values = c("#f2f0f7","#54278f")) +
    geom_hline(yintercept = 0.5 + 0:numseqs,
               color = "black",
               linewidth = 1.5) +
    theme_minimal() +
    theme(legend.position = "none",
          axis.title.y = element_blank(),
          panel.grid = element_blank(),
          axis.ticks.x = element_line(color = "black"),
          panel.background = element_blank())
}

heatmap_separate(msa1)
heatmap_separate(msa2)

```

#### Plot for fasta #1

#### Plot for fasta #2

The lines representing hits are all plotted as the same color (`#54278f`), but they appear to be different colors as an artifact of plotting lots of narrow tiles. This can be rectified by viewing the plots as svg files, though I can't display that file type here.

---
Tags: r, sequence-alignment, heatmap, motifs
---


title: DNA, molecular descriptors

text:
# Question
Title: DNA, molecular descriptors
Are there any programs (preferably command line tools) for calculating molecular descriptors of DNA? I am looking for something like Chemopy or RDKit but specifically for DNA.

Thanks in advance!

# Answer
Just found the right solution. iFeatureOmega does exactly what I was looking for. Sequence based and has CLI. https://ifeatureomega.erc.monash.edu/iFeature2/index.html Additionally it supports protein, DNA/RNA and also small molecule inputs.
> 1 votes

---
Tags: phylogenetics, phylogeny, rdkit, matrix
---


title: FlowSOM multi-step clustering

text:
# Question
Title: FlowSOM multi-step clustering
There are a few different commonly used clustering algorithms within the single-cell space, although Leiden seems to be the top choice these days. FlowSOM is a classic package for analyzing flow cytometry data. It has a two-step approach for clustering. First, it builds a self-organizing map (SOM) where cells are assigned to 100 grid points. Then, it uses ConsensusClusterPlus to meta-cluster the SOM codes. I have not seen other packages that take this approach. There is no publication for FlowSOM and the documentation does not explain this part. Why this two-step approach and why isn't this more common?

# Answer
The two step approach will be used in any cluster algorithm (except see "network theory" below), you just may not see it. Clustering methods are very common in genomics, they are popular because they are very visual and will stack against a heatmap nicely. They have been replaced by dimensionality reduction such as PCA type analysis. For data mining PCA is a much more robust method for a number of reasons for example they overcome the problems listed as "The Bad" and "The Ugly" below. PCA type analysis is difficult to visualise, particularly when it steps into multi-dimensions and thats why developers like clustering ... easy for the developer to understand and easy for the use to visually appreciate, particularly against a heat map.

**Any classical clustering** The first step is to build a matrix and the second is to reorganise the relationships within the matrix using a pre-defined algorithm, unweighted pairwise group method with arithmetic mean (UPGMA).

1. The matrix in this case will be the grid.
2. The second step is the clustering method. The method used ConsensusClusterPlus by Wilkerson and Hayes (2010) is, 0

> The consensus clustering (CC) method provides quantitative and visual stability evidence for estimating the number of unsupervised classes in a dataset.

All that means is it stacks a heat map onto a clustering tree. Here's the output of the paper and the 'workhorse' is the tree presented in "A",

Like all great movies, it has 3 faces,

**The Good** using `k-means clustering` and that is well respected method in supervised-learning and suspect it performs pretty well in unsupervised learning too. It's a good method.

**The bad** There could easily be artefacts in the deeper relationships of the tree and that would need a full empirical investigation. I know that because genetic diversity is often measured by trees and loads of attention has gone into modelling the process and CCP is assuming homogeneity in the rate process - which rarely hplds in biology. Thus, there will be situations where that falls over, however, what I strongly suspect is that the k-means methods is attempting to circumvent this situation, without doing so explicitly (if that makes any sense). I would need to study the algorithm to assess this and I would not use clustering in data mining/unsupervised learning.

**The Ugly** The basic problem with clustering and the reason it isn't used in genetic diversity *****, is because there's loads of clustering algorithms, they give different answers and there's no way to distinguish which is which. Having said all that `k-means clustering` is very respectable.

The reason it's not common is that there has never been a consensus on the cluster algorithm and modellers prefer to avoid them because they are only used for data exploration/mining and not useful for hypothesis testing. *However* I do stress that `k-means` supervised learning is a very good method, but that is not the same as clustering (which is unsupervised learning).

***** or at least if it is used its just a method to visual a matrix and a like *Warning clustering!* will appear somewhere in the legend.

---

**Parameters and conclusion** There is a parameter issue with the Figure above is tracking manually. It demonstrates *k* has a *HUGE* influence on the result. In real world supervised learning parameterisation is automated, never ever would it be done manually. The technique is called parameterisation and it's kinda important (so like the user can't use a parameter that gives an output that suits them!! \[I've not said that BTW\]). So you'd ask an algorithm above the *k*=*x* algorithm to find the optimal value of *k*.

*If* `FlowSOM` automates the value for *k* and then tells you ... thats good and I'd be happy with that. *If* it doesn't, you have to manually feed loads of values of *k* and provide an independent assessment on which was the best. That is not cool and you will have to perform a laborious manual calculation.

---

**Network theory** Leiden looks to be part of network theory e.g. Trag and van Eck (2019). That is a significant advantage over classical clustering, because network theory is being dumped over the top, i.e. the third step, which can yield real power over the data set. Personally I follow graph theory, but both have major potential and are heavily under-utilised. Its not really fair to compare a classical method with the power of network theory. Certainly when entering graph theory what you can leverage from the graph (network) is enormous. I personally, thought the whole area was being ignored (there's lots of applications), but maybe its time will yet come into the mainstream.

---

To @Pallie's comments of

**dimensionality and distance**. I do not disagree at all and is the reason why dimensionality reduction approaches such as PCA or multi-dimensional scaling are preferred.

**RAM overhead** Of course the calculation will be hard on RAM: if *k* was fixed at 2 or 3, it will run reasonably quickly, but might not be meaning. Any *k* parameterisation will be RAM intensive, but IMO thats the nature of big data calculation and cloud computing: you just pick a computer with loads of RAM or just fit more RAM. There will also be a risk of course of RAM bottlenecks on standard machines (where it just freezes). The key thing is whether the approaches are giving biologically significant different results and personally I would consider it empirically. These are different calculations and this is complex data.
> 1 votes

---
Tags: single-cell, clustering, flow-cytometry
---


title: Identify SNP from genetic position

text:
# Question
Title: Identify SNP from genetic position
I am looking for a way to obtain SNP names from their coordinates (Chromosome:BP) - very much the reverse of the question answered here:

(https://stackoverflow.com/questions/20251612/map-snp-ids-to-genome-coordinates)

*Example:* I have identified that a SNP at position 64377836 on Chromosome 17 is associated with a certain factor in my analysis and would like to find out whether it has a SNP ID and what this ID is.

Manually, I can search databases for Chr17:64377836 and find rs668 - but I would like to find a way to do this automatically, preferably in R.

Is this possible in R - and if so, how?

# Answer
You can do this in R using the library biomaRt from Bioconductor. Here an example starting from a table with 3 columns: Chr, Start, End. Positions are assumed to be in GRCh38.

```
library(biomaRt)
snpMart = useEnsembl(biomart = "snps",
             dataset = "hsapiens_snp")

# assuming your table is tab-separated, with column names
myData <- read.table("myTable.txt", sep="\t", header=True)

coords <- apply(myData, 1, paste, collapse = ":")

# submit the query
getBM(attributes = c('refsnp_id', 'chr_name', 'chrom_start', 'chrom_end', 'allele'),
  filters = c('chromosomal_region'),
  values = coords,
  mart = snpMart)

```

Source
> 2 votes

---
Tags: r, gwas
---


title: Manually order the set fails in UpsetR

text:
# Question
Title: Manually order the set fails in UpsetR
*This question was also asked on Biostars*

I have tried to use UpsetR to visualise the input file, which can be found here. How is it possible to make UpsetR accept an order defined in `sets=c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha")`?

**with keep.order = T**

**with keep.order = F**

```
> library("UpSetR")
> orthogroups_df<- read.table("orthogroups.GeneCount.tsv",  header=T, stringsAsFactors = F)
> #All species
> selected_species <- colnames(orthogroups_df)[2:(ncol(orthogroups_df) -1)]
> selected_species
 [1] "Atha" "Cann" "NQLD" "Natt" "Ngla" "Nlab" "Nsyl" "Ntab" "Ntom" "Slyc" "Stub" "Vvin"
> head(orthogroups_df)
  Orthogroup Atha Cann NQLD Natt Ngla Nlab Nsyl Ntab Ntom Slyc Stub Vvin Total
1  OG0000000    0    0  965    0    0    3    0    0    0    0    0    0   968
2  OG0000001    0    1    3    0    0  448    0    0    0    0    0    0   452
3  OG0000002    0    1  313    0    0  120    1    0    1    0    0    0   436
4  OG0000003    0   93   15   21   46   16   33   63   36   25   39   26   413
5  OG0000004    1   42    2   34  109    6    8  154   11    9    4    0   380
6  OG0000005    0    2   61    1   34   44   91   70   43   20    1    0   367
> ncol(orthogroups_df)
[1] 14
> orthogroups_df[orthogroups_df > 0] <- 1
> # we only show intersections of interest ,
> intersections=list(list(selected_species),
+                    list("NQLD", "Ngla", "Natt", "Nlab", "Nsyl", "Ntab", "Ntom"),
+                    list("Stub", "Slyc"),
+                    list("Atha", "Vvin"),
+                    list("Ntab", "Nsyl", "Ntom"),
+                    list("Nlab", "NQLD", "Ngla"),
+                    list("Nlab", "NQLD", "Nsyl"),
+                    list("Nlab", "Ngla", "Nsyl"),
+                    list("NQLD", "Nsyl", "Ngla"))
> upset(orthogroups_df,
+       text.scale = c(1.4),
+       sets=c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha"),
+       order.by = "freq",
+       keep.order=T,
+       intersections = intersections,
+       sets.x.label="Total number of orthogroups",
+       mainbar.y.label = "Number of orthogroups")

```

# Answer
> 2 votes
It's possible that `upset` is using factor levels for ordering, which is the default approach for ggplot2. If this is the case, then it can be fixed by converting the variable to a factor, and using an explicit level definition:

```
sets <- factor(sets, levels=sets)

```

For your specific example:

```
upSets <- c("Nlab", "NQLD", "Ngla", "Nsyl", "Ntom", "Ntab", "Natt", "Cann", "Stub", "Slyc", "Vvin", "Atha")
upSets <- factor(upSets, levels=upSets)
upset(orthogroups_df,
   text.scale = c(1.4),
   sets=upSets,
   order.by = "freq",
   keep.order=T,
   intersections = intersections,
   sets.x.label="Total number of orthogroups",
   mainbar.y.label = "Number of orthogroups")

```

---
Tags: r, ggplot2, upsetr
---


# Question
Title: CRISPR genome wide screen analysis using MAGeCK-- Loss of reads/sgRNA using MAGeCK count function?
Sorry if this question seems basic. I am trying to run the MAGeCK pipeline to analyze CRISPR knockout screen data produced by someone in my lab around 5 years ago. I was given the data as bam files and also have the sequencing sample statistics from when the samples were run including the number of total reads, aligned reads, and unaligned reads. I converted the bam files back to fastq files using

```
samtools view -h -F 2048 filename.bam > tmp.bam

bedtools bamtofastq -i tmp.bam -fq filename.fastq

```

After converting the bam files back to fastq files, I ran the MAGeCK count function. I saw that none of the total reads matched between the MAGeCK count summary output file and the original sequencing statistics.

More problematically, a few samples had significantly fewer total reads (like 50,000) shown on the MAGeCK count summary output files relative to the original sequencing statistics (~44 million). As a result, many of the sgRNAs are not represented at all in the sample (there are ~61,000 zero counts relative to the total 77441 sgRNAs) which is impacting my analysis.

Can anyone help me understand why I might be losing sgRNAs using MAGeCK? Is it a problem with my source files? I would appreciate any help!

# Answer
It's hard to say without looking more closely at the data, but you may be able to find a few guides that are drastically different in your counting, and see how prevalent they are within the original BAM file (before filtering). Perhaps the `-F 2048` being used to filter out chimeric reads is overfiltering, and causing the discrepancy? Are unaligned reads included in the BAM - if not, does this explain the discrepancy?

Since the total reads doesn't match the original, I'd look into this. Are the extracted FASTQs consistent read length? Maybe some trimming was performed that removed the guides?

Do you know what reference the BAM files were mapped to? Perhaps this reference did not include the guide sequences, and were mapped to the genome instead? You may be able to extract this information from the BAM like this:

```
 samtools view -H filename.bam > filename.header

```

Ideally you'd have the original FASTQs to work with, but it sounds like those are unavailable...
> 1 votes

---
Tags: samtools, fastq, crispr
---


# Question
Title: scRNA-Seq pseudobulk for marker detection
A standard approach for scRNA-Seq is to partition the single cells into individual clusters, then use a Wilcoxon test to find markers that characterize each cluster (or other statistical methods that consider single cells as replicates).

Recently, there has been a push towards using pseudobulk approaches for "differential state" (DS) analysis, which avoid the pseudoreplication bias of considering single cells as replicates. See e.g. Crowell (2020), Zimmerman (2021), Thurman (2021), Squair (2021), or the OSCA book.

However, all the above papers seem to use pseudobulk methods for DS of a single cluster between two conditions (e.g. drug treatment), or occasionally to compare a specific pair of clusters. I haven't seen any example of using pseudobulk to find markers distinguishing clusters (comparisons of one cluster against all others).

I can think of a few reasons, for example:

* these pseudobulk methods can be more computationally intensive,
* perhaps the Wilcoxon test is good enough to simply find markers,
* pseudobulk methods can only work when each cluster is represented in several samples, which is hard to guarantee for every single cluster in the dataset.

On the other hand it seems to me a pseudobulk approach should yield more relevant markers (genes that are robustly enriched across biological replicates). So, am I missing something? Is there a statistical reason not to use pseudobulk approaches for finding cluster-specific markers?

# Answer
The only reason I see for not using pseudobulks is lack of biological replication. Even though single-cell sequencing is common these days and standardized it is still expensive. If you have several conditions plus biological replicates and need many cells, maybe due to heterogeneity, the bill can **easily** reach the 5-digit realm of money. Biological replication is always desired but sometimes simply not possible. We recently did a study in which we were interested in organ heterogeneity in different conditions and both mouse work and library prep was so extensive that even hash tagging the mice we had to pool was simply not possible. So we were bound to using each cell as a replicate. I prefer limma-trend for this, not the Wilcox test. It is not computationally very demanding, neither is pseudobulks. But even if it was, I would happily wait 30min rather than 10sec to get the best possible rather than the fastest results. Pseudobulks conceptually have the advantage that you basically are forced to have biological replication (because only then they make sense) so findings are intrinsically more robust. Theoretically, you can of course have biological replicates and still treat each cell as replicate, and then blocking or adding the biological replication information into the DE testing design, but the statistical rigor of methods used for pseudobulks (DESeq2, edgeR and other) is for me the reason to always use pseudobulks if I can.

I personally use pseudobulks also for markers. After all marker detection comes downstream of DE testing by aggregating the results into lists that contain unique or enriched genes per cluster. I see no reason to use multiple testing strategies.
> 1 votes

# Answer
If you want a statistical reason not to use pseudobulk to do DE for marker finding, I think it is the same reason not to use any DE method for marker finding: the **double dipping problem** (see Gao and Witten 2020).

Briefly, if you cluster a dataset based on genes (we can safely say that clustering in PCA space is *still* based on genes), then asking which genes are differentially expressed between clusters is a statistical fallacy. We could just as easily divide apples in 2 groups according to their shade of red and then ask which RGB channels exhibit a statistically significant difference.

This means that even in the absence of real structure, you will identify very low p-values for your DE test. The OSCA book section you mention links to another section in the advanced book that tells you the same thing with a small practical demonstration generating a 200-dimensional Gaussian and clustering it, then testing for differences in clusters.

This does *not* mean, however, that markers found this way are meaningless; it simply means that you cannot use a p-value cutoff to define whether a gene is a marker in this scenario, or say that since you have statistically significant markers two clusters are different.

You can still rank these markers by fold change, difference in percentage expression, AUC, entropy, etc *as well as p-value*. Moreover, markers identified using the traditional methods tend to make biological sense, and statistically rigorous methods are still in development (there was a beautiful talk by Anna Neufeld on an approach involving count splitting, you can read the preprint here).

The test you suggest is interesting, but I don't think you would solve the double dipping problem; this being said it looks more robust and appropriate than using every single cell as a replicate.
> 1 votes

---
Tags: scrnaseq, differential-expression
---


# Question
Title: How to display novel genome assemblies or uncommon genome assemblies using the UCSC Genome Browser?
I want to display E.coli BW25113 (GenBank: CP009273.1) strain in UCSC browser. This strain is not listed in http://microbes.ucsc.edu/ browser. How can I display E.coli BW25113 assembly in the browser?

# Answer
> 4 votes
I found two methods:

1. You want to visit this page for instructions for novel assemblies: http://genomewiki.ucsc.edu/index.php/Assembly\_Hubs .
2. However, before doing that I highly recommend that you check Genome Assembly Hub first: https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/ . There is a good chance that your genome or strain is listed there. Here is an example of how you can do that. Example specifically shows how you can display E.coli BW25113 strain in the browser.

---

**Instructions:**

1. Go to NCBI Nucleotide: https://www.ncbi.nlm.nih.gov/nucleotide .
2. Enter your strain name: `E.coli BW25113 (GenBank: CP009273.1)`
3. Find Biosample ID on the page, located under DBLINK: `BioSample: SAMN03013572`
4. Go to GenBank assembly hubs: https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/ .
5. Select bacteria. This will take you to "Bacterial assembly hub" at https://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/bacteria/bacteria.html .
6. Look for *E. coli* hubs. At the time I looked #15-18 contain *E. coli* strains.
7. Click on each hub and search for `SAMN03013572` using browser search function (`Command`+`F` on macOS). I found my genome in hub #16.
8. Click on "common name" in the column-3. This will take you to https://genome-test.gi.ucsc.edu/cgi-bin/hgGateway?hubUrl=http://genome-test.gi.ucsc.edu/gbdb/hubs/genbank/bacteria/16/hub.txt&genome=GCA\_000750555.1\_ASM75055v1&position=lastDbPos , where `FindPosition` window will be populated with your genome Assembly.
9. Click on `Go`. This should take you to the UCSC Genome Browser with your genome assembly.

# Answer
> 1 votes
I am not sure if this qualifies as an official answer. I have contacted UCSC bioinformatics team. They were were very kind and told me that they would add any track that was requested by users. After my request, they kindly added the track for E.coli BW25113. If you want your reference tracks to be added you could submit a request at this website: Geneme assembly search and request page.

Then, at the upper right corner click on "Can't find your assembly?" button and submit your request. See the screenshot below from the webpage.

---
Tags: assembly, ucsc, bacteria, genome-browser
---


title: Reading a Fast5-file with Python

text:
# Question
Title: Reading a Fast5-file with Python
I am trying to extract data from fast5-file with python 3.9.13 in Ubuntu. I have found a library "fast5\_research"(This package comprises an API to HDF containers used by the research groups within Oxford Nanopore Technologies), installed on my PC and updated my software with "sudo apt update". According to tutorial I attempted to read fasta5-file with the next code:

```
from fast5_research import Fast5
filename='/home/master/Desktop/0.fast5'
with Fast5(filename) as fh:
    raw = fh.get_read(raw=True)
    summary = fh.summary()
print('Raw is {} samples long.'.format(len(raw)))
print('Summary {}.'.format(summary))

```

but as a result, I have an error:

```
KeyError                                  Traceback (most recent call last)
/tmp/ipykernel_410530/3650860803.py in <module>
      4 from fast5_research import Fast5
      5 filename='/home/master/Desktop/0.fast5'
----> 6 with Fast5(filename) as fh:
      7     raw = fh.get_read(raw=True)
      8     summary = fh.summary()

~/anaconda3/lib/python3.9/site-packages/fast5_research/fast5.py in __init__(self, fname, read)
     85
     86         # Attach channel_meta as attributes, slightly redundant
---> 87         for k, v in _clean_attrs(self[self.__channel_meta_path__].attrs).items():
     88             setattr(self, k, v)
     89         # Backward compat.

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

~/anaconda3/lib/python3.9/site-packages/h5py/_hl/group.py in __getitem__(self, name)
    175                 raise ValueError("Invalid HDF5 object reference")
    176         else:
--> 177             oid = h5o.open(self.id, self._e(name), lapl=self._lapl)
    178
    179         otype = h5i.get_type(oid)

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/_objects.pyx in h5py._objects.with_phil.wrapper()

h5py/h5o.pyx in h5py.h5o.open()

KeyError: 'Unable to open object (component not found)'

```

How to fix this error?

# Answer
The code looks fine `fast5` is simply `HDF5` format (rapid I/O).

The error,

> 'Unable to open object (component not found)'

You are not the only person to have this error and it could be a bug in the code with `hdf5` in Python. Just to note I use `hdf5` without a any problem.

1. Problem with the installation (?), this is the best guess to resolve it but I notice from your error log .. the errors are coming from Anaconda, I assume you made a `conda` installation, but the post implies `apt-get`?? Just check where the package is located.

```
conda install -c hcc fast5-research

```

Then perform a conda environment creation and activation of the same package as here. I think you've already done that. `anaconda3` is definitely your Python3 version ... and in that case it might be better to do the full `conda` installation because it stops things getting mixed up.

**NOTE** There is a Docker container

```
docker pull quay.io/biocontainers/fast5-research

```

I would recommend using this (it will almost certainly work), if you know how to leverage Docker (which is sort of complicated).

2. File corruption ... examine the file via the HDFView GUI, here https://www.hdfgroup.org/downloads/hdfview/ . If you can't see within the file via the viewer - I think the file has a problem. Perhaps examine via hdf5 ... wait no that gets complicated.
3. It can't find the file or a path error (below).

Then kindly try:

```
from fast5_research import Fast5
import pathlib as Path
import sys

filename=Path('/home/master/Desktop/0.fast5')
if filename.exists():
    print ('Okay {} is fine.'.format(filename))
    filename='/home/master/Desktop/0.fast5'
else:
    pathway = Path(__file__).absolute().parent
    if not Path(pathway, '0.fast5').is_file():
        raise ValueError ('Sorry filename has not been found')
        sys.exit()
    filename = str(Path(pathway, '0.fast5'))
try:
    if Path(filename).exists():
        with Fast5(filename) as fh:
            raw = fh.get_read(raw=True)
            summary = fh.summary()
        print('Raw is {} samples long.'.format(len(raw)))
        print('Summary {}.'.format(summary))
except Exception as err:
    print ("The error codes are {} and {}".format(err.message, err.args))

```

Please execute the code in the same `dir` as `'/home/master/Desktop/'` or maybe just `'/home/Desktop`

In `Fasta5` will almost certainly take a `path` from `pathlib` but I've erred in case not. `HDF5` will certainly accept a `pathlib` pathway.

I've not tested it but just post back what happens.
> 0 votes

---
Tags: python, nanopore, jupyter, base-calling, fast5
---


title: How to run how_are_we_stranded_here?

text:
# Question
Title: How to run how_are_we_stranded_here?
I'm trying to run the python library called how\_are\_we\_stranded\_here. I have installed with `pip install how_are_we_stranded_here`. I have also installed its dependencies i.e kallisto 0.44.0 and RSeQC. On how\_are\_we\_stranded\_here's pypi page, it gives the following command under the usage heading:

```
 check_strandedness --gtf Yeast.gtf --transcripts Yeast_cdna.fasta \
--reads_1 Sample_A_1.fq.gz --reads_2 Sample_A_2.fq.gz

```

So I tried running the program in the terminal using my own files, the command was:

```
    check_strandedness --gtf Glycine_max.Glycine_max_v2.1.55.gtf \
--transcripts Glycine_max.Glycine_max_v2.1.cdna.all.fa \
--reads_1 SRR2079644.fastq

```

But it keeps giving the error

> ./how\_stranded.sh: line 1: check\_strandedness: command not found

I will be very thankful for any help on how to run it.

# Answer
`pip` should just work, anyway... Two approaches, point 1 is a last resort approach (not recommended), I definitely recommend point 2.

First try,

```
python3 -m pip install --upgrade how_are_we_stranded_here

```

Try again does it work? Error logs?, then possibly ...

1. You need to find out where it is,

```
   pip show how_are_we_stranded_here
   pip show check_strandedness
   pip show how_stranded.sh

```

What you really want is the `check_strandedness`

```
cd ~
mkdir bin
./pathtoapp/check_stranded # make sure it works
ln -si /pathwaytoapp/how_stranded.sh /home/whatever/bin/
ln -si /pathwaytoapp/check_strandedness ~/bin/check_strandedness

```

Then write

```
export PATH="/home/whatever/bin:$PATH"

```

or modify `.profile` or `.bashrc`, then

```
source ~/.bash_profile

```

Try again ...

2. It is better to perform a conda installation

Install either,

```
    conda install -c bioconda how_are_we_stranded_here
    conda create --name stranded
    conda activate stranded
    check_strandedness -h

```
> 2 votes

---
Tags: rna-seq
---


# Question
Title: How do I determine read length from a user-inputted fastq file and then input that information into a tool's command line?
Can you guide me on how could take an input fastq file, read the length from read and then feed it into the findtail's command?

I have suggested concerns about documentation and program defaults to my prof and asked him to switch to cutadapt instead of using findtail that has no documentation. But he refuses so I have no option but to proceed with what he wants.

I want the user to input a RNA-seq data fastq file. Since the first line of each file contains the sequence length information, such as described:

> @SRR3003130.1 HWI-ST741:596:H5FMVBCXX:1:1101:1244:1952 length=51

How do I pick the *length=51*, the **51** from this and feed in into the following command's **--seqlength** parameter:

```
    ./findtail --input_file Rrp6_S38_all_R1_001_egal.fastq --seqlength 51 \
--endgap 5 --taillength 10 --identity 100 --ptype A --stype T \
--output_format fasta --output_type pr > pr.fasta

```

The tool's command that I'm using only requires us to enter --seqlength once. So I do not think that determining the length of each and every read in the fastq file is of any use. Just simply the length from the first line of first read would suffice.

*If something is lacking from my question, please let me know how can I improve.*

# Answer
I think that you cannot rely on that information being in a fastq header. **Note also that each read has its own header! This is not just the first line!** For example when I google "sample fastq file" the first examples do not have it. Spec doesn't seem to say anything about it.

For a one-off on a specific dataset this *might* be adequate:

```
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2
# prints "75" for each occurrence, so if you trust the first line you could pass to head:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2 | head -1
# or you could explicitly sort and grab the longest read len:
grep -Eo "length=([0-9]+)" test.fq | cut -d"=" -f2 | sort | tail -1

```

For a more general application I'd suggest using some existing/explicit tools for measuring read length like seqkit or awk, rather than trusting some information that may or may not be in the header.

I can't find any docs for `findtail`, depending on what your application is I'd suggest using cutadapt instead, which does not require a single sequence length input.
> 3 votes

# Answer
To get the length of the first sequence from a single-line-per-sequence fastq, the second line can be piped through `wc -c`, and stored in a shell variable:

```
seqLength=$( head -n 2 input_file_name.fastq | tail -n 1 | wc -c)

```

This variable can then be used for the program arguments using `${variable}`:

```
seqLength=$(head -n 2 input_file_name.fastq | tail -n 1 | wc -c)
./findtail --input_file input_file_name.fastq --seqlength ${seqLength} \
    --endgap 5 --taillength 10 --identity 100 --ptype A --stype T \
    --output_format fasta --output_type pr > pr.fasta

```

If the input file is gzipped \[I have no idea if findtail supports gzipped files\], then the length can be determined by piping through zcat:

```
seqLength=$( zcat input_file_name.fastq.gz | head -n 2 | tail -n 1 | wc -c)

```

As mentioned by many comments, note that some FASTQ files won't have the same length for each read, and it's unclear how `findtail` will behave for such a file.
> 3 votes

---
Tags: rna-seq
---


# Question
Title: Why does this modification site file contain duplicate sites and multiple related proteins?
The data below is coming from a **MaxQuant** post-translational protein modification (PTM) site file.

Why does it contain duplicate sites (yellow) and proteins (cyan)?

> *Source data*
>
> http://ftp.ebi.ac.uk/pride-archive/2020/12/PXD020517/MaxQuant\_BA39\_SOLUBLE.zip

# Answer
There are multiple protein isoforms that the gene *may* encode. The hyphen-number after the Uniprot ID represents the isoform —in a Uniprot you find the detail of these at the bottom. The canonical isoform in Uniprot is the most likely one to be biologically relevant, while in NCBI and Ensembl it is the longest. In a limited amount of case multiple isoforms are biologically relevant. For the sake of sanity most variant annotation tools have a `canonical` argument to use the canonical.
> 2 votes

# Answer
According to a proteomics researcher at the NIH:

* `proteins` column contains uniprot IDs.
* `position within proteins` column are the starting position of each ID.
* `position` is the *leading* protein
> 1 votes

---
Tags: proteins, mass-spectrometry, maxquant
---


# Question
Title: How to remove double '/' from file path, bash script
I'm using the following script to detect strandedness of my paired end RNA-seq data.

```
#!/usr/bin/env bash

cd ~/Desktop/fastq

for infile in *_1.fastq
do
   base=$(basename ${infile} _1.fastq)

    docker run quay.io/biocontainers/guessmylt:0.2.5--py_0 GUESSmyLT \
    --reads ${infile} ${base}_2.fastq

done

```

But when I run it, it keeps giving the error:

```
Building DAG of jobs... File path //SRR3475890_1.fastq contains double
'/'. This is likely unintended. It can also lead to inconsistent
results of the file-matching approach used by Snakemake. File path
//SRR3475890_2.fastq contains double '/'. This is likely unintended.
It can also lead to inconsistent results of the file-matching approach
used by Snakemake. MissingInputException in line 46 of
/usr/local/lib/python3.6/site-packages/GUESSmyLT/subsample: Missing
input files for rule subsample_fastq: //SRR3475890_2.fastq
//SRR3475890_1.fastq

```

Why is it appending double '/' to the fastq file's path? Nowhere in my script did I mention any path.

*I apologize for the very vague title, English is not my first language, so I did not know how else to phrase it.*

# Answer
The issue here is that the current working directory inside the container is `/`. You can see this with:

```
docker run --rm quay.io/biocontainers/guessmylt:0.2.5--py_0 python \
    -c 'import os; print(os.getcwd())'

```

When GUESSmyLT tries to get the absolute path of it's arguments (`--reads`, `--reference`, `--output`, etc), it actually appends a trailing slash to the current working directory, hence the double slash in the returned file names.

To pass the files into the container and ensure that the results are available when it exits, you'll need a volume mount anyway. For example:

```
docker run --rm -v "${PWD}:/working" quay.io/biocontainers/guessmylt:0.2.5--py_0 GUESSmyLT \
    --reads "/working/${base}_1.fastq" "/working/${base}_2.fastq" \
    --reference "/working/reference.fa" \
    --annotation "/working/annotation.gtf" \
    --output "/working"

```
> 4 votes

---
Tags: rna-seq, bash
---


# Question
Title: How to extract date of acquisiton from mzXML and mzML files with Python?
I would like to extract date of acquisiton from mzXML and mzML files, using Python. Is it possible e.g. with Pyteomics? pymzml is only for mzML files. Are there other libraries that can do that? I assume of course that the date of acquisition is recored in mzXML/mzML files.

# Answer
> 3 votes
You should be able to extract anything from mzML and mzXML with Pyteomics. You will have to look at your files and figure out what you need to do, though, because you need to know what element contains the information that you want. I don't often work with mzXML and in the files that I have I do not see an acquisition date at all. With mzML it appears that the `<run>` element has this information:

```
<run id="small_raw" defaultInstrumentConfigurationRef="IC1" startTimeStamp="2005-07-20T14:44:22">

```

So I can get it like this:

```
In [1]: from pyteomics import mzml

In [2]: with mzml.MzML('tests/test.mzML') as f:
   ...:     print(next(f.iterfind('run', recursive=False))['startTimeStamp'])
   ...:
2005-07-20T14:44:22

```

P.S. Disclaimer: I am the maintainer of Pyteomics.

# Answer
> 0 votes
`pymzml` certainly works via Bald *et al* 2012:

```
conda install -c bioconda pymzml

```

Whilst this only works with `mzML` you simply convert `mzXML` to `pymzml` via `ProteoWizard`

```
msconvert your_mzXML_file --mzML

```

It is described here with your previous question with @haci

I believe you managed this via

```
conda create -n openms -c bioconda openms
FileConverter -in file.mzXML -out file.mzML

```

I suspect the underlying question is direct parsing in `mzXML` without the need to perform a conversion from `mzXML` to `mzML` or without using `Bioconductor`

---

**Possibility** The following may be of use `massspec-toolbox` at:

http://code.google.com/p/massspec-toolbox/source/browse/#svn/trunk/mzxml

Note, I've not tested nor read this but the repository looks genuine.

---
Tags: python
---


# Question
Title: HTTP status code: 400 "Value must be an integer" error uploading JSON schema results to benchling
I am writing a benchling uploader script in python to upload data from a csv file in a JSON formatted schema to the benchling API platform. My function to create the JSON schema in the script looks like the following:

```
def create_antismash_schema(dataframe):
    """Create a json object for antismash output"""
    return {
        "assayResults": [
            {
                "fields": {
                    "entity": {"value": j['SBI ID']},
                    "sequence_id": {"value": j['sequence_id']},
                    "category": {"value": j['Category']},
                    "product": {"value": j['Product']},
                    "start": {"value": j['start']},
                    "end": {"value": j['end']},
                    "nucleotide_length": {"value": j['nucl_length']}
                },
                "schemaId": "assaysch_Z6B7TlWG",
                "projectId": 'src_tIUAcxDH'
            } for i, j in dataframe.iterrows()
        ]
    }

```

Now the data I have in benchling for the fields matches the keys in the schema and the values match the headers in the csv file I'm pulling the data from.

benchling results schema looks like:

csv file:

The error I am getting when I try to upload is the following:

```

Uploading Antismash data to Benchling...
  0%|                                                                                | 0/159 [00:00<?, ?it/s]{
    "assayResults": null,
    "errors": [
        {
            "fields": {
                "nucleotide_length": "Value must be an integer"
            },
            "global": []
        },
        {
            "fields": {
                "nucleotide_length": "Value must be an integer"
            },
            "global": []
        },
        {
            "fields": {
                "nucleotide_length": "Value must be an integer"
            },
            "global": []
        }
    ]
}
  0%|                                                                                | 0/159 [00:00<?, ?it/s]
Failed to upload summary.
HTTP status code: 400

```

The field is set to be an integer in benchling and all the values I'm pulling from in the csv file for the nucl\_length column are integer values so not sure why this error is beig thrown. Any guidance here appreciated, thanks.

# Answer
> 1 votes
There is a major issue about the `json` format required, however just converting to `json` is easy because `iterrows` has to be `pandas` and `pandas` has its own output `json` format ... For example,

```
dataframe.to_json('benchling.json'', orient='records', lines=True)

```

If that fails then,

```
jsonOut = dataframe.to_json(orient='records')[1:-1].replace('},{', '} {')
with open('benchling.json', 'w') as fout:
    fout.write(jsonOut)

```

However, I don't think it would work and at a minimum I think the headers of the `csv` will need to change.

---

The error

> "nucleotide\_length": "Value must be an integer"

Looking at your code ....

```
"nucleotide_length": {"value": j['nucl_length']}

```

Sure this is right?

Why not try ...

```
"nucleotide_length": int(j['nucl_length'])

```

Thats what the error seems to be saying.

---
Tags: python, dataframe, pandas, json
---


# Question
Title: How to include only the trimmed fastq files into a nextflow channel and exclude any secondary log files created?
I have the following simple nextflow script, that takes fastq paired reads and runs them through the atria adapter trimming.

```
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"
params.transcriptome_file = "$projectDir/data/Glycine_max.cdna.all.fasta"
params.gtf_file = "$projectDir/data/Glycine_max.gtf"
params.outdir = "results"

/*
 * define the `TRIM` process that creates trimmed.reads.fastq files
 * given the reads.fastq file.
 */

process TRIM {

    tag "$pair_id"

    input:
    tuple val(pair_id), path(reads)

    output:
    path(pair_id)

    script:
    """
    atria -r ${reads[0]} -R ${reads[1]} -t 1 -l 12 -q 30 -o $pair_id
    """
}

workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_ch=TRIM(read_pairs_ch)

}

```

Next what I want to achieve is input the atria.trimmed.fastq files generated by the trim process into a channel so that I could use it as an input for other process. But atria trimming tool also creates two additional log files along with atria.trimmed.fastq. How do i exclude those additional files and only input the trimmed files into a channel for later use.

Best wishes.

# Answer
With `output: path(pair_id)`, you're telling Nextflow to expect a file or directory with the value of *pair\_id* in the process working directory when the command exits (successfully). Atria outputs all files to a directory specified using `-o PATH`. This directory is created if it doesn't already exist. Note that the output directory defaults to the current working directory if `-o` is not specified. If you'd like to write the results to a sub-directory, you could try:

```
process TRIM {

    tag { pair_id }

    publishDir "${params.outdir}/TRIM", mode: 'copy'

    cpus 1

    input:
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path(trimmed_reads), emit: trimmed_fastqs
    path "${pair_id}/*.atria.log*", emit: logs

    script:
    def (r1, r2) = reads

    trimmed_reads = reads.collect {
        "${pair_id}/${it.baseName}.atria.fastq"
    }

    """
    atria \\
        -r "${r1}" \\
        -R "${r2}" \\
        -o "${pair_id}" \\
        -t ${task.cpus} \\
        -l 12 \\
        -q 30
    """
}

```

In the above, we used the emit option in the process output definition to assign a name identifier to each declaration. We can use this feature to reference the desired channel in the workflow block. For example:

```
workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    trim_out = TRIM(read_pairs_ch)

    trim_out.trimmed_fastqs.view()
}

```

Results:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [distraught_dalembert] DSL2 - revision: aea2d37092
executor >  local (3)
[5f/78908b] process > TRIM (A) [100%] 3 of 3 ✔
[B, [/path/to/work/6a/440a50e9a0725c18f9eede80549cb8/B/B_1.atria.fastq, /path/to/work/6a/440a50e9a0725c18f9eede80549cb8/B/B_2.atria.fastq]]
[A, [/path/to/work/5f/78908b96adee8b8d663d9b4eb4faba/A/A_1.atria.fastq, /path/to/work/5f/78908b96adee8b8d663d9b4eb4faba/A/A_2.atria.fastq]]
[C, [/path/to/work/b0/a3bb3024a4c05fffec555cce5097a7/C/C_1.atria.fastq, /path/to/work/b0/a3bb3024a4c05fffec555cce5097a7/C/C_2.atria.fastq]]

```

---

To input the trimmed reads into a downstream process, just ensure that the input tuple matches the input set cardinality declared by the process. If the input declaration is not the same, you may need to apply one or more operators to your channel to make it fit. For your alignment process, you might not need to make any changes. Also, if your alignment process needs to declare one or more reference files, just make sure these are value channels. Most of the time, what you want is one queue channel and one or more value channels when you require multiple input channels, otherwise you might observe unusual behavior (e.g. see this question). The section on multiple input channels in the docs is well worth taking the time to read in my opinion.

A dummy alignment process might look like:

```
process ALIGN {

    tag { sample }

    debug true

    input:
    tuple val(sample), path(reads)
    path fasta
    path gtf

    script:
    def (r1, r2) = reads

    """
    ls -g "${r1}"
    ls -g "${r2}"
    ls -g "${fasta}"
    ls -g "${gtf}"
    """
}

```

And the updated workflow:

```
workflow {

    read_pairs_ch = Channel.fromFilePairs( params.reads, checkIfExists:true )

    transcriptome = file( params.transcriptome_file )
    gtf = file( params.gtf_file )

    TRIM( read_pairs_ch )

    ALIGN( TRIM.out.trimmed_fastqs, transcriptome, gtf )
}

```

Results:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [extravagant_lattes] DSL2 - revision: 29b5315739
executor >  local (6)
[51/e4e261] process > TRIM (A)  [100%] 3 of 3 ✔
[cc/72c53d] process > ALIGN (A) [100%] 3 of 3 ✔
lrwxrwxrwx 1 users 57 Dec 14 00:21 C_1.atria.fastq -> ../../c1/9a85e19c601f2ab2c4af9c6fd203f0/C/C_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 C_2.atria.fastq -> ../../c1/9a85e19c601f2ab2c4af9c6fd203f0/C/C_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

lrwxrwxrwx 1 users 57 Dec 14 00:21 B_1.atria.fastq -> ../../3f/16747fbd2c9bc4c0dc945b24873cca/B/B_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 B_2.atria.fastq -> ../../3f/16747fbd2c9bc4c0dc945b24873cca/B/B_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

lrwxrwxrwx 1 users 57 Dec 14 00:21 A_1.atria.fastq -> ../../51/e4e261d019c65377d63826963b19b1/A/A_1.atria.fastq
lrwxrwxrwx 1 users 57 Dec 14 00:21 A_2.atria.fastq -> ../../51/e4e261d019c65377d63826963b19b1/A/A_2.atria.fastq
lrwxrwxrwx 1 users 40 Dec 14 00:21 Glycine_max.cdna.all.fasta -> ../../../data/Glycine_max.cdna.all.fasta
lrwxrwxrwx 1 users 29 Dec 14 00:21 Glycine_max.gtf -> ../../../data/Glycine_max.gtf

```
> 2 votes

---
Tags: nextflow
---


# Question
Title: How to keep the basename of an input file and change its extension in Nextflow?
The following Nextflow scripts takes the $gtf as an input and outputs a .bed file. What I want is to keep just the basename of the gtf file and add the .bed extension to it. I tried declaring a variable x and using basename from bash but it gives error.

```
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.gtf_file = "$projectDir/data/Glycine_max.gtf"

    /*
 * define the `GTF2BED` process that creates a BED file
 * given the params.gtf_file as input.
 */

process GTF2BED {

    input:
    path gtf

    output:
    "${gtf}/.bed"

    script:
    """
    gtf2bed --gtf $gtf --bed "${gtf}/.bed"
    """
}

```

Best wishes.

# Answer
> 3 votes
You can use the `getBaseName()` method on `file` objects to get the basename:

> The invocation of any method name starting with the `get` prefix can be shortcut by omitting the `get` prefix and `()` parentheses. Therefore, writing `myFile.getName()` is exactly the same as `myFile.name` and `myFile.getBaseName()` is the same as `myFile.baseName` and so on.

For example:

```
process GTF2BED {

    input:
    path gtf

    output:
    path "${gtf.baseName}.bed"

    script:
    """
    gtf2bed --gtf "${gtf}" --bed "${gtf.baseName}.bed"
    """
}

```

If you need to remove a double extension (.e.g. `.tar.gz`) and don't want to remove all extensions (e.g. using the `getSimpleName()` method), you can use the long form to supply the number of extensions to remove. For example:

```
tarball = file( 'my.file.tar.gz' )

println( tarball.getBaseName(2) )

```

Prints:

```
my.file

```

---
Tags: nextflow
---


title: SNPs location out side the gene

text:
# Question
Title: SNPs location out side the gene
I downloaded exome data from gnomad to get the allele frequency of some SNPs that I'm analyzing. The SNPs are from external resource. Some SNPs (in my data) are not present in the exome data (in gnomad):

> "No genes were found in this region"

As part of my report question for homework I need to answer this:

> "Are the SNPs in your data that are not in gnomad are outside genes?",

What does it mean SNP outside genes? Does it mean that they are not in the exome data but in the genome data on gnomad?

I found this in an article:

> "Single-nucleotide polymorphisms (SNPs) are single base-pair substitutions that occur within and outside genes"

Also found this:

> "SNP markers might be selected outside of exomes or the direct cause of a genetic mutation" does this means that it is found in the genome but not exome data in gnomad as it is outside the gene?

Basically, I can't seem to understand what is a SNP outside a gene and can I see this information in gnomad.

# Answer
> 0 votes
The answer is you can't see a SNP outside a gene for `gnomad` because you stated it is *exome data*. Thus, a minimal definition of the gene is it comprises: introns, exons and a promoter region.

*Exome* = *exon* so you downloaded a subsection the gene across the genome.

To reiterate, you don't have any external allele frequency data outside the gene. However, if your *in-house* data set is whole genome data, i.e. for making your SNP calls then most of these SNPs will be outside the exon (thats a result of comparative selection pressure). So the answer is yes. It isn't clear what the *in-house* data set comprises however and outside an *exon* does not necessarily mean its outside the gene. A lot of mutations will occur within introns (which is outside the exon but still within a gene), but in addition more mutation will occur in intergenic regions, i.e. completely outside the exon + intron + promoter. In this sense the question seems to possibly lack clarity.

To be honest, questions on bioinformatic practicals are often after-thoughts to the core activity, so I wouldn't worry too much about the semantics.

---
Tags: gene, snp, homework, exome, gnomad
---


title: Allele frequency gnomad data

text:
# Question
Title: Allele frequency gnomad data
I noticed that some gnomad SNPs have allele frequency =0 and some are just not found : No genes were found in this region. what is the difference between them : for example:Y-10010520-10010520.

Is all frq=0 means that it is very rare and therefore zero? then what about AF=NA? Also, does the gnomad genome data contain information about the exomes also, or is the information in the exome data is not included in the genome data in gnomad?

# Answer
I don't do eukaryotes, well humans anyway, but in population genetics an allele frequency of zero simply means it is present in no individuals against the reference allele. It doesn't mean the allele is heterozygous nor homozygous its simply absent. Whether that is informative or not is an interesting question. In addition, it does not necessarily mean the allele is rare, it is simply absent against the population in your analysis. The sickle cell anaemia allele is not rare, in certain parts of Africa it has high frequency, but it would be absent if all the individuals in the population of the analysis that did not have African ancestry. In this specific example, it would likely reflect inaccurate sampling of the population in the analysis.

If allele frequency was 1 it would mean complete homozygosity and identity, which is very informative.

I'd simply refer to the wiki which is pretty good for this. Allele frequency is the entire basis of population genetics.
> 1 votes

---
Tags: genome, exome, gnomad
---


# Question
Title: I have a partially described gene from a plant, which tool may I use for predicting regolatory and UTRs sequences from a single DNA fragment?
I have a partially described gene from a plant with CDS and 3'UTR described but no 5'UTR; if I know the locus and upstream and downstream regions, is there a tool I can use for predicting regulatory DNA and UTRs for a given gene?

My annotation data source is the assembly of *Phaseolus vulgaris* here.

# Answer
> 1 votes
Assuming that you know the gene name, you would be able to look it up in a transcriptome assembly presumably, or perhaps NCBI has a full length transcript (search `txid3885[Organism:exp]` with your gene name).

For nearby regulatory DNA, you could consult an ATAC-seq dataset or other similar, though I doubt that causal links of regulatory DNA to your gene will have been directly assayed.

---
Tags: sequence-annotation, functional-annotation, tools
---