# Question
Title: Convert Symmetrix Matrix into Lower Triangular Matrix
I'm using Scipy's pdist function to calculate a distance matrix which is working great but it gives me a symmetric matrix and I want to convert it into a lower-left triangular matrix.

This is what it prints out for me:

```
  a   b   c   d
a 0 ### ### ###
b ### 0 ### ###
c ### ### 0 ###
d ### ### ### 0

```

This is what I want it to be:

```
###
### ###
### ### ###

```

I know there's a function in R to do this (upper.tri) but it just replaces everything to the right of the 0 with more 0s. Is there a better/easier way to get it like this? I'm opening the file in excel and removing everything to the right with the 0 manually right now.

# Answer
NumPy has a convenient built-in function for accessing the lower-triangular or upper-trianglar entries of an array. To get the lower-triangular entries, you can use `numpy.tril_indices` to get the indices of these entries, and then use these indices to extract the entries.

```
import numpy as np

ar = np.array(
    [[0, 1, 2],
    [3, 4, 5],
    [6, 7, 8]]
)
idx = np.tril_indices(3, -1)
ar[idx]
# outputs array([3, 6, 7])

```

This doesn't automatically produce the spacing you want, but formatting the final display is likely context-dependent.
> 4 votes

# Answer
There seems to be a scipy method that retains the lower triangle but with the 0s. See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.tril.html#scipy.sparse.tril
> 1 votes

---
Tags: format-conversion, scipy, distance-matrix, megax, pdist
---


# Question
Title: BioPython SeqIO fails with Invalid Alphabet found
I am trying to store a GenBank file that could be read by Geneious later. I found an example how to create the GenBank file here

But now I run into `Alphabet` issues. I tried all protein classes from the documentation, but all fail with the same error:

`TypeError: Invalid alphabet found, <class 'Bio.Alphabet.ProteinAlphabet'>.`

I dont know what else to do here? Using `biopython==1.73`

```
    record = SeqRecord.SeqRecord(
        Seq.Seq("MIRQALAVAALLLAGTAQADGLIDN", alphabet=Alphabet.ProteinAlphabet),
        id="FooBar",
    )

    with NamedTemporaryFile(mode="w") as genbank:
        SeqIO.write(record, genbank, "genbank")
        pass

```

# Answer
> 1 votes
Quickly ...

```
from Bio.Alphabet import generic_protein
from Bio.Seq import Seq

```

then ..

```
myprot = Seq('EXAMPLEPROT', alphabet=generic_protein)

```

---
Tags: python, genbank, biopython
---


title: Retrieve ID ligand from PDB file

text:
# Question
Title: Retrieve ID ligand from PDB file
I have thousands of PDB files and I need to extract the following information from each of them: an ID or the name of the small molecule (ligand) present in the complex.

Is there a way to do that with BioPython or similar tools? Taking this PDB entry as an example, the tools should return 7FM.

# Answer
> 2 votes
I found out that the package PyPDB has the function `get_ligands` that does exactly what I'm looking for.

# Answer
> 1 votes
copying from How to download ligands for PDB structure one of the answers, one of mines, code :

```

import requests

import json

parsed_lig_dict = {}

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")
    parsed = json.loads(risposta.content.decode())
    # print(risposta.status_code)
    # print((parsed))
    # print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    # print('\n', type(parsed))

    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]
    # print('\n\n',parsed_lig, type(parsed_lig))

    # for i in parsed:
    #     print('\n', i  , parsed[i])

    for lig in parsed_lig:
        risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity/{entry_id}/{lig}")
        parsed = json.loads(risposta.content.decode())
        # print('\n',risposta.status_code, risposta.content)

        # print(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"],
              # parsed["pdbx_entity_nonpoly"]["comp_id"])

        parsed_lig_dict[lig] = [parsed["pdbx_entity_nonpoly"]["comp_id"] ,
                                dict.fromkeys(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"])]

    for lig in parsed_lig_dict:
        for chain in parsed_lig_dict[lig][1]:
            risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity_instance/{entry_id}/{chain}")
            # print(risposta.status_code)

            parsed = json.loads(risposta.content.decode())
            # print('\n',risposta.status_code, risposta.content)
            # print(chain , '--> ',parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"])

            parsed_lig_dict[lig][1][chain] = parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"]

    # print(parsed_lig_dict)

    # for i in parsed_lig_dict:
    #     print(parsed_lig_dict[i])

    # for i in parsed_lig_dict:
    #     for y in parsed_lig_dict[i][1]:
    #         print(i,parsed_lig_dict[i][0], y , parsed_lig_dict[i][1][y])

    cnt = 1
    for i in parsed_lig_dict:
        for y in parsed_lig_dict[i][1]:
            seq_id = parsed_lig_dict[i][1][y]
            chain = y
            comp_id = parsed_lig_dict[i][0]
            print(seq_id , y)

            risposta = requests.get(f"https://models.rcsb.org/v1/{pdb_file}/ligand?auth_seq_id={seq_id}&label_asym_id={chain}&encoding=sdf" ,
                              allow_redirects=True)
            print(risposta.status_code)

            filename = f"lig_{comp_id}_{chain}_{seq_id}_{cnt}.sdf"
            open(filename, 'wb').write(risposta.content)
            cnt += 1

if __name__ == '__main__':

   pdb_file = '5to8'

   main(pdb_file)

```

output file : `lig_7FM_B_701_1.sdf` :

```
7FM
  ModelServer 0.9.10

 39 43  0  0  1  0  0  0  0  0  0
   -2.0610  -10.3590   40.5700 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.1170  -11.7950   42.6720 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8640  -12.6100   43.6800 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6970  -10.3290   36.5750 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8210  -10.9080   35.9890 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.6600   -8.1940   35.9530 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.0300   -9.0340   41.4650 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.9390   -9.6680   41.0120 C   0  0  0  0  0  0  0  0  0  0  0  0
    0.7980  -10.0800   39.7250 N   0  0  0  0  0  0  0  0  0  0  0  0
    1.7720   -9.8430   38.8300 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.9200   -9.1930   39.2330 C   0  0  0  0  0  0  0  0  0  0  0  0
    2.9910   -8.8020   40.5580 C   0  0  0  0  0  0  0  0  0  0  0  0
    4.0340   -8.8690   38.3150 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0830   -9.8810   41.9510 N   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3380  -10.4670   41.7760 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.3020  -10.9670   40.4240 C   0  0  0  0  0  0  0  0  0  0  0  0
   -3.8290  -11.6870   41.4760 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.8800  -11.1950   42.8370 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.0400  -12.4130   41.5840 N   0  0  0  0  0  0  0  0  0  0  0  0
   -5.2460  -12.6990   43.0180 C   0  0  0  0  0  0  0  0  0  0  0  0
    1.5160  -10.2930   37.5490 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.4060  -11.1430   37.2090 C   0  0  0  0  0  0  0  0  0  0  0  0
   -2.8420  -10.2200   35.4050 N   0  0  0  0  0  0  0  0  0  0  0  0
   -2.7360   -8.8690   35.4000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.6240   -8.9330   36.5380 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.9670  -12.3410   35.9950 N   0  0  0  0  0  0  0  0  0  0  0  0
   -2.5680  -13.0870   37.1550 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.4690  -13.1930   34.5670 S   0  0  0  0  0  0  0  0  0  0  0  0
   -4.0440  -12.7310   37.3480 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.9670  -12.8650   40.6430 C   0  0  0  0  0  0  0  0  0  0  0  0
   -6.9720  -13.5050   40.9690 O   0  0  0  0  0  0  0  0  0  0  0  0
   -4.6200  -13.4130   38.5900 C   0  0  0  0  0  0  0  0  0  0  0  0
   -5.8010  -12.6200   39.1510 C   0  0  0  0  0  0  0  0  0  0  0  0
   -0.5520  -12.3690   33.7940 O   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1430  -14.5640   34.8900 O   0  0  0  0  0  0  0  0  0  0  0  0
   -3.0530  -13.2140   33.7590 C   0  0  0  0  0  0  0  0  0  0  0  0
    5.0490   -8.1440   38.8710 F   0  0  0  0  0  0  0  0  0  0  0  0
    4.6120   -9.9430   37.7390 F   0  0  0  0  0  0  0  0  0  0  0  0
    3.5350   -8.1080   37.3100 F   0  0  0  0  0  0  0  0  0  0  0  0
  1 16  2  0  0  0  0
  1 15  1  0  0  0  0
  2 17  2  0  0  0  0
  2 18  1  0  0  0  0
  2  3  1  0  0  0  0
  3 20  1  0  0  0  0
  4  5  2  0  0  0  0
  4 25  1  0  0  0  0
  4 22  1  0  0  0  0
  5 23  1  0  0  0  0
  5 26  1  0  0  0  0
  6 24  1  0  0  0  0
  6 25  2  0  0  0  0
  7 12  2  0  0  0  0
  7  8  1  0  0  0  0
  8  9  2  0  0  0  0
  8 14  1  0  0  0  0
  9 10  1  0  0  0  0
 10 21  1  0  0  0  0
 10 11  2  0  0  0  0
 11 13  1  0  0  0  0
 11 12  1  0  0  0  0
 13 39  1  0  0  0  0
 13 38  1  0  0  0  0
 13 37  1  0  0  0  0
 14 15  1  0  0  0  0
 15 18  2  0  0  0  0
 16 17  1  0  0  0  0
 17 19  1  0  0  0  0
 19 30  1  0  0  0  0
 19 20  1  0  0  0  0
 21 22  1  0  0  0  0
 23 24  2  0  0  0  0
 26 28  1  0  0  0  0
 26 27  1  0  0  0  0
 27 29  1  0  0  0  0
 28 36  1  0  0  0  0
 28 34  2  0  0  0  0
 28 35  2  0  0  0  0
 29 32  1  0  0  0  0
 30 33  1  0  0  0  0
 30 31  2  0  0  0  0
 32 33  1  0  0  0  0
M  END
> <model_server_result.job_id>
RFDalM1AmREvgkv0PqdfQw

> <model_server_result.datetime_utc>
2023-10-27 22:13:14

> <model_server_result.server_version>
0.9.10

> <model_server_result.query_name>
ligand

> <model_server_result.source_id>
pdb-bcif

> <model_server_result.entry_id>
5to8

> <model_server_params.name>
atom_site

> <model_server_params.value>
{"label_asym_id":"B","auth_seq_id":701}

> <model_server_stats.io_time_ms>
5

> <model_server_stats.parse_time_ms>
76

> <model_server_stats.create_model_time_ms>
2

> <model_server_stats.query_time_ms>
244

> <model_server_stats.encode_time_ms>
0

> <model_server_stats.element_count>
39

$$$$

```

picture of .sdf file :

# Answer
> 0 votes
my attempt, will get only non water (H2O) or non aa (aminoacid) present in the PDB .... tricky if the ligand is an aa :

```
from Bio.PDB import PDBParser

from Bio.PDB import Selection

d = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
     'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
     'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
     'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

pdb_file = '5to8'

struct = PDBParser(QUIET = True).get_structure(pdb_file, pdb_file+'.pdb')

res_list = Selection.unfold_entities(struct , 'R')

for i in res_list :

    if i.resname != 'HOH' and i.resname not in d :

        print(i, i.resname)

```

output :

```
<Residue 7FM het=H_7FM resseq=701 icode= > 7FM

```

---
Tags: python, pdb
---


title: how to use BLAST+ in an optimal way

text:
# Question
Title: how to use BLAST+ in an optimal way
I have a 32-bp query that I want to perform BLAST+. my query represents different RNA viruses and I do expect that my hits get more than 500,000. but in practice, when I perform it I get around only 700 hits. would you please kindly help me how to increase the efficacy of my BLAST+? I leave my syntax here.

```
blastn -task megablast -db nr -outfmt 7 -remote -max_target_seqs 99999999 -perc_identity 80 -query D:/Query.fas  -out D:/result.txt

```

# Answer
> 1 votes
Megablast is designed to find very similar sequences. It is highly specific but has a very low sensitivity. This is because it has a very large word size. By default, I think it will use a word size of 28 nt. That's almost your entire sequence, so if your target sequence doesn't have 28 contiguous nucleotides that are *identical* to 28 contiguous nucleotides of your query, no hit will be found.

If you are looking for distant homologous sequences, you should use simple `blastn` which defaults to a word size of 11. That will still likely be too small, you might want to change it to something lower (your blast tool should have an option for that), but it will give you more results than megablast.

Finally, if these 32nt are translated into protein, you can consider using tblastn instead.

---
Tags: blast, blastn
---


title: What happened to the PhyLoTA browser?

text:
# Question
Title: What happened to the PhyLoTA browser?
A paper I'm reading from 2021<sup>\[1\]</sup> mentions doing some analysis on data from the PhyLoTA database, with this reference citation:

This paper abstract claims to have an online interface at http://loco.biosci.arizona.edu/pb, but this link is not currently responsive. Does anyone know what happened to this project, and if it is currently accessible at some other location?

---

\[1\] S. M. Harrington et. al. 2021 Properties of Markov chain Monte Carlo performance across many empirial alignments MBE

# Answer
> 1 votes
Phylota was implemented in R and is present on GitHub here, described by Bennett *et al* (2018).

Here's the code...

```
library(phylotaR)
wd <- '[FILEPATH TO WORKING DIRECTORY]'
ncbi_dr <- '[FILEPATH TO COMPILED BLAST+ TOOLS]'
txid <- 9443  # primates ID
setup(wd = wd, txid = txid, ncbi_dr = ncbi_dr)
run(wd = wd)

```

Sanderson *et al* was 15 years ago, albeit Sys. Biol. is a very reputable journal in phylogenetics. I've deployed solutions published there with outstanding results.

---
Tags: phylogenetics, phylogeny, public-databases
---


title: How can I resolve `bind_rows()` error?

text:
# Question
Title: How can I resolve `bind_rows()` error?
So I wanted to download massive amounts of COX1 (or COI) sequences from any available database. I found this script in https://bioinformatics.stackexchange.com/a/13187/16657 this comment that achieves just that using the `bold` R package. I changed the script to fit my data preference:

```
# load packages
library(tidyverse)
library(rentrez)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output

set_entrez_key()

# get class-level taxa within "Mollusca" from NCBI taxonomy
taxa <- downstream("Mollusca", db = "ncbi", downto = "order")
# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Mollusca$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# Download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold, bold_seq, marker = 'COI-5P') %>%
  flatten() %>%
  bind_rows()

# Convert the list of sequences to a data frame
sequences_df <- do.call(rbind, sequences)

# Write sequences to a file
write.fasta(
  sequences = as.list(sequences_df$sequence),
  names = as.list(sequences_df$id),
  nbchar = 80,
  file.out = 'coi5p.fasta'
)

```

In the aforementioned thread, this script worked but now that I ran it I got the following error:

```
Error in `bind_rows()`:
! Argument 1 must be a data frame or a named atomic vector.
Run `rlang::last_trace()` to see where the error occurred.

```

I am not sure what to do for this. I commented on the original thread but since it is kind of old I decided to post about it here.

Note: the `set_entrez_key()` requires a NCBI API key which you can get by registering on NCBI.

Edit #1: The `sequences` object before the `flatten()` command is a large list containing 53 elements. After the command, the object becomes a large list with 257 elements.

# Answer
I managed to find the solution. After @gringer mentioned that they had not seen the `flatten` command I decided to dig deeper. Turns out the command has been superseded in purrr 1.0.0. So I will be posting the new script for anyone who seeks an updated version of the one in the old thread.

```
# load packages
library(tidyverse)
library(rentrez)
library(bold)    # API interface to BOLD
library(taxize)  # for NCBI taxonomy lookup
library(seqinr)  # for FASTA output
library(rBLAST)

set_entrez_key("d00fee673c4499df4c2479ba5fd71850b308") #my ncbi api key, this is essential for pulling the taxonomy through ncbi

# get class-level taxa within "Mollusca" from NCBI taxonomy
taxa <- downstream("Mollusca", db = "ncbi", downto = "class") #setting the taxonomic level of the search. it will pull sequences for all mollusca but seperately for each class. if the search is large change to a lower taxonomic level

# check if taxa present in BOLD
checks <- bold_tax_name(taxa$Mollusca$childtaxa_name)
taxa_bold <- checks[!is.na(checks$taxon),]$taxon

# Download sequences from BOLD for each class-level taxon
sequences <- map(taxa_bold, bold_seq, marker = 'COI-5P')
# Convert the list of sequences to a data frame
sequences_df <- do.call(rbind, sequences)

# Filter sequences based on the "COI-5P" marker. Despite the markes parameter, some files may have integrated multiple loci. The other sequences are removed from the R object
coi5p_sequences <- sequences_df %>%
  filter(marker == "COI-5P")

# Write sequences to a file
write.fasta(
  sequences = as.list(coi5p_sequences$sequence),
  names = as.list(coi5p_sequences$id),
  nbchar = 60,
  file.out = 'coi5p_Mollusca.fasta'
)

```
> 1 votes

---
Tags: database
---


# Question
Title: How do I select a subset of genes for functional enrichment(GO/KEGG) analysis from WGCNA results?
I am working on using WGCNA for a bulk RNA-seq experiment. I have three experimental conditions (control, treatment 1, and treatment 2) and have a total sample size of 18 (6 per group). Using the online WGCNA tutorials, I could analyze my datasets after variance stabilizing transformation (vst). I got the summary of network analysis results as a CSV with Gene significance (GS), p.GS, module membership(MM), and associated p-value (p.MM).

I want to extract a set of significant genes from the module with the highest positive correlation to my treatment conditions. As I am just getting acquainted with WGCNA interpretation, which of these parameters in the CSV should I use for filtering further (GS, p.GS, and module membership and associated p-value (p.MM))?

Also, how would I choose the limits of the criteria? Should it be using the scatter plot (MM vs GS)?

As I am pretty new to WGCNA, any help is appreciated. Thanks in advance!

# Answer
Consider filtering on every parameter: Gene Significance (GS), p-value for Gene Significance (p.GS), Module Membership (MM), and p-value for Module Membership (p.MM).

This doesn't mean you need to necessarily filter on every parameter, but it depends on what exactly you want to do, and you might reasonably choose any/all.

1. For GS, you may want to ensure you're capturing genes strongly associated with your experimental conditions.
2. For p.GS, you may want to ensure statistical relevance, (so filter out genes with high p-values). Similar for p.MM.
3. For MM, you may want it to be high to make sure the genes you select are strongly co-expressed and are relevant members of the module.

Plot MM vs GS to see if you can see an obvious filter threshold.
> 1 votes

---
Tags: networks, wgcna, go-enrichment, kegg
---


title: Standardizing Variant IDs best practise

text:
# Question
Title: Standardizing Variant IDs best practise
I'm currently pondering on the best approach to handle and standardize variant IDs within our department to mitigate the challenges associated with rsids and the potential ambiguity of just using chromosomal positions.

The GWAS glue (https://mrcieu.github.io/gwasglue2/articles/Strategy.html) packagae, state oa method where a variant ID is crafted by integrating the chromosomal position (chr), effect allele (ea), and non-effect allele (nea). It assumes the format: chr:pos\_ea\_nea. For alleles that exceed 10 characters, like certain indels, they apply the murmur32 hashing algorithm from the R digest package.

Given our objectives:

We likely prefer not to have excessively lengthy IDs. Yet, using hashed IDs for certain long variants might be perplexing for some who aren't informed about the approach. Hashing every ID could yield consistency, but perhaps might not be as intuitive for regular use. I'm keen on collecting insights:

Thoughts on the described method from Bristol? What might be the maximum length for an indel we could encounter? Is there a broader effort underway to address this concern? Would adopting the Bristol approach set a potential standard, or should we explore or develop alternative methods? For example should the genome build be built into the id? or species? I appreciate the community's feedback as we strive for clarity and efficiency in our variant naming conventions.

For clarity: I want to store data on variants. For example GWAS summary statistics. I want to be able to match variants from differnt studies so I want to have an identifier for the variants. I could use rs numbers I could use chr:pos I could use chr:pos:ref:alt I could use chr:pos:effect:noneffect I could use chr:pos:population:effect:noneffect or studyid:chr:pos:populatoin:effect:noneffect or etc

There are pros and cons for each of these but I would like to here what people on here think these are. I am also intrested in is there work out there to do this already. e.g. to update the idea of rsids to consider multiple alleles and effects in populations?

# Answer
> Thoughts on the described method from Bristol?

I don't immediately see the benefit to this over rsIDs.

> What might be the maximum length for an indel we could encounter?

Potentially very long. However, why does it matter? This will be rare.

> Is there a broader effort underway to address this concern?

dbSNP comes to mind.

> Would adopting the Bristol approach set a potential standard, or should we explore or develop alternative methods? For example should the genome build be built into the id? or species?

It depends on the type of analysis you're planning on doing. Are you doing cross-species analyses? And if so, what types? Do you expect to be comparing genomes with different genome builds?

> For example GWAS summary statistics. I want to be able to match variants from different studies so I want to have an identifier for the variants.

Both chr:pos as well as rsID work for this, but just make sure that the genome build matches if you go with chr:pos.

> I could use rs numbers I could use chr:pos I could use chr:pos:ref:alt I could use chr:pos:effect:noneffect I could use chr:pos:population:effect:noneffect or studyid:chr:pos:populatoin:effect:noneffect or etc

If you're sure your input data has rs numbers and always will, that might be easiest. However, you may find that sometimes it doesn't have rs numbers (e.g. some Affymetrix data), and it's surpassingly annoying to add them. I'd choose chr:pos, personally, and then consider the effect/non-effect/ref/alt/ on a per-analysis basis, because I would want to still identify the same position if there was a different ref/alt for example.

Edit: Based on the comment (OP wants to match effect direction), I would suggest chr:pos:effect\_allele, but of course this means there is an effect allele (i.e. the position has been studied).
> 1 votes

---
Tags: identifiers, variants
---


# Question
Title: How to get cytoband and gene level copy number from genome wide SNP array copy number data?
I have (human) Illumina genome wide SNP array copy number data. For each SNP genomewide, I have Log R Ratio (LRR) and B Allele Frequency (BAF).

What tool(s) can I use to get the integer copy numbers (either -3 to +3 or 0 to inf) for each cytoband and for each gene in the genome?

I assume that people routinely do this? Are there any R packages/pipelines/code to do this?

Do people use `DNAcopy` and/or `GISTIC` for this?

Thank you.

# Answer
> 2 votes
To obtain cytoband and gene-level copy numbers from genome-wide SNP array copy number data, you might want to consider using the ASCAT tool, as it makes allele-specific copy number segment data and derived integer gene-level copy number values from SNP array data. The SNP6 copy number analysis pipeline (ASCAT2) I believe is specifically for this purpose.

See also:

* CNVkit, which has a method for copy number detection via using targeted reads and the non-specifically captured off-target reads to infer copy number evenly across the genome.
* Array Comparative Genomic Hybridization (aCGH) and SNP Genotyping Technologies
* CNVannotator

---
Tags: r, genome, illumina, copy-number, snp-chip
---


# Question
Title: Generating simulated gold standard NGS dataset/ Any publicly available smal gold standard NGS dataset
I am trying to check my variant generation pipeline for humans. The pipeline is for NGS. I am searching for any dataset(.fastq+.bam+.vcf) for any specific sample. Now, GIAB has a nice dataset for that. However, the problem is that it's not a tiny dataset. And it's not feasible to test my pipeline with this data every time. I tried to generate a synthetic dataset with some read simulators(NEAT, MASON2). However, their aligned read files(.bam) files seem not very good. You can check my Issue for NEAT in link. So, my question is, **Is there any publicly available tiny dataset(.fastq,.bam, .vcf) with which I can validate my pipeline fast? Or, is there any simulator with which I can simulate .bam, .fastq files from a given input .vcf? I am using fastqc,bwa,samtools,GATK4 mainly and GRCH38 reference in my pipeline.**

# Answer
I would suggest subsampling an existing dataset such as the GIAB resource that you've already indicated.

Use a utility such as seqtk to subsample the appropriate number of reads from the dataset, and use that.

Example with `seqtk` to get 1 million reads each from paired fastqs while preserving pairing:

```
  # note -s100 random seed which should be used for both to keep pairs
  seqtk sample -s100 giab_read1.fq 1000000 > giab_sub_r1.fq
  seqtk sample -s100 giab_read2.fq 1000000 > giab_sub_r2.fq

```
> 2 votes

---
Tags: ngs, public-dataset, simulated-data, validation
---


# Question
Title: Current 'most comprehensive' whole genome alignment resource for human?
I found the human genome alignment sets here: https://www.ensembl.org/info/genome/compara/multiple\_genome\_alignments.html

but they doesn't include, for example, ancient human genomes, e.g. those here: https://reich.hms.harvard.edu/allen-ancient-dna-resource-aadr-downloadable-genotypes-present-day-and-ancient-dna-data

or, for example, any of the other 'pan genome' whole genome alignments here: https://humanpangenome.org/

However, those latter two resources don't include whole genome alignments...

About 10 years ago I remember reading a paper that used 'lots' of genomes to uncover novel features using different patterns of conservation, demonstrating the power of the approach.

Have I missed a new, more up to date, whole genome alignment resource?

Many thanks,

# Answer
> 2 votes
The obvious thing missing from your list is the telomere-to-telomere genome assemblies that have been coming out in the last couple years.

Assembly link here.

There are now a few others as well that are similar.

It's been demonstrated that this T2T assembly improves mapping rates, in some cases pretty dramatically.

I am not sure entirely what you want from whole genome alignments, but UCSC if not others has incorporated the T2T assembly. Here is a related commentary (which I haven't read in any detail).

---
Tags: ensembl, multiple-sequence-alignment, pangenome, whole-genome, ancient-dna
---


# Question
Title: Does using a microarray chip that matches the chip used in training lead to higher explained variance of a polygenic score?
I'm trying to replicate an existing polygenic score (i.e. test the accuracy in a new sample), and want to know if matching the original study's microarray chip will improve the accuracy (that is, variance explained in a new test set) of the polygenic score.

The labs I have access to have different multiple microarray chips to choose from (e.g., GSA, GDA, UK Biobank Axiom Array, UK BiLEVE Axiom array, GeneChip 2.0 array, Infinium Core-24 Kit, etc.). Will using a chip that is different than the chip used in the study that produced the polygenic score impact the explained variance (accuracy) of the polygenic score? If so, by how much?

Example: a study on heart disease runs a GWAS on 900,000 individuals who are genotyped with the UK BiLEVE Axiom array. This study uses these summary statistics to produce weights for a polygenic score. Lab A, B, and C are conducting polygenic score validation/replication studies, and they each have a random sample of 10,000 individuals drawn from the general population, and they all use the same summary statistics from the original study. After imputing the genotype data of the 10,000 people, they calculate the polygenic score for each person, and then evaluate how much variance in the phenotype they are able to explain. The only difference is that Lab A genotypes their 10,000 people with the UK BiLEVE Axiom array, whereas lab B uses the Illumina GSA, and lab C uses the Illumina GDA.

Do we have any reason to suspect that there will be a substantial difference in variance explained (r^2) by the polygenic score between labs? Are there any studies which have analyzed something that could be informative to this question?

I've been unable to find research on this question. This should not be dependent on which polygenic score is being studied; there is no reason to specify a particular phenotype.

Note: it is standard to impute both in the original study and the replication.

# Answer
> 0 votes
Genotyping can be validated by whole-genome sequencing, and this has been done on the International Genome Sample Resource. You should be less concerned about the *reproducibility* of genotyping results, and more concerned about the *correctness* or *validity* of genotyping results. Companies will often revisit old SNPchips to mark invalid spots, update rs numbers, and add new genotyping probes that are better on a population-wide scale. In that sense, choosing a SNPchip platform that has an established track record of QC and information about version updates would be a better idea than trying to match a previous study as precisely as possible.

However, with regards to replicating existing studies, a bigger concern for you should be the variation in sampled populations, which plays a big role in the outcome of a polygenic test.

This variation is so extreme that directions of associations can flip for the same variant, especially for low-frequency variants with extremely low p-values (i.e. the ones that often show up as manhattan plot peaks).

The issues around population and study replication are explained well in the following paper from 2001:

https://doi.org/10.1006/tpbi.2001.1543

This is not something that can be fixed by attempting to replicate studies in more populations. As Cecile Janssens says, larger sample size only affects random error, not systematic error.

https://twitter.com/cecilejanssens/status/1119235476567801856

---
Tags: snp, microarray, illumina, imputation, polygenic-risk-score
---


# Question
Title: Blastn through NCBIWWW not giving any results, as opposed to the online BLAST
I am trying to create a simple pipeline that takes a file containing various short sequences, and run them through blast. I have an example sequence;

```
>some_random_sequence
AAGGTTCGGTCCAAATTGAA

```

which returns a hit for FUT3 (using the RefSeq\_Gene DB and using grch38 as an entrez\_query). However, when I try to run this using NCBIWWW (Biopython) I don't get any results and I would like to figure out why.

A snippet that contains the call that I am using, is listed here;

```
from Bio.Blast import NCBIWWW
from io import StringIO

sequence=['>some_random_sequence\n','AAGGTTCGGTCCAAATTGAA']

lines = "".join(sequence)
file = StringIO(lines)
results = NCBIWWW.qblast(
    program="blastn",
    database="RefSeq_Gene",
    entrez_query="grch38",
    sequence=file.read(),
    # Manually setting the parameters because short_query checks the length
    # of the sequence, but for some reason includes the name (resulting in
    # a length of 42, and thus defeating the entire purpose of using
    # short_query (which does a if len(sequence) < 31).
    expect=1000,
    word_size=7,
    nucl_reward=1,
    filter=None,
    lcase_mask=None,
    # short_query=True,
)
with open('test_output.xml', 'w') as save_file:
    blast_results = results.read()
    save_file.write(blast_results)

```

# Answer
This is easy. If you set

```
expect=1000,

```

to

```
expect=100,

```

It will work. NCBI servers impose a tighter restriction on the CLI whilst the web restriction is 5000.

If you require a `1000` restriction then out of (US) hours blasting is required, e.g. weekends.

---

The Perl Blast API is here https://blast.ncbi.nlm.nih.gov/docs/web\_blast.pl

The documentation is here https://ncbi.github.io/blast-cloud/dev/api.html

Normally you need to modify it, that would be a separate question. This will work in default mode.
> 2 votes

---
Tags: blast, biopython
---


# Question
Title: error in downloading the annotation of sequences on NCBI
I have developed an R function to download the "organism" name of the numerous NCBI accession numbers which are already saved in an object. The function:

```
get_organism_names <- function(accession_numbers) {
  organism_names <- sapply(accession_numbers, function(accession) {
    record <- entrez_summary(db = "nucleotide", id = accession)
    return(record$organism)
  }, USE.NAMES = FALSE)
  names(organism_names) <- accession_numbers
  return(organism_names)
}

```

Although it works for some accession numbers, for others it encounters an error stating for example:

```
Warning: Esummary includes error message: Invalid uid KU253616.1 at position= 0Error: No esummary records found in file

```

How can I either solve this error or ask R to skip the problematic accession and keep downloading the rest?

# Answer
> 1 votes
You could use the `tryCatch` function to handle inputs that throw an error. A good example on how to write them can be found here. Though the message indicates it's a warning, so does the function stop running when you get this message?

---
Tags: r, ncbi
---


# Question
Title: How to analyse qualitatively the penetration ability of particles in spheroids using fluorescent z-stacks?
0

I have to establish the penetration profile of particles in tumour spheroids.

For this I have spheroids composed of cells which were exposed to particles, both of which are fluorescently labeled. The spheroids were then imaged using a confocal microscope by taking z-stacks of all spheroids. For analysis, I determined the intensity profile of each focal plane using an ImageJ plugin, namely "distribution\_shapes".

What should be done next?

I currently did not find any step-by-step guide for this analysis for manual assessment and any type of help would be greatly appreciated.

# Answer
> 1 votes
I am not sure if this is exactly what you're looking for, but I was able to find this ImageJ script that estimates a penetration profile in tumor spheroids.

I am not an ImageJ person but I suspect that with some minimal changes this could be used to answer your question.

# Answer
> 1 votes
The next step depends heavily on what your question is. Some general ideas, not specific to tumor particle penetration.

1. Your current visualization seems OK to me, if you want to make some features more obvious you can try density estimation.
2. If you have a generative model (for example you believe this is a mixture of a "shallow" and "deep" Gaussian distribution, or you have reason to believe the particles follow an exponential distribution), you can fit your data to this distribution and estimate its parameters.
3. If you want to compare conditions without too strong assumptions on the distributions, you could compute the cumulative distribution of the particles. There are statistical approaches to compare such curves.
4. To compare conditions, you might even be able to come up with a simpler metric, for example "total fluorescence in the core of the spheroid" normalized by "total fluorescence in its periphery".

Ultimately it all depends on your research question (are there treatments to compare to each others, biological replicates whose consistency you want to check, known structure within the tumor spheroid, a time course, ...). My suggestions above would be easier to implement in R/Matlab/Python than pure ImageJ.

---
Tags: 3d-structure
---


title: HISAT2: RNA strandedness

text:
# Question
Title: HISAT2: RNA strandedness
My library is unstranded and the code I'm trying to use is this:

```
# Run alignment with unstranded method
hisat2 -q --rna-strandness none -x HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq | samtools sort -o HISAT2/Kolf-1_R1_001.bam

echo "HISAT2 finished running!"

```

However, I get this error:

```
(base) ~ % RNASeq_pipline/data/script/RNASeqpipline.sh

Error: should be one of F, R, FR, or RF
Error: Encountered internal HISAT2 exception (#1)
Command: /usr/local/bin/hisat2-align-s --wrapper basic-0 -q --rna-strandness none -x
HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq
(ERR): hisat2-align exited with value 1
[W::hts_set_opt] Cannot change block size for this format
samtools sort: failed to read header from "-"
HISAT2 finished running!
0 minutes and 0 seconds elapsed.

```

Looking at the HISAT2 Manual:

For single-end reads:

> Use 'F' if the read corresponds to a transcript.
> Use 'R' if the read corresponds to the reverse complemented counterpart of a transcript.

For paired-end reads:

> Use 'FR' if the first read in the pair corresponds to a transcript on the forward strand, and the second read corresponds to the reverse strand.
> Use 'RF' if the first read in the pair corresponds to a transcript on the reverse strand, and the second read corresponds to the forward strand.
>
> When you use the `--rna-strandness` option with either 'FR' or 'RF' for paired-end reads, HISAT2 will assign an XS attribute tag to each read alignment, indicating whether the read belongs to a transcript on the '+' (plus) or '-' (minus) strand of the genome.

It doesn't say what to input in the case of an unstranded library. Any help would be greatly appreciated!

# Answer
According to the help produced when I type `hisat2 --help`, the default strandedness when running HISAT2 is unstranded:

```
  --rna-strandness <string>          specify strand-specific information (unstranded)

```

So taking out the `rna-strandedness` option should map the reads as if they were unstranded:

```
# Run alignment with unstranded method
hisat2 -q -x HISAT2/grch38/genome_index_grch38 -U hipsci/Kolf-1_R1_001.fastq | samtools sort -o HISAT2/Kolf-1_R1_001.bam

```

If this is a problem for you, I recommend you create an issue about it on the HISAT2 github page:

https://github.com/DaehwanKimLab/hisat2/issues/new
> 2 votes

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


# Question
Title: New to DNABERT confused about kmers input
I'm trying to use DNABERT to predict promoter regions. Right now, I have a df with 500k rows, 1st column is the DNA sequence (all different lengths), 2nd column says "promoter" vs. "non-promoter". I'm trying to understand the GitHub, read the paper, as well as other background reading, but I'm still confused.

If I want to get started, they said it's optional to pre-train (so you can skip to step 3). This is where I got tripped up: "Note that the sequences are in kmer format, so you will need to convert your sequences into that." From what I understand, you need to do this so that all of the sequences are the same length? So kmer=6 means all of the sequences are length 6? Someone suggested that I take the first nucleotide in the promoter and grab 3 nucleotides before and 3 nucleotides after (+/-3 bases). I don't think that's how the kmer thing works though? I tried replicating how I think it works down below (I got confused on the last row of the 'after' df). Please correct me if I'm wrong! Also, if there is a good package that does this for me, please let me know!

Before:

```
seq                          label
AAAGAGGT                      1 [promoter]

```

After:

```
seq                          label
AAAGAG                         1
AAGAGG                         1
AGAGGT                         1
??[GAGGT]??                    1

```

# Answer
The next line in the manual speaks of a python function they provide that can get k-mers from your sequence, so you don't need to worry about doing this a different way. The function slices the sequence starting from the first base with a sliding window of k to obtain k-sized pieces of the sequence and then space-separates these pieces into a single string.

```
seq="ATGCCGTA"
k=2
[seq[x:x+k] for x in range(len(seq)+1-k)]
['AT', 'TG', 'GC', 'CC', 'CG', 'GT', 'TA']

" ".join([seq[x:x+k] for x in range(len(seq)+1-k)])
'AT TG GC CC CG GT TA'

k=3
" ".join([seq[x:x+k] for x in range(len(seq)+1-k)])
'ATG TGC GCC CCG CGT GTA'

```
> 0 votes

---
Tags: python
---


# Question
Title: Getting ALL Protein entries from Entrez for a particular taxonomic id
I'm trying to get ALL Proteins for a Taxonomic ID from Entrez using Esearch/Efetch and the following command:

```
esearch -db "protein" -query "txid374666[Organism]" | efetch -format fasta > txid_374666.fasta

```

Can somebody tell if this command will get ALL proteins or not? Besides, this appears to be a bit slow.

Can somebody point out any Biopython workflow, particularly how to get ALL Proteins

Thanks

# Answer
> 3 votes
As suggested by somebody, the datasets tools from NCBI is blazing fast:

```
datasets download genome taxon 187492  --include protein

```

NCBI Datasets Tool

---
Tags: sequence, entrez, biopython, efetch
---


# Question
Title: How to reverse complement a genbank file?
I am visualizing the alignment of a particular genomic region across different genomes with mauve. In certain genomes the region of interest is in the negative strand and in others the positive strand. What tool can I use to reverse complement my genbank files without losing the features and annotation?

# Answer
> 2 votes
I have not really used Mauve, so I don't know whether it might have an in-built setting for that.

SnapGene Viewer should be able to do that - it can open genbank files and export them, and View \> Flip Sequence gets you the reverse complement without affecting the annotations. SnapGene can also do that; I only specified Viewer because it's free.

# Answer
> 0 votes
From the biology side of bioinformatics and the easier option, I would use Mauve. If you want to make the project more specific and tailored to your needs, and you have python skills, I would write Python code to complete this task. I would specifically use the bioinformatics Python library called Biopython: https://biopython.org pip install biopython

The second option is more involved but it gives you more control on your project/research since you're the one creating it. I personally would go with the 2nd option if you have coding skills. If you don't have coding skills, stick to the 1st option. Im sure there are also tools you can use that are found on Galaxy: https://usegalaxy.org

I am not sure what tools there are that would be good for your task and needs, but from using Galaxy in the past, im sure you can find something! :)

---
Tags: genomics, genbank, mauve
---


# Question
Title: How can I limit blast results to exclude a specific taxonomic ID?
I would like to perform blast+, using the command line, here is the code:

```
-task blastn -db nr -outfmt 7 -query myQuery.fas -out result.txt -remote.

```

I have to exclude SARS-CoV-2 hits from my results, so I must include `-negative_taxids` mentioning SARS-CoV-2 NCBI taxid number. However, the problem is, `-negative_taxids` is NOT compatible with `-remote` command and my command line encounters an error. I have to use `-remote` as I cannot download entire +400 GB NR database to run it locally. Could you please kindly help me out with this contradiction in my work?

My output txt file for this command line is:

```
# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 600483 hits found
AstV-VA1    MN464146.1  100.000 32  0   0   1   32  6362    6393    7.14e-06    59.0
AstV-VA1    MN148428.1  100.000 32  0   0   1   32  6685    6716    7.14e-06    59.0
AstV-VA1    MK671313.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671312.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671311.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0

```

# Answer
> 1 votes
This isn't particularly elegant, but it should work. First, extract the accessions into a file. Next, get the taxid for each of the accessions. Once you have that, you can identify the accessions that are for SARS-CoV-2 sequences, and remove them from your output. Something like this:

1. Extract the accessions into a new file

```
awk '!/^#/{print $2}' blast.out > names

```
2. Map those accessions to their taxids.

```
while read -r accession; do
   printf '%s\t%s\n' "$accession"  \
       $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
          awk -F'[<>]' '/Name="TaxId"/{print $3}');
done < names > taxids

```

With your example input, that produces:

```
$ while read -r accession; do
>    printf '%s\t%s\n' "$accession"  \
>        $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
>           awk -F'[<>]' '/Name="TaxId"/{print $3}');
> done < names
MN464146.1   1435485
MN148428.1   2661902
MK671313.1   1239566
MK671312.1   1239566
MK671311.1   1239566

```

Now, this will probably be way too slow, doing them one at a time. So if you have many IDs, try and combine multiple queries into one. Like this:

```
perl -ne 'chomp; print $. % 50 == 0 ? "$_\n" : "$_," ' names |
 | sed 's/,$/\n/' |
    while read accession; do
      curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" |
        awk -F'[<>]' '/Name="Caption"/{ acc=$3 }
                   /Name="TaxId"/{printf "%s\t%s\n", acc, $3}'
    done  < names > taxids

```

Here, the first Perl command will read through the list of accessions and print the accession followed by a comma and, once every 50 lines (`$. % 50 == 0`), it will instead print it with a newline. The `sed` command will replace a trailing comma with a newline for the very last line of output where, unless the number of accessions is an exact multiple of 50, you will have a line ending with a comma instead of a newline.

The result is a list of 50 comma-separated identifiers per line (fewer, on the last line). This is then fed to the `curl` command and then a naive awk that parses it and prints out the accession and the taxid.

**Important**: this approach loses the version information. So `MK671311.1` becomes `MK671311`. This shouldn't be a problem, and we will take it into account later, but it should be understood. This is what the output looks like:

```
MN464146    1435485
MN148428    2661902
MK671313    1239566
MK671312    1239566
MK671311    1239566

```
3. Remove the SARS-CoV-2 sequences.

The taxid for SARS-CoV-2 is 2697049, so you can get the list of its accessions with:

```
awk -F'\t' '$2=="2697049"{print $1}' taxids

```

You can now filter your file to remove them using:

```
awk '(NR==FNR){
       if($2=="2697049"){
         sars[$1]
       }
       next
      }
      {
        acc=$2;
        sub(/\.[0-9]*$/, "", acc);
        if(acc in sars){ next }
          print
      }' taxids blast.out > blast.filtered.out

```

This awk script will read the first input file (`NR` is the current input line number while `FNR` the line number of the file being read; the two are only equal while reading the first file) and if its second field is the desired taxid, it saves the first field, the accession, in the array `sars`. The `next` ensures we don't do anything else while reading the first file.

When we move to the second file, we make a copy of the second field, the accession, and save it as `acc`, and then we remove the version information ( `sub(/\.[0-9]*$/, "", acc)`. If the resulting accession is found in the `sars` array, we skip the line with `next` and if not, we print.

This assumes you are using the faster, batch retrieval I gave above, the one that removes the versioning information. If you instead use the very slow first one, you can filter with

```
 awk '(NR==FNR){ if($2=="2697049"){ sars[$1]} next }
      { if($2 in sars){next} print}' taxids blast.out > blast.filtered.out

```

That should give you a filtered blast output file that no loner has sequences from taxid 2697049. **However**, please make sure to test this, count the number of input and output lines and make sure everything makes sense. I don't have access to your actual data and while this does work with the example you provided, I cannot guarantee it will work out of the box for your real data. As always, you need to be careful!

# Answer
> 0 votes
I would personally look at is the version number of the blast you are using.

Bugs around blast+ versions \>= 2.9.0 (April 2019) for `-taxids` and I would assume `-negative_taxids` flag should have been fixed as described in the the Blast+ release notes here.

> Fix taxID filtering combined with mask-based alias BLAST databases.

Bug fix descriptions are usually short, hence the lack of further detail.

Blast `taxids` can be implemented via lists here https://www.ncbi.nlm.nih.gov/books/NBK569846/

---
Tags: blast, sars-cov-2, ncbi, command-line
---


# Question
Title: Downloading a list of all RCSB-PDB entries
I'm trying to download a list of all of the non-redundant protein chain IDs in the RSCB-PDB (not the actual FASTA entries, just the IDs), filtered by some % identity.

Does anyone know how to do this? I haven't been able to figure it out with the advanced search function .

# Answer
Did you try RCSB sequence similarity search? These two announcement highlight the options you ask about:

Alternatively, if you'd like to calculate identity on your own, you can download sequences of all the PDB entries from: https://files.wwpdb.org/pub/pdb/derived\_data/pdb\_seqres.txt (270MB).
> 1 votes

---
Tags: proteins, protein-structure, pdb
---


# Question
Title: Annotation result for phylogenetic analysis shows no common evolutionary gene in contigs, can i change the contigs?
I'm currently doing my thesis with the topic of phylogenetic analysis and is taking references from the previous person in my university (who have done the same topics but different species).

They used the largest contigs comprising well known markers such as matK, rbcL, and ITS. However my largest contigs doesn't have them.

Is it possible for me to use the 2nd/3rd largest contigs that have those markers?

# Answer
ITS: what is needed is the aim of the phylogenetic investigation. ITS evolves very quickly so it's good for a few closely related species and sub-species. Understanding the behaviour of a locus is as important as the volume of the locus analysed. Thus, the answer is yes.

I have a feeling these are insects/arthropods simply because of the fragmented nature of the genome: they're hard to assemble and genetically diverse so there's often no reference genome. The issue of genome vs locus was addressed by David Penny years ago: genome trees can be misleading and that would be justification for using the contigs to derive select loci. The reason why genome trees can be misleading is complex and a separate question.

Moreover, if you can map the loci to a chromosome and demonstrate the contigs are on different chromosomes you can use differences/similarities in trees between chromosomes as the first stage of disequilibrium analysis, i.e. selection acting via Mendelian exchange.

To reiterate there is good phylogenetic justification for switching contigs because the loci are more important than the size of the contig (you'll need supporting references however).
> 1 votes

---
Tags: phylogenetics, phylogeny, genome, sequence-annotation, scaffold
---


# Question
Title: installing and using differential methylation DSS software
I am attempting to install and use the Bioconductor package DSS on Ubuntu 22.04.

Following the instructions here: https://bioconductor.org/packages/release/bioc/html/DSS.html

I get to this error:

```
hfile_libcurl.c:45:10: fatal error: curl/curl.h: No such file or directory
   45 | #include <curl/curl.h>
      |          ^~~~~~~~~~~~~
compilation terminated.
make[1]: *** [Makefile.Rhtslib:128: hfile_libcurl.o] Error 1
make[1]: Leaving directory '/tmp/RtmpNry6uq/R.INSTALL10673560786c68/Rhtslib/src/htslib-1.7'
make: *** [Makevars.common:24: htslib] Error 2

```

# Answer
In Ubuntu 22.04, install:

```
sudo apt install libcurl4-openssl-dev

```

Then reattempt:

```
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("DSS")

```

Then upon completion:

```
library(DSS)
Loading required package: Biobase
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs
...

```
> 1 votes

---
Tags: dss, differential-methylation
---


# Question
Title: VEP annotation INFO field Ensembl IDs and locations
I have a vcf file that I annoteted with VEP, for human data. I have run VEP to annotate my files with some additional parameters (as shown below in the ##VEP-command-line).

However, my output is rather strange (mainly the INFO column).

```
##VEP="v108" time="2023-04-27 15:13:08" cache="workflow/resources/variants/cache_vep/homo_sapiens/108_GRCh38" ensembl-funcgen=108.56bb136 ensembl-variation=108.a885ada ensembl-io=108.58d13c1 ensembl=108.d8a9c80 1000genomes="phase3" COSMIC="96" ClinVar="202205" HGMD-PUBLIC="20204" assembly="GRCh38.p13" dbSNP="154" gencode="GENCODE 42" genebuild="2014-07" gnomADe="r2.1.1" gnomADg="v3.1.2" polyphen="2.2.2" regbuild="1.0" sift="sift5.2.2"
##INFO=<ID=CSQ,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|Existing_variation|DISTANCE|STRAND|FLAGS|SYMBOL_SOURCE|HGNC_ID|GENE_PHENO|SIFT|PolyPhen|AF|AFR_AF|AMR_AF|EAS_AF|EUR_AF|SAS_AF|gnomADe_AF|gnomADe_AFR_AF|gnomADe_AMR_AF|gnomADe_ASJ_AF|gnomADe_EAS_AF|gnomADe_FIN_AF|gnomADe_NFE_AF|gnomADe_OTH_AF|gnomADe_SAS_AF|CLIN_SIG|SOMATIC|PHENO">
##VEP-command-line='vep --dir_cache workflow/resources/variants/cache_vep --vcf -i WESP_run/results/variant_calling/VP_DBS_25_7_S7.vcf -o WESP_run/results/variant_annotation/VP_DBS_25_7_S7.annotated.vcf.gz --af --af_1kg --af_gnomade --sift b --polyphen b --gene_phenotype --fork 24 --offline --compress_output gzip'
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  VP_DBS_25_7_S7
chr1    13273   .   G   C   25.3    PASS    CSQ=C|non_coding_transcript_exon_variant|MODIFIER|DDX11L1|ENSG00000223972|Transcript|ENST00000450305|transcribed_unprocessed_pseudogene|5/6||||313|||||rs531730856||1||HGNC|HGNC:37102||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|non_coding_transcript_exon_variant|MODIFIER|DDX11L2|ENSG00000290825|Transcript|ENST00000456328|lncRNA|3/3||||521|||||rs531730856||1||EntrezGene|||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|downstream_gene_variant|MODIFIER|WASH7P|ENSG00000227232|Transcript|ENST00000488147|unprocessed_pseudogene||||||||||rs531730856|1131|-1||HGNC|HGNC:38034||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401||||||||||||,C|downstream_gene_variant|MODIFIER|MIR6859-1|ENSG00000278267|Transcript|ENST00000619216|miRNA||||||||||rs531730856|4096|-1||HGNC|HGNC:50039||||0.0950|0.0204|0.1455|0.0625|0.1471|0.1401|||||||||||| GT:GQ:DP:AD:VAF:PL  0/1:23:55:28,27:0.490909:25,0,27

```

For instance, at the first location/variant I look at the Ensembl gene id and see multiple genes on the same locations and strand. This is very strange to me, as I would understand different transcripts, but different genes on the same strand? Also, apart from that, the corresponding locations/regions of genes on the Ensembl website do not agree with the location in the vcf file. For instance, the location for this SNP is 13273, but the annotations are (corresponding ensembl archive version 108):

What am I missing here?

Kind regards, Eva

# Answer
> 1 votes
You need to set your expectations clearly instead of just saying "Does not agree". I'm taking a guess here - your statement is that the 2/4 of the gene co-ordinates do not overlap the variant location.

This is because you're viewing just part of the picture - the co-ordinates - and not seeing another important part, which is the kind of effect. For the latter two genes, VEP clearly says the effect is `downstream_gene_variant`, which Sequence Ontology defines as "A sequence variant located 3' of a gene."

Putting those two pieces together, the entire picture makes sense. Variants affect different transcripts (of possibly different genes) in different manners. Looking at just 1 of those 2 components will of course make things seem odd.

---
Tags: vcf, ensembl, identifiers, vep
---


# Question
Title: How to run a GATK Docker Image with local files?
I'm trying to use the HaplotypeCaller from the GATK toolkit but I keep getting an error. I pulled GATK through Docker and am using this command:

```
docker run -v /Users/rimo/ -it broadinstitute/gatk:latest gatk HaplotypeCaller -R /Users/rimo/reference.fasta -I /Users/rimo/sample1.bam -O /Users/rimo/sample1.g.vcf.gz -ERC GVCF

```

`/Users/rimo` is my home directory it's where the files are currently located on my computer and it's where I want GATK to look for files to use in the command.

but then I get this error:

```
A USER ERROR has occurred: The specified fasta file (file:///Users/rimo/reference.fasta) does not exist.

```

I don't fully understand how Docker commands work but I found a similar problem listed with this issue. I tried to use their solution but it's still giving me the same error as listed above.

# Answer
Read the manual please: https://docs.docker.com/storage/bind-mounts/

> `-v` or `--volume`: Consists of three fields, separated by colon characters (:). The fields must be in the correct order, and the meaning of each field is not immediately obvious.
>
> * In the case of bind mounts, the first field is the path to the file or directory on the host machine.
> * The second field is the path where the file or directory is mounted in the container.
> * The third field is optional, and is a comma-separated list of options, such as ro, z, and Z. These options are discussed below.

Try

```
docker run -v /Users/rimo/:/myhome ...

```

Also, why are you combining a `-it` run with an actual command? I thought the `-it` is to start an interactive shell, which is the opposite of running a single command.
> 3 votes

# Answer
I went back to the documentation listed on the GATK website instead of Docker because the Docker documentation was a bit confusing. And was able to figure it out with this command:

```
docker run -v ~/data:/gatk/my_data -it broadinstitute/gatk:latest

```

And for those who might use this as a resource because they don't understand how to use local files on a Docker image the breakdown is as follows:
> 0 votes

---
Tags: variant-calling, gatk, docker, command-line, troubleshooting
---


title: Mendelian randomization

text:
# Question
Title: Mendelian randomization
I would be thankful to you if you can help with how I can use MR from my case-control study to link DNA methylation (epic array) with toxic element exposure (arsenic) and health outcome (Cardiovascular diseases).

I do have data from 125 Cardiovascular diseases and 125 non-Cardiovascular diseases individuals matched with age, sex, and smoking.I do have phenotypes files in a csv format having sex, age, smoking, disease, samples id, iAs (arsenic) values in log format and then a separate files for methylation data. we also have SNP data for these samples.

Many thanks,

# Answer
I understand linkage disequilibrium (LD) versus panmixia, which the technical term for Mendelian randomisation and for humans PLINK is a population package. The best outcome is a correlation between LD and cardiovascular disease. It is however an unlikely outcome.

This really is a machine learning question and you would need to understand "feature selection" and will require very high accuracy scores \>0.95 - however you could triple the size of our control group quite easily (recommended). This is a separate question and I'm happy to answer this there.

---

Just to address the comments. Inheritance probability, if that is the intended statistics, I personally don't think it would work because a 125 sample is low and the patients would need to be in cohorts of familial generations. LD is unlikely the explanation with humans but other non-humans species that might work (selection advantage).

My personal view is this is a straight machine learning calculation, but it's a separate question. If your sample is heavily stratified across each familial generation then Hardy-Weinberg might work, but my experience of patient cohorts is that is really unusual and there could still be confounders. Supervised learning (subset of ML) is 100% the way forward.

.... possibly use GWAS, but you've multiple features/training targets. GWAS is usually when there's one phenotype, ML would account for interaction in a way GWAS could not (it couldn't account for non-sequence features, i.e. mixed data).
> 1 votes

---
Tags: r, statistics, epigenetics
---


# Question
Title: How many Ns and ns in GRCh37 / GRCh38 per 'canonical' chromosome?
This is kind of pedantic, but I'm not sure where to look...

For GRCh38 (and a lot of work...) I have the following...

```
Chr          Length         Ns   ns
chr1    248,956,422 18,475,229  181
chr2    242,193,529   1645,291   10
chr3    198,295,559    195,420    4
chr4    190,214,555    461,888    0
chr5    181,538,259    272,881    0
chr6    170,805,979    727,255  202
chr7    159,345,973    375,841    1
chr8    145,138,636    370,500    0
chr9    138,394,717 16,604,127   40
chr10   133,797,422    534,287  173
chr11   135,086,622    552,880    0
chr12   133,275,309    137,270  223
chr13   114,364,328 163,81,200    3
chr14   107,043,718 164,75,569    0
chr15   101,991,189 173,49,841   23
chr16    90,338,345   8532,401    1
chr17    83,257,441    336,367  870
chr18    80,373,285    283,627   53
chr19    58,617,616    176,858    0
chr20    64,444,167    499,586  324
chr21    46,709,983  6,621,300   64
chr22    50,818,468 11,658,691    0
chrX    156,040,895  1,147,765  101
chrY     57,227,415 30,812,232  140
chrM         16,569          1    0

```

(Where) could I have looked this up?

What are the counts for GRCh37?

# Answer
> 1 votes
I don't know where you could look this up and of course this will depend on the versions of the relevant genomes and on which exact assemblies you are working with. I doubt you would have found this specific number, how many Ns anywhere, but it's really not very hard to get them. For instance, I have hg19 and hg38 versions locally:

```
$ for genome in hg19 hg38 ; do ls $genome/$genome.fa; done
hg19/hg19.fa
hg38/hg38.fa

```

This means I can get the Ns for both using:

```
for genome in hg19 hg38 ; do
  FastaToTbl "$genome"/$genome.fa |
    perl -ne 'print if /^(chr)?[0-9MXY]*\s/' |
      awk -F"\t" '{ print $1, split($2,a,/[nN]/)}' > $genome.ns
done

```

And then combine into a nice table with:

```
$ printf 'chr\thg19\thg38\n-----------------------\n';
  join -j1 -t ' ' <( sed 's/chr//' hg19 | sort -gk1,1 ) \
                  <( sed 's/chr//'  hg38 | sort -gk1,1) |
                   column -t | sed 's/^/chr/'
chr hg19    hg38
-----------------------
chrM   1         2
chrX   4170001   1147867
chrY   33720001  30812373
chr1   23970001  18475411
chr2   4994856   1645302
chr3   3225296   195425
chr4   3492601   461889
chr5   3220001   272882
chr6   3720002   727458
chr7   3785001   375843
chr8   3475101   370501
chr9   21070001  16604168
chr10  4220010   534461
chr11  3877001   552881
chr12  3370503   137494
chr13  19580001  16381204
chr14  19060001  16475570
chr15  20836627  17349865
chr16  11470001  8532403
chr17  3400001   337238
chr18  3420020   283681
chr19  3320001   176859
chr20  3520001   499911
chr21  13023254  6621365
chr22  16410022  11658692

```

The fancy sed stuff is just to ensure that sorting works correctly by removing and then re-adding the `chr` string. Note that this is counting both `N` and `n`. The FastaToTbl script can be found here.

---
Tags: human-genome, repeat, grch37
---


# Question
Title: What is the exact definition of "breakend"?
I'm reading the Manta User Guide → it is a method to discover structural variants and indels from next generation sequencing data.

The word 'breakend' is used in several contexts

I have searched for a definition in `duckduckgo`, `qwant` and `google` but have not found a concise definition.

Thanks for you help

# Answer
A structural variant is an inversion, deletion, duplication, or translocation. A "breakend" is one endpoint of a structural variant. This paper (Alkan et al. *Nature Review Genetics* 2011) may help you understand: https://www.nature.com/articles/nrg2958

The Manta paper and documentation seem sparse on diagrams. DELLY, another structural variant caller, has diagrams in their publication. Their paper may be more of use: https://academic.oup.com/bioinformatics/article/28/18/i333/245403

**Specific example:** Say you're sequencing a chronic myeloid leukemia (CML) tumor. You will likely find reads that half map to the ABL1 gene on chr9 and half map to the BCR gene on chr22 (https://en.wikipedia.org/wiki/Philadelphia\_chromosome). The genomic position of the location where chr9 turns to chr22 is called the breakpoint. The positions on chr9 and chr22 are the breakends.
> 3 votes

# Answer
I think the term "breakend" is derived specifically from the VCF file spec:

https://samtools.github.io/hts-specs/VCFv4.2.pdf

Starting on page 12:

"5.4 Specifying complex rearrangements with breakends

An arbitrary rearrangement event can be summarized as a set of novel adjacencies. Each adjacency ties together 2 breakends. The two breakends at either end of a novel adjacency are called mates.

There is one line of VCF (i.e. one record) for each of the two breakends in a novel adjacency. A breakend record is identified with the tag “SVTYPE=BND” in the INFO field. The REF field of a breakend record indicates a base or sequence s of bases beginning at position POS, as in all VCF records. The ALT field of a breakend record indicates a replacement for s. This “breakend replacement” has three parts:

1. The string t that replaces places s. The string t may be an extended version of s if some novel bases are inserted during the formation of the novel adjacency.
2. The position p of the mate breakend, indicated by a string of the form “chr:pos”. This is the location of the first mapped base in the piece being joined at this novel adjacency.
3. The direction that the joined sequence continues in, starting from p. This is indicated by the orientation of square brackets surrounding p. These 3 elements are combined in 4 possible ways to create the ALT. In each of the 4 cases, the assertion is that s is replaced with t, and then some piece starting at position p is joined to t. The cases are:

REF ALT Meaning

s t\[p\[ piece extending to the right of p is joined after t

s t\]p\] reverse comp piece extending left of p is joined after t

s \]p\]t piece extending to the left of p is joined before t

s \[p\[t reverse comp piece extending right of p is joined before t

The example in Figure 1 shows a 3-break operation involving 6 breakends. It exemplifies all possible orientations of breakends in adjacencies. Notice how the ALT field expresses the orientation of the breakends.

#CHROM POS ID REF ALT QUAL FILTER INFO

2 321681 bnd W G G\]17:198982\] 6 PASS SVTYPE=BND

2 321682 bnd V T \]13:123456\]T 6 PASS SVTYPE=BND

13 123456 bnd U C C\[2:321682\[ 6 PASS SVTYPE=BND

13 123457 bnd X A \[17:198983\[A 6 PASS SVTYPE=BND

17 198982 bnd Y A A\]2:321681\] 6 PASS SVTYPE=BND

17 198983 bnd Z C \[13:123457\[C 6 PASS SVTYPE=BND"
> 1 votes

---
Tags: structural-variation
---


# Question
Title: How does the BWA-MEM algorithm assign its mapping qualities?
Is there any resource (paper, blogpost, Github gist, etc.) describing the BWA-MEM algorithm for assigning mapping qualities? I vaguely remember that I have somewhere seen a formula for SE reads, which looked like

$C * (s\_1 - s\_2) / s\_1,$

where $s\_1$ and $s\_2$ denoted the alignment scores of two best alignments and `C` was some constant.

I believe that a reimplementation of this algorithm in some scripting language could be very useful for the bioinfo community. For instance, I sometimes test various mapping methods and some of them tend to find good alignments, but fail in assigning appropriate qualities. Therefore, I would like to re-assign all the mapping qualities in a SAM file with the BWA-MEM algorithm.

Btw. This algorithm must already have been implemented outside BWA, see the BWA-MEM paper:

> GEM does not compute mapping quality. Its mapping quality is estimated with a BWA-like algorithm with suboptimal alignments available.

Unfortunately, the BWA-MEM paper repo contains only the resulting `.eval` files.

**Update:** The question is *not* about the algorithm for computing alignment scores. Mapping qualities and alignment scores are two different things:

* Alignment score quantifies the similarity between two sequences (e.g., a read and a reference sequence)
* Mapping quality (MAQ) quantifies the probability that a read is aligned to a wrong position.

Even alignments with high scores can have a very low mapping quality.

# Answer
> 3 votes
Yes, there bwa-mem was published as a preprint

> BWA-MEM’s seed extension differs from the standard seed extension in two aspects. Firstly, suppose at a certain extension step we come to reference position x with the best extension score achieved at query position y.
>
> ...
>
> Secondly, while extending a seed, BWA-MEM tries to keep track of the best extension score reaching the end of the query sequence

And there is a description of the scoring algorithm directly in the source code of bwa-mem (lines 22 - 44), but maybe the only solution is really to go though the source code.

# Answer
> 1 votes
The MAPQ score appears to be calculated in the mem\_approx\_mapq\_se function. I tried porting it to Python and inserting comments using ChatGPT 4, though I haven't tested it to ensure the output is identical to the original.

import math

# Constants
MEM_MAPQ_COEFFICIENT = 30.0
MEM_MAPQ_MAXIMUM = 60

def calculate_mapping_quality(alignment_options, alignment_region):
    # Determine the suboptimal alignment score, which is the second-best SW alignment score if available,
    # otherwise calculated from minimum seed length and match score.
    suboptimal_alignment_score = (alignment_region.suboptimal_score if alignment_region.suboptimal_score
                                  else alignment_options.minimum_seed_length * alignment_options.match_score)

    # Update the suboptimal alignment score to the higher of the two: current suboptimal score or tandem hit score.
    suboptimal_alignment_score = (alignment_region.tandem_hit_score if alignment_region.tandem_hit_score > suboptimal_alignment_score
                                  else suboptimal_alignment_score)

    # Early exit with zero mapping quality if the suboptimal score is greater than or equal to the alignment score.
    if suboptimal_alignment_score >= alignment_region.best_alignment_score:
        return 0

    # Calculate the maximum alignment length between the query and reference sequences.
    alignment_length = max(alignment_region.query_end - alignment_region.query_start,
                           alignment_region.reference_end - alignment_region.reference_start)

    # Calculate the identity of the alignment.
    identity = (1.0 - (alignment_length * alignment_options.match_score - alignment_region.best_alignment_score)
                / (alignment_options.match_score + alignment_options.mismatch_penalty) / alignment_length)

    # Compute the mapping quality based on the alignment score, identity, and other factors.
    if alignment_region.best_alignment_score == 0:
        mapping_quality = 0
    elif alignment_options.mapQ_coef_length > 0:
        # Adjust mapping quality based on alignment length and identity.
        length_factor = (1.0 if alignment_length < alignment_options.mapQ_coef_length
                         else alignment_options.mapQ_coef_fac / math.log(alignment_length))
        length_factor *= identity * identity
        mapping_quality = int(6.02 * (alignment_region.best_alignment_score - suboptimal_alignment_score)
                              / alignment_options.match_score * length_factor * length_factor + 0.499)
    else:
        # Use the constant MEM_MAPQ_COEFFICIENT for mapping quality calculation.
        mapping_quality = int(MEM_MAPQ_COEFFICIENT * (1.0 - suboptimal_alignment_score / alignment_region.best_alignment_score)
                              * math.log(alignment_region.seed_coverage) + 0.499)
        mapping_quality = int(mapping_quality * identity * identity + 0.499) if identity < 0.95 else mapping_quality

    # Penalize the mapping quality for the presence of suboptimal hits.
    if alignment_region.suboptimal_hit_count > 0:
        mapping_quality -= int(4.343 * math.log(alignment_region.suboptimal_hit_count + 1) + 0.499)

    # Ensure the mapping quality does not exceed the maximum allowed value.
    mapping_quality = min(max(mapping_quality, 0), MEM_MAPQ_MAXIMUM)

    # Adjust the mapping quality based on the fraction of repetitive regions in the alignment.
    mapping_quality = int(mapping_quality * (1.0 - alignment_region.repetitive_region_fraction) + 0.499)

    # Return the final mapping quality value.
    return mapping_quality

```

---
Tags: sequence-alignment, bwa, read-mapping
---


# Question
Title: Intepreting and applying ordinal logistic regression coefficients to calculate probabilities?
Can someone help hint me how I can interpret ordinal logistic regression coefficients and how I can use the .L, .Q and .C terms to calculate probabilities? I am analysing a dataset, where people complain of chest or not pain when they receive external supplements for a vitamin X or not. so I performed a logistic regression, using categorized and ordered (so it is an ordinal categorical variable) values of vitamin X from these observations. So I have grouped these vitamin X values within my data as very low, sub-optimal, optimal and high values and ordered that way.

The data are as such:

chest\_pain\_status: binary, 0 means not having chest pain, 1 means having chest pain

suppl = binary, 0 means not supplementing and 1 means supplementing

vitx\_levels: very\_low \< suboptimal \< optimal\< high

My approach is as follows:

```
mod <- glm(chest_pain_status ~ suppl*vitx_levels, dataset = data,family = binomial())

summary(mod)
                        Estimate       std.Error     z value          Pr( <|z|)
(intercept)             -2.7816        0.2993        -9.2927          0.000
suppl1                   0.9822        0.3209         3.0803          0.0021
vitx_levels.L           -0.1751        0.0634        -2.7391          0.0062
vitx_levels.Q           -0.1537        0.0552        -2.7689          0.0053
vitx_levels.C            0.0164        0.0460         0.3579          0.7203
suppl1:vitx_levels.L     0.16821       0.0739         2.7472          0.0229
suppl1:vitx_levels.Q     0.2433        0.0641         3.7932          0.00015
suppl1:vitx_levels.C    -0.0096        0.0537        -0.1787          0.8582

```

**Question1**: Please how do I interpret coefficients of `vitx_levels.L`, `vitx_levels.Q` and `vitx_levels.C`? From my googling I can understand that R tries to test linear, quadratic and cubic relationships between `vitx_levels` and the log odds ratios of developing chest pain. It therefore means that there is a significant linear and negative relationship between `vitx_levels` and `chest_pain`. Same for the quadratic term. I hope I interpreted it right? Any additional ideas?

**Question2**: What values could I substitute in to `vitx_levels.L`, `vitx_levels.Q` and `vitx_levels.C` If I wanted to create a formula like the one below where I will later transform y to get probabilities:

```
y <- -2.7816 + 0.9822 * suppl1 - 0.1751 * vitx_levels.L - 0.1537 * vitx_levels.Q + 0.0164 * vitx_levels.C +  0.16821 * suppl1:vitx_levels.L + 0.2433 * suppl1:vitx_levels.Q - 0.0096 * suppl1:vitx_levels.C

```

So for the function above, what values could I fit into `vitx_levels.L`, `vitx_levels.Q` and `vitx_levels.C` to be able to compute the value of y?

# Answer
> 1 votes
Its a while since I've done GLM because ML is the method of choice.

**Question1** Yes multi-linear regression is an approximate description of GLM.

*Generally* the weights simply mean whether a feature is positively or negatively associated with the phenomena. It's important to be sure it's the right way around (negatives and positives can flip, if that makes sense).

`suppl1` is significantly and relatively strongly associated with symptoms, whilst `L` and `Q` `vitamin` levels appear to have some *protective* effects. `Suppl1` has a dose dependent effect with `L` and `Q` having a lower correlation with the observed symptoms.

The dose dependency result therefore confirm suppl1 as associated with the symptoms. In all cases levels `C` are not significant, I assume this is a low level supplement.

**Summary** `vitx_levels` appear to have some protective effect.

I do stress it is important to ensure the "negative versus positive" is the right way round. The residual is important and I suspect you've omitted this in the GLM. The lower the residual the better the fit of the model.

BTW LD is different from GLM although they are related (title versus body of the question).

**Question 2** Yes that is correct `0.9822 * 1` I assume you mean `0.9822 * suppl1` ... However all `C` weights are not significant, so you could omit them.

---

residual versus intercept ... I'm pretty certain that these are the same thing (its the *c* variable) and the intercept here appears quite large.

---
Tags: r, statistics, ggplot2, machine-learning
---


# Question
Title: What is the best Query to retrieve DNA from NCBI?
I want to retrieve a sequence for many species from the Nucleotide database in NCBI.

I'm using a command line approach and I have to figure out what is the best query that will return exclusively the region of DNA I am interested in and filter out unwanted noisy sequences.

I am using a wide range of species (that I have stored as a list of TaxIDs of about 2000 species): it includes small crustaceans, invertebrates, some algae and small vertebrates (reptiles, amphibians and few mammals).

My final goal is to obtain a phylogenetic tree for all or most of the species.

I have been suggested to use these genes to create such phylogeny

* CO1 mithocondrial DNA
* 16S rRNA
* 18S rRNA

I want to formulate a specific query that will return exclusively those sequences.

I'm using the GenBank query builder to visually check for the accuracy of my search and when I find a good query I will use it in the API.

So far I came up with the following queries:

* `(COX1[Title] OR CO1[Title]) AND complete[Title]` --\> 63/2000 Species
* `16S[title] AND complete[title] AND rRNA[title] NOT partial[title]` --\> 9/2000 species
* `18S[title] AND complete[title] AND rRNA[title] NOT partial[title]` --\> 15/2000 species

As you can see the number of species that I get sequence for is very low compared to the initial 2000 species. I doubt that we have so little available sequences (especially for COX1 that is used for barcoding)

Can you help me understanding whether my queries are good or not? And if possible suggest a better alternative

### More Info

A subset of 10 of my species of interest is

```
Rasbora heteromorpha
Elasmopus rapax
Gasterosteus aculeatus
Palaemonetes pugio
Catostomus commersoni
Daphnia magna
Oryzias latipes
Xenopus laevis
Tigriopus japonicus
Oncorhynchus mykiss

```

Out of these 10 only 2 have sequences available using the COX1 query

| Name | SeqID |
| --- | --- |
| Oncorhynchus mykiss | EU186789.1 |
| Xenopus laevis | AB278691.1 |

But Daphnia magna is one of the most commonly used organisms in lab and I found this paper regarding the complete mitogenome of a specific strain of D. magna thus implying that the complete mitogenome of the species is already known. That means that there must be a way to retrieve all the other species that have full mitogenome but no COX1 gene mentioned in their GenBank title.

# Answer
> 1 votes
The search term you need is **\[gene name\]** rather than *\[title\]*. For example,

```
(Gasterosteus aculeatus[Organism] AND cox1[Gene Name]) OR (Gasterosteus aculeatus[Organism] AND coi[Gene Name])

```

Results in **266** hits which all contain COX1 COI

Using Daphnia magna ...

```
(Daphnia magna[Organism] AND cox1[Gene Name]) OR (Daphnia magna[Organism] AND coi[Gene Name])

```

This results in **775** hits.

*Daphina magna* comprises **61** mtDNA genomes, you perform the above search and sort by sequence length and set the sequences per page to 200(easy).

Hope that resolves the issue.

# Answer
> 2 votes
There are two main problems I can see in your approach. First, gene names are not standardized across species and the "same" gene (what that means is a whole different discussion) can have different names. For example, human "Cytochrome c oxidase assembly factor 1" has the symbol "COA1".

Next, I don't quite understand why you expect the word "complete" to be in the title. That is actually why you don't find the *D. magna* sequence: it isn't complete, the start codon hasn't been correctly annotated as you can see on the relevant NCBI gene page:

I would instead try something like this:

```
((chrM[Chromosome]) OR (chromosome M[Chromosome]) OR (MT[Chromosome]) OR (chromosome MT[Chromosome])) AND ((CO1[Gene Name]) OR (COX1[Gene Name]) OR ((cytochrome[Gene Name]) AND (oxidase[Gene Name]))

```

or

```
((16S[Gene/Protein Name]) and (ribosomal[Gene/Protein Name]))

```

All that said, trying to rely on gene annotations is rarely the best approach. There just isn't enough standardization. I would instead look at databases that provide lists of homologous genes. For instance, NCBI's Clusters of Orthologous Genes (COGs), which provides grouped homologous sequences.

---
Tags: phylogenetics, phylogeny, sequence-annotation, genbank, barcode
---


# Question
Title: Filtering columns from large file based on a column list
I have list of samples as a txt file and want to only keep these samples from this large file; column name is samples name in the large file.

```
library(readr)
library(data.table)

df1 <- read_csv("control.csv", col_types = cols(Chromosome = col_skip(), Start = col_skip(), End = col_skip(), Strand = col_skip()))

df1[1:5,1:5]
# A tibble: 5 × 5
  ID          DN54 DC146 DDN44 DC179
  <chr>      <dbl> <dbl> <dbl> <dbl>
1 cg26928153 0.848 0.879 0.872 0.897
2 cg16269199 0.731 0.732 0.759 0.775
3 cg13869341 0.842 0.863 0.819 0.896
4 cg24669183 0.824 0.866 0.822 0.910
5 cg26679879 0.425 0.363 0.420 0.358

```

# Answer
You can save time and memory by preprocessing your file before loading it into R.

Use `miller` toolkit to manipulate csv or tab-delimited files using column names.

Below is a one-liner that reads csv file, and writes columns specified in option `-f`, in the order (`-o`) listed on the command line:

```
mlr --csv cut -o -f col3,col1 in.csv > out.csv

```

To read the list of columns from a file, 1 column name per line:

```
mlr --csv cut -o -f $(paste -s -d, columns.txt) in.csv > out.csv

```

**See also:**

Note that `miller` can be easily installed, for example, using `conda`, specifically `miniconda`, like so:

```
conda create --channel conda-forge --name miller miller

```
> 3 votes

# Answer
It's really just `count[,col$V1]`. Please learn about basic subsetting of columns and rows, it's simple yet powerful.
> 3 votes

# Answer
Below code works perfectly of me.

```
library(readr)
library(data.table)

df1 <- read_csv("control.csv", col_types = cols(Chromosome = col_skip(), Start = col_skip(), End = col_skip(), Strand = col_skip()))

df1[1:5,1:5]
# A tibble: 5 × 5
  ID          DN54 DC146 DDN44 DC179
  <chr>      <dbl> <dbl> <dbl> <dbl>
1 cg26928153 0.848 0.879 0.872 0.897
2 cg16269199 0.731 0.732 0.759 0.775
3 cg13869341 0.842 0.863 0.819 0.896
4 cg24669183 0.824 0.866 0.822 0.910
5 cg26679879 0.425 0.363 0.420 0.358

col <- unlist(fread("samples.csv", header=F))
col <- c("ID", col)
extracted <- df1[,col, with=F]

write.table(extracted, file="control.beta.txt",sep=",", quote=F, col.names=TRUE,row.names=FALSE)

df <- read_csv("control.beta.txt")
df[1:5,1:5]
# A tibble: 5 × 5
  ID         DC103R DC110R DC115 DC120R
  <chr>       <dbl>  <dbl> <dbl>  <dbl>
1 cg26928153  0.929  0.885 0.862  0.877
2 cg16269199  0.836  0.826 0.604  0.739
3 cg13869341  0.826  0.903 0.830  0.841
4 cg24669183  0.767  0.719 0.808  0.812
5 cg26679879  0.388  0.366 0.428  0.360

```
> 1 votes

---
Tags: r, bash
---


# Question
Title: retrieving past BLAST requests in bioPython
1. is there any way to retrieve a BLAST request made using the Biopython library using the request's query\_id?
2. Alternatively is it possible to retrieve a query\_id before the blast result is complete, and checking on it later?
3. (bit off-topic) I seem to only retrieve 60 results even at higher  `hitlist_size` although i get much more results when blasting the same entry on uniprot

I'm working on a script which processes a BLAST result and having to re-request the search every time is quite temporally cumbersome since i'm still learning the library. currently i'm using:  `blastResult = NCBIWWW.qblast("blastp", "swissprot", targetAminoAcids,format_type="Tabular", hitlist_size=200)`

# Answer
1 and 2 —\> Yes, but you need to have the query ID in order to do so: having the description would be a bonus, but not required. If you go to 8.1.1 Query Result of the Biopython manual, it will tell you how to get a query id, if you don’t already have it! Go to section 7.1 as a reminder of how to run BLAST via the internet. Go to section 7.2 to. Run BLAST locally. Weigh your options out with running via Internet or locally. Running BLAST locally is faster usually and gives you more user freedom to create your own database with your own parameters, etcetera., but certain tools will need to be downloaded via the command line/terminal. Follow this link for the BLAST command line user manual: https://www.ncbi.nlm.nih.gov/books/NBK279690/ I will paste the link to the Biopython user manual if you do not already have it: http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec123

3 —\> Try the Protein Data Bank (PDB) instead of UniProt or you can combine both UniProt and PDB (I recommend this route as it’s a combination of the 2 databases). You can use these two databases to process a BLAST result! I have found in my experience, combining data from more than 1 database has been beneficial in my case!
> 1 votes

---
Tags: python, blast, biopython
---


# Question
Title: If I modify a PDB file with a specific mutation, how to minimize energy?
I have a series of proteins that are all phylogenetically related. Only one of this proteins is currently X-ray crystallized.

Is it valid to modify the reference PDB file (I mean by hand) to match the mutated sequence?

And if so, is there any way to minimize the protein energy to find the most likely structure? (Trying to match as closely as possible the reference structure)

# Answer
> 3 votes
I would recommend you try the FoldX mutation engine. It is very easy to use and state of the art (it is the main competitor to expensive commercial software). It can do both of what you want - the mutation and energy optimisation. It may not work as well as molecular dynamics software, but for your simple task, it is a reasonable solution.

A simple command line can do it:

First, repair your structure (since, usually PDBs are not energetically optimised, so you could get better energy anyway without doing this):

```
FoldX --command=repairpdb --pdb=1AJY.pdb

```

Then, either use AlaScan or PositionScan to mutate. AlaScan only mutates to Alanine. PositionScan can mutate to any amino acid. More advanced mutations can be done with the BuildModel command.

```
FoldX --command=PositionScan --pdb=1AJY_Repair.pdb --output-file=1AJY_Repair_MA30d --out-pdb=false --positions=MA30d

```

Note: 'd' means try mutate the residue to 24 different amino acids; or, you can specify any of those 24, including the standard 20.

Which produces two files:

1. PS\_1AJY\_Repair\_MA30d\_scanning\_output.txt
2. energies\_30\_1AJY\_Repair

(It can also output the mutated pdb file by setting --out-pdb=true.)

Either of those files can be used to determine change in free energy upon mutation (ddg ~ delta delta g). The second file has additional information (the individual values used to calculate total change in free energy).

Note: The FoldX website has very good documentation for each FoldX command. You can find any additional information you require there.

# Answer
> 5 votes
Yes you can modify the reference PDB file and look for the changes and for this purpose you need visualizers. One among these is Chimera. You can easily carry out the energy minimization steps using Chimera, first mutate your PDB file manually and then provide it as an input for Chimera and the steps for doing energy minimization are fairly simple as given in Chimera docs here.

# Answer
> 2 votes
The previous answers suggest Chimera and FoldX, which require a license for commercial use.

Commercial users can use pdbfixer or a similar tool to perform mutations (using e.g. applyMutation in the API).

You can then use OpenMM to refine the mutated structure.

Granted, this requires some comfort writing Python scripts to use the APIs.

---
Tags: proteins, protein-structure, pdb
---


# Question
Title: Confirming Identity of clone in Phylogeny
This might be something very trivial to ask, but I am new to this and was confused as in how to identify or differentiate between very closely related isolates from clonal isolate from the whole genome/core-genome phylogeny?

For instace, I am drawing a core-genome phylogeny of Staphylococcus aureus genomes taken from NCBI as well as from my own sequenced ones. Now, my isolates have grouped with several isolates from NCBI (Bottom four greens with four red ones). I want to know, weather these isolates are actually closely related to NCBI ones or they are just those strains themselves?

Thank you all in advance!

# Answer
> 0 votes
Your isolates are both robustly part of the NCBI lineage of 4, i.e. now making a lineage of 8 taxa and genetically very close to it, but slightly distinct.

The robustness is given by the bootstrap values (100) and the genetic proximity is given by the genetic distances for each branch. What you want to do is look at the epidemiological information for your strains and those in Genbank, there should be some sort of recognised association, e.g. country, possibly locality, anything else?

The association has epidemiological significance, thus it's not random, if that is your question. There will definitely something connecting this lineage of 8.

---
Tags: phylogeny, genomics
---


title: biomarRT conversion of symbols

text:
# Question
Title: biomarRT conversion of symbols
I wanted to use the code to convert mouse to human symbol, but I have problem with the mirror, I've tried to establish another one, but continue to have the same problem.

```
musGenes <- c("Hmmr", "Tlx3", "Cpeb4")

# Basic function to convert mouse to human gene names
convertMouseGeneList <- function(x){
    require("biomaRt")
    human = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
    mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl")
    genesV2 = getLDS(attributes = c("mgi_symbol"), filters = "mgi_symbol", values = x , mart = mouse, attributesL = c("hgnc_symbol"), martL = human, uniqueRows=T)
    humanx <- unique(genesV2[, 2])

    # Print the first 6 genes found to the screen
    print(head(humanx))
    return(humanx)
}

```

# Answer
My background in coding is more extensive in Python. If you have Python skills, I can say for sure, use the Bioservices library:

```
from bioservices import BioMart

```

Install Bioservices via pip:

```
pip install bioservices

```

You can also use the web-based tool on Ensembl called BioMart. It does the same thing, just easier to use if you continue to have difficulties with Python or in your case, R. I looked up the R package for BioMart, it is called: biomaRt

Here's the link for the Ensembl BioMart web-based tool if you need to use it and are in a crunch for time on your project:

https://useast.ensembl.org/info/data/biomart/how\_to\_use\_biomart.html#introduction

I hope this helps! :)
> 1 votes

---
Tags: r, biomart
---


# Question
Title: Mitochondrial genes - TPM calculation bulk RNA-Seq
*This question was also asked on Biostars*

I was wondering if any of you have encountered a situation for bulk RNA-Seq where, possibly due to low sample quality or presence of dead cells, mitochondrial genes are expressed to a very large degree relative to other genes, thus skewing TPM values of all nuclear genes (by effectively scaling them down).

In such a situation, what could I potentially do to alleviate this beyond preparing fresh samples? Could I, for example, exclude 'MT-' genes from samples and then recalculate TPM based on this filtered set of genes?

My aim was to get some idea of within-samples expression levels, hence why going for TPMs as I know that counts normalised by DESeq2, edgeR and the like do not correct for gene length and are inappropriate for within-samples comparisons.

Many thanks in advance for any insight

# Answer
> 1 votes
Yes, this is a common situation. Not just for mitochondrial genes, but for any gene that may be an outlier (e.g. the B-cell receptor gene in a contaminant/unexpected B-cell population). The way to get around this is to not use TPM for differential expression analysis.

DESeq2 by default uses a quantile-based normalisation method for bulk differential expression calculations, which works well when the majority of genes are not differentially expressed. When this assumption is broken (e.g. for a targeted panel), a control gene list can be supplied, or custom size factors can be applied to each sample.

More information about using DESeq2 for differential expression can be found in the DESeq2 manual:

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

For comparing expression from different genes, I normalise variance-stabilised data using the longest annotated gene length from transcripts.

See this answer and this paper for more details.

---
Tags: rna-seq, tpm
---


title: getting NCBI TAXname by R

text:
# Question
Title: getting NCBI TAXname by R
I have a data frame on R, called qIns8, with 14 columns and 204 rows. the column called "subject tax id" has NCBI taxIDs for the corresponding row, I need a script to add another column to the data frame called Organism for which the script has to get the organism name (NCBI taxname) from the column subject tax id. I already tried this script using Taxize package, but it did not work:

```
getTAXname <- function(tax_id) {
  name <- tryCatch({
    tax_name <- sci_name(tax_id)
    return(tax_name)
  }, error = function(e) {
    return(NA)
  })
  return(name)
}

# Apply the function to each row of the dataframe
qIns8$Organism <- sapply(qIns8$subject_tax_id, getTAXname)

# Check the updated dataframe
head(qIns8)

```

I also tried this from Taxonomizr package and it did not work either:

```
getTaxonomy(qIns8$`subject tax id`, sqlFile = qIns8$`Organism`)

```

# Answer
> 1 votes
I developed the function and it worked perfectly,

```
Converst_taxID_to_taxName <- function(tax_id) {
  tax_info <- entrez_summary(db="taxonomy", id=tax_id)
  return(tax_info$scientificname)
}

# usage example
myDataFrame <- myDataFrame %>%
mutate(Organism = sapply(`subject tax id`, Converst_taxID_to_taxName))

# in case one tax ID has a problem or cannot be proceeded, it puts "NA" and automatically goes the next ID
Converst_taxID_to_taxName2 <- function(tax_id) {
  tryCatch({
    tax_info <- entrez_summary(db="taxonomy", id=tax_id)
    return(tax_info$scientificname)
  }, error = function(e) {
    return(NA)
  })
}

```

---
Tags: r, taxonomy, ncbi
---


# Question
Title: Sequence Alignment for sequences with the same length
I am doing research on a new method of optimizing sequence alignment process (Needleman - Wunsch algorithm) but the idea would only work with sequences that have the same length. I am wondering if that will have any real-life utility. Specifically, how often do we have those situations where we are aligning same-length sequences (or at least sub-sequences) and what are those situations?

# Answer
> 1 votes
I see ... the application is if two genes had the same number of indels\* at different homologous positions within the gene. This happens but equally one gene might have an indel and another gene might lack it, in which case these genes would be different lengths.

So just to reiterate your method would find genes with an identical number of indels at different homologous positions in a gene. Its application is protein modelling. This is because indels cause a large change in the entropy of a protein, so even if these two genes had identical sequence (barring the indels obviously) they would be structurally different.

Beyond this however the application is limited because apart from the known changes to protein structure, little is understood about quantifying the function of indels. Certainly in phylogenetics we specifically avoid them, because we struggle to fit them into a known model.

---

\*, indel = insertion deletion event. In a protein they occur in triplets.

# Answer
> 0 votes
My 2 cents: Sequence alignment is a well studied problem in bioinformatics with some great existing methods that do handle uneven length sequences (i.e. indels). With tons of genomic data the trend these days is towards faster methods with more heuristics like pseudo-alignment or k-mer based approaches rather than theoretically optimal but slow and unscalable approaches like Needleman Wunsch. There might be some niche applications on DNA motifs or short peptide analyses but I think these are quite limited.

---
Tags: sequence-alignment, rna, sequence, dna
---