# Question
Title: How can I find the maximum percent identity between two sets of (unaligned) sequences?
As the title states, I've got two sets of unaligned amino acid sequences (~25k sequences in one, ~3k in the other). I want to find the minimum distance between a sequence in the first set and a sequence in the second, expressed in terms of percent identity (so, maximum percent identity). So far I've only been able to find R packages that will compute percent identity for two already aligned sequences, which would then allow me to do a double loop over both sets of sequences to find the shortest distance between them. But this seems woefully inefficient - is there an easier/faster way? I also tried playing around with mmseqs2 but couldn't tell if there was a way to use it for anything other than clustering.

To clarify - there is a large amount of variation between the sequences in question, including in overall length (ranging from 148 residues to 872). They don't all belong to one family of homologs or anything. I'm expecting many of them to have quite low identity with each other. But the constraints of the problem I'm working on require me to put a number on that identity, or at least an upper bound.

I'm not constrained to using R, it's just what I'm most familiar with. Solutions with command line tools or python would also be fine - I'm sure there are better tools out there that I'm not aware of.

# Answer
> 1 votes
When I've been in 'possibly' comparable situations I've just used `local Blast` and parsed out the top hit using Python. It may not be the most optimal package (actually its not), but its easy and works, whilst NCBI's Blast support desk is efficient. Each query in turn is tested via a loop against the whole local database.

# Answer
> 1 votes
Try mmseqs2. It is largely a faster version of blastp. I haven't used mmseqs2 before but its documentation seems simple enough.

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


# Question
Title: Gene not found in Affymetrix expression profiles
I am studying the ABA network in A. thaliana, consisting of HB7, ABI1 and AREB2. The AGI code I was given are, respectively: AT2G46680, AT4G26080 and AT1G45249.

I downloaded the following file in order to convert the array element name to the AGI one, but I am not able to find AT1G45249 in the above file. How can I check whether AT1G45249 has another name? Or even better, how can I check what is the array element name in Affymetrix (GPL198) corresponding to AT1G45249. Thank you.

# Answer
> 1 votes
You can use `org.At.tair.db`: Genome wide annotation for Arabidopsis, which is an R package available on Bioconductor to extract gene symbols.

---
Tags: microarray, gene-expression, identifiers
---


# Question
Title: What is the most accurate approach for de Novo sequencing?
I'm trying to decide between PacBio HiFi or Illumina sequencing platforms for sequencing the genome of aChrysina scarab. We want to identify the color pattern loci and also perform a phylogenetic analysis as we'll be sequencing several species of that genus.

When considering a hybrid approach I was told by someone at a sequencing facility that adding short reads to polish the HiFi reads assembly lowers the quality of the assembly and introduces more errors that corrections. Why is this? I haven't been able to find this statement in the literature.

# Answer
HiFi reads can resolve most repeats in a genome. Illumina reads from these repetitive regions are likely to get mismapped and these mapping errors may become consensus errors. Most short-read polishing tools developed in the old days (e.g. pilon) will reduce the contig consensus quality. Furthermore, when you have a non-inbred diploid genome, HiFi reads can locally phase a large portion of the genome, while Illumina reads are likely to cancel/discard most HiFi phasing information, which would lead to a worse and less informative assembly.

In theory, if you use Illumina reads the right way, you can further improve HiFi contig consensus like what the telomere-to-telomere team is doing here. However, this pipeline is for haploid/inbred genomes only. There are no tools for generic genomes. Developing your own is beyond the capability of most people in the field. To this end, the person at your sequencing facility is mostly right: don't polish HiFi assemblies with Illumina reads.
> 4 votes

# Answer
Very interesting. Personally, I think the most important thing than anything I can input is it is worth asking about the HiFi quality index being considered and whether this was compatible with your biological objective.

Just to highlight that insect genomes are very different to genomics here - traditionally and even now - there has been a disproportionate focus on dipteran genomes. Okay certain Lepidoptera, namely the silk moth (Bombyx), is an advanced area of study. Beetles genomes are even more different still, agreed they holometabolous insects but still along way from dipterans. In addition, its worth pointing out beetles were late on the genome sequencing scene.

I don't know the precise biological objective - but I'm going to guess it. I'm assuming your wanting SNP variation to associate with phenotype colour variation, for the target genes (just guessing).

> HiFi reads

With regarding mixing the genome sequencing approaches: what is the quality scale being used? If its lowering the chromosomal assembly (traditionally a key difficulty in insects), or weakening contig assembly then ok I could imagine. This is very different DNA. If the quality index is relating to lowering SNP accuracy I would be surprised, for others to comment, but just can't see it.

I am aware that a key issue is that insect genome, specifically with reference to dipteran genomes, has caused extreme difficult assembling due to the high AT content and repeat regions, with chromosomal identification a traditiona nightmare and requiring wet lab chromosomal mapping (Diptera). I think beetles are an outlier even on this spectrum - against other insects. The advice you might have been given could be with this issue in mind, particularly if it was from another entomologist.

I don't know if this genome has been sequenced, the implication of the question (title) is that is hasn't been sequenced yet and which leads to the second issue,

JBS Haldane

> God has an inordinate fondness for stars and beetles.

There are an awful lot of beetles and a lot of those genomes not sequenced yet.

Personally I think you don't want accurate (really to maximise) genome assembly. In that scenario its got to be Illumina short-read, because this is a SNP study, where high accuracy calling from read depth is needed, given my guessed objective.

If the genome hasn't been sequenced there could be pressure to enhance genome assembly rather than target your biological objective. This could be why there has been mixed advice, between what's good for the community and what's needed for your biological objective, this being colour variation genes - for very pretty beetles.

I think, personally, its not necessarily a technical issue its either a win for others, or win for you. The only caveat is if the genes of interest are some extensive, physically linked gene family, then of course short-read may not be cool. I am not aware of that is a key feature of insect genomes.

---

In insect genomes there has been misplaced excitement ... new insect genome ... hang on thats just unordered contigs. This will happen using a short-read only approach.
> 0 votes

---
Tags: ngs, illumina, pacbio
---


# Question
Title: How to merge sam files together with adding read groups
I have three sequencing libraries of single individual mapped to a reference using `bwa-mem`. I would like to merge the three unsorted `.sam` files I have so, I can call variants and heterozygosity estimates using atlas. Atlas requires one input mapping file (`bam`) with defined read groups because it estimates the error profiles of different libraries separately.

How can I merge multiple sam files? Preferably avoiding java (Picard tools).

I tried to figure out a solution using `samtools 1.3`. I sorted individual files using `samtools sort`, then I used `samtools merge -r merged.bam s1.sort.sam s2.sort.sam s3.sort.sam` to merge the sorted files. However, the read group didn't make it to the header (and the variant caller I use is complaining about it), also the read group is stupidly the file name.

I tried to define meaningful readgroup names using procedure described on BioStars, but I found that this will just change the header, it does not adjust the names of read groups defined by `samtools merge` (the file names).

Following this related thread on SeqAnswers, I tried to define the correct header with read groups corresponding to merged file names:

```
samtools -rh rg.txt merged.bam s1.sort.sam s2.sort.sam s3.sort.sam

```

where `rg.txt` is a file with header

```
@RG s1.sort
@RG s2.sort
@RG s3.sort
...
output of samtools view -H s1.sort

```

However, the header still had not the read group, I guess because sam header accepts only tagged items to be specified (something like `@RG XY:s1.sort`). So, I looked into the merged `bam` and I found that the tag of RG is `Z:`. So I tried to just rename header of the merged file using `samtools reheader`, but then samtools complain about the fact that a tag needs to be of length 2:

```
Malformed key:value pair at line 123: "@RG    Z:s1.sort"
Segmentation fault (core dumped)

```

I have opened an issue to report this strange incompatibility of read groups generated by `samtools merge` with `samtools reheader`.

I would like solution to :

* create a bit more standardized read group names (`SM:Sample\tLB:library` format)
* avoid pointless writing to disk like in sam -\> sorted.sam -\> merged.bam case (can be probably achieved through "pipes and tees", thanks @bli)

I also know that I can specify RG to bwa, so the sam files will have read groups defined in the first place. But I do not like the idea of remapping three libraries just to create correct formatting of read groups.

# Answer
> 5 votes
I have used https://github.com/ekg/bamaddrg to add read groups quickly to multiple sam files. And then you can do a samtools merge of the tagged files.

# Answer
> 3 votes
**Not so elegant but working solution**

I found a solution that satisfy several of my conditions, basically I just have to assign read groups to individual mapping files, which can be just added to sorting process and then merge the sorted sam files with read groups.

```
function sort_and_assign_RG {
    # $1 input file
    # $2 read group id
    # $3 library

    OFILE=$(basename $1 .sam).sort.RG.sam
    HEADER=$(basename $1 .sam).header.sam

    # since I take header from unsorted file, so I need to add this information to header as well as the readgroup
    echo -e "@HD\tVN:1.3\tSO:coordinate" > $HEADER
    # add read group to header
    echo -e "@RG\tID:"$2"\tLB:""$3" >> $HEADER
    # and the rest of the header
    samtools view -H $1 >> $HEADER

    # now sort the input sam, remove header,
    # attach to each alignment read group ID (RGID)
    # and cat header and alignment together
    samtools sort -- $1 | samtools view - | \
      awk -v RGID="$2" '{ printf "%s\tRG:Z:%s\n", $0, RGID; }' | \
      cat $HEADER - > $OFILE

    # remove temp header
    rm $HEADER
}

sort_and_assign_RG s1.sam s1 is180
sort_and_assign_RG s2.sam s2 is350
sort_and_assign_RG s3.sam s3 is550

samtools merge merged.bam s{1,2,3}.sort.RG.sam

```

It would be nice to avoid creation of three `sam` files on the way, but at least I get standardized read groups in a merged bam tile.

*--- edit ---*

I also found where my confusion is coming from. The read group in header is in format

```
@RG    ID:foo    LB:lib_foo    ...`
@RG    ID:bar    LB:lib_bar    ...`

```

where `ID` is read group ID (not sample ID as I thought) and other tags are just specifications of the read group. Then individual mapping have read group assigned in different format

```
{aliment1}    RG:Z:foo
{aliment2}    RG:Z:bar

```

where `RG` is just a tag for read group and `Z` is just a mark that says that this tag is just a "printable string". Therefore I think that another solution would just to merge these sorted files and then just adding three **correct** lines to the header.

I have to say, `.sam` is a unbelievable intuitive format.

# Answer
> 0 votes
If the files already have defined read groups inside them, then doing a merge without the `-r` option should work:

```
samtools merge merged.bam s1.sort.sam s2.sort.sam s3.sort.sam

```

If there are multiple input files that share the same read group, then by default they will have random strings appended to make the read groups unique. This can be stopped by using the `-c` option, as mentioned in `man samtools merge`:

```
   -c      When  several  input files contain @RG headers with the same ID, emit only
           one of them (namely, the header line from the first file we find  that  ID
           in) to the merged output file.  Combining these similar headers is usually
           the right thing to do when the files being merged originated from the same
           file.

```

---
Tags: sam
---


# Question
Title: How can Iso-Seq reverse transcriptase artifacts be avoided?
My end goal is annotate a de-novo assembled genome. When trying to select the best method for transcriptome assembly I read Iso-seq was the preferred method. However other people suggested Nanopore as reverse transcriptase artifacts can be avoided. Is there any other way one can avoid these artifacts without the need for Nanopore? Also, what is the recommended depth? Is using more SMRTcells worth the cost?

# Answer
It is possible to filter out iso-seq reads with RT artifacts, though I don't know how effective these informatics approaches are. With nanopore, you will get more false junctions in alignment or transcript assembly due to sequencing errors. I more often see iso-seq used for new species.
> 2 votes

---
Tags: rna-seq, ngs, nanopore, pacbio
---


# Question
Title: Error of duplicated rownames although there are no duplicates
I have a data frame that I want to switch its row names from EnsembleID (GENEID) to gene symbol (SYMBOL). When I try to switch, I get this error:

```
Error in `.rowNamesDF<-`(x, value = value) :
  duplicate 'row.names' are not allowed
In addition: Warning message:
 Error in `.rowNamesDF<-`(x, value = value) :
  duplicate 'row.names' are not allowed

```

The code I'm using works for other similar data frames. Only this one is giving me trouble. There are no duplicates in the row names, I checked, so the source of this error is strange.

This is the code:

```
GetSymbols <- function(df) {   # Switch from geneID format to Symbol format

  ensembl.genes = as.vector(rownames(df))
  geneIDs = ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype = "GENEID", columns = c("SYMBOL","GENEID"))
  dupli = duplicated(geneIDs$GENEID)
  geneIDs = geneIDs[!dupli,]
  rownames(geneIDs) = geneIDs$GENEID
  rownames(df) = geneIDs$SYMBOL
  return(df)

}

#Calling the function:
 df = GetSymbols(df)

```

#### A sample of the dataframe:

```
structure(list(logFC = c(0.362414024114299, 0.270660335317552,
0.250197714648158, 0.355106754020497, 0.860436247608723, 0.331838240496886,
0.642000001979217, -0.506238094622956, 0.199130048525963, 0.226010606161764,
-0.300273572714138, 0.763447672404514, -0.282220097838676, 0.686284440644756,
0.209266959343866, 0.201216572039244, 0.467207891400503, 0.207476905368996,
0.548263917432251, 0.582095719576479, 0.167369145307915, -0.239886539805618,
-0.133088162158502, 0.432528267703316, -0.232843401197013), AveExpr = c(2.63663927914016,
6.23730774047258, 6.74130527553022, 6.51256660562955, 3.51386518575641,
0.962682352871946, 3.60110717340049, -1.07537043066235, 4.55672701732791,
3.96125159652725, 1.20473881401599, 1.72955742088622, 2.10495435034315,
2.95446094011986, 4.31205976523475, 4.85412008885978, 4.09171246258306,
3.94681319329962, 6.27193781389662, 0.936821016422818, 7.97270280545117,
2.4470079966345, 5.75991260844846, 3.84225682530045, 3.35515997418382
), t = c(5.15993300980101, 4.74644046948475, 4.67895886939685,
4.65270476519835, 4.39369295809775, 4.51083389626156, 4.20489542862464,
-4.49987927160777, 4.1096966271885, 4.03273843222089, -4.09143780518073,
4.07713970281633, -3.99157698713235, 3.97814493972727, 3.94402861982903,
3.92673590933521, 3.92436825669156, 3.90538427413924, 3.89364498564491,
3.93360332510695, 3.80233390398449, -3.76935782400013, -3.7013905266149,
3.68313224298669, -3.68114307289328), P.Value = c(4.50582660244125e-07,
3.21245260466636e-06, 4.37260426507432e-06, 4.92523952290176e-06,
1.54869545509172e-05, 9.2842399715271e-06, 3.45373717288344e-05,
9.74370585569957e-06, 5.12014322967072e-05, 7.00192255628614e-05,
5.51721193906694e-05, 5.84844087520547e-05, 8.26184479488647e-05,
8.71763749953859e-05, 9.9850615705057e-05, 0.000106924308803713,
0.000107929018751166, 0.000116314737750083, 0.000121805424073131,
0.000104060063751203, 0.000173710084683786, 0.000197140244520663,
0.000255164374909145, 0.000273300166932535, 0.000275347732934516
), adj.P.Val = c(0.00727510763230165, 0.0198807293341929, 0.0198807293341929,
0.0198807293341929, 0.0357217668827299, 0.0262203124576875, 0.06970505049172,
0.0262203124576875, 0.0858444785191522, 0.09421086799483, 0.0858444785191522,
0.0858444785191522, 0.0968123298197961, 0.0968123298197961, 0.0968123298197961,
0.0968123298197961, 0.0968123298197961, 0.098333518854239, 0.098333518854239,
0.0968123298197961, 0.130401412989609, 0.138392451653505, 0.150726778618736,
0.150726778618736, 0.150726778618736), B = c(5.59649208873277,
4.23844568248881, 3.95926511142815, 3.85011503575047, 2.70095619489428,
2.52151007226572, 2.02920447371096, 1.82077681362955, 1.73104697407381,
1.44381690648743, 1.32264026808754, 1.26349133643967, 1.17476645994761,
1.17423988785926, 1.13844103561967, 1.07813935334313, 1.06670248721183,
1.0004861785984, 0.949862540556422, 0.634079467108277, 0.620590447127368,
0.489071423128816, 0.29054426889472, 0.254833988789798, 0.246174538031329
)), row.names = c("ENSG00000188938", "ENSG00000182768", "ENSG00000174579",
"ENSG00000119523", "ENSG00000271798", "ENSG00000204923", "ENSG00000281910",
"ENSG00000271092", "ENSG00000153767", "ENSG00000116903", "ENSG00000168101",
"ENSG00000206834", "ENSG00000180035", "ENSG00000251733", "ENSG00000058729",
"ENSG00000132467", "ENSG00000212464", "ENSG00000115561", "ENSG00000230669",
"ENSG00000206612", "ENSG00000275272", "ENSG00000162139", "ENSG00000143889",
"ENSG00000207233", "ENSG00000119878"), class = "data.frame")

```

# Answer
> 4 votes
I found an answer shortly after posting this question:

```
ensembl.genes <- rownames(df)
geneIDs <- ensembldb::select(EnsDb.Hsapiens.v86, keys= ensembl.genes, keytype = "GENEID", columns = c("SYMBOL","GENEID"))
dupli = duplicated(geneIDs$SYMBOL)
geneIDs = geneIDs[!dupli,]
ind = intersect(rownames(df), geneIDs$GENEID)
df = df[ind,]
rownames(df) = geneIDs$SYMBOL

```

---
Tags: r, gene, ensembl
---


title: Cellranger results have too many cells

text:
# Question
Title: Cellranger results have too many cells
I have the results of cellranger analysis for single-nucleus RNA seq data that was done by someone else. So I do not know which parameters were used for cellranger. In the results, there are 77000 estimated cells (which is way too high) and around 800 median reads per cell. Those are filtered counts. Raw `feature_bc_matrix` has almost 2 million cells.

How to tune cellranger count to get lower number of cells?

# Answer
> 2 votes
The Barcode Rank plot shows very high background in your sample.
The cell estimation from `Cell Ranger` will look for a sudden drop in the number of UMIs and use this to call cells. In your case, the drop is too much to the right because the empty droplet plateau is almost completely merged with the real cells. If you have already looked at CellBender your case best corresponds to the rightmost example image.

With high background of free mRNA, even empty droplets will contain a significant number of mRNA molecules and become hard to distinguish from real cells.

This likely is a consequence of the single nulcei preparation. To isolate the nuclei, the cells necessarily need to be broken and all the cytoplasmic mRNA is freed contributing to the background. You could try removing this background in future experiments by increasing the number of washing steps after cell lysis.

Using software such as CellBender, SoupX or DropletUtils may help with identifying and reducing the background. CellBender also will rerun the cell estimation step, automatically adjusting the number of called cells.

In any case you can manually reduce the number of cells by filtering based on the minimum number of UMIs or genes. However, the higher the background the more arbitrary such a thresholding strategy will become.

---
Tags: scrnaseq, single-cell, 10x-genomics, cellranger
---


# Question
Title: Are the doublets in the tutorial dataset in Scanpy/ Seurat already filtered?
I noticed the tutorials that Scanpy and Seurat use do not demonstrate doublet removal in their down stream analysis. Is the dataset output of cellranger count already doublet removed or do I need to incorporate doublet removal before doing quality control etc?

# Answer
No in general cellranger does not remove doublets. You can apply a tool like Scrublet in scanpy or Doubletfinder in Seurat to infer doublets based on a mixed-transcriptional signal.

I say in general because cellranger does detect certain types of doublets in two contexts:

1. Multi-species experiments, e.g. the so-called barnyard experiment, cellranger detects cross-species doublets
2. With CellPlex sample multiplexing, cellranger detects cross-tag (cross-CMO) doublets
> 2 votes

---
Tags: scrnaseq, transcriptome, 10x-genomics
---


title: RNA strings and amino acid string

text:
# Question
Title: RNA strings and amino acid string
Which of the following RNA strings could translate into the amino acid string PRTEIN? (Select all that apply.)

I am not sure if the following mass of amino acid is relevant to the question.

```
CCAAGUACAGAGAUUAAC

CCAAGAACAGAUAUCAAU

CCACGUACUGAAAUUAAC

CCGAGGACCGAAAUCAAC

```

# Answer
> 3 votes
I don't quite get the question but here's the answer and the code

```
#!/usr/bin/env python3
from Bio.Seq import Seq

for dna in ['CCAAGUACAGAGAUUAAC','CCAAGAACAGAUAUCAAU', 'CCACGUACUGAAAUUAAC','CCGAGGACCGAAAUCAAC']:
    print(Seq(dna).translate())

```

> PSTEIN PRTDIN PRTEIN PRTEIN

# Answer
> 0 votes
Here's how you can do the translation yourself, using the genetic code in form of a dictionary from Biopython:

```
import Bio

# Load the genetic code (a mapping from RNA triplet to amino acid)
d = Bio.Data.CodonTable.standard_rna_table.forward_table
print('Genetic code:', d)
print()

# Use it to translate an RNA sequence into an amino acid sequence
s = 'CCAAGUACAGAGAUUAAC'
for triplet in (s[i:i+3] for i in range(0, len(s), 3)):
    print(triplet, '->', d[triplet])

```

Printed result:

```
Genetic code: {'UUU': 'F', 'UUC': 'F', 'UUA': 'L', 'UUG': 'L', 'UCU': 'S', 'UCC': 'S', 'UCA': 'S', 'UCG': 'S', 'UAU': 'Y', 'UAC': 'Y', 'UGU': 'C', 'UGC': 'C', 'UGG': 'W', 'CUU': 'L', 'CUC': 'L', 'CUA': 'L', 'CUG': 'L', 'CCU': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P', 'CAU': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q', 'CGU': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'AUU': 'I', 'AUC': 'I', 'AUA': 'I', 'AUG': 'M', 'ACU': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T', 'AAU': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K', 'AGU': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R', 'GUU': 'V', 'GUC': 'V', 'GUA': 'V', 'GUG': 'V', 'GCU': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A', 'GAU': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E', 'GGU': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G'}

CCA -> P
AGU -> S
ACA -> T
GAG -> E
AUU -> I
AAC -> N

```

Note that the same RNA sequence can also be translated into another sequence of amino acids by shifting the reading frame 1 or 2 positions. For example, the RNA sequence `'CCAAGUACAGAGAUUAAC'` can be split into different triplets:

1. `'CCA', 'AGU', 'ACA', ...`
2. `'CAA', 'GUA', 'CAG', ...`
3. `'AAG', 'UAC', 'AGA', ...`

In biological systems, an mRNA would have a start codon that marks the beginning of an open reading frame, which is often the triplet `AUG`.

---
Tags: python, phylogenetics, phylogeny, rna
---


# Question
Title: Sequence format: grouping sequences within a flat file
Is there a sequence format allowing multiple sequences to be grouped across a flat file? The specific application is stacking separate nucleotide sequences against a given RNA secondary structure.

To rephrase the same question I'm looking to stack viral nucleotide data in different groups across a single flat file. FastaQ provides an example:

Traditional fastQ,

```
@sequence_name
AGCTGAT
+
...(quality score)

@sequence2
AGCTGAC
+
...(quality score)

```

What I'm looking for is a format,

```
@Gp1_virus1-2-3
AGACGTAG
+
AGAGACAA
+
AGAGACAT

@Gp2_virus1-2
AGATGTAG
+
AGATACAA

```

In the above example there are two groups comprising 3 (group 1) and 2 (group 2) individual sequences (viruses). Here 'group' represents different structures.

Stockholme format has potential. My understanding of this format, which is basic, is it stacks a sequence, together with structural information. What is sought is stacking multiple sequences and then a structure format. Any clarification would help.

There are multiple XML formats, which may be relevant.

Just to mention a 'flat file' is needed to allow manual inspection. Dumping JSON is a solution, but isn't friendly for optical inspection.

---

**Application** The immediate downstream application is structure/sequence viewers. I think at least one. viewer will accept Stockholme 1.0 format. rFAM is very similar (@terdon's comment)

**Summarised Question** The question could be rephrased as to whether multiple RNA families can be stacked inside one Stockholme 1.0 or else rFAM format?

Some of the members here are likely familiar with **pFAM**. So could multiple protein structures be stored and retrieved inside one **pFAM** file.

@gringer's contribution is very welcome.

# Answer
An unmapped SAM file format allows for read groups. You can use `samtools import` to create unmapped SAM files from FASTQ files:

```
samtools import -R virus-1-2 in12.fq > 12_reads.sam
samtools import -R virus-1-2-3 in123.fq > 123_reads.sam
samtools merge out.sam 12_reads.sam 123_reads.sam

```
> 2 votes

# Answer
I'm awarding 'solved' to @gringer, because thats the context of the question.

The best solution I could find is stacking `Stockholm (sth)` formats together. In this example the RNA families are associated with one structure:

```
# STOCKHOLM 1.0
#=GF ID    UPSK
AF035635.1/619-641             UGAGUUCUCGAUCUCUAAAAUCG
M24804.1/82-104                UGAGUUCUCUAUCUCUAAAAUCG
J04373.1/6212-6234             UAAGUUCUCGAUCUUUAAAAUCG
M24803.1/1-23                  UAAGUUCUCGAUCUCUAAAAUCG
#=GC SS_cons                   .AAA....<<<<aaa....>>>>
//

```

`Biopython` allows `sth` format to be input.

```
from Bio import AlignIO
from pathlib import Path

file = 'test.sth'
path = '/Volumes/mypath/testSth'
alignment = AlignIO.read(open(Path(path, file)), "stockholm")
print("Alignment length %i" % alignment.get_alignment_length())
print (alignment.column_annotations['secondary_structure'])

```

Output,

```
Alignment length 23
.AAA....<<<<aaa....>>>> # this is the RNA structure

```

If we look the `alignment` object via

```
from pprint import pprint
pprint(vars(alignment)) # or via print(alignment.__dict__)

```

The output

> {'\_annotations': {}, '**\_per\_col\_annotations**': {**'secondary\_structure'**: '.AAA....\<\<\<\<aaa....\>\>\>\>'},
> '\_records': \[SeqRecord(seq=Seq('UGAGUUCUCGAUCUCUAAAAUCG'), id='AF035635.1/619-641', name='AF035635.1', description='AF035635.1/619-641', dbxrefs=\[\]), SeqRecord(seq=Seq('UGAGUUCUCUAUCUCUAAAAUCG'), id='M24804.1/82-104', name='M24804.1', description='M24804.1/82-104', dbxrefs=\[\]), SeqRecord(seq=Seq('UAAGUUCUCGAUCUUUAAAAUCG'), id='J04373.1/6212-6234', name='J04373.1', description='J04373.1/6212-6234', dbxrefs=\[\]), SeqRecord(seq=Seq('UAAGUUCUCGAUCUCUAAAAUCG'), id='M24803.1/1-23', name='M24803.1', description='M24803.1/1-23', dbxrefs=\[\])\], 'annotations': {}} {**'secondary\_structure'**: '.AAA....\<\<\<\<aaa....\>\>\>\>'}

The reason `vars(alignment)` is required is to internally construct an alignment object in Stockholm format, here's the first part of the skeleton 'pseudo-code',

```
from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord

outAlign = MultipleSeqAlignment([])
record in alignment:
    name = record.id
    myseq = record.seq
    sr = SeqRecord(Seq(name), myseq, '', '')
    outAlign.append(sr)

```

The above code constructs a basic alignment, the bit thats missing is to add the following key/value to the `alignment` object of `{'secondary_structure': '.AAA....<<<<aaa....>>>>'}`.

Its complicated, but you get use to it.

---

https://github.com/biopython/biopython/blob/master/Bio/AlignIO/StockholmIO.py
> 0 votes

---
Tags: fastq, file-formats, format-conversion, secondary-structure
---


# Question
Title: Lower case vs. upper case nucleotids in sequence vs. dots at the end
What is the difference between lower case and upper case nucleotides in a sequence? My other question is what are the dots at the end of the sequence? Some examples are shown below:

* GGgG,GGgG,GGgG,GGgG,GGGG,GGGG…
* GGGG,GGGG,GGGG,GGGG,GGGG,GGGG…
* GGGG,gggg,GGGG,GGGG,UUUU
* ggGg,ggGg,ggGg,GGgG,GGgG,GGgG…

These sequences are real experimental samples. I obtained them from QuadBase. Here is the link below:

https://quadbase.igib.res.in/

Hopefully you will provide some literature within the thread in order to understand this concept in bioinformatics. BTW I'm a beginner in this field and hopefully this Bioinformatics Stack Exchange site will help me find out the answers.

# Answer
* GGGG,gggg,GGGG,GGGG,UUUU

Quadbase appears to be about DNA secondary structures, albeit the above example is RNA. The website states,

> G-Quadruplexes or G4s are secondary DNA structures that are now established to have functionally crucial roles in various cellular processes.

Two guesses, firstly its all just a bit weird, secondly its a representation of hypoxanthine. Its more likely to be the second guess than the first.

1. The nomenclature is a bit weird: gggg with UUUU would infer non-canonical binding/bounding in WUSS (Washington University Secondary Structure) format for a psuedoknot (type of secondary structure). If that wasn't the interpretation it would be kinda weird to share nomenclature with a major secondary structure format (WUSS) for secondary structure notation.

The 'dots', i.e. '.' appear to represent continuation, although in dot-bracket format (secondary structure format) this represents a nucleotide that doesn't hybridise.

---

2. Its worth reading the paper for hypoxanthine representation (normally its 'I'), because 'I' will behave like a guanine residue and would then make sense if it was being represented as 'g'. Thus, the following means no hypoxanthine,

* GGGG,GGGG,GGGG,GGGG,GGGG,GGGG…

My bet is its point 2, otherwise it doesn't make sense.
> 0 votes

---
Tags: sequence-alignment, genome, ngs, database
---


# Question
Title: Read PDB file, extract dihedral angles, modify dihedral angles, reconstruct Cartesian coordinates, and write PDB file
As the title summarizes, I am trying to:

1. Read a PDB file (for example, `1enh.pdb`).
2. Extract the backbone dihedral angles (phi, psi).
3. Modify the dihedral angles (phi, psi) arbitrarily (for example, replace them with `np.random.uniform(-np.pi, np.pi, size = (56, 2))` in 1enh, since there are 56 pairs of dihedrals).
4. Reconstruct the Cartesian coordinates of the modified backbone.
5. Write a PDB file based on the modified backbone.

From this question it seems the tools for what I want are available in BioPython. There is also a `set_angle()` method. But I only have a very limited knowledge of Python and BioPython. I am stuck at Steps 3 and 4, and dubious on Step 5.

This is what I have tried so far:

```
#!/usr/bin/python

import numpy as np
from Bio.PDB import PDBParser, PICIO, PDBIO

## Step 1

parser = PDBParser()
structure = parser.get_structure("", "1enh.pdb")

## Step 2

structure.atom_to_internal_coordinates()
chain = list(structure.get_chains())[0]
ic_chain = chain.internal_coord
dihedrals = ic_chain.dihedra

## Step 3 - ??

# Modifiy (phi,psi). For example, replace (phi,psi) with
newdihedrals = np.random.uniform(-np.pi, np.pi, size = (56, 2))

## Step 4 - ??

modified = PICIO.read_PIC('internal')
modified.internal_to_atom_coordinates()

## Step 5

io = PDBIO()
io.set_structure(modified)
io.save('1enh_modified.pdb', preserve_atom_numbering = True)

```

Many thanks in advance.

**EDIT:** I do not care about possible clashes if the side-chain is included in the modified PDB. In any case, writing a modified PDB file that only includes the backbone is also fine.

# Answer
ok, here my input pdb, `example_short.pdb` :

```
ATOM      1  N   MET A   1      14.067  31.664  -1.639  1.00 80.87           N
ATOM      2  HN1 MET A   1      14.250  31.259  -0.699  1.00  0.00           H
ATOM      3  HN2 MET A   1      14.891  32.220  -1.945  1.00  0.00           H
ATOM      4  HN3 MET A   1      13.229  32.279  -1.596  1.00  0.00           H
ATOM      5  CA  MET A   1      13.824  30.564  -2.612  1.00 79.17           C
ATOM      6  HA  MET A   1      14.702  29.921  -2.671  1.00  0.00           H
ATOM      7  C   MET A   1      13.556  31.158  -3.989  1.00 78.52           C
ATOM      8  O   MET A   1      13.301  32.365  -4.085  1.00 79.40           O
ATOM      9  CB  MET A   1      12.634  29.724  -2.148  1.00 83.38           C
ATOM     10  HB1 MET A   1      12.600  29.755  -1.059  1.00  0.00           H
ATOM     11  HB2 MET A   1      11.725  30.170  -2.552  1.00  0.00           H
ATOM     12  CG  MET A   1      12.680  28.264  -2.580  1.00 85.71           C
ATOM     13  HG1 MET A   1      12.066  28.145  -3.473  1.00  0.00           H
ATOM     14  HG2 MET A   1      12.274  27.650  -1.776  1.00  0.00           H
ATOM     15  SD  MET A   1      14.354  27.697  -2.951  1.00100.24           S
ATOM     16  CE  MET A   1      15.137  27.675  -1.338  1.00 97.10           C
ATOM     17  HE1 MET A   1      15.118  28.678  -0.913  1.00  0.00           H
ATOM     18  HE2 MET A   1      14.597  26.991  -0.683  1.00  0.00           H
ATOM     19  HE3 MET A   1      16.170  27.342  -1.441  1.00  0.00           H
ATOM     20  N   GLU A   2      13.604  30.306  -5.018  1.00 75.32           N
ATOM     21  H   GLU A   2      13.628  29.291  -4.794  1.00  0.00           H
ATOM     22  CA  GLU A   2      13.626  30.684  -6.444  1.00 72.66           C
ATOM     23  HA  GLU A   2      13.508  29.732  -6.961  1.00  0.00           H
ATOM     24  C   GLU A   2      14.903  31.384  -6.843  1.00 65.44           C
ATOM     25  O   GLU A   2      15.489  31.127  -7.892  1.00 67.46           O
ATOM     26  CB  GLU A   2      12.516  31.672  -6.810  1.00 72.71           C
ATOM     27  HB1 GLU A   2      12.094  32.068  -5.886  1.00  0.00           H
ATOM     28  HB2 GLU A   2      12.960  32.486  -7.383  1.00  0.00           H
ATOM     29  CG  GLU A   2      11.397  31.081  -7.622  1.00 77.72           C
ATOM     30  HG1 GLU A   2      10.901  31.877  -8.178  1.00  0.00           H
ATOM     31  HG2 GLU A   2      11.809  30.354  -8.321  1.00  0.00           H
ATOM     32  CD  GLU A   2      10.385  30.393  -6.731  1.00 79.89           C
ATOM     33  OE1 GLU A   2       9.171  30.574  -6.987  1.00 82.69           O
ATOM     34  OE2 GLU A   2      10.824  29.702  -5.779  1.00 72.73           O
ATOM     35  N   ASN A   3      15.194  32.399  -6.043  1.00 60.12           N
ATOM     36  H   ASN A   3      14.712  32.444  -5.123  1.00  0.00           H
ATOM     37  CA  ASN A   3      16.140  33.452  -6.363  1.00 56.21           C
ATOM     38  HA  ASN A   3      16.072  33.715  -7.419  1.00  0.00           H
ATOM     39  C   ASN A   3      17.548  32.968  -6.047  1.00 50.40           C
ATOM     40  O   ASN A   3      18.513  33.666  -6.331  1.00 50.95           O
ATOM     41  CB  ASN A   3      15.821  34.667  -5.493  1.00 57.10           C
ATOM     42  HB1 ASN A   3      15.913  34.381  -4.445  1.00  0.00           H
ATOM     43  HB2 ASN A   3      16.538  35.457  -5.717  1.00  0.00           H
ATOM     44  CG  ASN A   3      14.425  35.196  -5.732  1.00 62.07           C
ATOM     45  OD1 ASN A   3      14.220  36.001  -6.641  1.00 69.12           O
ATOM     46  ND2 ASN A   3      13.464  34.748  -4.923  1.00 57.21           N
ATOM     47 1HD2 ASN A   3      12.486  35.079  -5.045  1.00  0.00           H
ATOM     48 2HD2 ASN A   3      13.693  34.068  -4.170  1.00  0.00           H
ATOM     49  N   PHE A   4      17.643  31.764  -5.492  1.00 47.14           N
ATOM     50  H   PHE A   4      16.781  31.193  -5.378  1.00  0.00           H
ATOM     51  CA  PHE A   4      18.914  31.219  -5.038  1.00 46.77           C
ATOM     52  HA  PHE A   4      19.696  31.915  -5.340  1.00  0.00           H
ATOM     53  C   PHE A   4      19.239  29.888  -5.695  1.00 48.73           C
ATOM     54  O   PHE A   4      18.490  28.919  -5.594  1.00 49.50           O
ATOM     55  CB  PHE A   4      18.922  31.095  -3.507  1.00 42.32           C
ATOM     56  HB1 PHE A   4      18.091  30.456  -3.207  1.00  0.00           H
ATOM     57  HB2 PHE A   4      19.862  30.635  -3.201  1.00  0.00           H
ATOM     58  CG  PHE A   4      18.787  32.415  -2.806  1.00 39.42           C
ATOM     59  CD1 PHE A   4      17.545  32.878  -2.408  1.00 40.76           C
ATOM     60  HD1 PHE A   4      16.673  32.232  -2.509  1.00  0.00           H
ATOM     61  CD2 PHE A   4      19.881  33.258  -2.648  1.00 35.13           C
ATOM     62  HD2 PHE A   4      20.857  32.941  -3.015  1.00  0.00           H
ATOM     63  CE1 PHE A   4      17.397  34.138  -1.889  1.00 36.49           C
ATOM     64  HE1 PHE A   4      16.400  34.520  -1.669  1.00  0.00           H
ATOM     65  CE2 PHE A   4      19.757  34.474  -2.043  1.00 32.88           C
ATOM     66  HE2 PHE A   4      20.640  35.090  -1.873  1.00  0.00           H
ATOM     67  CZ  PHE A   4      18.507  34.929  -1.643  1.00 43.95           C
ATOM     68  HZ  PHE A   4      18.401  35.892  -1.144  1.00  0.00           H
END

```

here my code :

```
from Bio import PDB

from Bio.PDB import PICIO, PDBIO

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser(PERMISSIVE=1, QUIET=1)

structure: PDB.Structure.Structure = parser.get_structure("prova", "example_short.pdb")

structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
d:  Dict[Tuple[PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey],
                 PDB.internal_coords.Dihedron] = ic_chain.dihedra

cnt = 1
for key in d:

    if key[0].akl[3] == 'N':
        if key[1].akl[3] == 'CA':
            if key[2].akl[3] == 'C':
                if key[3].akl[3] == 'N':

                    print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

                    d[key].angle += 45

                    cnt += 1

cnt = 1
for key in d:

    if key[0].akl[3] == 'N':
        if key[1].akl[3] == 'CA':
            if key[2].akl[3] == 'C':
                if key[3].akl[3] == 'N':

                    print ('\n',cnt,' :   ',  [x.akl[3] for x in key], d[key].angle)

                    cnt += 1

structure.internal_to_atom_coordinates(verbose = True)

io = PDBIO()

io.set_structure(structure)

io.save('atom_coord.pdb',  preserve_atom_numbering=True)

```

here an image of the original vs new quadripeptid (cartoon repr, just backbone) :

original is rainbowed one (blue to red), magenta is the new one

saved as `atom_coord.pdb`

I believe I changed φ (Phi) adding 45° since

> Now, a φ angle is N-CA-C-N

It is worth mentioning that changing

```
 `d[key].angle += 45` to `d[key].angle += 6645`

```

works too. Apparently, need confirmation from someone more knowledgeable,

```
`structure.internal_to_atom_coordinates`

```

can cope with more than 360° angle values.
> 5 votes

# Answer
I do not know if there is a function in Biopython.PDB for dihedral space transformations, but due to the lack of answers I would say no. Biopython.PDB is nice and clean, but I would use something with more molecular mechanics power... or should I say potential (:drum::laughing:).

Personally, I use PyRosetta, but it is not open source (code available to academics), so others may opt for OpenForceField. PyRosetta requires a licence (free for academia), which is different from the Rosetta proper licence.

The code below uses a helper module I wrote, so I'll use the installer therein (which is shorthand for the official pip wheel, which provides different options etc., while the official shorthand is actually `pyrosettacolabsetup`)

```
pip3 install pyrosetta-help
install_pyrosetta -u  -p

```

## Ligands?

First lets make sure there's nothing interesting ligand wise.

NB. RCSB has a new shiny API which may be better... Not tried it: currently using European.

```
import requests

def retrieve_nonpolymer_data(entry_id: str) -> dict:
    """cf. https://www.ebi.ac.uk/pdbe/api/doc/"""
    response: requests.Response = requests.get(f' https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{entry_id.lower()}')
    response.raise_for_status()
    return response.json()[entry_id.lower()]

info = retrieve_nonpolymer_data('1enh')
chem_resns = [c for entity in info for c in entity.get('chem_comp_ids', [])]
print(chem_resns)  # ['HOH']

```

If is were something that is not water, I made a list here, which is based off Peter Curran's repo for the HotSpots algorithm, which allows the elimination of worthless ligands like glycol.

I mention this because the route differs if there is a ligand.

## Init and read

```
import pyrosetta
import pyrosetta_help as ph
from types import ModuleType

logger = ph.configure_logger()
pyrosetta.init(extra_options=ph.make_option_string(no_optH=False,
                                                ex1=None,
                                                ex2=None,
                                                #mute='all',
                                                ignore_unrecognized_res=True,
                                                load_PDB_components=False,
                                                ignore_waters=True)
                               )

# shorthand for sanity
prc: ModuleType = pyrosetta.rosetta.core

# this route strips ligands
pose: Pyrosetta.Pose = pyrosetta.toolbox.pose_from_rcsb('1enh')

```

Note `ignore_unrecognized_res`, `load_PDB_components` and `ignore_waters`.

## Angles

This will not happen to an non-Rosetta user, but I get tripped out by the fact that ψ & φ data is not stored in residue, but in the pose.

```
for residue in pose.residues:  #: prc.conformation.Residue
    # residue.get_phi() # ---> Error!
    pass

```

For the residue index in pose, Rosetta was written in Fortran so is one-indexed.

```
import numpy as np
import numpy.typing as npt
phi:npt.NDArray[float] = np.fromiter(map(pose.phi, range(1, pose.total_residue()+1)), dtype=float)
psi:npt.NDArray[float] = np.fromiter(map(pose.psi, range(1, pose.total_residue()+1)), dtype=float)

```

## Change φ/ψ angles

Time for some fun in dihedral space, say:

```
r = 2
phi = 180
pose.set_phi(r, phi)

```

For modelling, I wrote some code for setting spans with the archetypical angles for different SS: here. I do not recall whence the data came (maybe Wikipedia?), but it may be handy anyway.

## Save

To save a PDB file:

```
pose.dump_pdb('.pdb')

```

There is also `pose.dump_scored_pdb`, but that is only for friends of J. Willard Gibbs.

## Minimisation

It is said in the comments that energy minimisation is not required. But here is a snippet that minimises the sidechains only

```
scorefxn = pyrosetta.get_fa_scorefxn()  # ref2015 will do.
before: float = scorefxn(pose)
cycles = 5
movemap = pyrosetta.MoveMap()
movemap.set_bb(False)  # backbone
movemap.set_chi(True)  # sidechain
relax = pyrosetta.rosetta.protocols.relax.FastRelax(scorefxn, cycles)
relax.set_movemap(movemap)
relax.apply(pose)
after: float = scorefxn(pose)
print(after - before)  # Willard says "Remember free energy is a potential so it should be negative"

```

If the BB tweaks are minimal, there's also the option of scoring/minimising against the electron density etc.
> 4 votes

# Answer
I answered about doing it with PyRosetta as I assumed that as BioPython PDB module is basic and the fact it was being asked made me think it was not possible. Turns out I was wrong. So I am making a separate answer —as I do not actually believe doing it with Bio.PDB is the best solution.

The provided code from above with type hints is:

```
import numpy as np
from Bio import PDB

from typing import TypedDict, Dict, Tuple

parser = PDB.PDBParser()
structure: PDB.Structure.Structure = parser.get_structure("", "/Users/matteo/Downloads/1enh.pdb")
structure.atom_to_internal_coordinates()
chain: PDB.Chain.Chain = list(structure.get_chains())[0]
ic_chain: PDB.internal_coords.IC_Chain = chain.internal_coord
dihedrals:  Dict[Tuple[PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey,
                       PDB.internal_coords.AtomKey],
                 PDB.internal_coords.Dihedron] = ic_chain.dihedra

```

I got the types simply by using `type` and the attributes of note with `dir`. The latter helps me explore the `Bio.PDB.internal_coords.Dihedron` telling me of the `angle` attribute:

```
d: PDB.internal_coords.Dihedron = list(dihedrals.values())[0]
d.angle: np.float64

```

A change here appears it the main `PDB.Structure.Structure` instance, i.e. it is not a disconnected copy.

The class `PDB.internal_coords.AtomKey` appears to be hard to instantiate by humans and may not work as expected —`help` is gibberish so it is a case for `inspect.getsource`. So iterating across the dihedrals for the Ramachandran dihedrals is the solution. Now, a φ angle is N-CA-C-N and a ψ is C-N-CA-C (in PDB atom names the carboxylic carbon is `C` ). When the `__str__` method of one of these `PDB.internal_coords.AtomKey` instances is called some values are given, which are available in the tuple attribute `akl`

```
atom_key: PDB.internal_coords.AtomKey = list(ic_chain.dihedra.keys())[0][0]
print(atom_key)  #3_R_N
print(atom_key.akl) #('3', None, 'R', 'N', None, None)

```

The first entry is the residue number and is a string because of insertion codes. The fourth entry is the atom name. So iterating across the keys for the correct sequence is the trick —I know there's bounty asking for a complete solution, but it's simple, so anyone care to fill out the blanks and get the meaningless points, please do!
> 4 votes

---
Tags: python, protein-structure, pdb, 3d-structure, biopython
---


# Question
Title: Why can't AlphaFold predict the consequences of point-mutations?
In the literature, it specifically states that AlphaFold has "Has not been trained to predict structural consequences of **point mutations**". See : https://alphafold.com/faq

> AlphaFold has not been validated for predicting the effect of mutations. In particular, AlphaFold is not expected to produce an unfolded protein structure given a sequence containing a destabilising point mutation

What I want to understand, is why?

If I change **one amino acid**, and it's **not destabilising**, why can't AlphaFold use its training and form a moderately confident structure? I mean, to a computer, what's the difference between a wild-type and a mutated protein? I know what you may be thinking, why don't you just **try it** yourself, and see if AlphaFold works. But if I don't have the **experimental structure** of the mutant, how on earth am I supposed to believe that this is even close to the theoretical real thing?

Maybe I'm not understanding the **main premise**, but it's something that's been bugging me for a while now. Any insight into this would be much appreciated.

# Answer
Other have addressed why AlphaFold2 cannot really get the necessary signal from MSA (Multiple Sequence Alignment) and memory to properly model a variant, but I thought I'd add my observations/rants.

I have modelled *in silico* both engineering and clinical variants in the past, using PyRosetta, MD and other tools. For me, AlphaFold2 has really stood out but in a different way than the hype.

## Variants can have complicated effects

I developed a tool Venus, https://venus.sgc.ox.ac.uk/, which collates different pieces of information to help figure out the effect a variant was having beyond a simple ∆∆G calculation. This is because often a variant is having a variety of effects, for example a stable but catalytically dead variant can still form a complex and sequester a critical component (e.g. in ref POLR2E-sequestering POR2A is more pathogenic than a simple truncation). This is not only true with human clinical variants, but also with engineering variants: for one protein I have worked on, DogCatcher, lowering the isoelectric point away from neutral pH was critical, not the stabilisation of the transition state (ref). Another example is that the entrance loops in cytochromes drive selectivity more so than the active site configuration.

## The unknowns: alt conformers and complexes

I will be the first to admit that most folk want a deleterious ∆∆G value and will p-value shop for it and not catering for this is a bad idea for a publication. But in terms of mechanism, the biggest issue I face are dynamic processes involving a residue where the other active/inactive conformer is unknown or when the residue is a surface residue possible at the interface in a complex, which is known simply to exist from a coIP western. No nearby gnomAD variants or presence of post-translational sites or sequence conservation do show that something may be important but not why.

## Enter AlphaFold2

AlphaFold2, via the more user-friendly ColabFold, allow the creation of wild type monomers, homo-oligomers, but also hetero-oligomers (complexes) and alternative conformations. These days I seen a cartoon of an interaction and I model it in quintuplicate. Sometimes nothing binds, but sometimes a complex is found 5 out of 5. The weirdest part is that it reveals some novel mechanism that makes sense that has probably never been seen before, but is beyond the scope of the experiment so one has the power to just put them aside —interfaces which are stronger phosphorylated, inhibitory domains binding in the same site but in a different way, unexpected complex arrangements etc. So the true revolutionary power of AlphaFold2 is predicting complexes and alternative conformations, which can be used to predict the mechanism of a variant.

## WIP

Complexes do not always resolve and there is a lot of discussion in terms of massaging the MSA —e.g. boney-fish+tetrapods only etc. And there is very little cross-linking mass-spec data available for validation and/or constraining. But hopefully things will get better.

## Footnote

The last paragraph sounds like an exaggeration, so I best give an example. Here is an example for the drawer from today. A target protein (catalytic domain = pale) has a regulatory terminal which inhibits in isolation (terminus = turquoise), however, when modelled with another protein (mustard) this terminal (terminus = coral) is flipped away. Thus revealing the mechanism by which the mustard protein sequesters the regulatory helix (coral), which would otherwise (turquoise) inhibit the catalytic domain (pale cyan). In the literature, this is a western blot. Where a variant of interest involved in this mechanism one state may be more energetically favourable by scoring the ∆G of the different "snapshots" (=static conformers).
> 3 votes

# Answer
This answer is from common knowledge rather than specific infield expertise. The question could be more expertly answered by other members.

If you really want to pursue this type of calculation, it's possible that something like `INPS-MD` *de novo* -\> protein structure stability from variants will provide an answer. So maybe compare the outputs with `AlphaFold`?

Intuitively, I don't think it's cool, just my opinion. In my opinion, an experimentally derived structure is required as the basis for your calculation, where you want to consider the impact of an amino acid mutation.

It's important to separate *de novo* prediction from variant prediction. Alphafold is *de novo* prediction. What you are looking at is structure stability from a variant (see Summary). That's my core answer, these different calculations. To explain this more clearly ...

Good variant stability calculators are DUET (2014) and INPS-MD (2016), amongst others. Just to reiterate, a good starting point is an experimentally derived protein structure.

The webservers are here:

I personally don't know `Alphafold` but understand its predecessors which accommodate convolutions by leveraging deep-learning. The basis of the predecessor's calculations firstly are *de novo* prediction and secondly can miss the impact of variants if they are not directly involved in the core-algorithm. It's about identifying which amino acids are interacting considering all available homologous amino acids within a database. The calculation can predict local interactions with extreme accuracy, but it's a prediction. It will also depend on whether there are single amino acid changes or two- or more amino acid changes (they could compensate each other, particularly for a natural variant).

In contrast, DUET and INPS-MD, amongst a lot of other similar programs, consider the property shift of an individual amino acid in relation to its immediate local stereo-chemical environment measured by a shift of free-energy.

Essentially, these are two sets of very different calculations.

**Summary** The problem I intuitively see is if *de novo* prediction in any form is used, okay they are very good, but they're just predictions. Variant structural predictions are also just predictions which are very widely respected. Stacking *de novo* prediction onto predicting protein structure stability resulting from variants has scope for error, especially if there was some hidden non-linear effect.

There could be additional issues because the *de novo* calculations immediately prior to `Alphafold` don't predict the precise local distances.

---

`AlphaFold` would simply respond

> 'it's not been trained'

They'd want a wet-lab experimental dataset with wild-type vs. variants, if these exist with sufficient abundance to prevent over-training. It would be a very skewed dataset so of limited application in my personal opinion.
> 1 votes

# Answer
As marcin points out in the comments: by "AphaFold (AF) has not been trained to predict structural consequences of point mutations" it is meant "AF is not able to tell you whether your point mutation is destabilizing or not, and if it is destabilizing then the structure it predicts will most probably be wrong".

It is totally possible that AF predicts the structure of a protein well when only non destabilizing mutations are introduced. As pointed out in the paragraph marcin quotes, however: "AlphaFold has not been validated for predicting the effect of mutations." And as long as we don't know, we should be warry of any of AlphaFold predictions with point mutations. Maybe it works, maybe it doesn't, or maybe it's case dependent.
> 1 votes

---
Tags: protein-structure, mutations, modelling, deep-learning, point-mutation
---


# Question
Title: drawing multiple sequence alignment as tree in R
I am very new here. Please help me with my codes as am trying to draw multiple sequence alignment as tree suing ggtree package in R and i have following error.

```
library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("get_result.nwk")
msaplot(p=ggtree(tree), fasta = "get_results", window=c(150, 175))

Error in inds:inde : NA/NaN argument
> msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "get_result.fasta", offset = 5)
Error in inds:inde : NA/NaN argument

```

My data looks like

```
LN681202.1_Ralstonia_solanacearum_isolate RsBe2
TGCAAGTCGAACGGCAGCGGGGGTAGCTTGCTACCTGCCGGCGAGTGGCGAACGGGTGAGTAATACATCGGAACGTGCCCTGTAGTGGGGGATAACTAGTCGAAAGACTAGCTAATACCGCATACGACCTGAGGGTGAAAGTGGGGGACCGCAAGGCCTCATGCTATAGGAGCGGCCGATGTCTGATTAGCTAGTTGGTGGGGTAAAGGCCCACCAAGGCGACGATCAGTAGCTGGTCTGAGAGGACGATCAGCCACACTGGGACTGAGACACGGCCCAGACTCCTACGGGAGGCAGCAGTGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCTTCGGTTAATACCTGGAGTGGATGACGGTACCGGAAGAATAAGGACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTCCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTGTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAATTGCATTGGTGACTGCACAGCTAGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGATAACACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCCTAAACGATGTCAACTAGTTGTTGGGGATTCATTTCCTTAGTAACGTAGCTAACGCGTGAAGTTGA
MF179616.1_Kv-1
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTCCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MN725749.1_SA006
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTCCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF183971.1_Klebsiella_sp._strain_Kv-2
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MT507217.1_Uyi23r
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF183969.1_Kv-3
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MK929064.1_NM43TS4-1
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF183973.1_Kv-4
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MN082125.1_IILSFPR619
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF183974.1_Kv-5
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MK671348.1_1
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF183970.1_Kv-6
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF682957.1_BT3409
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF183972.1_Kv-7
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MZ475068.1_Klebsiella_variicola_EMBL-1
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
MF373424.1_KV8
GCTCTCGGGTGACGAGCGGCGGACGGGTGAGTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACCTGGAAACGGTAGCTAATACCGCATAACGTCGCAAGACCAAAGTGGGGGACCTTCGGCCTCATGCCATCAGATGTGCCCAGATGGGATTAGCTGGTAGGTGGGGTAACGGCTCCCTAGGCGACGATCCCTAGCTGGTCTGAGAGGATGACCAGCCAACTGGAACTGAGACACGGTCCAGACTCCTACGGGAGGCAGCAGTGGGGAAATTGCACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGACTTTCAGCGGGGAGGAAGGCGGTGAGGTTAATAACCTCATCGATTGACGTTACCCGCAGAAGAAGCACCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCGTCAAGTCGGATGTGAAATCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGATAACCGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGTTAAAACTCAAATG
KC853296.1_JM16
GCTCTCGGGTGACGAGCGGCGGACGGGTGAGTAATGTCTGGGAAACTGCCTGATGGAGGGGGATAACCTGGAAACGGTAGCTAATACCGCATAACGTCGCAAGACCAAAGTGGGGGACCTTCGGCCTCATGCCATCAGATGTGCCCAGATGGGATTAGCTGGTAGGTGGGGTAACGGCTCCCTAGGCGACGATCCCTAGCTGGTCTGAGAGGATGACCAGCCAACTGGAACTGAGACACGGTCCAGACTCCTACGGGAGGCAGCAGTGGGGAAATTGCACAATGGGCGCAAGCCTGATGCAGCCATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGACTTTCAGCGGGGAGGAAGGCGGTGAGGTTAATAACCTTATCGATTGACGTTACCCGCAGAAGAAGCACCGGCTAACTCCGTGCCAGCAGCCGCGGTAATACGGAGGGTGCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCGTCAAGTCGGATGTGAAATCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCTCAGGTGCGAAAGCGTGGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGATAACCGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGTTAAAACTCAAATG

```

The tree info looks like this

```
(((((MF183974.1_Kv-5:0.0000000001,
(MN082125.1_IILSFPR619:0.0000000001,
(MF183973.1_Kv-4:0.0000000001,
(MK929064.1_NM43TS4-1:0.0000000001,
(MF183969.1_Kv-3:0.0000000001,
(MT507217.1_Uyi23r:0.0000000001,
MF183971.1_Klebsiella_sp._strain_Kv-2:0.0000000001)0.0000000000:0.0000000001)
0.0000000000:0.0000000001)0.0000000000:0.0000000001)0.0000000000:0.0000000001)
0.0000000000:0.0000000001)0.0000000000:0.0000000001,
MK671348.1_1:0.0000000001)0.8360000000:0.0012113715,
(MZ475068.1_Klebsiella_variicola_EMBL-1:0.0000000001,
(MF183972.1_Kv-7:0.0000000001,
(MF682957.1_BT3409:0.0000000001,
MF183970.1_Kv-6:0.0000000001)
0.0000000000:0.0000000001)0.0000000000:0.0000000001)0.0000000000:0.0000000001)
0.0000000000:0.0012061821,
((MN725749.1_SA006:0.0000000001,
MF179616.1_Kv-1:0.0000000001)
0.0000000000:0.0000024449,
LN681202.1_Ralstonia_solanacearum_isolate_RsBe2:99.9997264606)
0.0000000000:0.0000002190)1.0000000000:4.9050549661,
KC853296.1_JM16:0.0012459932,
MF373424.1_KV8:0.0000028057);

```

# Answer
> 1 votes
Firstly, do check `tree` is not an empty object, i.e. it uploads the tree file. Secondly, check the path to the alignment is correct `get_results`. Normally an alignment file would be suffixed with the alignment format, e.g. `get_results.fa`.

The immediate thing that strikes me, if all the paths are correct and the file names is the alignment format supplied isn't `fasta` format.

If you replace the contents of the file `get_results` (not `get_results.nwk` \- because that's your tree file, with the following (wait ... see **Note** and **Note2** ) ...

```
    >LN681202.1_Ralstonia_solanacearum_isolate_RsBe2
TGCAAGTCGAACGGCAGCGGGGGTAGCTTGCTACCTGCCGGCGAGTGGCGAACGGGTGAGTAATACATCGGAACGTGCCCTGTAGTGGGGGATAACTAGTCGAAAGACTAGCTAATACCGCATACGACCTGAGGGTGAAAGTGGGGGACCGCAAGGCCTCATGCTATAGGAGCGGCCGATGTCTGATTAGCTAGTTGGTGGGGTAAAGGCCCACCAAGGCGACGATCAGTAGCTGGTCTGAGAGGACGATCAGCCACACTGGGACTGAGACACGGCCCAGACTCCTACGGGAGGCAGCAGTGGGGAATTTTGGACAATGGGGGCAACCCTGATCCAGCAATGCCGCGTGTGTGAAGAAGGCCTTCGGGTTGTAAAGCACTTTTGTCCGGAAAGAAATCGCTTCGGTTAATACCTGGAGTGGATGACGGTACCGGAAGAATAAGGACCGGCTAACTACGTGCCAGCAGCCGCGGTAATACGTAGGGTCCAAGCGTTAATCGGAATTACTGGGCGTAAAGCGTGCGCAGGCGGTTGTGCAAGACCGATGTGAAATCCCCGGGCTTAACCTGGGAATTGCATTGGTGACTGCACAGCTAGAGTGTGTCAGAGGGAGGTAGAATTCCACGTGTAGCAGTGAAATGCGTAGAGATGTGGAGGAATACCGATGGCGAAGGCAGCCTCCTGGGATAACACTGACGCTCATGCACGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACGCCCTAAACGATGTCAACTAGTTGTTGGGGATTCATTTCCTTAGTAACGTAGCTAACGCGTGAAGTTGA
>MF179616.1_Kv-1
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTCCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
>MN725749.1_SA006
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTCCGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
>MF183971.1_Klebsiella_sp._strain_Kv-2
 GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC
>MT507217.1_Uyi23r
GGAATTACTGGGCGTAAAGCGCACGCAGGCGGTCTGTCAAGTCGGATGTGAAATCCCCGGGCTCAACCTGGGAACTGCATTCGAAACTGCAGGCTAGAGTCTTGTAGAGGGGGGTAGAATTCCAGGTGTAGCGGTGAAATGCGTAGAGATCTGGAGGAATACCGGTGGCGAAGGCGGCCCCCTGGACAAAGACTGACGCCAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCTGGTAGTCCACCTGTAAACGATGTCGATTTGGAGGTTGTGCCCTTGAGGCGTGGCTTCCGGAGCTACGTTAAATCGACCGCCTGGGGAGTACGGCCGCAAGGAAACTCAAATGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATGTGGTTTAATTCGATGAACGCGAAGAACCTTACCTGGTCTTGACATCCACAGAACTTTCCAGAGATGGATTGGTGCCTTCGGGAACTGTGAGACAGGTGCTGCATGGCTGTCGTCAGCTCGGTTGTGAAATGTTGGGTTAAGTCCCGCAACCGCAACCCTTATCCTTTGTTGCCAGCGGTTAGGCCGGGAACTCAAAGGAGACTGCAGTGATAAACTGGAGGAAGGTGGGGATGACGTCAAGTCATCATGGCCCTTACCCAGGGCTACACACGTGCTCAATGGCATATACAAAGAGAAGCGACCAGAGCAAGCGGACCTCATAAAGTATGTCGTAGTCGATTGGAGTCTGCAACTCGACTCCATGAAGTCGGAATCGCTAGTAATCGTAGATCAGAATGCTCGGTGAATACGTTCCCGGGCCTTGTACACACCGC

```

Also ....

* Please remove ALL space after the taxa name both in the alignment file and the treefile (I've done it for you in the alignment file)
* `LN681202.1_Ralstonia_solanacearum_isolate RsBe2` \- the space is a problem between the word `isolate` and `RsBe2` ... please make sure in the treefile the space is replaced with a `_` character.

The names in the treefile must exactly match the names in the alignment file. Just to confirm, the last alignment file you forwarded had more taxa (that space is still an issue BTW), however ...

* Fasta format must start the taxa name with a `>` otherwise it will not recognise it.

I think thats the basic error and why its saying `null` because it can't parse the alignment file. Getting the treefile names precisely aligned with the alignment taxa names might be an issue.

---

**Note** My fasta alignment format still looks like a dogs dinner, I can't get the fasta format on the markdown language working ....

1. fasta format there must be a carriage return between the taxa name and the sequence.
2. the `>` must not have any spaces between the start of the line (unlike my version) until the `>` ... nor any spaces anywhere on the line ...
3. Spaces within the taxa name could cause problems in `Newick` tree format (depends what program you've used).

**Note2** The alignment I edited and the new alignment you just supplied (presumably complete) obviously can't be simply cut and paste as I originally thought, but you get the gist. Writing code is a much easier way to do phylogenetics because there's very little manual editing and parsers do everything for you.

---

In the tree file you've correctly named `LN681202.1_Ralstonia_solanacearum_isolate_RsBe2`, but the same name must in the alignment file. Generally taxa names don't use spaces and are replaced with `_`, in `phylip` format for example this causes mayhem with the parser.

---

**Finally** If thats all thats wrong it would be cool, I think they'll be more to get it working correctly.

---
Tags: phylogenetics, phylogeny, multiple-sequence-alignment, ggtree
---


# Question
Title: CLI access to genomes/proteomes by BioSample ID
*This question was also asked on Biostars*

I'd like to download multiple genome assemblies or proteomes using a set of BioSample IDs from NCBI.

I'm able to find the assemblies belonging to the BioSample IDS in a browser (in the search field of https://www.ncbi.nlm.nih.gov/), but couldn't find a commandline interface solution (CLI).

E.g. for BioSample `SAMN09405588` the assembly id is `PDT000806148.1`, and from https://www.ncbi.nlm.nih.gov/assembly/GCA\_014136285.1/ I can download the proteome: `GCA_014136285.1_PDT000806148.1_protein.faa.gz`

Thank you for your help!

# Answer
> 0 votes
*Answer from @m, converted from comment:*

Biostars are recommending esearch ... its not the approach I'd have taken but presume will do the job. Its available via SRA, so I'd have gone for the SRA toolkit. Is SRA more appropriate in context? (I liked Biostars BeautifulSoup solution)

*\[Please update or improve this answer if you can\]*

---
Tags: python, phylogenetics, identifiers
---


# Question
Title: Expression analysis of miRNAs with normal RNA-seq data without small RNA-seq data
I am looking to perform expression analysis of miRNAs with normal RNA-seq data lacking small RNA-seq data?

Which path should I choose for known miRNAs and unknown new miRNAs?

**Data set:** rna-seq data that analyzes the expression differences of a plant under stress.

**Objective** To assess known and new miRNAs.

# Answer
> 1 votes
You are going to have to look up the library prep protocol used, and determine if it will capture miRNA. My guess is, in the course of most RNASeq preps, small products are removed as contamination, so you won't have any real miRNA left behind.

I believe you need to use a special library prep designed for miRNA to study them.

---
Tags: rna-seq, ngs, rna
---


# Question
Title: How can to validate the presence of a certain type of cells in a single cell dataset?
I have a single cell dataset which I consider to be a reference , let say, for an human organ, I have identified some new clusters that correspond to cell types. I also have WGCNA modules that correspond to a identified cluster. Given this background, is there any possibilities to validate the presence of the cell types from my reference dataset inside another single cell dataset ?

I was thinking about GSEA, but all the gene in a set have the same importance. Could I use eigen vectors from first principal component for a given module expressed in the corresponding cell type, as a weight. Then calculate inside another dataset a sort of "first principal component" using the previous eigen vectors. Then I would simply calculate the correlation between the first principal component and the calculated one to validate whether my cell type is present or not. I was also thinking about a sort of weigthed GSEA, if anything like this ever existed.

# Answer
I think the typical way to do this would be to identify a small number of marker genes whose presence and absence is specific for that cell type.
> 1 votes

---
Tags: single-cell, transcriptome, wgcna, gsea
---


# Question
Title: How to convert fasta alignment to nexus or phylip format?
I generated a multiple sequence alignment using mafft, but I would like to convert it to nexus or phylip format (required by a downstream software). What's the easiest way to do the conversion?

# Answer
EMBO's bioinformatic suite EMBOSS has seqret program that can do exactly that. You can install it via conda (`conda install -c bioconda emboss`) and then the conversion to `nexus` is

```
seqret -sequence input.mafft.fasta -outseq output.mafft.nex -osformat nexus

```
> 1 votes

# Answer
```
from Bio import AlignIO
from pathlib import Path
from Bio.AlignIO.PhylipIO import SequentialPhylipWriter
import re

file = 'input.mafft.fasta'
fileoutnex = re.sub('fasta', 'nex', file)
pathway = '/mypath/'
alignment = AlignIO.read(Path(pathway, file), 'fasta')

# write out in nexus
AlignIO.write(alignment, Path(pathway, fileoutnex), 'nexus')

# write out in relaxed phylip
fileoutphy = re.sub('fasta', 'phy', file)
with open(Path(pathway, fileoutphy), 'w') as output_handle:
    SequentialPhylipWriter(output_handle).write_alignment(alignment, id_width=30)

```

The advantage here is you can select what phylip format you want. That is a BIG advantage, traditional phylip format is terrible, particularly with taxa ids lengths. Here relaxed phylip format is being used.

The other advantages,

* Easily parse, change, alter the taxa names;
* Phylip taxa id can be ANY length (`id_width=`);
* save it as a script;
* no need to write an outfile (ok minor);
* pipe the alignment straight into downstream stuff;
* very easy manipulate the alignment.
> 1 votes

---
Tags: sequence-alignment, format-conversion
---


# Question
Title: Select representetive sequence within MSA (protein)
I am trying to pick the best representative sequence from an amino acid MSA. I do not want to use the consensus, as it generates an artificial one. Here is my approach:

1. I calculate a Percentage Identity Matrix (PIM) from the MSA file
2. I check the highest value in each column and rank the respective sequences. What I mean by this, for example, the values in the first column would be the Percentage Identity of the first seq compared the rest of the dataset. If I pick the highest value in there, that would give me the sequences that are highly identical to the first seq.
3. I do this for the every column and rank the sequences. The top ranked sequence would yield me the best representetive sequence within my MSA data set.

The problem with this approach is that its computationally expensive to calculate the PIM as its doing a character match, comparing each sequence's position with all of the dataset. Is there any easier way to achieve this? Im not sure, if we could achieve something similar by looking at the tree. I think we can find the lowest common ancestor by the tree but I do not think that would yield a representetive sequence.

# Answer
> 2 votes
I have two suggestions that I used in the past (both using Python/BioPython though). They are arguably the same, but obtained through different paths.

1. Calculate the (artificial) consensus sequence of your alignment. Then, compare each sequence in your MSA to the consensus in a pairwise alignment, keeping the score. The highest scoring **real** sequence will be your representative.

```
from Bio import AlignIO
from Bio import pairwise2
from collections import Counter

aln_file = "MSA.aln"
alignment = AlignIO.read(aln_file, "clustal") # whatever format you have

for i in range(0,alignment.get_alignment_length()):  # iterate position by position
    aa_position = []

    # Get all residues in that position
    for record in alignment:
        aa_position.append(record.seq[i])

     # Use the most common residue as representative
    consensus_list.append(Counter(aa).most_common()[0][0])
consensus_seq = "".join(consensus_list)

ref_seq = None
highest_score = 0
for record in alignment:
    # globalxx finds the best global alignment between both sequences.
    # The first x indicates that identical characters are given 1 point, 0 otherwise
    # Second x indicates no gap penalties
    pairwise_score = pairwise2.align.globalxx(record.seq,
                                              consensus_seq,
                                              score_only=True
                                              )
    if pairwise_score > highest_score:
        highest_score = pairwise_score
        ref_seq = record
    elif pairwise_score == highest_score:
        # If they are equal, you can prioritize other things
        pass

```
2. You can again obtain the consensus sequence, and then calculate the Levenshtein (edit) distance between each record and the consensus. The sequence with the lowest distance will be the most representative. The easiest way to calculate it is via installing the python-Levenshtein module with `pip install python-Levenshtein`.

```
from Levenshtein import distance as edit_distance

# First, calculate the consensus as described above

ref_seq = None
best_distance = 999999999999
for record in alignment:
    distance = edit_distance(record.seq, consensus_seq)
    if distance < best_distance:
        best_distance = distance
        ref_seq = record
    elif distance == best_distance:
        # If they are equal, you can prioritize other things
        pass

```

# Answer
> 0 votes
Suggestions in increasing order of relevance:

1. The only program that is close is rootdigger. This finds the sequence most likely to be the root, as the name suggests. Available via GitHub, but not via conda.
2. You could just write your `PIM` approach as small program? It wouldn't be difficult.
3. Alternatively, you could just select the longest, shortest or median branch length and use that as a representative sequence. This wouldn't be computationally expensive at all (you'd need the correct root however and a rooted tree\*).
4. The shortest branch length from the root is a reasonable solution. You can do it via a GUI (if its e.g. a Bayesian tree, there's one program that does this), write you own code, or else manually via optical analysis (easy to spot the shortest branch in a rooted tree).

The shortest branch-length is the closest sequence to the ancestral sequence for the in-groups (non-root). You can obtain a reconstruction of the ancestral sequence - it is computationally (very) expensive and isn't a real sequence.

\*, Point 3 and 4 require a phylogram rather than a dendrogram of course.

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


# Question
Title: Extracting base sequences from ABI/AB1 sanger sequencing chromatogram
I am trying to understand the sanger sequencing ABI/AB1 file format better, and extract base calls from given signal intensities over time.

As I understand, reading in a raw AB1/ABI file into python, I can access the channel corrected values for the different bases in \['DATA9'\] to \['DATA12'\], and understand which bases belong to a specific channel through \['FWO\_1'\] i.e either ordered GATC or ATGC etc.

However, the file I am looking at has identified 412 bases as the output from SnapGene, yet has 4950 data records for each channel in the raw \['DATA9'\] to \['DATA12'\]. How do you convert those 4950 data records into the bases correctly? Is there some form of normalisation - e.g every 10 raw records gives data for 1 base position, or peak calling signal detection algorithms. And if so, do you take the average over 10 records, or the highest peak in the 10 records from the channels and the channel with the highest peak/average is the correct base? Do you start at the beginning or do you have some form of cutoff? Does it make sense to convert this figure into a quality score similar to NGS, and do you use the highest peaks in the 10 record interval or the average?

I've tried using scipy's signal.find\_peaks() function (python) which can get me various peaks with a minimum height and a minimum distance to other peaks, currently using a minimum peak height of 150 and distance to other peaks of 5, but still getting around 20-30 more peaks than the bases identified.

Could anyone suggest any logic, methods, algorithms or resources that could help with translating the raw signal intensities for the 4 channels into a sequence of accurate base calls?

I know software exists to identify the bases, but that's not what I am looking for. I intend to do this from scratch to understand the workings better so I can write some custom processing algorithms.

I hope this question makes sense.

Thanks

# Answer
> 1 votes
You can access the peak positions of Sanger Sequencing chromatograms from `["PLOC1"]`.

The position values in `["PLOC1"]` correspond to the index values in `["DATA9"]`, `["DATA10"]`, `["DATA11"]`, `["DATA12"]`, so you can extract the signal values of each base at the positions.

---
Tags: python, ngs, genomics, genetics, base-calling
---


# Question
Title: How should I interpret DGE results if only one HLA-A gene shows up as significant but not the others?
I have done a DGE recently and have been looking at the DGE list. One of the genes is HLA-A. However, when I dug deeper I realised there are hundreds of HLA-A genes with unique ENSEMBL number (of course it's HLA...). A few other HLA-A genes also show up on the table but are not significant at all.

Considering the way Salmon deals with multimapping reads, can I assume these expressing genes are different enough that they shouldn't be consider the same? Or should the expression of all HLA-A genes considered as one single gene for the purpose of identifying differentially expressed genes?

I do understand that this may be more of a biological question than bioinformatics, but I also want to see if it is a practice pouring the HLA genes (or other variable genes/loci of similar nature) with different sequences together when doing DGE, and whether the statistics would be valid if I do so?

# Answer
> 0 votes
It is bioinformatics IMO because HLA is a massive area of human population genetics and its complex. As a personal opinion I think you need further supporting evidence.

If I understand the scenario correctly,

1. you are confident the data is correct (see below),
2. a given HLA-A genotype is significant above other HLA-A and HLA-B genotypes

Thats publishable result, cutting-edge potential, you just need to find the given genotype and represent the association within a given disease association study. Yes your result is absolutely of interest particularly if you consider the wider output of association studies (HLA-B). It is trendy to identify a given HLA to a given association. If this is non-pathogen related they HLA exhibits very strong population skews.

I have put a large caveat below, because:

* you shouldn't underestimate the complexity of this research field and
* I don't really understand the specifics of your study.

---

If you want to ditch or include HLA it would require carefully defined biological aims. It is a specialist area of investigation IMO (below).

HLA diversity and its link to population shifts in infection susceptibility is subject to repeated studies. Different populations show very distinct allelic frequencies. For example, "HLA-B\*40 Allele Plays a Role in the Development of Acute Leukemia". Distinct link between African population and malaria .. These just being a very small sample in a hugely active area of investigation. The idea is that specific genotypes have very specific outcomes in human health.

**Rationale** The single gist is HLA this represents the cellular immune response and is considered essential for infection control, with the notable exception of pathogenic bacteria. Each genotype is considered to be uniquely functional and again is well known to exhibit population skews. It is beyond 'trendy' is huge and established field of investigation.

**Caveat** If you were looking for 'wiggle room' to avoid HLA-A diversity and allelic frequency and all the associated work it requires ... or you specifically want to make it a standout result ... then the angles that might work are ...

1. alot of examples of association focus more on HLA-B than HLA-A. This is a rule of thumb and might be disputed. So if you're going to ignore HLA-A you could perform an extensive literature review to see the extent that associations of HLA-B predominate. I think you might find that holds. So its either really cool cause its an outlier, or you can ditch it because its an outlier.
2. If pathogenic bacteria were the focus of your investigation, then most would accept HLA is not relevant.

I'm sure there's other angles of justification, e.g. "I recognise HLA-A .... but this will be subject to further investigation by (putative) collaboration".

**Advice** If I was tackling this data I would contact a *bioinformatics immunologist*, because they will have databases and pipelines dedicated to the exact task you require. Personally, I would be reluctant to take this on without specialist assistance, at a minimum someone with current experience. Essentially, what I'm trying to say, is you will need to understand the HLA field because crossing into it isn't trivial. The situations where I could be wrong are that its simply so complex that no-one really understands it, in which case your good to go.

# Answer
> 3 votes
You need specialized tools to study HLA genes. Generic tools won't work well. Before other analysis, you should get the correct genotype first. Try arcasHLA. It has been a popular HLA genotyper for RNA-seq data. When you know the genotype, I guess you may get rid of other HLA-A alleles and redo the differential analysis if you are interested in HLA-A. I haven't done this, though. Just a thought.

---
Tags: r, rna-seq, differential-expression
---


title: PyMol alignment script

text:
# Question
Title: PyMol alignment script
Is there a way write a script in python to use "super" to align proteins from a 2D array? Each row of the array has 2 columns, and for each row, I want to "super" the two proteins (one in each column)

# Answer
You can use Python "normally" and import the `pymol` module. In the script that holds the proteins in your matrix, you can add a function that does the super alignment.

To use it, in your console type `pymol -cq this_script_name.py`.

You need to have pymol installed and accessed via console. Easiest way I found to manage it is to create a new environment via conda and install pymol there.

```
import pymol

def super_prots(prot_1, prot_2):
    pymol.cmd.fetch(prot_1, 'prot_1')
    pymol.cmd.fetch(prot_2, 'prot_2')

    # If you are using already downloaded pdb files, use load instead:
    # pymol.cmd.load(prot_1, 'prot_1')
    # pymol.cmd.load(prot_2, 'prot_2')

    pymol.cmd.hide('all')
    pymol.cmd.set('ray_opaque_background', 0)

    # Show cartoon and paint them different
    pymol.cmd.show('cartoon', 'prot_1')
    pymol.cmd.color('orange', 'prot_1')
    pymol.cmd.show('cartoon', 'prot_2')
    pymol.cmd.color('lightblue', 'prot_2')

    # Align both proteins
    pymol.cmd.super('prot_1', 'prot_2')

    # Improve image quality, and get png.
    pymol.cmd.set('depth_cue', 0)
    pymol.cmd.set('spec_reflect', 0)
    pymol.cmd.set('cartoon_sampling', 15)
    pymol.cmd.set('ribbon_sampling', 15)
    pymol.cmd.set('antialias', 2)
    pymol.cmd.space('cmyk')

    # This creates high-quality images, but might take a long time to
    # process. Adjust parameters accordingly.
    pymol.cmd.png(f'{prot_1}_{prot_2}.png',
                  width=2000,
                  dpi=300,
                  ray=1
                  )
    pymol.cmd.reinitialize()

if __name__ == '__main__':
    # These are entries from PDB. It can work with already downloaded
    # files, but you have to change the function accordingly
    protein_matrix = [['1EWL', '4W5B'],
                      ['4E2D', '3ATZ'],
                      ]
    for prot_1, prot_2 in protein_matrix:
        super_prots(prot_1, prot_2)

```

This will generate .png images of the alignment, but you can also ask it to save the alignment into .pdb files adding `pymol.cmd.save(f'{prot_1}_{prot_2}.pdb')`

The commands' documentation can be found here.
> 1 votes

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


# Question
Title: Step-by-Step Construction of Gene Co-expression Networks from High-Throughput Arabidopsis RNA Sequencing Data
I am following the tutorial here to learn how to process raw RNA-seq data and get gene read counts. I've had no issues until I got to step 3.9, "Assignment of RNA-Seq Reads to Genes." Here, the authors run the command

```
fc0 <- featureCounts(sam.list, annot.ext =
"Ara-port11_GFF3_genes_transposons.201606.gtf", allowMultiOverlap = T, isPairedEnd = F,
nthreads = 2, strandSpecific = 0, isGTFAnnotationFile = T)

```

I ran the same command except I set isGTFAnnotationFile = F because I couldn't find the gtf file, so I used the gff file `TAIR10_GFF3_genes_transposons.gff` instead, downloaded from here (they mentioned that using gtf only speeds up the process but that you can still do this using gff). I'm not sure if this is what is causing the issue, but the error I get is:

```
Error: Line 1 contains a format error. The expected annotation format is SAF.
Error in featureCounts(sam.list, annot.ext = "TAIR10_GFF3_genes_transposons.gff",  :
  No counts were generated.

```

# Answer
*Answer from @an-ignorant-wanderer, converted from comment:*

I ended up using the annotation file from https://support.illumina.com/sequencing/sequencing\_software/igenome.html and it worked.

*\[Please update/improve this answer if possible\]*
> 1 votes

---
Tags: r, rna-seq
---