# Question
Title: Why is COI metabarcoding not used for prokaryotes?
I'm a beginner in bioinformatics working with environmental metabarcoding, and while reading up on different marker genes, I understand that COI is preferably used for animals thanks to its higher variability over 18S (although it depends on which taxonomic resolution fits best and which reference data is available). What I don't understand though, and can't find an answer to, is why COI metabarcoding isn't used for bacteria and archaea as well. Why is 16S preferred, although bacteria and archaea clearly have the COI gene? Is it not variable enough? Or are the primer binding sites too variable? Something else?

# Answer
The exact answer to your question is bacteria do not have a mitochondrion, therefore *cox1* is not a universal gene combining both Archaebacteria and bacteria. In contrast 16S is universal for these two kingdoms AND it is also single copy. Moreover, bacteria have diverse tropisms beyond the Krebs or Calvin cycles. Extensive studies have been performed on 16S and it provides an excellent taxonomic target in bacteria.

18S is a problematic target in eukaryotes because it can paralogous and that is not an ideal target in taxonomic studies, but has nicely conserved primer site: thus single copy genes are required/preferred. In addition, 18S is prone to indels (insertions/deletions) and these are difficult to model, align and understand the underlying molecular mechanism or selection advantage.

Every animal on earth has a mitochondrion and therefore *cox1*. Mitochondrion genes are guaranteed to be single copy in higher eukaryotes and *cox1* has been extensively sequenced throughout metazoans. Thereby barcoding is appropriate because there's an exact reference to identify a given species is usually available. It has been extremely successful in insects which are the most biodiverse and species rich taxonomic Class on earth: there are more insects than all other higher eukaryotic species combined.

---

In answer to the comments ...

* Whether mitochondrion were originally bacteria is a very difficult question and remains a theory (compartmentalisation theory - very old). I have not seen any rigorous phylogenetic support for this.
* *cox1* (COI) or *cox2* (COI) is not universal to bacteria, e.g. here. Its primary use is for aerobic organisms (oxygen consuming). Thus *cox1* has a strong Krebs cycle association but there are a multitude of anaerobic pathways.

Anaerobic bacteria don't necessarily need *cox1*/*cox2* and there are a lot of obligate anaerobic bacteria in addition to non-obligatory anaerobic bacteria (which can do both).

Marine metazoan are obligate aerobic organisms.

* There is no reason to suppose bacteria are strictly single copy *cox1*. mtDNA is tiny and hence copy number heavily constrained, bacterial genomes by contrast are much larger and could easily accommodate duplication events.
* There isn't much work on bacterial *cox1*, it has been pursued to attempt to augment antibiotic therapy. Degli Esposti (2020) is the closest I know to a comprehensive discussion on the bacterial *cox* genes.
> 2 votes

---
Tags: phylogenetics, phylogeny, barcode
---


# Question
Title: Automate the def function in pandas for correlation
```
def correlations(group):
    # Remove columns
    group = (
        group.reset_index()
        .drop(columns=["RNA1_Approved_Symbol"], level=0)
        .set_index("RNA1_ID")
    )
    # compute correlations
    correlation = group.T.corr()

    # rename pair column
    correlation = correlation.reset_index().rename(
        columns={"RNA1_ID": "RNA1_ID_new"}
    )

    # Remove NAs
    corr_result = (
        correlation.melt(value_name="Correlation", id_vars=["RNA1_ID_new"])
        .set_index(["RNA1_ID", "RNA1_ID_new"])
        .dropna()
    )
    return corr_result

def correlation_1(data):
    ## Create Pivot table
    pivot_library = data.pivot(
        index=["RNA1_Approved_Symbol", "RNA1_ID"], columns=["RNA2_ID"]
    )

    # Measure correlations between four id's
    return (
        pivot_library.groupby("RNA1_Approved_Symbol")
        .apply(correlations)
        .reset_index()
        .dropna()
    )

```

I am using the above function to calculate the correlation

```
data = df1
correlations = correlation_1(df1)

```

I want to automate the above function when I am using different data frame where there is `RNA2_Approved_Symbol` instead of `RNA1_Approved_Symbol` and `RNA2_ID` instead of `RNA1_ID`. Like using index 1 and index 2.

# Answer
If I have understood the question is to generalise the column names to permit a function to accept dataframes where you don't have the same column name, or it's different from the previous dataframe column names. Here the column could be called `RNA1_Approved_Symbol` or `RNA2_Approved_Symbol` in addition a different column could be called `RNA2_ID` or `RNA_ID`. No problem providing the column order is identical for the columns in question. If the orders are not identical see the code below the `----` to correct that.

If the order is identical the simple answer is `df.columns[i]`. Thus it's obtaining the column name on its position within the dataframe (`columns` command), where `i` is the column location.

Example, to rename `RNA1_ID` to `RNA1_ID_new`, or `RNA2_ID` to `RNA2_ID_new`

```
import re

df1 = pd.DataFrame({"perch": [5,6,7], "RNA1_ID": [1,2,4], "roach": [30,42,51]})
df2 = pd.DataFrame({"perch": [5,6,7], "RNA2_ID": [1,2,4], "roach": [30,42,51]})

def renamecolumn(df):
    myname = df.columns[1]
    newname = re.sub(r'$','_new', myname)
    df = df.rename(columns={myname: newname})
    return df

print (renamecolumn(df1))
print (renamecolumn(df2))

```

Output,

| - | perch | RNA1\_ID\_new | roach |
| --- | --- | --- | --- |
| 0 | 5 | 1 | 30 |
| 1 | 6 | 2 | 42 |
| 2 | 7 | 4 | 51 |

| - | perch | RNA2\_ID\_new | roach |
| --- | --- | --- | --- |
| 0 | 5 | 1 | 30 |
| 1 | 6 | 2 | 42 |
| 2 | 7 | 4 | 51 |

I noted there was a fair number of `set_index` commands so the earlier the column name can be trapped inside a variable the better, because `set_index` will change the column order.

---

Alternatively, the columns can be re-arranged to be in a given order if the original order between `df1` and `df2` is different. Here I'll move the column to the front.

```
mycols = df1.columns.tolist()
mycols = mycols[-2:] + mycols[:-2]
df1 = df1[mycols]
print (df1)

```

| - | RNA1\_ID | roach | perch |
| --- | --- | --- | --- |
| 0 | 1 | 30 | 5 |
| 1 | 2 | 42 | 6 |
| 2 | 4 | 51 | 7 |

Thus `RNA1_ID` has been moved from the second to the first column. If `df2` was in the last column repeat as above use `-1` instead of `-2` and both columns will now be at the start of the dataframe. Modify the core dataframe generalisation renaming accordingly, viz.

```
myname = df.columns[0]

```

... sorted.
> 1 votes

---
Tags: python, statistics, pandas
---


# Question
Title: Would duplicate climatic variables impact maxent?
I want to use maxent to make some future distribution predictions. I downloaded future climate data from worldclim. In worldclim's page there are three variables (max temp, min temp and precipitation) as well as a fourth titled "bioclimatic variables" which mentions that it contains the other three and some other new ones. My question is: if i download all four and use them for my analysis, will maxent produce reliable results? Or would the duplicate variables be a problem?

# Answer
Ok I understand better now so I'll write an answer.

Look at the bioclimatic data site. You can download them from the main page, but the suggestion to generate them from the time series data using an R package (`dismo`). The bioclimatic variables are thus various linear combinations of the 3 variable time series input.

Not having run this myself, I speculate, but it doesn't make sense to me that the bioclimatic variables are also time series, they appear to be summaries across the time series. So they're not *exactly* identical in information content, but they're close enough that I'm not sure it's worth the effort of figuring out how to munge the monthly data into a format that works with the summaries.

Someone already went to a lot of trouble to think about the biologically relevant representation of this data (e.g. what ANUCLIM does). I recommend relying on that work, or at the least looking at their arguments and deciding whether you agree enough to use their summaries.

It doesn't particularly matter what software you use, I'd argue (probably MaxEnt the species distribution package, but you haven't clarified). Most methods are going to be doing similar things with this data. They are going to find the variables with the best explanation and weigth them more heavily. Throwing in a couple extra correlated variables might get you a little margin, but it is unlikely to increase explanatory power.

If you want a formal methodology to compare the two, you could use something like AIC or BIC to compare the models formally with the 2 variable sets. But it hardly seems worth the effort of the data munging.

## Assumption:

I am assuming that you are not talking about adding ~identical variables also into your model, e.g. files here, but are instead talking about the less identical monthly data here. I'm *pretty sure* this is right, but I wanted to clarify.
> 1 votes

---
Tags: modelling, ecology
---


# Question
Title: Barrnap Bacterial rRNA Predictor script permission denied error when running Prokka
I am trying to run the **Prokka annotation tool** on my linux machine and have installed **Barrnap and added it to $PATH and given the folder chmod 777 permission**. However, when running Prokka, it gives the line `Could not parse version from 'sh: 1: barrnap: Permission denied'`.

Running Barrnap solely through the command line with a valid input works fine but does not seem to work when another script calls barrnap. I have tried checking path variables and giving max permissions to barrnap and Prokka but it does not seem to change this error. Here is the full log file:

```
[15:01:28] This is prokka 1.14.6
[15:01:28] Written by Torsten Seemann <torsten.seemann@gmail.com>
[15:01:28] Homepage is https://github.com/tseemann/prokka
[15:01:28] Local time is Thu Feb 23 15:01:28 2023
[15:01:28] You are saurab
[15:01:28] Operating system is linux
[15:01:28] You have BioPerl 1.7.8
[15:01:28] System has 16 cores.
[15:01:28] Will use maximum of 8 cores.
[15:01:28] Annotating as >>> Bacteria <<<
[15:01:28] Generating locus_tag from 'jobs/NC_004556.faa' contents.
[15:01:28] Setting --locustag ECJHPMJE from MD5 ec31963ed5d3d76e28501fa7e45ca273
[15:01:28] Creating new output folder: PROKKA_02232023
[15:01:28] Running: mkdir -p PROKKA_02232023
[15:01:28] Using filename prefix: PROKKA_02232023.XXX
[15:01:28] Setting HMMER_NCPU=1
[15:01:28] Writing log to: PROKKA_02232023/PROKKA_02232023.log
[15:01:28] Command: /home/saurab/BIOIN/basys2-app/sub_programs/prokka/bin/prokka jobs/NC_004556.faa
[15:01:28] Appending to PATH: /home/saurab/BIOIN/basys2-app/sub_programs/prokka/bin/../binaries/linux
[15:01:28] Appending to PATH: /home/saurab/BIOIN/basys2-app/sub_programs/prokka/bin/../binaries/linux/../common
[15:01:28] Appending to PATH: /home/saurab/BIOIN/basys2-app/sub_programs/prokka/bin
[15:01:28] Looking for 'aragorn' - found /home/saurab/BIOIN/phastest-app/sub_programs/aragorn
[15:01:28] Determined aragorn version is 001002 from 'ARAGORN v1.2   Dean Laslett'
[15:01:28] Looking for 'barrnap' - found /home/saurab/BIOIN/phastest-app/sub_programs/barrnap
[15:01:28] Coult not parse version from 'sh: 1: barrnap: Permission denied'

```

# Answer
> 3 votes
My suggestion is to run `prokka` from within `conda` by simply downloading miniconda then

```
conda create --name prokka
conda activate prokka
conda install -c conda-forge -c bioconda -c defaults prokka
prokka --help
prokka --version
prokka mycontigs.fa # running prokka

```

You could always perform a

```
chmod 700 /home/saurab/BIOIN/phastest-app/sub_programs/barrnap

```

or some variant like `chmod 500 ...`

or possibly using a recursive `chmod -R 500` that sort of thing. Eventually it will work.

**However**, it strikes me you are installing `prokka` outside any virtual environment such as `conda` or `brew` and I would not do that personally, because you are giving significant permissions to the algorithm.

`conda` (or `brew`) will sort out all the permission needed and there's no need to do anything with `chmod` stuff.

Good luck with the annotation.

---
Tags: sequence-annotation, perl, bacteria, functional-annotation
---


# Question
Title: Determine expression values for a defined list of genes across different clusters in scRNA data using seurat
I am analyzing publically available scRNA seq datasets on R using Seurat. I have created a seurat object, clustered and annotated different cell types. Now, I am interested in comparing the expression values and percentage of expressing cells within each cluster for a given set of genes I am interested in. I can visualize their expression using DotPlots, feature plots, etc... but how can I generate a table with the expression values (log2FC) for my defined set of genes across the different clusters/cell types? Thanks!

# Answer
> 1 votes
Assuming you're looking for *differential* expression values, the `FindMarkers` function within Seurat would be a good thing to try first:

```
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)
##                       p_val avg_log2FC pct.1 pct.2     p_val_adj
## FCGR3A        8.246578e-205   4.261495 0.975 0.040 1.130936e-200
## IFITM3        1.677613e-195   3.879339 0.975 0.049 2.300678e-191
## CFD           2.401156e-193   3.405492 0.938 0.038 3.292945e-189
## CD68          2.900384e-191   3.020484 0.926 0.035 3.977587e-187
## RP11-290F20.3 2.513244e-186   2.720057 0.840 0.017 3.446663e-182

```

\[Note the `avg_log2FC` column\]

More details can be found at the bottom of the Seurat tutorial pages for the PBMC 3K Guided Tutorial, *Finding differentially expressed features (cluster biomarkers)*:

https://satijalab.org/seurat/articles/pbmc3k\_tutorial.html

---
Tags: scrnaseq, seurat, gene-expression
---


title: Math on Pandas Columns

text:
# Question
Title: Math on Pandas Columns
I have a pandas dataframe that reads in a PAF file from minimap2. What I would like to do is take the first 5 columns of the data from to create a BED file.

I used this to extract the first 5 columns:

```
import pandas as pd
import numpy as np
df=pd.read_csv(tag1_paf, delimiter = "\t")
#print(df[4])

specific_columns=df.iloc[:,[0,1,2,3,4]]
print(specific_columns)

```

Which outputs this:

```
read1  47215  20591  21087 +
read2  14111  3478   3973  +
read3  30861  21367  21855 +
read4  4647   767    1257  -
read5  11706  9569   9990  +

```

But now I want to do math on the third and fourth columns depending on the strand its on. For example, if column 5 = "+" add 1000 to the third column.

I can't figure it out on my own - I've tried looking it up and using the series add/sub methods from pandas (https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.Series.add.html) but I'm still confused.

If anyone can help that would be very appreciated!

# Answer
The vectorised solution (fast) is,

```
import pandas as pd
filepath = '/pathtodir/test/bioinfo.csv'

df = pd.read_csv(filepath, delimiter=" ", header=None)

df.loc[df[4] == "+", 'new2'] = df[2] + 1000
df.loc[df[4] == "-", 'new2'] = df[2]
df = df.drop([2], axis=1)
df = df[[0,1,'new2',3,4]]
print (df)

```

| - | 0 | 1 | new2 | 3 | 4 |
| --- | --- | --- | --- | --- | --- |
| 0 | read1 | 47215 | 21591.0 | 21087 | + |
| 1 | read2 | 14111 | 4478.0 | 3973 | + |
| 2 | read3 | 30861 | 22367.0 | 21855 | + |
| 3 | read4 | 4647 | 767.0 | 1257 | - |
| 4 | read5 | 11706 | 10569.0 | 9990 | + |

An alternative method to re-order the columns then the answer is Automate the def function in pandas for correlation

Original data

|  | 0 | 1 | 2 | 3 | 4 |
| --- | --- | --- | --- | --- | --- |
| 0 | read1 | 47215 | 20591 | 21087 | + |
| 1 | read2 | 14111 | 3478 | 3973 | + |
| 2 | read3 | 30861 | 21367 | 21855 | + |
| 3 | read4 | 4647 | 767 | 1257 | - |
| 4 | read5 | 11706 | 9569 | 9990 | + |

---

Final solution

```
df = pd.read_csv(filepath, delimiter=" ", header=None)
df.loc[df[4] == "+", 'new2'] = df[2] + 1000
df.loc[df[4] == "-", 'new2'] = df[2] - 1000
df.loc[df['new2'] < 0, 'new2'] = 1
df = df.drop([2], axis=1)
df = df[[0,1,'new2',3,4]]
print (df)

```

| - | 0 | 1 | new2 | 3 | 4 |
| --- | --- | --- | --- | --- | --- |
| 0 | read1 | 47215 | 21591.0 | 21087 | + |
| 1 | read2 | 14111 | 4478.0 | 3973 | + |
| 2 | read3 | 30861 | 22367.0 | 21855 | + |
| 3 | read4 | 4647 | 1.0 | 1257 | - |
| 4 | read5 | 11706 | 10569.0 | 9990 | + |

... worked too
> 4 votes

---
Tags: python, phylogenetics, phylogeny, sequence-analysis, pandas
---


# Question
Title: Is there a graphical/interactive 16S rRNA clustering method?
I've been doing phylogenetics with large (hundreds) 16S rRNA sequences lately.

Usually I'm focusing on one order, and using a combination of trees and sequence similarity to assess stuff like 'is this a genus', 'is this a separate genus', 'is this sequence actually in the order of interest'. It feels like a lot of what I do is repetitive and could be automated more. In particular, it feels like there should be a simpler way for having an assemble view at sequence similarity rather that calculating a Percentage Identity Matrix and then putting it in Excel with conditional formatting. I'm imagining something like this clustering diagram for species (https://www.researchgate.net/publication/360738216/figure/fig4/AS:1182363421806594@1658908815491/The-t-SNE-clustering-of-bacterial-genomes-in-the-code-RA-space-at-different-taxonomic.png), but for sequences instead (using sequence similarity values, at least). Still, I can't seem to find something like that for 16S rRNA that's as visually compelling (or even as visually helpful as my Excel sheet). There are a number of programs for clustering OTUs, but the ones I've encountered don't seem to have a graphical / interactive interface.

Does something like this exist? Thank you for your time!

# Answer
My recommendation is MegaX. It's the closest thing in phylogenetics to an all encompassing GUI.

1. This has a nice visual alignment editor which highlights conserved, variable and parsimonious sites.
2. It will produce distance matrices based on any number of criteria, e.g. % homology, Jukes-Cantor, Kimura-2 parameter etc ...
3. Good for bootstrapping distance-calculations as well as can perform ML based tree building methods with full parameterisation.
4. Nice tree visualisation, which can be edited, switched to viewing bootstrapping (if performed), or to branch genetic distances (again under any distance criteria, e.g. %age homology etc ..). It does loads of other stuff, which in this case are not needed.

You can export the tree and import into FigTree to produce collapsed clades represented by coloured triangles. In summary for easy GUI stuff MegaX, then FigTree.

There are no PCA/MDA based analysis for 16S sequence data because eiganvalue-based calculations would override lots of theory involved in point mutation modelling and the resulting trees. I could replicated the PCA output in the question with sequence data, its not hard to do - but it's not publishable.
> 1 votes

---
Tags: phylogenetics, phylogeny, sequence-analysis, rna, software-recommendation
---


title: How to interpret fish plot in R

text:
# Question
Title: How to interpret fish plot in R
I am working on subclonal evolution of multiple myeloma. After literature survey, I found people use Fish plot to represent tumor evolution. The example of fish plot can be found here. But I am not able to interpret this plot. How this plot represent tumor evolution and subclones in tumor. Can anybody explain with example. Thanks.

# Answer
> 2 votes
Author of the package here, stumbled across this today. The existing explanation is close, but not entirely correct.

Our paper describes it like this: "We first created a fish plot of an AML patient with a chemotherapy-induced bottleneck that eliminated one subclone, while another survived and drove the relapse".

In other words, the grey represents a tumor, starting from a single cell (point at far left). The tumor grows and expands, acquiring additional mutations (red) that increase the relative growth rate of a subclone. When sampled at day 0 (when the patient showed up with cancer), the red subclone made up about a third of the tumor cells. Chemotherapy was applied, shrinking the tumor dramatically and completely eliminating the red subclone (patient went into remission).

Unfortunately, remission doesn't always mean "cured", as this tumor persisted, and picked up more mutations (orange subclone), and then within that subclone, picked up even more (yellow subclone). Each gave those cells sufficient growth advantage that at relapse, essentially all of the tumor cells were comprised of the orange subclone, with a subset of cells that contained both the orange and yellow mutations.

There are also subtleties about whether the orange and yellow mutations pre-existed at very low levels (below the sensitivity of the sequencing that was performed), or were actually brand new mutations, but that's outside the scope of this explanation!

# Answer
> 5 votes
Each color is a clone (set of mutations) or an individual mutation (depends on how the input data was generated). This is a cancer with the bottleneck likely being induced by a treatment (e.g. UV, chemo, drugs, etc...)

Based on the vignette this is likely an AML (leukemia).

It looks like grey are healthy cells. At **day 0** a specific clone (red) existed at a low frequency (~20-40%). The patient (or cells) was then exposed to a treatment and on **day 150** another measurement was taken. In this case two novel clones (subclones of red) have taken over the population. They seem to co-exist or are co-dominate clones with orange being 50-60% of the population and yellow being 30-40%.

**The example says "*Simple plot with multiple subclones*" so: the orange and yellow populations have additional (likely mutually exclusive) mutations that were not present in the red clone - but biologically also have all the mutations in the red clone. So my first assumption was right - each color is a proxy for a set of mutations.**

You could have generated the mutation data (i.e. allele frequencies) from a wide variety of sources (scRNA-Seq, exome-sequencing, PCR genotyping, etc...)

---
Tags: cancer, figure-explanation
---


# Question
Title: variant visualisation of VCF files for 30+ samples
I am relatively new to genome analysis and would like to compare 30+ vcf files, 1 sample per file. After filtering there are about ~20 variants per vcf file.

I would like to evaluate to mutational landscape between samples and for each sample as well.

One option that I found is a ***mutation heatmap*** where each row is a sample and each column a position in the genome. Is this the preferred way of displaying this type of information and how would I include gene information?

At this point I am interested in getting an overview first of the affected positions and then dig deeper using e.g. zygosity.

Thank you very much for helping out!

**2nd edit** I would also like to evaluate if there is a bias towards certain genes. For example are certain oncogenes mutated more frequently in this cohort or is the accumulation random and unique within each sample?

Best, M

# Answer
> 3 votes
Single sample PVCF files are generally not combinable in the way you describe. PVCF stands for "project" VCF. As far as I know, PVCF isn't defined anywhere but, it refers to a multi-sample VCF that does not contain reference blocks. It only contains lines for which at least one sample has at least one alternate allele. We distinguish this from a GVCF which is a single-sample VCF that does contain reference blocks. We can't combine two PVCFs in a principled way because some variant might only appear in one PVCF. We lack the reference block information necessary to decide if a sample in the other VCF is 0/0 or "no call".

Do all the VCF files have the same set of variants? If all the VCF files have the same set of variants, you might try using Hail (disclaimer: I'm a Hail maintainer):

```
import hail as hl

vcf_files = ['sample1.vcf.bgz', ...]
mts = [hl.import_vcf(f) for f in vcf_files]
mt = hl.MatrixTable.union_cols(mts)  # assumes each file has the same set of variants
mt.write('combined.mt')  # save to Hail's native format so the following code is fast

mt = hl.read_matrix_table('combined.mt')
mt.show(n_rows=5, n_cols=10)  # see a table of the first 5 rows and 10 cols

mt = hl.variant_qc(mt)  # adds a bunch of variant QC annotations
mt.describe()  # inspect available annotations

```

You mentioned "genomic distribution". Are you referring to the allele frequency in this sample? You can see a genome-wide scatter plot of reference allele frequency like this:

```
from hail.ggplot import *
fig = ggplot(mt) + geom_point(aes(x=mt.locus, y=mt.variant_qc.AF[0]))
fig.show()

```

---
Tags: vcf, variant-calling, visualization
---


# Question
Title: calculating mutation frequencies for every gene
I have a dataset for mutation data and I want to calculate mutation frequencies across all genes

**df** (This is only the small subset of data)

```
Gene name   Sample id   MUTATION_ID Mutation Description
ARID1B  2719660 171258500   Substitution - Missense
ARID1B  2719659
ARID1B  2719661 171258501   Substitution - Missense
ARID1B  2719662
ARID1B  2719663 171258501   Substitution - Nonsense
CD58    2878555 110346783   Substitution - Nonsense
CD58    2877956
CD58    2878557
CD58    2877958 110346784   Substitution - Nonsense
CD58    2878559 110346785   Substitution - Nonsense
CD58    2877960
MRE11   2861617 123320443   Substitution - coding silent
MRE12   2861617 123320444   Substitution - coding
MRE13   2861617
MRE14   2861617 123320445   Substitution - coding silent
MRE15   2861617
MRE16   2861617 123320446   Substitution - coding

```

The formula for calculating the mutation is

`Positives ÷ (Positives + Negatives) x 100`

where,

`Positives` = No of samples where `MUTATION_ID` is present

`Negative` = No of samples `MUTATION_ID` for the sample

I want to calculate mutation frequency for every gene in the `column_1:Gene name` with python script

I tried the following code

```
df = df.groupby("Gene name").count()
Positives = df["MUTATION_ID"]
Negatives = df["Sample id"] - df["MUTATION_ID"]
df['Mutation_Frequency'] = Positives / (Positives + Negatives) * 100

```

# Answer
The code is fine. The only issue is you'd put \s characters in-between the column header names. I've place an additional "mutation\_frequency" column because thats how I'd do it.

```
df = pd.read_csv('/pathtodir/test/bioinfo.csv',sep="\s",engine='python')
df = df.groupby("Genename").count()
df['Mutation_Frequency1'] = ((df["MUTATION_ID"] / (df["Sampleid"] + (df["Sampleid"] - df["MUTATION_ID"]))) * 100).round(2)
df['Mutation_Frequency2'] = (df["MUTATION_ID"] / df["Sampleid"])*100
df.drop(['Description', 'Mutation'], axis=1, inplace=True)
print(df)

```

| - | Sampleid | MUTATION\_ID | Mutation\_Frequency1 | Mutation\_Frequency2 |
| --- | --- | --- | --- | --- |
| ARID1B | 5 | 3 | 42.86 | 60.0 |
| CD58 | 6 | 3 | 33.33 | 50.0 |
| MRE11 | 1 | 1 | 100.00 | 100.0 |
| MRE12 | 1 | 1 | 100.00 | 100.0 |
| MRE13 | 1 | 0 | 0.00 | 0.0 |
| MRE14 | 1 | 1 | 100.00 | 100.0 |
| MRE15 | 1 | 0 | 0.00 | 0.0 |
| MRE16 | 1 | 1 | 100.00 | 100.0 |
> 2 votes

---
Tags: python, mutations, pandas
---


# Question
Title: Inconsistent replicate numbers in RNA-seq
I have 3 groups for the RNA-seq analysis (Control, treatment A and treatment B). There are 2 replicates for control and treatment A and 3 replicates for treatment B (lost 2 replicates due to a mistake in sample handling). May I ask how will the inconsistent replicate numbers impact the results (e.g. DEGs)?

# Answer
> 2 votes
It will not do any harm other than a reduction in statistical power compared to more samples per group. Common analysis tools such as DESeq2, edgeR and limma will work on unequal sample size just fine.

---
Tags: rna-seq, genomics, differential-expression, gene-expression
---


# Question
Title: Which basecaller for nanopore is the best to produce event tables with information about the block size/move table?
I previously used albacore version 2.3.1 to make initial move tables, but then I re-squiggle using Tombo version 1.5.1 to fix the errors.

Example of move table produced by albacore, in which all the moves have the same size:

When it is re-squiggled by Tombo (which produces variable length moves/blocks):

I am starting a new project, and I can repeat the same pipeline, however, my concern is that albacore version 2.3.1 is too old now and I was wondering to ask if there are better pipelines/software for base-calling that produces initial move tables so that I can build on top of.

# Answer
> 3 votes
ONT's newest basecaller, Dorado, can do this for SAM output with the `--emit-moves` argument. It was added in this commit.

Dorado can be downloaded from here:

https://github.com/nanoporetech/dorado

Here's an example command line:

```
$ ~/install/dorado/dorado-0.2.1-linux-x64/bin/dorado basecaller --emit-moves \
  ~/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0 DAE_pod5 | \
  samtools view -b - > called_DAE_gDNA_reads_all.bam

```

And a few lines of example output:

```
@HD     VN:1.6  SO:unknown
@PG     ID:basecaller   PN:dorado       VN:0.2.1+c70423e        CL:dorado basecaller --emit-moves /home/minion/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0 DAE_pod5
@PG     ID:samtools     PN:samtools     PP:basecaller   VN:1.10 CL:samtools view -b -
@PG     ID:samtools.1   PN:samtools     PP:samtools     VN:1.10 CL:samtools view -h called_DAE_gDNA_reads_all.bam
@RG     ID:460fc479-3b30-4ed5-801f-a82bbe556101_dna_r10.4.1_e8.2_260bps_sup@v4.1.0      PU:PAK83314     PM:     DT:2023-02-08T10:22:53.256+00:00        PL:ONT  DS:basecall_model=dna_r10.4.1_e8.2_260bps_sup@v4.1.0 runid=460fc479-3b30-4ed5-801f-a82bbe556101 LB:DAE_gDNA     SM:DAE_gDNA
@RG     ID:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0      PU:PAK83314     PM:     DT:2023-02-08T20:28:34.34+00:00 PL:ONT  DS:basecall_model=dna_r10.4.1_e8.2_260bps_sup@v4.1.0 runid=227379b1-c5fc-465d-9e24-eb771133f11e LB:no_sample    SM:no_sample
003d84eb-ce0f-47ea-9c7c-fe939fe8702e    4       *       0       0       *       *       0       110     GTGCCTCCTTGCTTTGGTTCAGTTCTATCCAGCCTTTAGATATCCCAGCCCTCACAGGGTTTGTGTGACCTTGTTGCATTATCATGCGAAACACTTTCGCGTTTTTAGTG  %%%%%(&%%&%%'((''))*.-,-''&&&'&&&(%$$$%&&&'%%%%&(*))'&&&&)-05>6;>>777887++'&&&&&&%%%&*.....0820001:;;;<<87''&&  qs:i:7  du:f:0.526      ns:i:2104       ts:i:10 mx:i:4  ch:i:1037       st:Z:2023-02-09T07:07:00.271+00:00      rn:i:96929      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:98.674     sd:f:22.395     sv:Z:quantile   RG:Z:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,1,0,0,1,0,0,1,0,0,0,0,1,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,1,0,0,0,0,0,1,0,1,0,1,1,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,0,0,1,0,0,0,1,0,0,1,0,0,1,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,1,0,0,0,1,0,1,0,1
002e48ca-84bf-4419-8370-943d1ca3ac64    4       *       0       0       *       *       0       107     TATGTCTAACTTACTTCGTTCAGTTGTGAAAACTGGGTGTTTAACCTTTGATTCTATCGTGTTTCCCTAGTTTTCGCATTTGACGTGAAACGCTTTCACATTTTGGT     &&'''&%%%%((((*+)***())((')&%%%%$%%%'--,,,,-6996/,+-5ABBBCCFBCE>88:1764424/--,,--&%$$%*08??;;;;<B<:((('()//     qs:i:9  du:f:0.92975    ns:i:3719       ts:i:10 mx:i:4  ch:i:1058       st:Z:2023-02-09T07:07:08.905+00:00      rn:i:62167      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:85.094     sd:f:23.944     sv:Z:quantile   RG:Z:227379b1-c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,0,0,0,0,0,1,1,1,0,1,0,1,0,0,1,0,0,1,0,1,1,0,0,1,0,1,0,0,0,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,1,0,1,0,0,1,0,1,0,1,1,0,1,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,1,0,0,1,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,0,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,1,0,1,0,0,1,0,1,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,0,1,0,1,0,0,1,0,0,0,1,0,0
0138d83a-da9f-4468-bafc-b974c4e9eece    4       *       0       0       *       *       0       234     ATGTTTATGTGTTCTTCGTTAGTAAGCCACCGATATTGACTCCAGGGTTTGTGTAACCTTGTTTTCGCATTTATCGTGAAACGCTTTCGCGTTTTTCGTCGCTTCAAATAGTTTTAATTCAGTTGTCTAAAAATGTTAACATATGGGTATTATTTCTCTAGAACATTTTAGTTGTGTAATTTCTACCATTTAAATGGGGTGAATAAATAACCAAATCATTGAACAACTGGAATG      &&&')'&%$$##$&&''+,-)''&'')**'&%&''))&&&)'((*+8<AEEFGHJIJIJIICA@;843344:254,,,,CEE>====7<<<>>@=?+***'&&''(''&%%%%$%%&)1//..-+,+*'&'(+2+++45<ACCC66677@75556+((((55>@AEFHFHGFEHHJHHFHIKJGFFFKLHZEGEGGIDFDDCDDERKIEGDG>=>5332/'&&&&)**,2242.      qs:i:10 du:f:0.69975    ns:i:2799       ts:i:10 mx:i:4  ch:i:1933       st:Z:2023-02-09T07:07:01.924+00:00      rn:i:77788      f5:Z:PAK83314_fail_unclassified_227379b1_ea81734c_853.pod5      sm:f:97.397     sd:f:27.276     sv:Z:quantile   RG:Z:227379b1c5fc-465d-9e24-eb771133f11e_dna_r10.4.1_e8.2_260bps_sup@v4.1.0    mv:B:c,5,1,0,0,0,1,0,1,0,1,0,1,0,0,1,1,0,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,0,0,1,0,0,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,1,0,1,1,0,0,1,1,0,0,0,1,0,0,0,1,0,0,1,0,0,1,0,1,0,0,0,0,1,1,0,0,1,0,1,0,0,1,0,1,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,1,0,0,1,0,0,1,0,0,1,0,0,0,1,0,1,0,1,0,1,1,0,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,1,0,0,0,0,0,0,0,1,1,0,0,0,1,0,1,0,1,0,1,0,0,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,1,0,1,1,0,1,0,1,1,0,1,0,1,0,0,0,0,0,0,1,0,1,0,1,0,1,0,1,1,1,0,1,0,1,0,1,1,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,1,0,1,1,0,1,0,1,1,1,1,0,1,0,1,1,1,0,0,1,1,0,0,0,0,1,0,1,0,1,0,0,1,0,0,0,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,1,1,1,0,0,1,0,1,1,0,1,0,1,0,0,0,0,0,1,0,1,0,0,0,0,0,0,0,0,1,0,1,0,1,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,0,0,1,0,1,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,0,1,0,1,0,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0,1,0,0,1,0,0,1,0,0,1,0,1,1,0,1,0,1,0,1,1,0,1,1,0,1,1,0,1,0,1,0,0,1,0,0,1,1,0,0,1,0,0,1,0,0,1,1,1,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1

```

Based on the commit, the move table is attached to the `mv` tag for each read, indicating an array of 8-bit signed integers (represented as comma-separated values in the SAM output). I'd guess (based on how the calling has worked in the past) that these represent strides comprising a fixed number of consecutive samples, so (as before) it'll need correction to convert the stride width to variable width events.

For the model I'm using here, my guess from the model toml file is that the stride width is 5:

$ grep -B 8 stride ~/install/dorado/dorado-0.2.1-linux-x64/dna_r10.4.1_e8.2_260bps_sup@v4.1.0/config.toml
[encoder]
type = "serial"
[[encoder.sublayers]]
type = "convolution"
insize = 1
size = 16
bias = true
winlen = 5
stride = 1
--
max = 3.5

[[encoder.sublayers]]
type = "convolution"
insize = 16
size = 16
bias = true
winlen = 5
stride = 1
--
max = 3.5

[[encoder.sublayers]]
type = "convolution"
insize = 16
size = 1024
bias = true
winlen = 19
stride = 5
```

---
Tags: nanopore, base-calling
---


title: Sanger Sequencing Knitting Error

text:
# Question
Title: Sanger Sequencing Knitting Error
I am doing a project where I am reproducing the analysis from the article "sangeranalyseR: Simple and Interactive Processing of Sanger Sequencing Data in R". Below is the example chunk for Sanger Alignment. Below are the package and datasets.

```
library(sangeranalyseR)
data("sangerAlignmentData")
data("sangerContigData")
data("sangerReadFData")

```

Below is the code for Sanger Alignment, similar to the code for Sanger Contig and Sanger Reading.

```
rawDataDir <- system.file("extdata", package = "sangeranalyseR")
fastaFN <- file.path(rawDataDir, "fasta", "SangerAlignment", "Sanger_all_reads.fa")
my_sangerAlignmentFa <- new("SangerAlignment",
                            inputSource          = "FASTA",
                            processMethod        = "REGEX",
                            FASTA_File           = fastaFN,
                            REGEX_SuffixForward  = "_[0-9]*_F$",
                        REGEX_SuffixReverse  = "_[0-9]*_R$",
                            refAminoAcidSeq      = "SRQWLFSTNHKDIGTLYFIFGAWAGMVGTSLSILIRAELGHPGALIGDDQIYNVIVTAHAFIMIFFMVMPIMIGGFGNWLVPLMLGAPDMAFPRMNNMSFWLLPPALSLLLVSSMVENGAGTGWTVYPPLSAGIAHGGASVDLAIFSLHLAGISSILGAVNFITTVINMRSTGISLDRMPLFVWSVVITALLLLLSLPVLAGAITMLLTDRNLNTSFFDPAGGGDPILYQHLFWFFGHPEVYILILPGFGMISHIISQESGKKETFGSLGMIYAMLAIGLLGFIVWAHHMFTVGMDVDTRAYFTSATMIIAVPTGIKIFSWLATLHGTQLSYSPAILWALGFVFLFTVGGLTGVVLANSSVDIILHDTYYVVAHFHYVLSMGAVFAIMAGFIHWYPLFTGLTLNNKWLKSHFIIMFIGVNLTFFPQHFLGLAGMPRRYSDYPDAYTTWNIVSTIGSTISLLGILFFFFIIWESLVSQRQVIYPIQLNSSIEWYQNTPPAEHSYSELPLLTN",
                            minReadsNum           = 2,
                            minReadLength         = 20,
                            minFractionCall       = 0.5,
                            maxFractionLost       = 0.5,
                            geneticCode           = GENETIC_CODE,
                            acceptStopCodons      = TRUE,
                            readingFrame          = 1,
                            processorsNum         = 1)

launchApp(my_sangerAlignmentFa)

writeFasta(my_sangerAlignmentFa,
           outputDir         = tempdir(),
           compress          = FALSE,
           compression_level = NA,
           selection         = "all")

generateReport(my_sangerAlignmentFa,
               outputDir           = tempdir(),
               includeSangerRead   = TRUE,
               includeSangerContig = TRUE)

```

There are two problems that occurred for me. One is that the launchapp() function created an error for me stating

```
INFO [2023-03-01 14:25:52] Your input is 'SangerAlignment' S4 instance
INFO [2023-03-01 14:25:52] SangerAlignment with 'FASTA' inputSource
cannot run Shiny app (You don't need to do trimming or base calling)
NULL

```

How can I fix this error?

Another error that occurred within knitting. I ran the same code without the launchapp() function and an error that I never seen before said this

```
Quitting from lines NA-219 (SangerRead_Report_fasta.Rmd)  Error in parse_block(g[-1], g[1], params.src, markdown_mode) :
Duplicate chunk label 'unnamed-chunk-1', which has been used for the chunk:
library(sangeranalyseR) library(rentrez) Calls: <Anonymous> ...
process_file -> split_file -> lapply -> FUN -> parse_block

Execution halted

```

Is it my RStudio or is the code itself. Is there anything I am missing?

sessionInfo(sangeranalysesR)
R version 4.2.2 (2022-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 22621)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.utf8
[2] LC_CTYPE=English_United States.utf8
[3] LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C
[5] LC_TIME=English_United States.utf8

attached base packages:
character(0)

other attached packages:
[1] sangeranalyseR_1.8.0

loaded via a namespace (and not attached):
  [1] utils_4.2.2            nlme_3.1-161
  [3] fs_1.6.1               bitops_1.0-7
  [5] bit64_4.0.5            httr_1.4.4
  [7] GenomeInfoDb_1.34.6    tools_4.2.2
  [9] bslib_0.4.2            utf8_1.2.2
 [11] R6_2.5.1               DT_0.27
 [13] DBI_1.1.3              lazyeval_0.2.2
 [15] BiocGenerics_0.44.0    colorspace_2.1-0
 [17] rhdf5filters_1.10.0    ade4_1.7-20
 [19] withr_2.5.0            tidyselect_1.2.0
 [21] gridExtra_2.3          phangorn_2.11.1
 [23] bit_4.0.5              compiler_4.2.2
 [25] cli_3.5.0              datasets_4.2.2
 [27] shinyjs_2.1.0          plotly_4.10.1
 [29] ggdendro_0.1.23        base_4.2.2
 [31] sass_0.4.5             scales_1.2.1
 [33] quadprog_1.5-8         stringr_1.5.0
 [35] digest_0.6.31          rmarkdown_2.20
 [37] sangerseqR_1.34.0      XVector_0.38.0
 [39] pkgconfig_2.0.3        htmltools_0.5.4
 [41] sessioninfo_1.2.2      fastmap_1.1.0
 [43] excelR_0.4.0           grDevices_4.2.2
 [45] htmlwidgets_1.6.1      rlang_1.0.6
 [47] rstudioapi_0.14        RSQLite_2.2.20
 [49] shiny_1.7.4            jquerylib_0.1.4
 [51] generics_0.1.3         jsonlite_1.8.4
 [53] BiocParallel_1.32.5    zip_2.2.2
 [55] dplyr_1.1.0            RCurl_1.98-1.10
 [57] magrittr_2.0.3         GenomeInfoDbData_1.2.9
 [59] biomformat_1.26.0      Matrix_1.5-3
 [61] Rhdf5lib_1.20.0        Rcpp_1.0.10
 [63] munsell_0.5.0          S4Vectors_0.36.1
 [65] fansi_1.0.4            logger_0.2.2
 [67] DECIPHER_2.26.0        ape_5.6-2
 [69] shinycssloaders_1.0.0  lifecycle_1.0.3
 [71] yaml_2.3.7             stringi_1.7.12
 [73] MASS_7.3-58.2          zlibbioc_1.44.0
 [75] rhdf5_2.42.0           plyr_1.8.8
 [77] grid_4.2.2             blob_1.2.3
 [79] parallel_4.2.2         promises_1.2.0.1
 [81] shinydashboard_0.7.2   crayon_1.5.2
 [83] methods_4.2.2          lattice_0.20-45
 [85] Biostrings_2.66.0      zeallot_0.1.0
 [87] knitr_1.42             pillar_1.8.1
 [89] igraph_1.3.5           seqinr_4.2-23
 [91] reshape2_1.4.4         codetools_0.2-19
 [93] stats4_4.2.2           fastmatch_1.1-3
 [95] reprex_2.0.2           BiocVersion_3.16.0
 [97] glue_1.6.2             evaluate_0.20
 [99] BiocManager_1.30.20    data.table_1.14.6
[101] vctrs_0.5.2            httpuv_1.6.8
[103] graphics_4.2.2         gtable_0.3.1
[105] purrr_1.0.1            tidyr_1.3.0
[107] cachem_1.0.6           ggplot2_3.4.0
[109] openxlsx_4.2.5.1       xfun_0.36
[111] mime_0.12              xtable_1.8-4
[113] later_1.3.0            viridisLite_0.4.1
[115] tibble_3.1.8           stats_4.2.2
[117] memoise_2.0.1          IRanges_2.32.0
[119] shinyWidgets_0.7.6     ellipsis_0.3.2
[121] BiocStyle_2.26.0
```
</code></pre>

# Answer
> 3 votes
Your code runs without error on my laptop (macOS 13.2.1; R v4.2.2; Rstudio 2022.07.2), so I have some further questions that might help with troubleshooting (too many to include in a comment):

* Does the example from the docs run without errors/warnings on your system?
* Did you install the sangeranalyseR package via Bioconductor (per the docs)?
* Do you have the most up-to-date Bioconducter installed?
* Have you tried running the example code in a 'fresh' session?
* Are you able to please run the command `sessionInfo()` and edit your question to include the output?
* If you run your example via the reprex package, are there any differences between your output and my output?

The code from your question run on my laptop:

```
# BiocManager::install("sangeranalyseR")
library(sangeranalyseR)
#> Loading required package: stringr
#> Loading required package: ape
#> Loading required package: Biostrings
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> Loading required package: stats4
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:base':
#>
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: XVector
#> Loading required package: GenomeInfoDb
#>
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:ape':
#>
#>     complement
#> The following object is masked from 'package:base':
#>
#>     strsplit
#> Loading required package: DECIPHER
#> Loading required package: RSQLite
#> Loading required package: parallel
#> Loading required package: reshape2
#> Loading required package: phangorn
#> Loading required package: sangerseqR
#> Loading required package: gridExtra
#>
#> Attaching package: 'gridExtra'
#> The following object is masked from 'package:BiocGenerics':
#>
#>     combine
#> Loading required package: shiny
#> Loading required package: shinydashboard
#>
#> Attaching package: 'shinydashboard'
#> The following object is masked from 'package:graphics':
#>
#>     box
#> Loading required package: shinyjs
#>
#> Attaching package: 'shinyjs'
#> The following object is masked from 'package:shiny':
#>
#>     runExample
#> The following object is masked from 'package:RSQLite':
#>
#>     show
#> The following object is masked from 'package:Biostrings':
#>
#>     show
#> The following object is masked from 'package:GenomeInfoDb':
#>
#>     show
#> The following object is masked from 'package:XVector':
#>
#>     show
#> The following object is masked from 'package:IRanges':
#>
#>     show
#> The following object is masked from 'package:S4Vectors':
#>
#>     show
#> The following object is masked from 'package:stats4':
#>
#>     show
#> The following objects are masked from 'package:methods':
#>
#>     removeClass, show
#> Loading required package: data.table
#>
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:reshape2':
#>
#>     dcast, melt
#> The following object is masked from 'package:IRanges':
#>
#>     shift
#> The following objects are masked from 'package:S4Vectors':
#>
#>     first, second
#> Loading required package: plotly
#> Loading required package: ggplot2
#>
#> Attaching package: 'plotly'
#> The following object is masked from 'package:ggplot2':
#>
#>     last_plot
#> The following object is masked from 'package:phangorn':
#>
#>     add_boxplot
#> The following object is masked from 'package:XVector':
#>
#>     slice
#> The following object is masked from 'package:IRanges':
#>
#>     slice
#> The following object is masked from 'package:S4Vectors':
#>
#>     rename
#> The following object is masked from 'package:stats':
#>
#>     filter
#> The following object is masked from 'package:graphics':
#>
#>     layout
#> Loading required package: DT
#>
#> Attaching package: 'DT'
#> The following objects are masked from 'package:shiny':
#>
#>     dataTableOutput, renderDataTable
#> Loading required package: zeallot
#> Loading required package: excelR
#> Loading required package: shinycssloaders
#> Loading required package: ggdendro
#> Loading required package: shinyWidgets
#>
#> Attaching package: 'shinyWidgets'
#> The following object is masked from 'package:shinyjs':
#>
#>     alert
#> Loading required package: openxlsx
#> Loading required package: tools
#> Loading required package: rmarkdown
#> Loading required package: knitr
#> Loading required package: seqinr
#>
#> Attaching package: 'seqinr'
#> The following object is masked from 'package:shiny':
#>
#>     a
#> The following object is masked from 'package:sangerseqR':
#>
#>     read.abif
#> The following object is masked from 'package:Biostrings':
#>
#>     translate
#> The following objects are masked from 'package:ape':
#>
#>     as.alignment, consensus
#> Loading required package: BiocStyle
#>
#> Attaching package: 'BiocStyle'
#> The following objects are masked from 'package:rmarkdown':
#>
#>     html_document, md_document, pdf_document
#> The following object is masked from 'package:shiny':
#>
#>     markdown
#> Loading required package: logger
#> Welcome to sangeranalyseR
data("sangerAlignmentData")
data("sangerContigData")
data("sangerReadFData")

rawDataDir <- system.file("extdata", package = "sangeranalyseR")
fastaFN <- file.path(rawDataDir, "fasta", "SangerAlignment", "Sanger_all_reads.fa")
my_sangerAlignmentFa <- new("SangerAlignment",
                            inputSource          = "FASTA",
                            processMethod        = "REGEX",
                            FASTA_File           = fastaFN,
                            REGEX_SuffixForward  = "_[0-9]*_F$",
                        REGEX_SuffixReverse  = "_[0-9]*_R$",
                            refAminoAcidSeq      = "SRQWLFSTNHKDIGTLYFIFGAWAGMVGTSLSILIRAELGHPGALIGDDQIYNVIVTAHAFIMIFFMVMPIMIGGFGNWLVPLMLGAPDMAFPRMNNMSFWLLPPALSLLLVSSMVENGAGTGWTVYPPLSAGIAHGGASVDLAIFSLHLAGISSILGAVNFITTVINMRSTGISLDRMPLFVWSVVITALLLLLSLPVLAGAITMLLTDRNLNTSFFDPAGGGDPILYQHLFWFFGHPEVYILILPGFGMISHIISQESGKKETFGSLGMIYAMLAIGLLGFIVWAHHMFTVGMDVDTRAYFTSATMIIAVPTGIKIFSWLATLHGTQLSYSPAILWALGFVFLFTVGGLTGVVLANSSVDIILHDTYYVVAHFHYVLSMGAVFAIMAGFIHWYPLFTGLTLNNKWLKSHFIIMFIGVNLTFFPQHFLGLAGMPRRYSDYPDAYTTWNIVSTIGSTISLLGILFFFFIIWESLVSQRQVIYPIQLNSSIEWYQNTPPAEHSYSELPLLTN",
                            minReadsNum           = 2,
                            minReadLength         = 20,
                            minFractionCall       = 0.5,
                            maxFractionLost       = 0.5,
                            geneticCode           = GENETIC_CODE,
                            acceptStopCodons      = TRUE,
                            readingFrame          = 1,
                            processorsNum         = 2)
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.05 secs
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.03 secs
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.03 secs
#> Assessing frameshifts in nucleotide sequences:
#> ================================================================================
#>
#> Time difference of 0.03 secs

writeFasta(my_sangerAlignmentFa,
           outputDir         = tempdir(),
           compress          = FALSE,
           compression_level = NA,
           selection         = "all")

generateReport(my_sangerAlignmentFa,
               outputDir           = tempdir(),
               includeSangerRead   = TRUE,
               includeSangerContig = TRUE)
#> processing file: SangerRead_Report_fasta.Rmd
#> output file: SangerRead_Report_fasta.knit.md
#> /usr/local/bin/pandoc +RTS -K512m -RTS SangerRead_Report_fasta.knit.md --to html4 --from markdown+autolink_bare_uris+tex_math_single_backslash --output /private/var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T/Rtmplj2Uy8/SangerAlignment/Achl_ACHLO006-09/Sanger_all_reads/SangerRead_Report_fasta.html --lua-filter /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rmarkdown/rmarkdown/lua/pagebreak.lua --lua-filter /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rmarkdown/rmarkdown/lua/latex-div.lua --embed-resources --standalone --variable bs3=TRUE --section-divs --table-of-contents --toc-depth 3 --variable toc_float=1 --variable toc_selectors=h1,h2,h3 --variable toc_collapsed=1 --variable toc_smooth_scroll=1 --variable toc_print=1 --template /Library/Frameworks/R.framework/Versions/4.2/Resources/library/rmarkdown/rmd/h/default.html --no-highlight --variable highlightjs=1 --variable theme=bootstrap --mathjax --variable 'mathjax-url=https://mathjax.rstudio.com/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML' --include-in-header /var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T//Rtmplj2Uy8/rmarkdown-str3926187706b9.html
#>
#> Output created: /private/var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T/RtmplI8nQX/SangerAlignment/SangerAlignment_Report.html
#> [1] "/private/var/folders/gf/3p_ynkts411bs238rtw3y0b40000gn/T/RtmplI8nQX/SangerAlignment/SangerAlignment_Report.html"

```

<sup>Created on 2023-03-02 with reprex v2.0.2</sup>

---

Example from the docs:

```
library(sangeranalyseR)
#> Loading required package: stringr
#> Loading required package: ape
#> Loading required package: Biostrings
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#>     anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#>     colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#>     get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#>     match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#>     Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#>     table, tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#> Loading required package: stats4
#>
#> Attaching package: 'S4Vectors'
#> The following objects are masked from 'package:base':
#>
#>     expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: XVector
#> Loading required package: GenomeInfoDb
#>
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:ape':
#>
#>     complement
#> The following object is masked from 'package:base':
#>
#>     strsplit
#> Loading required package: DECIPHER
#> Loading required package: RSQLite
#> Loading required package: parallel
#> Loading required package: reshape2
#> Loading required package: phangorn
#> Loading required package: sangerseqR
#> Loading required package: gridExtra
#>
#> Attaching package: 'gridExtra'
#> The following object is masked from 'package:BiocGenerics':
#>
#>     combine
#> Loading required package: shiny
#> Loading required package: shinydashboard
#>
#> Attaching package: 'shinydashboard'
#> The following object is masked from 'package:graphics':
#>
#>     box
#> Loading required package: shinyjs
#>
#> Attaching package: 'shinyjs'
#> The following object is masked from 'package:shiny':
#>
#>     runExample
#> The following object is masked from 'package:RSQLite':
#>
#>     show
#> The following object is masked from 'package:Biostrings':
#>
#>     show
#> The following object is masked from 'package:GenomeInfoDb':
#>
#>     show
#> The following object is masked from 'package:XVector':
#>
#>     show
#> The following object is masked from 'package:IRanges':
#>
#>     show
#> The following object is masked from 'package:S4Vectors':
#>
#>     show
#> The following object is masked from 'package:stats4':
#>
#>     show
#> The following objects are masked from 'package:methods':
#>
#>     removeClass, show
#> Loading required package: data.table
#>
#> Attaching package: 'data.table'
#> The following objects are masked from 'package:reshape2':
#>
#>     dcast, melt
#> The following object is masked from 'package:IRanges':
#>
#>     shift
#> The following objects are masked from 'package:S4Vectors':
#>
#>     first, second
#> Loading required package: plotly
#> Loading required package: ggplot2
#>
#> Attaching package: 'plotly'
#> The following object is masked from 'package:ggplot2':
#>
#>     last_plot
#> The following object is masked from 'package:phangorn':
#>
#>     add_boxplot
#> The following object is masked from 'package:XVector':
#>
#>     slice
#> The following object is masked from 'package:IRanges':
#>
#>     slice
#> The following object is masked from 'package:S4Vectors':
#>
#>     rename
#> The following object is masked from 'package:stats':
#>
#>     filter
#> The following object is masked from 'package:graphics':
#>
#>     layout
#> Loading required package: DT
#>
#> Attaching package: 'DT'
#> The following objects are masked from 'package:shiny':
#>
#>     dataTableOutput, renderDataTable
#> Loading required package: zeallot
#> Loading required package: excelR
#> Loading required package: shinycssloaders
#> Loading required package: ggdendro
#> Loading required package: shinyWidgets
#>
#> Attaching package: 'shinyWidgets'
#> The following object is masked from 'package:shinyjs':
#>
#>     alert
#> Loading required package: openxlsx
#> Loading required package: tools
#> Loading required package: rmarkdown
#> Loading required package: knitr
#> Loading required package: seqinr
#>
#> Attaching package: 'seqinr'
#> The following object is masked from 'package:shiny':
#>
#>     a
#> The following object is masked from 'package:sangerseqR':
#>
#>     read.abif
#> The following object is masked from 'package:Biostrings':
#>
#>     translate
#> The following objects are masked from 'package:ape':
#>
#>     as.alignment, consensus
#> Loading required package: BiocStyle
#>
#> Attaching package: 'BiocStyle'
#> The following objects are masked from 'package:rmarkdown':
#>
#>     html_document, md_document, pdf_document
#> The following object is masked from 'package:shiny':
#>
#>     markdown
#> Loading required package: logger
#> Welcome to sangeranalyseR
rawDataDir <- system.file("extdata", package = "sangeranalyseR")
parentDir <- file.path(rawDataDir, 'Allolobophora_chlorotica', 'ACHLO')

ACHLO_contigs <- SangerAlignment(ABIF_Directory     = parentDir,
                                 REGEX_SuffixForward = "_[0-9]*_F.ab1$",
                             REGEX_SuffixReverse = "_[0-9]*_R.ab1$")

launchApp(ACHLO_contigs)
writeFasta(ACHLO_contigs)
generateReport(ACHLO_contigs)

```

<sup>Created on 2023-03-02 with reprex v2.0.2</sup>

---

My sessionInfo output:

```
sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur ... 10.16
#>
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8
#>
#> attached base packages:
#>  [1] tools     parallel  stats4    stats     graphics  grDevices utils
#>  [8] datasets  methods   base
#>
#> other attached packages:
#>  [1] sangeranalyseR_1.8.0  logger_0.2.2          BiocStyle_2.26.0
#>  [4] seqinr_4.2-23         knitr_1.42            rmarkdown_2.20
#>  [7] openxlsx_4.2.5.2      shinyWidgets_0.7.6    ggdendro_0.1.23
#> [10] shinycssloaders_1.0.0 excelR_0.4.0          zeallot_0.1.0
#> [13] DT_0.27               plotly_4.10.1         ggplot2_3.4.1
#> [16] data.table_1.14.8     shinyjs_2.1.0         shinydashboard_0.7.2
#> [19] shiny_1.7.4           gridExtra_2.3         sangerseqR_1.34.0
#> [22] phangorn_2.11.1       reshape2_1.4.4        DECIPHER_2.26.0
#> [25] RSQLite_2.3.0         Biostrings_2.66.0     GenomeInfoDb_1.34.9
#> [28] XVector_0.38.0        IRanges_2.32.0        S4Vectors_0.36.2
#> [31] BiocGenerics_0.44.0   ape_5.7               stringr_1.5.0
#>
#> loaded via a namespace (and not attached):
#>  [1] nlme_3.1-162           bitops_1.0-7           fs_1.6.1
#>  [4] bit64_4.0.5            httr_1.4.5             R.cache_0.16.0
#>  [7] bslib_0.4.2            utf8_1.2.3             R6_2.5.1
#> [10] DBI_1.1.3              lazyeval_0.2.2         colorspace_2.1-0
#> [13] ade4_1.7-22            withr_2.5.0            tidyselect_1.2.0
#> [16] bit_4.0.5              compiler_4.2.2         cli_3.6.0
#> [19] sass_0.4.5             scales_1.2.1           quadprog_1.5-8
#> [22] digest_0.6.31          R.utils_2.12.2         pkgconfig_2.0.3
#> [25] htmltools_0.5.4        styler_1.9.0           fastmap_1.1.1
#> [28] htmlwidgets_1.6.1      rlang_1.0.6            rstudioapi_0.14
#> [31] jquerylib_0.1.4        generics_0.1.3         jsonlite_1.8.4
#> [34] zip_2.2.2              dplyr_1.1.0            R.oo_1.25.0
#> [37] RCurl_1.98-1.10        magrittr_2.0.3         GenomeInfoDbData_1.2.9
#> [40] Matrix_1.5-3           Rcpp_1.0.10            munsell_0.5.0
#> [43] fansi_1.0.4            lifecycle_1.0.3        R.methodsS3_1.8.2
#> [46] stringi_1.7.12         yaml_2.3.7             MASS_7.3-58.2
#> [49] zlibbioc_1.44.0        plyr_1.8.8             grid_4.2.2
#> [52] blob_1.2.3             promises_1.2.0.1       crayon_1.5.2
#> [55] lattice_0.20-45        pillar_1.8.1           igraph_1.4.1
#> [58] codetools_0.2-19       fastmatch_1.1-3        reprex_2.0.2
#> [61] glue_1.6.2             evaluate_0.20          BiocManager_1.30.20
#> [64] vctrs_0.5.2            httpuv_1.6.9           gtable_0.3.1
#> [67] purrr_1.0.1            tidyr_1.3.0            cachem_1.0.7
#> [70] xfun_0.37              mime_0.12              xtable_1.8-4
#> [73] later_1.3.0            viridisLite_0.4.1      tibble_3.1.8
#> [76] memoise_2.0.1          ellipsis_0.3.2

```

---
Tags: r, sequence-alignment, bioconductor
---


# Question
Title: How to align multiple contigs? To generate full length genome
I have MAG geneome sequence released in NCBI. In the released sequence there are 4915 contigs present. I want to get the full length genome sequence as FASTA file to proceed with further works like circular genome viewing and other molecular typing. Kindly let me know how to proceed in this regard. I am not much aware of bioinformatics so need even the basic information in this regard.

Thanks

# Answer
Simple answer is, it might not be possible to do it.

The large number of contigs suggests to me that a short-read sequencing method (such as Illumina sequencing) was used. Alone, short-read sequencing technologies cannot resolve repetitive sequences that are longer than the length of the reads. This is why you have a large number of contigs.\*

\*You might already have done this, but if you did not, it might be worth checking that your MAG is not actually a mixture of different genomes. This can also somewhat increase the contig number and, more importantly, can skew your interpretation if you assume it's only one organism. There are various ways to test this, but the simplest I know is the program CheckM. To run it you will need a way to use Unix/Linux (could be on a virtual machine). If you are not familiar with Unix, bash and conda, I would heartily suggest software carpentry lessons such as The Unix Shell (https://swcarpentry.github.io/shell-novice/) and possibly this conda lesson although it is in beta and I haven't tested it https://carpentries-incubator.github.io/introduction-to-conda-for-data-scientists/

I'm not familiar with the term 'molecular typing', but it looks like you're interested in what organisms / genes there are in your sample? You probably might not need a full genome for that.

Small probability of helping - but it might also be good to have a look at the sizes of your contigs. You might luck out and have a few large contigs (hundred thousand - million basepairs) which could be compared to existing database genomes.
> 1 votes

---
Tags: fasta, wgs
---


# Question
Title: To find amino acid substitution from fasta alignment
I have fasta multiple sequence alignment file of protein sequence. I want to look for mutation/substitution in amino acid between query with reference to one query. Is there any R package or any method I can go through.

# Answer
> 1 votes
User-friendly alternative:

Get the alignment into Jalview, set a colour scheme you like, and inspect alignment visually. You can move sequences up and down in the alignment with keyboard arrows, so you can focus on what you need.

Online, no-install-needed alternative: (same program as M\_'s method) blastn on NCBI website. Tick 'align two or more sequences'. Query in the top box, all the others in the bottom box. (Would not use this method if you had\> 20 seqences).

If you have many sequences, I reccomend going with M\_ 's advice instead.

# Answer
> 0 votes
Use local Blast

```
makeblastdb -in myproteins.fa -parse_seqids -blastdb_version 5  -out myfasta.db -title "protein db" -dbtype protein

blastn -query reference.fas -db myfasta.db -evalue 1e-6 -num_threads 4 -out blasthits.txt

```

`-db` made from the fasta alignment file (preferably without gaps) `-query` is the sequence your reference sequence

Output there will be a pairwise comparison for every sequence against your reference.

---
Tags: r, proteins, mutations
---


# Question
Title: eggNOG: DNA-DNA or protein-protein orthologue identification?
I'm trying to work out is how reliable/accurate my eggNOG metagenomic outputs are. Critical for this is how its identifying orthologues.

**Background** `eggNOG` is commonly used in metagenomics gene identification. Take a sample, randomly sequence all the DNA therein, then "annotate" via eggNOG. It identifies the putative gene function against lots of databases e.g. NCBI, KEGG, Brite, Go ...

eggNOG uses is own database for gene identification via 'seed orthology' each seed has its own KEGG ID.

**Question** Does the comparison between the query (metagenomics data) and the *seeds* take place using DNA-DNA comparisons? Alternatively, is it performing protein-protein comparisons? Algorithmically, by the easiest is DNA-DNA, however DNA-DNA would result in limitations of interpretation ('cause of complex evolutionary stuff). It could be using an optimised a protein database instead and returning the output as DNA.

# Answer
> 1 votes
I can answer this now. I've been 'under the hood' of eggNOG and its using "hmm" (hidden Markov model) very likely, almost certain, the Eddie (sp?) group's `hmmer` to cluster proteins associated. It represents the cluster members strictly via seed orthologue code, which then are aligned to form a 'raw' and 'trimmed' alignment. Each alignment is assigned a code which corresponds to its know function. This is a database of alignments. Once it gets a protein - protein hit against one of its orthologue (I suspect it will be using a nested system) it will then assign the query to that alignment - each of which is assigned to a protein function.

It representative database not an all encompassing database, i.e. sequences are being omitted.

Taxonomically it is simply declaring a the familial name. There will be an internal database mapping the `seed orthologue` onto the species, but that is not being publicly disclosed.

I make a better explanation than that if anyone wants it.

---

There is some description in the original paper

---
Tags: metagenome, kegg, go, eggnog
---


# Question
Title: Extracting all reads from bam file which match read IDs in another file
I have a long list of read IDs of interest to me in a file called read\_names.txt. it is simply in the format:

```
m54197_180831_211346/4981510/ccs
m54197_180831_211346/6226723/ccs
m54197_180831_211346/6619607/ccs
...

```

etc where these are the actual read names from a fastq file. I am then trying to find those in a bam file for which the fastq has been mapped, and I can accomplish this like:

```
for ID in `cat read_names.txt`
do
samtools view inbam.bam | grep $ID >> read_locs.sam
done

```

However, this method is obscenely slow because it is rerunning samtools view for every ID iteration (several hours now for 600 read IDs), and I was hoping to do this for several `read_names.txt` files.

I tried sort of flipping the script a bit and running samtools view first but it only returned the first read ID present in the file and stopped:

```
samtools view inbam.bam | for ID in `cat read_names.txt`; do grep $ID >> read_locs.txt; done

```

Am I missing something or is this the only way to accomplish this task?

# Answer
It is still slow but grep has a `-f` option to take in a file

`samtools view inbam.bam | grep -f read_names.txt > read_locs.txt`
> 9 votes

# Answer
`samtools` can do this natively too using `-N, --qname-file` CLI option:

```
samtools view -N read_names.txt in.bam > read_locs.sam

```

```
❯ samtools --version
samtools 1.16.1
Using htslib 1.16

```
> 2 votes

---
Tags: sam, samtools, reads
---


# Question
Title: How to extract a list of hydrogen bonds from a PDB?
I am analyzing the 7VF2 structure, from the RCSB-database. I noticed that some residues belonging to one of its chains interact with the other chains that are present within the assembly. My question is if there is a method to extract a list of interacting residues that are present in a PDB.

EDIT:

Just as an example, we have Lys 1616 belonging to ZC3H13 chain that interacts with Glu 621 of the protein virilizer homolog (the first chain)

Mol\* 3D visualizer can be used for major information in the 3D View in the link page I provided

# Answer
You can do that with GROMACS. You probably need to check the manual for this, but roughly these are the steps:

```
gmx pdb2gmx -f input.pdb -o output.gro -water spc
gmx solvate -cp output.gro -cs spc216.gro -o output_solv.gro -p output.top
gmx pdb2gmx -f output_solv.gro -o output_processed.gro -p output.top -ignh
gmx make_ndx -f output_processed.gro -o index.ndx
gmx hbond -f output_processed.gro -n index.ndx -s output.tpr -num hbond.xvg -hbn hbond.ndx

```

First, you convert the PDB file into a GROMACS topology file using the pdb2gmx command. Then you can solvate the system with water molecules and add hydrogen atoms. Then you cCreate an index file containing the groups of atoms you want to analyze. Then you can use the hbond command to calculate the hydrogen bonds.

I have not used GROMACS in a while. Hope that provides a starting point. With GROMACS you can even run a simulation and monitor the h-bonds over time.

You actually might need to run an energy minimziation first, as the PDB might not in an optimal state for the h-bonds detection.
> 1 votes

---
Tags: pdb, structural-biology
---


# Question
Title: Transcript vs Primary transcript on phytozome
Could someone help me understand what the difference between transcript and primary transcript on phytozme is? For example, this dataset of *A.thaliana* has "primary transcript CDS" vs CDS.

Off the top of my head primary transcript represents the initial strand of RNA made from DNA and transcripts represent genes post-processing of the initial strand.

The thing that is confusing me is, why is the transcript file larger than the primary transcript file? The transcript file has about 21k more headers than the "primary transcript". The only explanation I can think of is that due to alternate splicing you could end up with a bunch of isoforms that lend more headers to the transcript file. Is that it?

# Answer
> 1 votes
According to a team lead at Phytozome:

> Multi-exon genes can produce different transcripts via process of alternative splicing: inclusion/exclusion of particular exons in the spliced transcript.Designation of one of these transcripts as “primary” will vary based on the genome under consideration: the primary might be the one with the highest expression under “normal” conditions, it might be the one whose mutations was first studied, or it might simply be the transcript with the longest CDS.

> In all cases, Number\_of\_transcripts \>= Number\_of\_primary\_transcripts Number\_of\_primary\_transcripts == Number\_of\_genes

---
Tags: genomics, public-databases, transcriptome
---


title: Non-ACGTN alternate alleles in VCF

text:
# Question
Title: Non-ACGTN alternate alleles in VCF
I was trying to left normalize my vcf file using following command:

`bcftools norm -f reference.fasta -m -any input.vcf > output.norm.vcf`

However this gives the following error:

`Non-ACGTN alternate allele at 1:278165 .. VCF:'-3AGG'`

It looks like the variant caller returned all deletions and insertions like this:

```
1   5387285 .   C   +2TG
1   5387316 .   G   -8TGTGTATC

```

This file format was produced by SCCaller 2.0.0, there are several unconventional output fields produced by this caller but but the repository isn't maintained.

Is this an issue of the vcf format or does this means that there is a non-ACGTN character on position 1:278165 in the fasta reference file? I don't see anyting out of the ordinary at this position using IGV.

What would be the best way forward? I tried the following command to extract non-ACGTN allels but it outputs all alleles:

`bcftools view -i 'ALT !~ "^[ACGTN]+$"' input.vcf.gz > errors.vcf`

EDIT: including several lines of vcf file

CHROM   POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  SAMPLE1
1   89946   .   A   T   68  .   NS=1    GT:SO:AD:BI:GQ:PL   0/1:False:3,18:0.973:68:273,164,49,117
1   286158  .   A   +2TG    78  .   NS=1    GT:SO:AD:BI:GQ:PL   0/0:NA:11,10:0.982:78:154,97,175,456
1   900119  .   A   G   150 .   NS=1    GT:SO:AD:BI:GQ:PL   0/0:False:17,14:0.995:150:214,136,321,657
1   932613  .   G   -4TTTC  62  .   NS=1    GT:SO:AD:BI:GQ:PL   0/1:NA:4,18:0.998:62:274,165,103,173
1   1168310 .   T   +2CA    35  .   NS=1    GT:SO:AD:BI:GQ:PL   0/0:NA:18,5:0.985:35:77,55,91,712```

</code></pre>

# Answer
That isn't a valid ALT field according to the VCF standard:

> ALT — alternate base(s): Comma-separated list of alternate non-reference alleles. These alleles do not have to be called in any of the samples. Options are base Strings made up of the bases A,C,G,T,N (case insensitive) or the ‘\*’ symbol (allele missing due to overlapping deletion) or a MISSING value ‘.’ (no variant) or an angle-bracketed ID String (“”) or a breakend replacement string as described in Section 5.4. If there are no alternative alleles, then the MISSING value must be used. In other words, the ALT field must be a symbolic allele, or a breakend replacement string, or match the regular expression ^(\[ACGTNacgtn\]+|\*|.)$. Tools processing VCF files are not required to preserve case in the allele String, except for IDs, which are case sensitive. (String; no whitespace, commas, or angle-brackets are permitted in the ID String itself)

So yes, you can have things other than ACTG, `N` most commonly or `.` or `*`. And you can also have things like `<DEL>` and `<DUP>` for structural variants. But a `+` or a digit make no sense and are not allowed so this is not a valid VCF file.

So, filter them out. You can use something like this:

```
zcat input.vcf.gz | awk -F'\t' '/^#/ || $5 ~ /^[ACGTNacgtn]+|\*|\.)$/' |
    bgzip > fixed.vcf.gz

```

I would also report this as a bug to the developers of your variant caller.
> 1 votes

# Answer
I would personally use plink2 for this. Nice and simple and unerstandable by all.

```
plink2 --vcf in.vcf --snps-only just-acgt --export vcf --out out.vcf

```
> 0 votes

---
Tags: vcf, variant-calling, bcftools
---


# Question
Title: How to calculate cell type percentage in every sample
I have a Seurat object (metadata) with the single R samples consisting of cell types and sample types columns. I am trying to make a table that has a sample and percentage of cell types for each sample so that I can run a t-test across the samples that came out of 3 prime vs multi-omics.

I have a seuratobject meta data with singler\_labels (cell labels like chondrocytes, astrocytes, macrophage...), sample names (S0030, S0031, S0032.....), cell id/barcodes, and everything a seurat object has. I want to make a table with samples and percent of cell types per sample just to see what sample has what percent types of cells.

How do I do this?

Thank you!

# Answer
If your cell types and sample names are in separate metadata variables attached to the Seurat object, then you can use `table` to count up the pairings:

```
table(pbmc.seurat$singler_labels, pbmc.seurat$sample_name)

```

This will output a matrix containing these counts
> 1 votes

# Answer
The following package performs this type of analysis and can be directly used on a Seurat object:

paper: propeller: testing for differences in cell type proportions in single cell data

github: Speckle

You mention that you would like to do some statistics so keep in mind that this requires replicated measurements to handle the large variability in proportions between samples. Your metadata should contain:

***celltype, sample, group\_to\_test***

once the package is installed you can simply run the propeller function and set group to technology (in your post you mention you want to evaluate 3' vs multi-omics):

```
propel <- propeller(clusters = object@meta.data$celltype,
                    sample = object@meta.data$sample,
                    group = object@meta.data$technology)

```

This will return the statistics: t-test for two groups or ANOVA for multiple testing.

The following function will return the transformed proportions that can then be used in limma for example to test for more complex designs or simple sample-wise plotting using ggplot. For `transform` you will offcourse have to decide which transformation is applicable to your data.

```
getTransformedProps(clusters = object@meta.data$celltype,
                             sample = object@meta.data$sample, transform="asin")

```

There are a variety of other packages that allow for more complex analysis and a recent benchmark can be found here
> 1 votes

---
Tags: rna-seq, scrnaseq, seurat, single-cell, multi-omics
---


# Question
Title: Which non-matching DNA bases clash the least?
I am building primers for LAMP which I know will not match all my target sequences perfectly.

To some extent, having primers with multiple variants at variable sites (degenerate primers) and longer sequences (tagged degenerate primers) will help. However, I can't just increase degeneracy (variability) forever because then I will have no amplification efficiency.

Therefore, especially at the 5' end which is less important from an amplification perspective, I will have some calculated mismatches. So I have this puzzle:

If those positions of interest (closer to the 5' end) are likely to be any base with about the same frequency, which is either a) the base or b) the pair of bases, which are likely to clash the least?

My instinctive reaction is that for a single base, a pyrimidine would fit better, because even if it doesn't match, it's smaller? And for a pair of bases, probably a purine and a pyrimidine? I would appreciate your thoughts!

# Answer
I would just use Primer(3)-Blast here: https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?GROUP\_TARGET=on

I would avoid the approach you outline because:

1. Primer design algorithms incorporating degeneracy are very good;
2. I'm not sure you're right (but could be wrong).

The algorithm behind Primer-Blast is Primer3, and is legendary in the primer design business and of course Blast.

I agree that strong conservation at 3' is key to the design, but based on homology rather than chemistry.

The thermodynamic calculation ... I would suggest that a mismatch with AT is going to less severe than GC. GC is a stronger bond, but I'd let the algorithm sort this out and prioritise using personal judgement:

1. Conservation by homology at the 3' end.
2. Conservation by homology for GC regions.

Ultimately, primers are cheap so ordering loads of algorithmic designs is usually the best approach.

---

Just to add an explanation on the subject of why is this bioinformatics?

One of the issues of primer selection is self-hybridisation within the primer or annealing between the primer pair. This is bad news for a primer - because it might not work.

Algorithm primer design programs are good at spotting and avoiding unwanted hybridisation. When degeneracy is involved this becomes complicated.

What Primer-Blast does is assess all combinations of degeneracy against issues such as self-hybridisation and primer-annealing. Humanely it is not possible to assess all those combinations and that is why an algorithm is needed.
> 2 votes

---
Tags: sequence-analysis, primer
---


# Question
Title: Can you recommend me some programs or R packages for species distribution modelling?
I am working with a whole genus of land snails, trying to predict its distribution in the coming decades and how it might get affected by climate change. I have already worked with MaxENT and the R package Biomod2 (I have tried out all the important algorithms within it as well). I was wondering if any of you can recommend another option thats reliable.

# Answer
Spatial species distribution always follow a negative binomial distribution at any point of *most* other factors e.g. climate. When a species enters a human monoculture type environment - then deviation from the negative binomial can occur within that monoculture. The clumping will be approximated via trigonometry: 'cause the clumping can be approximated via a circle, very clumped - small circle, more dispersed bigger circle.

The parameters of the negative binomial are likely to change in relationship to climate and thats what I'd model if I was doing this work. The 'clumping' (which is what the negative binomial models) is likely to tightened as global warming increases. In practical terms the snails start drying out so they have to clump in the remaining pockets of moisture (this might be more true of slugs than snails ... at a guess). Alot of invertebrates that cannot or don't really fly - are prone to desiccation.

Personally, I go to empirical stats rather than use a package. However, the packages you are using will follow these two principles. So you can do it from first principles. Its fairly simple ecology its not like phylogenetics that have been through decades of development.

If you've DNA from the snails then the phylogenetic Bayesian spatial models around Beast provide a powerful basis for the investigation. That would be very strongly recommended if you've DNA sequence from a good proportion of the snails.
> 0 votes

---
Tags: modelling, ecology
---


# Question
Title: How to calculate frequency of a category with respect to the value in another column?
I was trying to calculate the frequency of disease\_present (`yes`) when smoking status is `y` (yes) for each group (A, B, C, D)

`The formula is = Number of samples with disease_present (yes) / No of samples with smoking status 'y' * 100`

**df**

```
Sample_code Group   smoking disease_present
1   A   y   yes
2   B   n
3   C   n   yes
4   D   y
5   B   n   yes
6   A   y   yes
7   C   n
8   D   y
9   A   n   yes
10  B   y
11  D   n   yes
12  C   n
13  A   y   yes
14  B   y   yes

```

I tried this

```
df.groupby("Group").count()

```

This gave the same number in every category

1. Pandas by default turns empty cells into NaN
2. I tried to replace NaN with empty strings with `df1 = df.fillna("")`. I still got the same results

# Answer
I am not able to reproduce the issue.

However I think this is a crosstab problem:

```
>>> new_df = df.fillna("no")
>>> new_df
# prints:
    Sample_code Group smoking disease_present
0             1     A       y             yes
1             2     B       n              no
2             3     C       n             yes
3             4     D       y              no
4             5     B       n             yes
5             6     A       y             yes
6             7     C       n              no
7             8     D       y              no
8             9     A       n             yes
9            10     B       y              no
10           11     D       n             yes
11           12     C       n              no
12           13     A       y             yes
13           14     B       y             yes

>>> pd.crosstab([new_df.disease_present, new_df.smoking], new_df.Group)
# prints
Group                    A  B  C  D
disease_present smoking
no              n        0  1  2  0
                y        0  1  0  2
yes             n        1  1  1  1
                y        3  1  0  0

```
> 2 votes

---
Tags: python, pandas
---


# Question
Title: Where can I find lncRNA expression data for different cell types?
Are there any publicly available databases providing expression data for long non-coding RNAs (lncRNAs) across cell types of multicellular organisms? Alternatively, are there lesser known UCSC tracks for this?

* **Example 1:** I want to compare tissue-specific expression of the human lncRNA *Xist* across HeLa, GM12878, and HEK293 cell types \[or more\].
* **Example 2:** I want to compare tissue-specific expression of the plant lncRNA *COLDAIR* across cortex, phloem, and trichoblast cell types \[or more\].

# Answer
As long as the lncRNA is annotated you can use https://gtexportal.org/home/gene.
> 2 votes

---
Tags: public-databases, ucsc
---


# Question
Title: In cellranger, can I change the minimum number of reads for autodetecting chemistry?
I ran cellranger multi on a TotalSeq B dataset and got this error

```
[error] Pipestance failed. Error log at:
sampleUnderscore2/SC_MULTI_CS/SC_MULTI_CORE/MULTI_CHEMISTRY_DETECTOR/_GEM_WELL_CHEMISTRY_DETECTOR/DETECT_COUNT_CHEMISTRY/fork0/chnk0-u0e6a05f224/_errors

Log message:
There were not enough reads to auto detect the chemistry: Sample CITE_2 in "/scratch/palmer/E222532_PR_Run01_051222/outs/fastq_path/AAAVFT3HV"
Note that you can avoid auto detection by specifying the specific chemistry type and version.
- Minimum number of required reads = 10000
- Number of reads available = 1408

```

I'm not sure what to specify the chemistry as. Does anyone know what it might be for CITEseq as part of TotalSeq B? Alternatively, can I reduce the minimum number of required reads?

Thank you

# Answer
> 2 votes
I don't think you can reduce the number of reads required for autodetection, and it would be a bad idea anyway, because the guess could go wrong.

From 10X Q&A:

> TotalSeq B is specifically designed to be compatible with the Single Cell 3’ v3 and 3.1 workflows

And from CellRanger support:

```
--chemistry     Optional. Assay configuration. NOTE: by default the assay configuration is detected automatically, which is the recommended mode. You should only specify chemistry if there is an error in automatic detection. Select one of:

    auto for auto-detection (default),
    threeprime for Single Cell 3′,
    fiveprime for Single Cell 5′,
    SC3Pv2 for Single Cell 3′ v2,
    SC3Pv3 for Single Cell 3′ v3,
    SC3Pv3LT for Single Cell 3′ v3 LT,
    SC3Pv3HT for Single Cell 3′ v3 HT,
    SC5P-PE for Single Cell 5′ paired-end (both R1 and R2 are used for alignment),
    SC5P-R2 for Single Cell 5′ R2-only (where only R2 is used for alignment).
    SC3Pv1 for Single Cell 3′ v1. NOTE: this mode cannot be auto-detected. It must be set explicitly with this option.
    ARC-v1 for analyzing the GEX portion of multiome data. NOTE: this mode cannot be auto-detected.

```

That leaves more than one possible chemistry option that could be applicable to your data. Check with the people who prepared the library if you don't know which library kit was used.

---
Tags: scrnaseq, cellranger, 10x
---


# Question
Title: Run nextflow pipeline multiple times depending on length of a list
*This question was also asked on Stack Overflow*

I have a Nextflow pipeline that runs a BUSCO search on input data, returns the output and then creates a report based on the output. One of my parameters is the lineage dataset. I would like to have my pipeline run multiple times (preferably at the same time) based on how long my list of lineage datasets is. For example, this is my nextflow.config file. The goal is that my pipeline would run twice in this case. If there are n lineages specified I would like it to run n times.

```
params {
    inputFilePath           = "$PWD/input_data/test_data.fasta"
    lineage_datasets        = [
        "fungi_odb10",
        "microsporidia_odb10"
        ]
    augustus_species        = "fumigatus"
    mode                    = "geno"

    threads                 = 16
    help                    = false
    outdir {
        main                = "$PWD/results"
    intermediate        = "$PWD/busco_outputs"
    }
}

```

The first part of my busco process looks like this:

```
  publishDir params.outdir.intermediate, mode: 'copy'

  input:
    path inputFile
    string lineage_dataset

  output:
    path 'hap1/*hap1.txt'

```

The idea is that each of the summary.txt files will be collected inside of my outdir.intermediate folder which is where my createReport process will target. For context, this is my createReport process:

```
process createReport {

  publishDir params.outdir.main, mode: 'copy'

  input:
    path 'busco_outputs/*'

  output:
    path 'report.html'

  script:
    """
    source "$launchDir/venv/bin/activate"
    report.py > report.html
    """
}

```

I have tried the following in my workflow process but it does not seem to work:

```
workflow {
  runBuscoChannel = Channel.fromPath(params.inputFilePath).cross(Channel.of(params.lineage_datasets))
  | mix
  | collect
  | set { outputs }

  createReport(outputs)
}

```

When I try to run the pipeline it does not run properly

Does anyone have any ideas of how I could change my file to allow for this? Thanks very much!

# Answer
> 2 votes
I think the answer to your question is the `each` input repeaters. You can also check Nextflow FAQ for the question on how to iterate over a process or a file.

Not in front of a PC now but you could try something like `each lineage_dataset`.

---
Tags: bash, nextflow
---


# Question
Title: Which human reference genome version does VARAdb use?
I am interested in annotating a list of variants using VARAdb. My coordinates use the GRCh38 (hg19) reference genome. However, I do not know which coordinate system VARAdb uses.

A section of the VARAdb publication reads: "GWAS have provided a large amount of data associating genetic variants with common phenotypes. We collected risk SNPs fromfive sources, including theNHGRI GWAS Catalog (29), GWASdb v2.0 (30), GAD (31), Johnson and O’Donnell (32), and GRASP v2.0 (33). Then, we filtered risk SNPs ofwhich ‘Variant ID’which are not ‘rsID’ and **transformed locations of risk SNPs from the GRCh38 assembly version into the hg19 version in the GWAS Catalog**. Finally, we obtained 1 515 001 risk SNPs associated with diseases, traits and phenotypes."

which makes me think that they may be using the hg19 genome, but I am not sure.

So, my question is:

Which human reference genome version does VARAdb use?

VARAdb database: http://www.licpathway.net/VARAdb/search.php

VARAdb publication: https://pubmed.ncbi.nlm.nih.gov/33095866/

# Answer
They are using hg19.

I found this by simply going to "Browse" then copying taking the first rsID in the list (rs1891805) and comparing the chromosomal position they give in VARAdb to what I see in dbSNP. VARAdb gives it as position 154904840 on chromosome 1 and dbSNP has:

```
Sequence name       Change
GRCh37.p13 chr 1    NC_000001.10:g.154904840T>C
GRCh38.p14 chr 1    NC_000001.11:g.154932364T>C
PMVK RefSeqGene     NG_053028.1:g.15287A>G

```

GRCh37 is hg19, and GRCh38 is hg38. So, as you can see above, the position 154904840 is from hg19 and the hg38 position for this rsID is 154932364.
> 3 votes

---
Tags: sequence-annotation, database, reference-genome
---


# Question
Title: Getting nearest gene from SNP data using SNP ID, CHR, and BP
I have data from a GWAS that provides me with the SNP, Chromosome, and base-pair. My data set has thousands of SNPs. What is the easiest way to find the nearest gene for each SNP using this information? My data looks like this:

# Answer
> 2 votes
This answer assumes that you have a terminal that runs linux commands and python, and that you can install and run bedtools.

**step 1:** convert your input spreadsheet to BED format. There is probably an awk one liner that does this but something like this python script will do it based on csv input (assuming that your coordinates are 1-indexed):

```
# writes BED format to standard output, capture it in input_file.bed
# not actually run
with open("input_file.csv") as infile:
  for line in infile:
    fields = line.split(",")
    coord = int(fields[2])
    outline = f"chr{fields[0]}\t{coord - 1}\t{coord}\t{fields[1]}"
    print(outline)

```

**step 2:** sort the BED.

```
sort -k1,1 -k2,2n input_file.bed > input_file.sorted.bed

```

**step 3:** Get genome annotations. I am assuming that you are using the grch38 build of the human genome.

```
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_43/gencode.v43.annotation.gff3.gz

```

**step 4:** Get nearest gene for each snp.

```
bedtools closest -a input_file.sorted.bed -b gencode.v43.annotation.gff3.gz > input_annotated_with_genes.bed

```

---
Tags: gene, snp, gwas
---


# Question
Title: Is it possible to create SMILES representation out of DNA/RNA-Oligonucleotides (~50 nts)?
I want to use physico chemical descriptors of single stranded DNA/RNA-Oligonucleotides in order to model certain properties. To generate these descriptors i need a computer understandable representation of these oligonucleotides, e.g. the SMILES representation. So far i could not find of any method to generate these SMILES representations out of the primary structure of these oligonucleotides. I am wondering if someone knows of a software, that takes the primary sequence of a single stranded DNA-Oligonucleotide with e.g. 50 nts as input and gives the SMILES representation as output. So far i could find nothing that would accomplish this task, although, at least to me, it seems like a straight forward task.

I will appreciate every constructive thought on that.

Thanks!

**EDIT:** Add some form of input and output file

As requested here an example of input and output pair:

input should be primary sequence: TCA

output should be the corresponding SMILES string: O=C(C(C)=CN1\[C@@\]2(\[H\])C(\[H\])(\[H\])C@@(\[H\])C@(\[H\])O2)NC1=O

This example was generated with chemdraw 20.1.1.

For my use case i need a software, that convert many (like 50k - 1000k) primary sequences of single stranded DNA/RNA with a length of about 50 nts into the corresponding SMILES string. Ideally would be some sort of python package but anything else would help also a lot!

# Answer
It looks like various chemical software toolkits can export SMILES, for example OpenBabel, using the babel command line tool. I am not readily seeing the best way to read e.g. a FASTA, though. Another option would be ChemDraw.

This researchgate thread has a lot of information on peptide-\>SMILES conversion, I don't see nucleotides though.

One option, though it would be a little finicky, would be to look up the nucleotides SMILES strings on PubChem and then write a script to string them together (with a monophosphate) based on an input sequence. It seems like there is some related tooling already in rdkit, but a convenient API doesn't seem to exist yet. The developer does have one suggestion for peptides, which would translate into something like this for nucleotides:

```
from rdkit import Chem

nt_mol = Chem.MolFromFASTA("ACATATAT", flavor=9)  # flavor param notes as DNA
Chem.MolToSMILES(nt_mol)

```

To do this across many nucleotides you could hook it up to biopython to iterate across a FASTA file.

The nucleotide `flavor` option exists, but there is a note that only peptide is supported. I don't know whether that means the feature is experimental or is not implemented.
> 3 votes

---
Tags: rna
---