# Question
Title: Convert ensemble genes to gene names in my sigle cell signature matrix annotated with seurat in R?
I have single cell data which I want to annotate with the seurat package, so I ran this code in R:

```
lung.data <- Read10X(data.dir ="C:/Users/Desktop/SINGLE-CELL/DATA",gene.column = 1)

lung <- CreateSeuratObject(counts = lung.data, project = "lung", min.cells = 3, min.features = 200)

```

And the output of the features of this seurat object is:

```
    >rownames(lung)
    [1] "ENSG00000000003" "ENSG00000000419"
 [3] "ENSG00000000457" "ENSG00000000460"
 [5] "ENSG00000000938" "ENSG00000000971"
 [7] "ENSG00000001036" "ENSG00000001084"
 [9] "ENSG00000001167" "ENSG00000001460"
[11] "ENSG00000001461" "ENSG00000001497"
[13] "ENSG00000001561" "ENSG00000001617"
[15] "ENSG00000001626" "ENSG00000001629"
[17] "ENSG00000001630" "ENSG00000001631"
[19] "ENSG00000002016" "ENSG00000002079"

```

How can I change the row names of this Seurat Object from Ensembl to gene name, knowing that not all the ensemble genes have conversion to gene name, in order to obtain the seurat object with the gene names as row names and cell names as column names?

# Answer
I think BioMart will do the gene ID -\> gene name conversion you're after: https://www.ensembl.org/info/data/biomart/index.html

To use BioMart to convert the gene IDs to gene names, you can use the BioMaRt bioconductor R package. Installation and running instructions can be found here: ensembl.org/info/data/biomart/biomart\_r\_package.html

For this query, you will need to use the human genes dataset:

```
ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")

```

Then use your list of gene IDs as the 'Filter', with list of IDs as the values:

```
filters = 'ensembl_gene_id'

```

Most human protein-coding genes have an associated HGNC symbol, and these are used in Ensembl. Therefore, you can select the 'gene\_name' attribute to retrieve the human gene names, and this will usually correspond to the HGNC gene names, as stored in the Ensembl database.

More information on gene names can be found in the Ensembl documentation: https://ensembl.org/info/genome/genebuild/gene\_names.html
> 1 votes

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


# Question
Title: BioPython internal_coords module returns different dihedral angles for the (seemingly) same protein structure
I am using the `internal_coords` module from BioPython in order to compare dihedral angles of two different conformations of the same protein. The conformations are processed by me in order to make sure that their residues are aligned and they always contain exactly the same atoms, so that in the end I expect to have the same number of dihedral angles in the two structures.

However, despite the match in the residues and atoms, I obtain a different number of dihedral angles for the two conformations.

You can find a reproducible code sample here.

```
from Bio.PDB import *

def extract_torsion_angles(structure):
    ic_chain = internal_coords.IC_Chain(structure)
    ic_chain.atom_to_internal_coordinates()

    d = ic_chain.dihedra

    """
    Update the keys of a dihedra dictionary to make searching easier.

    The keys of the dihedra dictionary (as obtained from
    extract_torsion_angles()) are 4-tuples where each element is
    an AtomKey containing 5 fields (respos, icode, resname, atm,
    altloc, occ). We use the residue position, residue name and
    atom name in order to create a new representation for each
    AtomKey. The final dihedra dictionary will have 4-tuple keys
    where each element of the tuple is one such string.
    """

    new_d = {}
    for key, value in sorted(d.items(), key=lambda x: x[0]):
        new_key = []
        for atom_key in key:
            # Keep: residue position + resiude name + atom name.
            atom_key_str = (atom_key.akl[0] + '_' +
                atom_key.akl[2] + '_' + atom_key.akl[3])

            new_key.append(atom_key_str)

        new_d[tuple(new_key)] = value

    return new_d

f1 = 'protein1.pdb'
f2 = 'protein2.pdb'

parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('', f1)
s2 = parser.get_structure('', f2)

c1 = list(s1.get_chains())
c2 = list(s2.get_chains())

print('Chains: ', len(c1), ' ', len(c2))

r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())

print('Residues: ', len(r1), ' ', len(r2))

for i in range(len(r1)):
    resname1 = r1[i].get_resname()
    resname2 = r2[i].get_resname()

    if resname1 != resname2:
        print('Residues not matching')
        break

    atoms1 = list(r1[i].get_atoms())
    atoms2 = list(r2[i].get_atoms())

    atom_names1 = sorted([atom.get_id() for atom in atoms1])
    atom_names2 = sorted([atom.get_id() for atom in atoms2])

    if atom_names1 != atom_names2:
        print('Different atoms!')

d1 = extract_torsion_angles(s1)
d2 = extract_torsion_angles(s2)

print('Torsion angles: ', len(d1), ' ', len(d2))

diff_keys = set(d2.keys()) - set(d1.keys())
print(diff_keys)

```

These are the two PDB files protein1.pdb and protein2.pdb used.

As you can see, the number of chains, the number of residues, the residue names and the atom names in each residue are the same, yet there are 4 extra dihedral angles in the second protein.

My question is, am I missing something from my comparison? Is it unreasonable to expect the same dihedral angles for these two structures? (And if yes, why?). Is this the expected behaviour of the `atom_to_internal_coordinates()` function?

And extra: do you have any recommendations for other tools that can easily compute all dihedral angles (including the ones in side chains) for a protein?

# Answer
> 7 votes
The problem got solved after posting on the official BioPython repo.

Because of the processing (where some residues were removed in order to align the proteins), some consecutive residues were too far away from each other and thus the bond length was too large. The solution was to increase the maximum peptide bond allowed, as follows:

```
internal_coords.IC_Chain.MaxPeptideBond = 5.0
ic_chain = internal_coords.IC_Chain(structure, verbose=True)

```

# Answer
> 0 votes
OK , visualizing both proteins in **PyMOL** , using `set cartoon_gap_cutoff = 0`  ,

I get :

Apparently the two conformations , have more than one difference in various loops.

I modify the code above , thanks to your hint :

```
internal_coords.IC_Chain.MaxPeptideBond = 5.0
ic_chain = internal_coords.IC_Chain(structure, verbose=True)

```

My code :

```
import Bio

print('########## Biopython  VERSION ##########################################')
print('\n                BIOPYTHON_VERSION : ', Bio.__version__)
print('\n########################################################################\n\n')

from Bio.PDB import *

def extract_torsion_angles(structure):

    print('\nSTRUCTURE ---> ' , structure)

    ic_chain = internal_coords.IC_Chain(structure, verbose=True)

    ic_chain.atom_to_internal_coordinates()

    d = ic_chain.dihedra

    """
    Update the keys of a dihedra dictionary to make searching easier.

    The keys of the dihedra dictionary (as obtained from
    extract_torsion_angles()) are 4-tuples where each element is
    an AtomKey containing 5 fields (respos, icode, resname, atm,
    altloc, occ). We use the residue position, residue name and
    atom name in order to create a new representation for each
    AtomKey. The final dihedra dictionary will have 4-tuple keys
    where each element of the tuple is one such string.
    """

    new_d = {}
    for key, value in sorted(d.items(), key=lambda x: x[0]):
        new_key = []
        for atom_key in key:
            # Keep: residue position + resiude name + atom name.
            atom_key_str = (atom_key.akl[0] + '_' +
                atom_key.akl[2] + '_' + atom_key.akl[3])

            new_key.append(atom_key_str)

        # new_d[tuple(new_key)] = value

        new_d[tuple(new_key)] = tuple((value , structure.get_full_id()))

    return new_d

f1 = '18760-BIOINFORMATICS-protein1.pdb'
f2 = '18760-BIOINFORMATICS-protein2.pdb'

parser = PDBParser(PERMISSIVE=1, QUIET=0)
s1 = parser.get_structure('protein1', f1)
s2 = parser.get_structure('protein2', f2)

c1 = list(s1.get_chains())
c2 = list(s2.get_chains())

print('Chains: ', len(c1), ' ', len(c2))

r1 = list(c1[0].get_residues())
r2 = list(c2[0].get_residues())

print('Residues: ', len(r1), ' ', len(r2))

for i in range(len(r1)):
    resname1 = r1[i].get_resname()
    resname2 = r2[i].get_resname()

    if resname1 != resname2:
        print('Residues not matching')
        break

    atoms1 = list(r1[i].get_atoms())
    atoms2 = list(r2[i].get_atoms())

    atom_names1 = sorted([atom.get_id() for atom in atoms1])
    atom_names2 = sorted([atom.get_id() for atom in atoms2])

    if atom_names1 != atom_names2:
        print('Different atoms!')

def get_no_tors_diff(structure1 , structure2) :

    pept_bond = 1.4

    while pept_bond < 10.0 :

        pept_bond = round(pept_bond , 2)

        print('\n\n Peptide Bond Lenght = ' , pept_bond)

        internal_coords.IC_Chain.MaxPeptideBond = pept_bond

        d1 = extract_torsion_angles(s1)
        d2 = extract_torsion_angles(s2)

        print('Torsion angles: ', len(d1), ' ', len(d2))

        diff_keys = sorted((set(d1.keys()) ^ set(d2.keys())) , key = lambda tup : tup[0].split('_')[0])

        print('diff_keys lenght : ', len(diff_keys))

        print(diff_keys , type(diff_keys))

        for i in diff_keys :

            print(i , i[0].split('_')[0] , type(i))

        diff_keys_values = {}

        for i in diff_keys :

            if i in d1.keys() :

                diff_keys_values[i] = d1[i]

            if i in d2.keys()  :

                diff_keys_values[i] = d2[i]

            else:

                pass

        print('diff_keys_values lenght : ', len(diff_keys_values))

        for i in diff_keys_values :

            print(i , diff_keys_values[i][1])

        print('\n\n')

        if len(diff_keys) == 0 :

            print('\n\n Peptide Bond Lenght without differences = ' , pept_bond)

            return pept_bond

        pept_bond += 0.1

    print('\n\n Cannot get rid of differences up to 10.0A in Peptide Bond ')

    return None

print('get_no_tors_diff(structure1 , structure2)' , get_no_tors_diff(s1 , s2))

```

First and last bit of output :

```
Chains:  1   1
Residues:  313   313

 Peptide Bond Lenght =  1.4

STRUCTURE --->  <Structure id=protein1>
chain break at GLY  6  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  37  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  38  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  75  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  84  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  92  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PRO  163  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  164  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at HIS  198  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  224  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  246  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  251  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  279  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  285  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (1.4 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at PHE  1  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at GLY  6  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  37  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  38  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  75  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  84  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  92  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PRO  163  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at VAL  164  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at HIS  198  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  224  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LYS  246  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at LEU  251  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  278  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ALA  285  due to MaxPeptideBond (1.4 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (1.4 angstroms) exceeded
Torsion angles:  2504   2501
diff_keys lenght :  9
[('0_N_C', '1_F_N', '1_F_CA', '1_F_C'), ('0_N_CA', '0_N_C', '1_F_N', '1_F_CA'), ('0_N_N', '0_N_CA', '0_N_C', '1_F_N'), ('277_D_C', '278_A_N', '278_A_CA', '278_A_C'), ('277_D_N', '277_D_CA', '277_D_C', '278_A_N'), ('277_D_CA', '277_D_C', '278_A_N', '278_A_CA'), ('278_A_N', '278_A_CA', '278_A_C', '279_K_N'), ('278_A_C', '279_K_N', '279_K_CA', '279_K_C'), ('278_A_CA', '278_A_C', '279_K_N', '279_K_CA')] <class 'list'>
('0_N_C', '1_F_N', '1_F_CA', '1_F_C') 0 <class 'tuple'>
('0_N_CA', '0_N_C', '1_F_N', '1_F_CA') 0 <class 'tuple'>
('0_N_N', '0_N_CA', '0_N_C', '1_F_N') 0 <class 'tuple'>
('277_D_C', '278_A_N', '278_A_CA', '278_A_C') 277 <class 'tuple'>
('277_D_N', '277_D_CA', '277_D_C', '278_A_N') 277 <class 'tuple'>
('277_D_CA', '277_D_C', '278_A_N', '278_A_CA') 277 <class 'tuple'>
('278_A_N', '278_A_CA', '278_A_C', '279_K_N') 278 <class 'tuple'>
('278_A_C', '279_K_N', '279_K_CA', '279_K_C') 278 <class 'tuple'>
('278_A_CA', '278_A_C', '279_K_N', '279_K_CA') 278 <class 'tuple'>
diff_keys_values lenght :  9
('0_N_C', '1_F_N', '1_F_CA', '1_F_C') ('protein1',)
('0_N_CA', '0_N_C', '1_F_N', '1_F_CA') ('protein1',)
('0_N_N', '0_N_CA', '0_N_C', '1_F_N') ('protein1',)
('277_D_C', '278_A_N', '278_A_CA', '278_A_C') ('protein1',)
('277_D_N', '277_D_CA', '277_D_C', '278_A_N') ('protein1',)
('277_D_CA', '277_D_C', '278_A_N', '278_A_CA') ('protein1',)
('278_A_N', '278_A_CA', '278_A_C', '279_K_N') ('protein2',)
('278_A_C', '279_K_N', '279_K_CA', '279_K_C') ('protein2',)
('278_A_CA', '278_A_C', '279_K_N', '279_K_CA') ('protein2',)

.........................
.................
..........
......
 Peptide Bond Lenght =  4.6

STRUCTURE --->  <Structure id=protein1>
chain break at PHE  160  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.6 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at LEU  84  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  160  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.6 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.6 angstroms) exceeded
Torsion angles:  2546   2543
diff_keys lenght :  3
[('83_K_CA', '83_K_C', '84_L_N', '84_L_CA'), ('83_K_N', '83_K_CA', '83_K_C', '84_L_N'), ('83_K_C', '84_L_N', '84_L_CA', '84_L_C')] <class 'list'>
('83_K_CA', '83_K_C', '84_L_N', '84_L_CA') 83 <class 'tuple'>
('83_K_N', '83_K_CA', '83_K_C', '84_L_N') 83 <class 'tuple'>
('83_K_C', '84_L_N', '84_L_CA', '84_L_C') 83 <class 'tuple'>
diff_keys_values lenght :  3
('83_K_CA', '83_K_C', '84_L_N', '84_L_CA') ('protein1',)
('83_K_N', '83_K_CA', '83_K_C', '84_L_N') ('protein1',)
('83_K_C', '84_L_N', '84_L_CA', '84_L_C') ('protein1',)

 Peptide Bond Lenght =  4.7

STRUCTURE --->  <Structure id=protein1>
chain break at PHE  160  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.7 angstroms) exceeded

STRUCTURE --->  <Structure id=protein2>
chain break at PHE  160  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at PHE  245  due to MaxPeptideBond (4.7 angstroms) exceeded
chain break at ARG  296  due to MaxPeptideBond (4.7 angstroms) exceeded
Torsion angles:  2546   2546
diff_keys lenght :  0
[] <class 'list'>
diff_keys_values lenght :  0

 Peptide Bond Lenght without differences =  4.7
get_no_tors_diff(structure1 , structure2) 4.7

```

Using `internal_coords.IC_Chain.MaxPeptideBond = 4.7` both models have same Chain Breaks at resi 160, 245 e 296 and are described by 2546 dihedral angles.

My question now is how, if feasible , to force PyMOL to paint petide bonds of such lenght both in cartoons and lines/sticks mode.

Second question if Biopython has a tool to detect chain breaks in Structures other than just printing it out with `ic_chain = internal_coords.IC_Chain(structure, verbose=True)` , asking this because setting `internal_coords.IC_Chain.MaxPeptideBond = 50` gets rid of all the chainbreaks but don't have a way to monitor their disappearing. Pics below show chain breaks at 4.7.

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


# Question
Title: Statistical approach to link DNA methylation with toxic element exposure and health outcome
I would be thankful to you if you can help with the statistical approach for 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.

---

**Notes**

1. The sample size is small.
2. Maybe structural equation modelling could be used to find CpG sites associated with arsenic exposure and disease condition as I have very small size samples
3. R is a preferred language.

# Answer
> 1 votes
Having given some thought to this, I'd still use machine learning. I would attempt to augment the negative controls however. Is there a suitable published data set? R program is caret but it high barrier to entry without appropriate support.

250 samples is just about enough. If you can augment that to 350 with negatives ... should be enough. What you do is dump all your data ("features") into a single record, including the methylation data and use toxic element exposure as the training target. The output you want is called "feature selection", what you are wanting is methylation sites to carry more "weight" than any other feature. Then you start other analytics to see what it's associating with. Its not trivial getting it working and a sceptic would complain about the lack of significance in the results, but its a good starting point.

GWAS might work, just considering methylation vs toxic element exposure. SEM, I've never used it so I don't know.

---
Tags: r, statistics, linear-regression, epigenetics, epidemiology
---


# Question
Title: Get a certain gene sequence from bam/vcf and reference
I need to get a fasta sequence of a certain gene for a certain worm strain that is different from reference. I have a reference genome, BAM for the strain of interest, and coordinates of the gene. I know that vcftools can convert bam to fasta, but I do not see if I can do it only for a certain coordinates (like, for example, chr1:10100500-10200500). I converted the whole bam to fasta and then I tried using bedtools getfasta, but my system ran out of space. Is there a more rational way to do it?

# Answer
> 4 votes
The first approach that comes to mind is to cut the bam file and then convert the already cut bam to fasta:

```
samtools view -hb chr1:10100500-10200500 > small.bam
samtools fasta small.bam > small.fa

```

(note that I am not certain about the `samtools fasta` command, I have never used it and I don't know if it will produce many sequences, one for each haplotype supported by reads; please test this.)

---
Tags: fasta, samtools, bedtools
---


# Question
Title: Are there applications of semi-supervised learning to genomics?
Semi-supervised learning is when a portion of the training data set is labeled, and the rest is unlabeled; the unlabeled data is used to improved learning performance / accuracy.

Are there (i.e., have there been) applications of semi-supervised learning to genomics? If so, what are they? If not, is there a reason why not?

# Answer
**Personal view** Semi-supervised learning is basically unsupervised learning with weights, so the "label" or "category" can be quantified rather than approximately clustered as occurs in unsupervised learning.

Unsupervised learning is PCA and tSNE and other multi-variate approaches. The advantages of semi-supervised learning is unlike unsupervised learning the feature in the data that shows greatest proportionate variance can be precisely identified and all features ranked; thus there aren't loads of components where the result could be anywhere.

Sure I've interpreted semi-supervised learning and found it was not useful for any comparative purposes. Basically the same variance Eigen-vector weight observed in one gene in one unique environment, is the strongest Eignen-vector weight for the same gene in a completely different environment. However, these are two independent semi-supervised learning analyses. Thus, it's not measuring comparative difference but relative variance between all genes in that given experiment.

Thus I personally didn't find it useful for my metric, HOWEVER such a result is good for some people where the overall structure of the variance identified by Eigen vectors is very similar between contrasting environments. Basically it's a comparison to all other genes in the analysis. To get it working you'd need to combine all genes in all environments into a single data set.

That would work, but ensure that the label for each feature (e.g. gene expression) comprised the environment and the gene. What the statistical theoretical validity of such a result would be - I haven't a clue - but it should give some interesting results. To be honest, if I'm not sure about the theoretical validity, I will not be alone.
> 1 votes

---
Tags: genomics, machine-learning
---


# Question
Title: Does AlphaFold actually calculates the fold?
After reading through papers and the code, I'm still not clear if AlphaFold and the likes do actually calculate the 3d structure. My current feeling they are more like information retrieval engines. They find statistical 'similarity' features between the primary structures of the query protein and a database protein (or a family of those) with a known 3d structure, and that known fold is the 'calculated' 3d structure. This would mean that if there is no 'similar' protein with a known 3d fold then the program cannot 'predict' anything. I realize this could be a stupid layman question but here it goes.

# Answer
Alphafold calculates the "folds" throughout a protein, thus its protein structural prediction. It uses convolution-based calculation, which is a modification of co-variance. The reason it doesn't use co-variance directly is complicated, but can result in artefacts.

Thus it looks for "simultaneous" amino acid mutations in protein alignments and believes those represent close spatial proximity - based on previous supervised learning of proteins, comprising an alignment and a known structure.

It has won two CASP competitions (world competition for protein folding prediction). The training is likely deeper both in data sets and ANN/CNN/RNN "layers" (deep learning jargon) than previous protein structural prediction models and hence its success. "Deepmind" the company behind it will be Tensorflow experts, rather than structural biologists.

I don't use it but have seen its structural accuracy first-hand. The accuracy statistic RMSD was impressive.

It is important to note it is built on the back of a long term effort by this community as to calculating protein structural prediction. So it's not "radical technology", they've generally increased the accuracy of established methods. Again they've likely done this with more extensive training.
> 0 votes

---
Tags: protein-structure, alphafold
---


# Question
Title: Where to access the WGCNA tutorial documents: Horvath lab site down
I am currently using the WGCNA package for some analysis and it seems the Horvath lab site is down. Does anyone know of anywhere else I can access the tutorial documents?

# Answer
> 10 votes
Apologies for the web site trouble, we are working with UCLA to restore the web pages. In the meantime, I uploaded the tutorials to Dropbox at https://www.dropbox.com/scl/fo/4vqfiysan6rlurfo2pbnk/h?rlkey=thqg8wlpdn4spu3ihjuc1kmlu&dl=0. I suggest downloading the entire folder, unzipping and opening the file index.html in a browser. You should see a stripped down but usable version of the web page and all tutorial sections should be accessible.

---
Tags: wgcna
---


title: PANTHER FUNCTIONAL CLASSIFICATION in R

text:
# Question
Title: PANTHER FUNCTIONAL CLASSIFICATION in R
Does anyone know of a way to programmatically get panther protein class terms and IDs for a gene list, in a way that the classification structure is preserved (PARENT CLASS + CHILD CLASS)?

i.e. For a given gene, I would want to get the list of classes, in away that I can then use the top classification for downstream analyses.

GENEX, PC00244 zinc finger transcription factor \> PC00218 DNA-binding transcription factor \> PC00264 gene-specific transcriptional regulator

# Answer
> 1 votes
One approach is to get the gene's sequence (DNA/AA) and use the interproscan API (e.g., https://www.ebi.ac.uk/Tools/common/tools/help/index.html?tool=iprscan5) with the httr2 library in R to POST gene sequence(s) and GET the resulting Panther annotation. Make sure to include `Panther` in the `appl` parameter for the post request.

Here's a function that does this for a single sequence.

```
submit_fasta <- function(
    fasta,
    path_outfolder,
    email,
    stype = "n",
    appl = c("Panther"),
    title = 'stackoverflow'
    ) {
  stopifnot("stype must be 'n' for nucleotide seqs or 'p' for protein seqs" = stype %in% c('n', 'p'))
  url_post <- "https://www.ebi.ac.uk/Tools/services/rest/iprscan5/run"
  chr_seq <- as.character(fasta)
  req <- httr2::request(url_post)
  req <- req |> httr2::req_body_form(
    "email" = email,
    "title" = title,
    "goterms" = "false",
    "pathways" = "false",
    "stype" = stype,
    "appl" = paste0(appl, collapse = ","),
    "sequence" = chr_seq
  )
  req <- req |> httr2::req_headers("Accept" = "text/plain", "Content-type" = "application/x-www-form-urlencoded")
  resp_post <- req |> httr2::req_perform()
  chr_retrieval_code <- resp_post$body |> rawToChar()
  url_get <- "https://www.ebi.ac.uk/Tools/services/rest/iprscan5/result"
  req <- httr2::request(url_get)
  req <- req |> httr2::req_url_path_append(chr_retrieval_code)
  req <- req |> httr2::req_url_path_append("tsv")
  req <- req |> httr2::req_retry(
    max_tries = 15L, backoff = ~ 30L,
    is_transient = function(resp) {httr2::resp_status(resp) == 400L}
  )
  req <- req |> httr2::req_headers(Accept = NULL)
  req |> httr2::req_dry_run(quiet = FALSE)
  resp_get <- req |> httr2::req_perform()
  cat('status code:', resp_get$status_code, '\n')
  result <- resp_get$body |> rawToChar()
  result <- read.table(text = result, sep = '\t')
  return(result)
}

```

Importantly, the fasta input parameter is expected to be XString object from the Biostrings package. However, the code could easily be adjusted. The sequence eventually needs to be coerced to a string for the POST request.

```
#fasta <- Biostrings::readDNAStringSet('seq.fa')
#resp <- submit_fasta(fasta = fasta_single, path_outfolder = 'x', email = 'test@gmail.com', title = 'stackoverflow')

```

**edit**: set the `is_transient` param to handle the actual transient error case where the data is not yet available (response code = 400). Also, it's important to note there's a "status" endpoint https://www.ebi.ac.uk/Tools/common/tools/help/index.html?tool=iprscan5#/Status) which is probably much more appropriate than pinging the data endpoint until either success or reaching max tries. However, I hope this answer suffices!

---
Tags: r, public-databases, functional-annotation
---


# Question
Title: How long should it take to cluster Uniref?
I have the Uniref90 dataset containing around 76M sequences. I have made modifications to all sequences based on a reduced amino acid alphabet.

I am trying to cluster it again at 90% identity. I am using `mmseq cluster` on a 64 core, 500GB RAM machine. It has been more than 12h now.

```
Starting prefiltering scores calculation (step 1 of 1)
Query db start 1 to 74390789
Target db start 1 to 74390789
[=

```

How long does it normally take to cluster such a database? I know I can use `linclust` instead to cluster in linear time, but it might not be what I need in this case.

I am using `mmseq cluster` as follows:

```
mmseqs createdb "uniref90-$1.fasta" "uniref90-DB-$1"
mmseqs cluster "uniref90-DB-$1"  "uniref90-DB-$1_clu" tmp/ --cov-mode 2 -c 0.8 --min-seq-id 0.9

```

# Answer
> 0 votes
When I want to estimate the time taken to do something big, I subset my input data to something that takes a few minutes or less, and calculate forward from there. If there's alignment involved, it's usually a good idea to work it out for multiple dataset sizes, as there could be a linear relationship between input size and time taken, square relationship, or worse.

Have you looked at how long it takes for, say, 1000 sequences, 2000 sequences, and 5000 sequences?

# Answer
> 0 votes
The **mmseqs2** code should be specifically compiled and optimised for your system. **mmseqs2** optimisation is automated and might be available via Spack\* (I'm not sure its actually done for mmseqs2), it is definitely available for **EasyBuild**. I wouldn't advise attempting a manual optimisation unless you've prior experience.

How much time this will save is more difficult to assess but 50% faster is a conservative assessment. This approach is good practice for a job like this size and its good for the associated report, its standard HPC practice. If you're doing cloud computation (i.e. $$$), HPC optimisation is essential. AWS make serious money from suboptimal HPC practice (in fairness to them they also do a lot of development work to minimise costs).

Finally, if you're an expert in optimisation you can save **serious** computational time, but this is a major investment of your personal development time and technique. The rationale of how and why is complicated.

\*, its very likely available on AWS.

---

You might simply consider its been running for a few days already so HPC technique isn't needed for this specific job. However, you'll know for the future.

---
Tags: proteins, clustering, uniprot
---


title: How to add charges to a PDB file?

text:
# Question
Title: How to add charges to a PDB file?
I did some calculations using a PDB file as input. As a result, I got a file with the atomic partial charges (only the values in one column, no other information) and the output PDB file structure without the charge information.

My question: Is there a tool to add the charges to the output PDB file using the separate file?

My PDB file:

```
    TITLE     sistema05
    REMARK   4      COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
    REMARK 888
    REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
    HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C
    HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C
    HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C
    HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C
    HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C
    HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C
    HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C
    HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C
    HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C
    HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C
    HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C
    HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C
    HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C
    HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C
    HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C
    HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C
    HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C
    HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C
    HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C
    HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C
...
...
...
    CONECT    1  212    2   20
    CONECT    1   20
    CONECT    2    1    3   21
    CONECT    2   21
    CONECT    3    2  213    4
    CONECT    3    4
    CONECT    4    3    5   23
    CONECT    4    3
    CONECT    5    4  214    6
    CONECT    5    6
    CONECT    6    5    7   25
    CONECT    6    5
    CONECT    7    6  215    8
    CONECT    7    8
    CONECT    8    7    9   27
    CONECT    8    7
    CONECT    9    8  216   10
    CONECT   10    9   11   29
    CONECT   11   10  217   12
    CONECT   12   11   13   31
    CONECT   13   12  218   14
    CONECT   14   13   15   33
    CONECT   14   15
    CONECT   15   14  219   16
    CONECT   15   14
    CONECT   16   15   17   35
    CONECT   16   17
    CONECT   17   16  220   18
    CONECT   17   16
    CONECT   18   17   19   37
    CONECT   18   19
    CONECT   19   18  221   20
    CONECT   19   18
    CONECT   20    1   19   39
    CONECT   20    1
    CONECT   21    2   22   40
    END

```

My CHARGE file:

```
   -0.02474906
    0.02291960
   -0.02789302
    0.02227110
   -0.02831805
    0.02230034
   -0.02876149
    0.02225527
   -0.02774076
    0.02216971
   -0.02074528
    0.02275309
   -0.02662750
    0.02177531
   -0.03007329
    0.02283631
   -0.02779055
    0.02312780
   -0.02312680
    0.02315866

```

# Answer
> 3 votes
OK given , `input_connect.pdb` , using **Biopython** library

```
TITLE     sistema05
REMARK   4      COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
REMARK 888
REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C
HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C
HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C
HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C
HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C
HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C
HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C
HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C
HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C
HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C
HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C
HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C
HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C
HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C
HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C
HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C
HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C
HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C
HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C
HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C
CONECT    1  212    2   20
CONECT    1   20
CONECT    2    1    3   21
CONECT    2   21
CONECT    3    2  213    4
CONECT    3    4
CONECT    4    3    5   23
CONECT    4    3
CONECT    5    4  214    6
CONECT    5    6
CONECT    6    5    7   25
CONECT    6    5
CONECT    7    6  215    8
CONECT    7    8
CONECT    8    7    9   27
CONECT    8    7
CONECT    9    8  216   10
CONECT   10    9   11   29
CONECT   11   10  217   12
CONECT   12   11   13   31
CONECT   13   12  218   14
CONECT   14   13   15   33
CONECT   14   15
CONECT   15   14  219   16
CONECT   15   14
CONECT   16   15   17   35
CONECT   16   17
CONECT   17   16  220   18
CONECT   17   16
CONECT   18   17   19   37
CONECT   18   19
CONECT   19   18  221   20
CONECT   19   18
CONECT   20    1   19   39
CONECT   21    2   22   40

```

and charges list `input_charge.txt` :

```
-0.02474906
 0.02291960
-0.02789302
 0.02227110
-0.02831805
 0.02230034
-0.02876149
 0.02225527
-0.02774076
 0.02216971
-0.02074528
 0.02275309
-0.02662750
 0.02177531
-0.03007329
 0.02283631
-0.02779055
 0.02312780
-0.02312680
 0.02315866

```

with code :

```
import Bio
print('########## Biopython  VERSION ##########################################')
print('\n\n BIOPYTHON_VERSION : ', Bio.__version__)
print('###################################################################')

from Bio.PDB import (
                    PDBParser,
                    Atom,

                    )

                    # PDBIO)

import PDBIOcopy

# from Bio.PDB import  PDBIO

pdb = 'input'

charges = 'input_charge.txt'

charges_list = []

with open(charges , 'r') as handle :

    for i in handle.readlines() :

        charges_list.append(i.strip())

print(charges_list)

parser = PDBParser(PERMISSIVE=True , QUIET = False)

structure = parser.get_structure(pdb, pdb+'.pdb')

atoms = structure.get_atoms()

idx = 0

for atom in atoms :

    print(atom)

    atom.set_charge(charges_list[idx])

    idx += 1

atoms = structure.get_atoms()

for atom in atoms :

    print(atom , atom.get_charge())

struct_to_save = PDBIOcopy.PDBIO()

# struct_to_save = PDBIO()

struct_to_save.set_structure(structure)

struct_to_save.save('result_input.pdb')

```

and modifying `Bio.PDB.PDBIO` as :

line 22 :

```
_ATOM_FORMAT_STRING = (
    "%s%5i %-4s%c%3s %c%4i%c   %8.3f%8.3f%8.3f%s%6.2f      %4s%2s%12s\n"
)

```

line 164 :

```
# charge="  ", #deleted
pqr_charge     # insertes

```

line 213:

```
charge = atom.pqr_charge  #inserted

```

line 389:

```
atom.pqr_charge,  #inserted

```

output `result_input.pdb`:

```
HETATM    1  C1  BAS A   1      10.151  -3.071  -1.881  1.00  0.00           C -0.02474906
HETATM    2  C2  BAS A   1       9.389  -3.683  -0.874  1.00  0.00           C  0.02291960
HETATM    3  C3  BAS A   1       9.976  -3.936   0.382  1.00  0.00           C -0.02789302
HETATM    4  C4  BAS A   1       9.194  -3.742   1.530  1.00  0.00           C  0.02227110
HETATM    5  C5  BAS A   1       9.781  -3.201   2.695  1.00  0.00           C -0.02831805
HETATM    6  C6  BAS A   1       9.029  -2.357   3.510  1.00  0.00           C  0.02230034
HETATM    7  C7  BAS A   1       9.651  -1.313   4.236  1.00  0.00           C -0.02876149
HETATM    8  C8  BAS A   1       8.948  -0.139   4.488  1.00  0.00           C  0.02225527
HETATM    9  C9  BAS A   1       9.635   1.100   4.502  1.00  0.00           C -0.02774076
HETATM   10  C10 BAS A   1       9.001   2.259   4.094  1.00  0.00           C  0.02216971
HETATM   11  C11 BAS A   1       9.724   3.177   3.329  1.00  0.00           C -0.02074528
HETATM   12  C12 BAS A   1       9.148   3.815   2.228  1.00  0.00           C  0.02275309
HETATM   13  C13 BAS A   1       9.909   3.945   1.081  1.00  0.00           C -0.02662750
HETATM   14  C14 BAS A   1       9.332   3.754  -0.198  1.00  0.00           C  0.02177531
HETATM   15  C15 BAS A   1      10.098   3.266  -1.245  1.00  0.00           C -0.03007329
HETATM   16  C16 BAS A   1       9.501   2.445  -2.232  1.00  0.00           C  0.02283631
HETATM   17  C17 BAS A   1      10.234   1.431  -2.847  1.00  0.00           C -0.02779055
HETATM   18  C18 BAS A   1       9.575   0.238  -3.215  1.00  0.00           C  0.02312780
HETATM   19  C19 BAS A   1      10.249  -0.986  -3.121  1.00  0.00           C -0.02312680
HETATM   20  C20 BAS A   1       9.529  -2.124  -2.713  1.00  0.00           C  0.02315866
TER      21      BAS A   1
END

```

That unfortunately is not a pdb file format anymore:

PDB Format uses for RECORD ATOM CHARGE columns 79-80 , see Protein Data Bank Contents Guide: Atomic Coordinate Entry Format Description

At the end of all of this I realized that using python and just adding the charge values at the end of the pdb file read line by line would be a lot easier & faster but anyway I had the chance to explore `PDBIO`

If you want to add the charges in the `bfactor` or `occupancy` columns just use :`atom.set_bfactor(charges_list[idx])` or `atom.set_occupancy(charges_list[idx])` be sure to modify `PDBIO` to make space for your values

There are other examples on SO about same problem too :

How to modify columns contents in a protein PDB file by keys and values stored in a text file using only bash tools

Just Python code :

```
pdb = 'input_connect.pdb'

charges = 'input_charge.txt'

charges_list = []

with open(charges , 'r') as handle :

    for i in handle.readlines() :

        charges_list.append(i.strip())

print(charges_list)

with open('input_connect.pdb' , 'r') as handle :

    pdb_file = handle.readlines()

print(pdb_file)

idx = 0

with open('result_python.pdb' , 'w') as handle :

    for i in pdb_file :

        i = i.strip()

        if i[0:6] == 'HETATM' :

            i = i +'  '+  charges_list[idx] +'\n'

            idx += 1

            handle.writelines(i)

        else :

            i = i + '\n'

            handle.writelines(i)

```

Am sure there is a bash or perl one liner for the code above.

... How to merge two files consistently line by line

# Answer
> 1 votes
The existing answer works I'm sure, however there are existing tools that do this or something quite similar.

1. OpenBabel is able to add charges using the option `--partialcharge`, though the docs are not impressive.
2. Autodock vina has a utility script that adds charges and extends the pdb to pdbqt, which is admittedly not exactly what you want.

From the comments: @pippo1980 suggests this page which shows how to use the openbabel Python API to set charges

Openbabel models **(note that it includes `fromfile` meaning you can supply your own)**:

```
obabel -L charges
eem    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Bultinck B3LYP/6-31G*/MPA
eem2015ba    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/AIM
eem2015bm    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/MPA
eem2015bn    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf B3LYP/6-311G/NPA
eem2015ha    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/AIM
eem2015hm    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/MPA
eem2015hn    Assign Electronegativity Equilization Method (EEM) atomic partial charges. Cheminf HF/6-311G/NPA
eqeq    Assign EQEq (charge equilibration) partial charges.
fromfile    Assign charges from file containing {'atom-name', charge} pairs
gasteiger    Assign Gasteiger-Marsili sigma partial charges
mmff94       Assign MMFF94 partial charges
none    Clear all partial charges
qeq    Assign QEq (charge equilibration) partial charges (Rappe and Goddard, 1991)
qtpie    Assign QTPIE (charge transfer, polarization and equilibration) partial charges (Chen and Martinez, 2007)

```

---
Tags: pdb
---


title: issue with rna seq analysis

text:
# Question
Title: issue with rna seq analysis
I am working on RNA seq analysis and I would like to know the following things:

**Current Methods:**

* downloaded genome fasta file for non-coding rna from ensembl and got the gtf file for hg38
* performed hist2 and got 17% alignment for my sample against ncrna genome
* sorted the file and then ran feature counts with the hg38 gtf file but got zero counts

**Question:**

* what could be the issue that is causing zero counts?

My bam file after sorting has the SRR ID at first whereas the GTF file has chromosome no: . Could this be the issue?

# Answer
Several issues:

1. There is no non-coding RNA genome. There are non-coding genes but the genome fasta files from Ensembl are always the entire genome and this is what you should align against. Probably you downloaded a transcriptome fata file.
2. If you aligned against a transcriptome fasta file (like ncRNAs) then it makes no sense to use a GTF file since GTF files have genome coordinates but you mapped against a selection of a transcriptome fasta, hence you have (specific) transcriptome coordinates. GTF and these coordinates are not compatible.

Solution: Align your data against the entire genome (hg38), then use the GTF to get counts for all annotated genes, then do subsetting to specific genes if that is what you need.
> 3 votes

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


title: stacked barplot with error bars

text:
# Question
Title: stacked barplot with error bars
I need to create a barplot with error bars, my data consists of 2 groups called "Genot", with 3 sub-groups "Suberin" for each, associated with a value "Percent".

```
summary_data <- FY %>%
  group_by(Genot, Suberin) %>%
  summarise(
    Mean_Percent = mean(Percent),
    SD_Percent = sd(Percent))

# plot
ggplot(summary_data, aes(x = Genot, y = Mean_Percent, fill=Suberin)) +
  geom_bar(stat = 'identity', aes(fill = Suberin)) +
  geom_errorbar(data = summary_data,
                aes(x = Genot, ymax = Mean_Percent + SD_Percent, ymin = Mean_Percent - SD_Percent),
                , position = "identity",width = 0.2)

```

My main problem is that the error bars are placed from 0, but I have to order my data as the following to get the graph in the order I want.

```
FY$Suberin <- ordered(FY$Suberin,
                    levels = c("Con", "Patchy","Non" ))

```

I can't figure out a way to describe the positions

# Answer
> 1 votes
(Following from comments)

I am 99% sure that, in the second code block, you are reordering only that one column and the rest of the columns are remaining unordered. That means that that column will be permuted with respect to the rest of the data, which I doubt is what you intend. I'd suggest instead running this, which will reorder the whole dataframe:

```
# not actually run cause I don't have the data, check it!
FY <- FY[
         order(factor(FY$Suberin,
               levels = c("Con", "Patchy","Non" )
         ),
]

```

Consider this analogous example:

```
> test = data.frame(
  "one"=c(1,2,3),
  "two"=c(6,5,4),
  "three"=c("this", "that", "which")
)
> test
  one two three
1   1   6  this
2   2   5  that
3   3   4 which
> test$three = ordered(test$three, levels=c("that", "this", "which"))
# only column "three" is reordered!!!
> test
  one two three
1   1   6  this
2   2   5  that
3   3   4 which

# THIS reorders the whole dataframe
> test[order(test$three),]
  one two three
2   2   5  that
1   1   6  this
3   3   4 which

# impose a specific ordering by making it a factor, maybe there's a better way
> test[order(factor(test$three, levels=c("that", "which", "this"))),]
  one two three
2   2   5  that
3   3   4 which
1   1   6  this

```

I don't know that that impacts your error bar issue. But the underlying data will certainly be affected and possibly made inaccurate by permuting only one column.

Another thing that's a little weird is that you're computing `Mean_Percent` in two places: first, in your `summarise` expression, and then again when you actually make the plot with `geom_bar(stat="identity"...`, which will fill the axes you've set up with newly estimated values based on the relative frequencies of the `Suberin` column. So you aren't plotting the `Mean_Percent` values that you are getting from `summarise`, as far as I can tell.

I'm a little scared of the way that you are computing mean/sd percentages and stacking them, especially not knowing what the data look like. I'd recommend confirming that your `Mean_Percent` that you compute with `summarise` adds up to 100%; for the actual bars in the barplot you are letting `geom_bar` recompute the proportions directly from the data, so I don't think the fact that it apparently adds up to 100% in the plot is trustworthy. To summarize: it'd make sense to confirm that what's in your dataframe matches what's in the plot since they are separately computed.

I don't know exactly that either of these issues is leading to your error bars having the wrong central tendencies, but investigating a couple of these lines may lead you to discover something.

I'm not a `ggplot` person so I can't diagnose just by looking at your code, and obviously I don't have your data so I can't tinker with it. Someone with more experience there might be a good idea. However, there are some related questions on SO that may be worth reading over.

Note again that none other than Hadley Wickham finds this kind of plot confusing.

---
Tags: r, ggplot2
---


# Question
Title: Kallisto RNA-seq errors because of incompatible indices when running <kallisto quant> using downloaded mouse index
I'm trying to run RNA-seq on single-end mouse RNA data (50bp reads) but have very limited RAM (\<4 gb) - I don't think I really have a way to work around this memory restriction, and wanted to try Kallisto's pseudoalignment (I conda installed it) as indexing and/or aligning with other tools like STAR and HISAT2 have resulted in the processes being killed. Building the index with `kallisto index` still results in the process being killed, so I downloaded the index files from https://github.com/pachterlab/kallisto-transcriptome-indices/releases/tag/v1 instead. I ran `kallisto quant -i index.idx -o output --single -l 180 -s 20 sample1.fastq` (I don't entirely know what the mean and sd arguments should be for my data), but this resulted in the following:

```
[quant] fragment length distribution is truncated gaussian with mean = 180, sd = 20
Error: incompatible indices. Found version 13, expected version 10
Rerun with index to regenerate

```

I've seen that other people have come across this error before when quantifying with kallisto, but it seems that they are building the index themselves while I was trying to use the pre-generated index from the unpacked `mouse_index_standard.tar.xz` file downloaded from the link above.

How can I get around this error? And in general, are there any suggestions on how to get raw/normalized count data from RNA-seq data with such memory limitations? I'm rather new to this.

# Answer
It means that the version of the tool used for indexing and your current version are different and not compatible. Either downgrade your version to the version used for indexing, or build a new index with your current version. As by the GitHub repo, these indices that you link `(supports kallisto versions >= 0.50.1)`, so use this version if you want to use the index.

If that does not work then I would open an issue at their GitHub.
> 3 votes

---
Tags: rna-seq, kallisto
---


# Question
Title: How can I convert a BED file to GTF/GFF with gene_ids?
Given a .bed (BED12), how can I convert it to GTF/GFF formats with gene\_id attributes? What is the fastest way or available tools to do it?

For example, given an input like this:

```
chr27 17266469 17281218 ENST00000541931.8 1000 + 17266469 17281218 0,0,200 2 103,74, 0,14675,

```

How can I convert it to:

```
chr27 bed2gtf gene 17266470 17285418 . + . gene_id "ENSG00000151743";
chr27 bed2gtf transcript 17266470 17281218 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8";
chr27 bed2gtf exon 17266470 17266572 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8"; exon_number "1"; exon_id "ENST00000541931.8.1";

```

# Answer
Here I provide an ordered list of options (note that I am the author of `bed2gtf` and `bed2gff`):

## bed2gtf

A high-performance BED-to-GTF converter written in Rust from https://github.com/alejandrogzi/bed2gtf.

```
Usage: bed2gtf[EXE] --bed/-b <BED> --isoforms/-i <ISOFORMS> --output/-o <OUTPUT> --threads/-t <THREADS>

where:
--bed <BED>: a .bed file
--isoforms <ISOFORMS>: a tab-delimited file
--output <OUTPUT>: path to output file (*.gtf)

```

The isoforms file specification:

> a tab-delimited .txt/.tsv/.csv/... file with genes/isoforms (all the transcripts in .bed file should appear in the isoforms file):

```
> cat isoforms.txt

ENSG00000198888 ENST00000361390
ENSG00000198763 ENST00000361453
ENSG00000198804 ENST00000361624
ENSG00000188868 ENST00000595977

```

Converts

* Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 3.25 seconds.
* Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 1.99 seconds.
* Canis lupus familiaris ROS\_Cfam\_1.0 Ensembl 110 (55,335 transcripts) in 1.20 seconds.
* Gallus galus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.36 seconds.

## bed2gff

A Rust BED-to-GFF3 translator that runs in parallel from https://github.com/alejandrogzi/bed2gff.

```
Usage: bed2gff[EXE] --bed/-b <BED> --isoforms/-i <ISOFORMS> --output/-o <OUTPUT> --threads/-t <THREADS>

where:
--bed <BED>: a .bed file
--isoforms <ISOFORMS>: a tab-delimited file
--output <OUTPUT>: path to output file (*.gff)

```

The isoforms file specification:

> a tab-delimited .txt/.tsv/.csv/... file with genes/isoforms (all the transcripts in .bed file should appear in the isoforms file):

```
> cat isoforms.txt

ENSG00000198888 ENST00000361390
ENSG00000198763 ENST00000361453
ENSG00000198804 ENST00000361624
ENSG00000188868 ENST00000595977

```

Convert

* Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 4.16 seconds.
* Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 2.15 seconds.
* Canis lupus familiaris ROS\_Cfam\_1.0 Ensembl 110 (55,335 transcripts) in 1.30 seconds.
* Gallus gallus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.51 seconds.

## bedToGenePred + genePredToGtf + refTable

UCSC offers a fast way to convert BED into GTF files through KentUtils or specific binaries using:

```
bedToGenePred in.bed /dev/stdout | genePredToGtf file /dev/stdin out.gtf

```

You can install these tools with bioconda, or download them here. The gene\_id is only achieved when using refTables (a format specified in UCSC's web browser), you can see a more elaborate answer here Obtaining Ucsc Tables Via Ftp And Converting Them To Proper Gff3 Via Genepredtogtf?.

## Other options

Other scripts/tools **That DO NOT produce a complete GTF file (lacking gene\_id attributes)** are:

* gtf2bed

```
gtf2bed < foo.gtf | sort-bed - > foo.bed
awk '{print $1"\t"$7"\t"$8"\t"($2+1)"\t"$3"\t"$5"\t"$6"\t"$9"\t"(substr($0, index($0,$10)))}' foo.bed > foo_from_gtf2bed.gtf

```

-kscript from https://github.com/holgerbrandl/kscript:

```
kscript https://git.io/vbJ4B my.bed > my.gtf

```
* pfurio/bed2gtf

from https://github.com/pfurio/bed2gtf:

```
python bed2gtf [options] <mandatory>

```
* AGAT

AGAT

Considering only the options that produce gene\_ids attributes, bed2gtf and bed2gff are faster by ~3-4 seconds than UCSC's C binaries. More detailed instructions of this tools are explained in the sources linked.
> 6 votes

---
Tags: bed, format-conversion, gtf, gff
---


title: MMSeqs taxonomy running for over a day

text:
# Question
Title: MMSeqs taxonomy running for over a day
I've been trying to run mmseqs2 on a few metagenomic assemblies and despite my best efforts in reading the wiki and playing with parameters, the process is taking over a day.

In their paper they claim to be 2 - 18x faster than CAT. I ran CAT on the same datasets and they took at most a few hours. So I can only assume that the long run time for mmseqs is down to user error (i.e. me).

I am currently running the commands (general format):

```
mmseqs createdb metagenome1_scaffolds.fasta metagenome1_query_DB
mmseqs taxonomy metagenome1_query_DB /home/mmseqs2/NR/NR metagenome1_mmseqs_taxonomy.txt metagenome1_temp --tax-lineage 1 --orf-filter 1 --threads 24

```

Could someone please let me know where I'm going wrong or what parameters I can tweak to get this running in the expected time frame?

# Answer
After much googling about I realised that I was missing a step! Here are the commands I used to get it running:

```
mmseqs createdb metagenome1_scaffolds.fasta metagenome1_query_DB
mmseqs taxonomy metagenome1_query_DB /home/mmseqs2/NR/NR metagenome1_mmseqs_taxonomy.txt metagenome1_temp --tax-lineage 1 --orf-filter 1 --threads 24 --split-memory-limit 500G
mmseqs createtsv metagenome1_DB metagenome1_mmseqs_taxonomy.txt metagenome1_mmseqs_taxonomy_final.txt

```

The split-memory-limit flag sped up the process and the final line put all the split files back together
> 1 votes

---
Tags: assembly, metagenome, taxonomy
---


# Question
Title: mrbayes, Sumt ; Error when setting parameter "Filename" (1)
## question 1

I wanted to use sumt command after completion of mcmc run. Whenever I use, `sumt Filename=infile.nex`, it gives me the error:

```
A matrix must be specified before sumt can be used
Error when setting parameter "Filename" (1)

```

But for the same analysis, I can use sump without any problem. What should I do now?

## question 2

Another issue, I used burn in as 25% where as in Tracer, it is always showing 10%. No matter how much burn in I use, Tracer always shows 10% as burn in.

### system stats

```
Operating system : Centos7
Version of MrBayes: 3.2.6

```

Thanks in advance, Best Regards, Zillur

# Answer
> 2 votes
Having thought this through, I think you've made a mistake. The file it wants is the outfile, you have given it the infile only. My original answer is below.

====================

From memory (my personal memory, not RAM, SSD or Stata), in the run command, you need to specify,

```
nruns=2

```

sumt will then combine both runs. If the file name is mystuff.run and nruns=2, there will be an output of mystuff.run1.t and mystuff.run2.t.

If nruns isn't set then sumt will not work, or if the MCMCMC output ends up with different names again sumt will not work.

With the burnin, you simply specify this in MrBayes. It is easy to write a script that will perform a 25% burnin. Tracer is nicely written, presumably C, its fast, but it has been around for a long time.

# Answer
> 1 votes
This is a bit late. However, I am hoping someone else will find this useful since I had the same problem.

The data matrix is your infile.nex file (the nex file you supplied for the analysis)

1. Bracket out the mcmc command in your infile.nex
    * Prevents the file from being run again.
    * May not be present.
2. In Mr. Bayes: Execute infile.nex
    * Loads you infile as the data matrix
3. In Mr. Bayes: sumt Filename=infile.nex nruns=2
    * If you have nruns set in you Mr. Bayes block you can skip this.

# Answer
> 0 votes
This question was posted and answered as a Github Issue in the MrBayes repository; I am copying the answer here for better visibility.

---

> Dear zillurbmb51,
>
> Thank you for your question. You need to load the data matrix in memory before issuing the `sumt` command, or in other words: "A matrix must be specified before sumt can be used".
>
> **Example 1** (where you haven't yet run a MCMC):
>
> Start mb, then type:
>
> ```
> exe infile.nex
> mcmc
> sumt
>
> ```
>
> **Example 2** (where you already did a MCMC, and have `.t` files in your current working directory, and they are named `infile.nex.run1.t`, `infile.nex.run2.t`, etc):
>
> Start mb, then type:
>
> ```
> exe infile.nex
> sumt filename=infile.nex
>
> ```
>
> Cheers
>
> Johan

---
Tags: phylogenetics, phylogeny, mrbayes
---


title: BioPython bootstrap is not reliable?

text:
# Question
Title: BioPython bootstrap is not reliable?
Here i will show you a minimal working example of code and as you can see the support values for the tree is always 100.

I am using synthetic sequences of 100bp for 6 elements. The sequences have been generated at random choosing from ATCG for each position with equal probabilities.

The input fasta file is

```
>S1
CTCAGGCACTAGGAGCTTCTCCAGGGCAAAGTTGTCTACGAATATGCCGACTCAGAAGGTTATCAATACGGTTTACTTATCTGCACGCCATTTTCCTATG
>S2
TAACAAATGTTTCTCGCTCGAACCCCGTGGTGCGAGGGGTCACGATAAAAGGTCCCTCCTTCCACGGCATAATTGCTCCCTTTCTTTTCCGTGGGCAGGA
>S3
TTCCTAAAGGACGACTGGAAGCCGGGTACACGTCAACGGAGCTTCTAGCCCGAGGTCTATAGACCGGTTAATGAACTAGACTTAATGTGGAGCTCGTAGA
>S4
ACTTCGTGCCCAACCACTCTATGAGCAGAGTGCTCGAATGAACCTCAAAAGGATATCCGCATTTACTCTATATAACAAACGGCCGTCCCCCCATCTGTCA
>S5
TGACACTGGGTCATTTTACCCGCACTGATCGCTGGGCAGGTCGGCAATTTCGTCAGAATGCCGTGGCCGCACTGAAAATATCATTACCGGACTGAGTATG
>S6
ATAAGCCGAGGGGTAGCCCTCTATTTCGCACCGTATAGACGAAGTGATAAACTTTCTAACCACGTGTCGCCTATCTCTACCTAGCCACATTTGAGGTGCG

```

The following python code makes the multiple alignment using MUSCLE algorithm

```
### alignment
from Bio.Align.Applications import MuscleCommandline
from Bio import AlignIO, SeqIO
import os

muscle_exe = "C:/Scripts/muscle3.exe"
path = os.getcwd()
with open('tmp_in.fa', 'w') as fp:
    fp.write(fasta)

## replacing \\ with /
inp_fp = path.replace('\\', '/') + '/tmp_in.fa'
out_fp = path.replace('\\', '/') + '/tmp_out.afa'

muscle_cline = MuscleCommandline(muscle_exe, input=inp_fp, out=out_fp)

```

The following code block takes the previous alignment `out_fp` and creates a tree using UPGMA algorithm.
The number of bootstrap trees is set to 50 and the consensus is majority based.

```
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio.Phylo.Consensus import *

# read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')

calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)

## get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)

print(consensus_tree)

```

The final consensus tree looks like this

```
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.350925925925926, name='S5')
        Clade(branch_length=0.05833333333333341, confidence=100.0)
            Clade(branch_length=0.05879629629629617, confidence=100.0)
                Clade(branch_length=0.26481481481481484, name='S3')
                Clade(branch_length=0.020370370370370372, confidence=100.0)
                    Clade(branch_length=0.24444444444444446, name='S1')
                    Clade(branch_length=0.24444444444444446, name='S4')
            Clade(branch_length=0.031018518518518404, confidence=100.0)
                Clade(branch_length=0.2925925925925926, name='S6')
                Clade(branch_length=0.2925925925925926, name='S2')

```

My question is: why does the tree have all support values equal to 100? I am sure there must be an error either on my side or in the package module itself.

# Answer
> 2 votes
I think this is a bug.

It seems to work if you do this, creating an equivalent Alignment object instead of a MultipleSeqAlignment to give the bootstrap step:

```
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)

```

---

`bootstrap_consensus` calls `bootstrap_trees`, which makes however many randomly shuffled alignments you asked for. But if the input is a MultipleSeqAlignment object like you're using it builds a tree over and over with *that* same input to the function as the tree function input, so it always gets the same tree over and over. The other branch of `bootstrap_trees`, which runs if you have an Aignment object, builds trees based on the randomly-selected columns of the input, like you'd expect. Pretty sure that's just a variable naming screw-up one of us should report to Biopython.

Self-contained version of your tree part:

```
#!/usr/bin/env python
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio import AlignIO
from Bio.Phylo.Consensus import *

out_fp = "tmp_out.fa"

# read the multiple sequence alignment (msa)
with open(out_fp,'r') as afa:
    alignment = AlignIO.read(afa, 'fasta')

calculator = DistanceCalculator('identity')
constructor = DistanceTreeConstructor(method='upgma', distance_calculator=calculator)

## get the list of bootstrap trees
consensus_tree = bootstrap_consensus(alignment=alignment, times=50, tree_constructor=constructor, consensus=majority_consensus)

print(consensus_tree)

# Give bootstrap_consensus an Alignment object instead of MultipleSeqAlignment
# object
from Bio.Align import Alignment
alignment2 = Alignment(list(alignment))
consensus_tree = bootstrap_consensus(alignment=alignment2, times=50, tree_constructor=constructor, consensus=majority_consensus)
print(consensus_tree)

```

Example outputs for the two different input alignment objects:

```
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.012239583333333349, confidence=100.0)
            Clade(branch_length=0.31666666666666665, name='S3')
            Clade(branch_length=0.033333333333333284, confidence=100.0)
                Clade(branch_length=0.2833333333333333, name='S2')
                Clade(branch_length=0.2833333333333333, name='S4')
        Clade(branch_length=0.033072916666666716, confidence=100.0)
            Clade(branch_length=0.29583333333333334, name='S5')
            Clade(branch_length=0.03750000000000007, confidence=100.0)
                Clade(branch_length=0.2583333333333333, name='S1')
                Clade(branch_length=0.2583333333333333, name='S6')
Tree(rooted=True)
    Clade()
        Clade(branch_length=0.2583333333333333, name='S2')
        Clade(branch_length=0.03405145202020203, confidence=66.0)
            Clade(branch_length=0.275, name='S5')
            Clade(branch_length=0.049892241379310344, confidence=58.0)
                Clade(branch_length=0.21250000000000002, name='S1')
                Clade(branch_length=0.21250000000000002, name='S6')
        Clade(branch_length=0.039416666666666676, confidence=50.0)
            Clade(branch_length=0.2583333333333333, name='S4')
            Clade(branch_length=0.321875, name='S3')

```

# Answer
> 0 votes
Initially, 100 looks like a good value. However, if you lower to bootstrap value, let's say maybe to 50 as you did, you will probably get an erroneous support value. Bootstrapping (you probably already know this) is how many times a branch in a phylogenetic tree is observed. A good bootstrapping value would be anywhere between 80 and 100. 50 is not a good bootstrapping value. This means the branches are only being observed 50 times as opposed to 85, 90, etc... If I were you, try changing the bootstrapping value of 50 to at least 80. See what that gets you. You may be satisfied with it, but you may want to increase it depending on your output. By doing this, you will see a more realistic and phylogenetically accurate output of your tree. Most bioinformaticians use 100 as a bootstrap value. The only issue, I wouldn't really call it an issue, but it will take longer for an output to occur due to computing time and computing memory available. BUT you will still get an output, you just may need to wait a bit. My bioinformatics background is rooted in phylogenetic and horizontal gene transfer research involving phages and bacteria, so I have seen this issue a lot and have experienced myself a lot! I hope this information provides some insight to your project and bioinformatics journey! :)

---
Tags: sequence-alignment, multiple-sequence-alignment, biopython
---


# Question
Title: What are the applications of Tries(data structure) of an ordered sequence of strings in bioinformatics?
*This question was also asked on reddit*

Tries are a data structure that can be used to efficiently store and search for strings. Tries created from an *ordered sequence of strings* differ from the regular Tries in the following way:

If there is a string of a certain character $w$ at a character position $x$ and is between two other strings that have both the same character $p$ where $p \neq w$ at position $x$(and all their characters before position $x$ are the same in all three strings) then the Trie must have an edge for the first occurrence of $p$, an edge for $w$ and an edge for the second occurrence of $p$.

I have recently written a paper regarding a novel data structure that is related to tries and I wish to connect it to some applications in bioinformatics since the paper itself presents only one application related to string algorithms.

Since bioinformatics makes some use of string algorithms and since the structure is related to Tries I wish to ask you about the possible applications of Tries so I may find some application to my related data structure.

What are the applications of Tries(data structure) of an *ordered sequence of strings* in bioinformatics?

# Answer
> 1 votes
(*Should be a comment not an answer, but I weary of SE comment formatting- feel free to express disapproval in whichever way*)

When in doubt, google it and look at wikipedia:

> Tries are used in Bioinformatics, notably in sequence alignment software applications such as BLAST, which indexes all the different substring of length k (called k-mers) of a text by storing the positions of their occurrences in a compressed trie sequence databases.\[25\]

---
Tags: phylogenetics, sequence-alignment, genomics, algorithms
---


# Question
Title: MEGA11 gives the error "stop codons are found in the translated sequences", even though there are not!
I am using MEGA11 MEGA11: Molecular Evolutionary Genetics Analysis Version 11 , Mol. Biol. Evol. 38(7):3022–3027 to align protein coding sequences that I downloaded from KEGG database. I am doing MUSCLE alignment and I choose the "align codons" option. The problem that occurs has to do with the error "stop codons are found in the translated sequences" that pops as an warning during the process. When I check the sequences I cannot find any stop codons except the one at the end, as expect. I don't know how to interpret the message. Does anyone has any idea what is going wrong? Thank you in advance!

# Answer
> 1 votes
The solution is very simply, MegaX will not allow any stop codons within the "translated" sequence. The reason is simple

* MEGAX's function is phylogenetic models with nucleotide or protein data
* All phylogenetic protein calculations require a 20x20 matrix of amino acid frequencies, e.g. BLOSUM, Dayhoff, mtREV24 etc ... No matrix has ever calculated stop codon frequency (well NCBI have and sort of obviously always '1' or transformed '1' \[the definition of a gene is it has one stop codeon\])

Thus MEGAX will not permit a single stop codon, even the terminal codon of the gene, because it can't calculate phylogenetic information from them (doesn't fit the model). Thus the "getter" to pick up bugs occurs right at that beginning.

Two solutions:

1. Don't translate the nucleotide sequence on import - just stick with nucleotides (there are no codon models in MEGAX).
2. Remove the terminal 3 resides

```
     perl -p -i -e 's/(^[AGCT-NRY]{3}.*)[AGCT-NRY]{3}$/$1/g' alignment.fa

```

Note the above code is not tested and intended only for fasta sequence or MEGA input format - don't use it on Clustal or Phylip it will be disastrous, i.e. interleaved format. IN ADDITION make sure there are no carriage returns within the contiguous sequence (again the result is disastrous). Finally, note I forget the `{3}` stuff (its ages since I wrote perl) and should include all the degenerate nucleotide code beyond NRY, e.g. MB etc ...

Really trimming the 3' last 3 nucleotides is a good use of Biopython - no risk of errors regardless of alignment. Trimal can probably do it BTW or seqkit.

---

I'm still interested in how you interrogated KEGG gene DB BTW

---
Tags: errors, kegg, codon, megax
---


# Question
Title: single_to_multi_fast5 does not collect all the single files if the input folder contains mixed types of fast5 files
I have a dataset that contains thousands of mixed multiple and single fast5 files in a non-homogenous folder structure. The reads and fast5 files are coming from different labs and some labs gave single fast5 files and some multi fast5 files, without clear mappings. My final goal is to convert all the fast5 files to multi fast5 files.

There are commands in ont\_fast5\_api which is an interface to HDF5 files of the Oxford Nanopore. Using this API we can convert single fast5 files and multi fast5 files to each other. this is the link to this API github page: https://github.com/nanoporetech/ont\_fast5\_api

So I made a test mixed folder to check if I can use the API without missing reads as a result of conversion. The test folder structure is similar to my original folder structure:

* mixed
    + Subfolder\_A
        - single\_fast5 x 4000
    + Subfolder\_B
        - multiple\_fast5 (contains 4000 single reads fast5)
    + single\_fast5 x 2
    + multiple\_fast5\_1 (contains 4000 single reads fast5)
    + multiple\_fast5\_2 (contains 4000 single reads fast5)
    + multiple\_fast5\_3 (contains 100 single reads fast5)

Here is a screenshot of the structure with a few files from each folder:

and I created two other folders outside of the mixed folder as follows:

* single
* multi

Here my question is about *single\_to\_multi\_fast5* and *multi\_to\_single\_fast5* commands.

```
test_path=mixed
multi_path=multi
single_path=single
multi_to_single_fast5 -i $test_path/ -s $single_path/ --recursive

```

in the $single\_path, I expect to see, all of the reads in multi files converted to single files (12100 reads in 12100 single fast5 files), but I see 8100 single files.

So I have 4000 missing reads.

Now I convert single read fast5 file to multi:

```
single_to_multi_fast5 -i $test_path/ -s $multi_path/ --filename_base $output_name --batch_size 1000 --recursive

```

in the $multi\_path folder, I expect to see that generated multi files contain all my single-read files (4000 + 2), but I see 3006 in the mapping summary.

So again I missed 996 single reads in my multi files. Since I made a small subsample of data, I am sure that the files are not overwritten.

Because my final goal is **converting everything to multi** I was going to do the following **pipeline**:

```
orig_path=original_folder
intermediate_path=intermediate
final_path=final_multi
# converts/collects all multi files in the original folder
multi_to_single_fast5 -i $orig_path/ -s $intermediate_path/ --recursive
single_to_multi_fast5 -i $intermediate_path/ -s $final_path/ --filename_base $output_name --batch_size 1000 --recursive

# converts/collects all single files in the original folder
single_to_multi_fast5 -i $orig_path/ -s $final_path/ --filename_base $output_name --batch_size 1000 --recursive

```

However because it seems I get missing reads during the process, I cannot use this pipeline.

At the same time, my original dataset is super huge, I cannot check if individual files are single or multi. It would take ages (I guess). Is there a solution for solving this problem besides checking every file? Maybe my pipeline is not quite suitable for what I want to do.

# Answer
It looks like you're creating too much additional work for yourself. What I'm trying to lead to with my comments is that the file names should help in working out if a file represents single reads or multiple reads, in which case the single reads can be shifted into a separate folder.

Based on the file names that I would expect, a search for files containing '\_read\_\[0-9\]' should be sufficient to fish out the single fast5 files:

```
find . -name '*_read_[0-9]*' > single_reads.txt
mkdir -p single_reads
cat single_reads.txt | while read fname; do mv -i ${fname} single_reads; done

```

*\[That `single_reads.txt` file could be split into multiple files if directory entries get too large.\]*

After which you should have a single folder that contains thousands of single-read fast5 files, and other folders that contain multi-fast5 files, removing the problem of needing to deal with mixed file types.
> 1 votes

---
Tags: nanopore, fast5
---


title: Fetching uniprot seq using bash command

text:
# Question
Title: Fetching uniprot seq using bash command
I am trying to get a protein sequence using this command:

```
curl https://www.uniprot.org/uniprot/A2Z669.fasta

```

But I don't get any output:

```
$ curl https://www.uniprot.org/uniprot/A2Z669.fasta
$

```

And no file is created.

# Answer
The URL https://www.uniprot.org/uniprot/A2Z669.fasta redirects to https://rest.uniprot.org/uniprotkb/A2Z669.fasta. You can see this if you open it in your browser: clicking on the first link will take you to the address of the second.

Now, `curl` doesn't follow redirects by default, you need the `-L` option for that (see `man curl`):

> -L, --location
>
> (HTTP) If the server reports that the requested page has moved to a different location (indicated with a Location: header and a 3XX response code), this option makes curl redo the request on the new place. If used together with -i, --include or -I, --head, headers from all requested pages are shown.

So you need `-L` for `curl` to work here:

```
$ curl -L https://www.uniprot.org/uniprot/A2Z669.fasta
>sp|A2Z669|CSPLT_ORYSI CASP-like protein 5A2 OS=Oryza sativa subsp. indica OX=39946 GN=OsI_33147 PE=3 SV=1
MRASRPVVHPVEAPPPAALAVAAAAVAVEAGVGAGGGAAAHGGENAQPRGVRMKDPPGAP
GTPGGLGLRLVQAFFAAAALAVMASTDDFPSVSAFCYLVAAAILQCLWSLSLAVVDIYAL
LVKRSLRNPQAVCIFTIGDGITGTLTLGAACASAGITVLIGNDLNICANNHCASFETATA
MAFISWFALAPSCVLNFWSMASR

```

Or, alternatively, you can use `wget` which does follow redirects by default:

```
$ wget "https://www.uniprot.org/uniprot/A2Z669.fasta"
--2023-12-10 16:21:46--  https://www.uniprot.org/uniprot/A2Z669.fasta
Loaded CA certificate '/etc/ssl/certs/ca-certificates.crt'
Resolving www.uniprot.org (www.uniprot.org)... 193.62.193.81
Connecting to www.uniprot.org (www.uniprot.org)|193.62.193.81|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://rest.uniprot.org/uniprot/A2Z669.fasta [following]
--2023-12-10 16:21:46--  https://rest.uniprot.org/uniprot/A2Z669.fasta
Resolving rest.uniprot.org (rest.uniprot.org)... 193.62.193.81
Connecting to rest.uniprot.org (rest.uniprot.org)|193.62.193.81|:443... connected.
HTTP request sent, awaiting response... 301 Moved Permanently
Location: https://rest.uniprot.org/uniprotkb/A2Z669.fasta [following]
--2023-12-10 16:21:46--  https://rest.uniprot.org/uniprotkb/A2Z669.fasta
Reusing existing connection to rest.uniprot.org:443.
HTTP request sent, awaiting response... 200 OK
Length: 314 [text/plain]
Saving to: ‘A2Z669.fasta.1’

A2Z669.fasta             100%[=======================================>]     314  --.-KB/s    in 0s

2023-12-10 16:21:47 (8.57 MB/s) - ‘A2Z669.fasta.1’ saved [314/314]

```

That creates a `A2Z669.fasta` file in the current directory which will contain the sequence.
> 1 votes

---
Tags: uniprot
---


# Question
Title: How to download Sequencing data on Windows using SRA toolkit?
I have downloaded and installed the SRA toolkit, but there seems to be no online material on how to download SRR files using the SRA toolkit on Windows.

All online materials are Ubuntu or Linux based.

Is there some obvious resource that I am missing?

# Answer
I found a solution to this problem.

After having installed SRA toolkit:

`prefetch -v SRR9....1`

this helps download the SRA file.

Then to extract it as FASTAQ

`fasterq-dump SRR9....1`
> 1 votes

---
Tags: sratoolkit, sra
---


# Question
Title: How to remove third codon positions from a charset in iqtree?
I need to build a phylogenetic tree using IQ-TREE, starting from a sequence alignment in CODON format of several invertebrate mitochondrial genes. These are my charsets:

```
nexus
BEGIN SETS;

charset atp6.fasta = 1-651;
charset atp8.fasta = 652-729;
charset cob.fasta = 730-1749;
charset cox1.fasta = 1750-3237;
charset cox2.fasta = 3238-3858;
charset cox3.fasta = 3859-4632;
charset nad1.fasta = 4633-5484;
charset nad2.fasta = 5485-6156;
charset nad3.fasta = 6157-6459;
charset nad4.fasta = 6460-7683;
charset nad4L.fasta = 7684-7899;
charset nad5.fasta = 7900-9294;
charset nad6.fasta = 9295-9663;
END;

```

I would need to remove the third bases, how should I write them? These are the commands with which I would launch IQ-TREE:

```
iqtree -s "fasta_file" -bb 1000 -alrt 1000 -m MFP+MERGE -spp "charset.nexus" -st CODON5

```

Please, any suggestions would be very helpful.

# Answer
Yes, it's easy.

```
#nexus
BEGIN SETS;

charset atp6.fasta = 1-651\3, 2-651\3;
charset atp8.fasta = 652-729\3, 653-729\3;
charset cob.fasta = 730-1749\3, 731-1749\3;
# ... so on for every gene
END;

```

So for atp6 you will get 1,2,4,5,7,8,10,11 .... up to 651, and so on for each gene. Thus the residues at position 3, 6, 9, 12 are excluded.
> 0 votes

---
Tags: phylogenetics, phylogeny, fasta, mtdna, iqtree
---


# Question
Title: Featureplot by co-expression of some genes
I want to identify `MDSCs` cells which don't not express `HLA-DRs` (HLA-DR negative) but express `CD14`, `CD33`, and `ITGAM`

For this, I firstly plotted a feature plot to locate `HLA-DR` by

```
cd.features <- list(grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE))

pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')

FeaturePlot(object = pbmc_small, features = 'CDFeatures1')&
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

```

Which returned

For `CD14`, `CD33`, and `ITGAM` I done

```
> cd.features=list(c("CD14", "CD33",  "ITGAM"))
> pbmc_small <- AddModuleScore(object = pbmc_small, features = cd.features, ctrl = 5, name = 'CDFeatures')
> FeaturePlot(object = pbmc_small, features = 'CDFeatures1')&
+     scale_colour_gradientn(colours = rev(brewer.pal(n = 11, name = "RdBu")))

>

```

Which returned

By looking at these two plots I really can not say which part of plot is HLA-DR negative but expressing `CD14`, `CD33`, and `ITGAM`

Do you know a way to do this?

Plotting options don't help because looking at the below plots you can not say which cells are HLA negative and CD14 positive

# Answer
Very old question, but recently modified, which means that someone might still be looking into this.

I usually prefer to use Ucell as the signatures are calculated within cells, which makes comparisons a bit easier. In addition is allows the formulation of signatures with both `positive` and `negative` genes.

```
hladr_genes <- grep(pattern = '^HLA-DR', x = rownames(x = pbmc_small), value = TRUE)
signatures_list <- list(MDSCs = c("CD14+","CD33+","ITGAM+",paste0(hladr_genes,"-")))

> signatures_list
$MDSCs
[1] "CD14+"     "CD33+"     "ITGAM+"    "HLA-DRA-"  "HLA-DRB5-" "HLA-DRB1-"

# calculate signature
pbmc_small <- AddModuleScore_UCell(pbmc_small,signatures_list,ncores = 8,name = "")
# perform smoothing
pbmc_small <- SmoothKNN(pbmc_small, signature.names = names(signatures_list), reduction="pca")

```
> 2 votes

# Answer
You need to make a transformed variable to the gene sets of interest. I would rescale the data on a range of 0-1. Then you can make a log2 ratio of the sum(set1)/sum(set2). You could explore other transformations, but this will give contrasting values for high set1 and high set2, but neutral if the expression is similar between the sets.
> 1 votes

---
Tags: r, scrnaseq, seurat
---


# Question
Title: PCR primers for detecting hetero vs homozygous and WT mice in our colony
My PI has asked me for advice on making a PCR screen to identify mice from a particular Cre line. We have multiple Cre lines in the colony so it can't just be for mice with a copy of Cre in their genome the primer. The insertion site is known for the mutation, is there a way to use blast or some other search tool to pull the sequences we should base our primers on?

We also want to identify which mice are hetero vs homozygous. Are there particular kinds of PCR we should use or will a well designed primer enable this?

---

**Details**

The strain: https://www.jax.org/strain/004146

I was thinking of using this tool: https://www.ncbi.nlm.nih.gov/tools/primer-blast/index.cgi?LINK\_LOC=BlastHome

The strain reference from Jax, not sure the recommended **cre** primers are going to be sufficient for what we want. So I'm trying to find the FASTA file for this strain in the region of interest, then I can throw that into the BLAST tool for primer design I listed above. Not sure where to search for the FASTA for this strain though. Strain information: https://www.informatics.jax.org/reference/J:66884

Jax recommends this protocol so I suggested it: https://www.jax.org/Protocol?stockNumber=004146&protocolID=20627

# Answer
Primer-blast is okay, but I presume these are lab mice so the variation is low, thus you can just use Primer3. SNPs are detected via TaqMan here, i.e. a separate probe that sits over the mutation site, it will detect heterozygosity via qPCR signal (half strength than homozygous).
> 1 votes

---
Tags: primer, pcr, qpcr, mice
---


# Question
Title: Using Diffusion Map on normalized - or PCA reduced data?
*This question was also asked on GitHub*

In the R package destiny, a Diffusion Map is often performed on the first 50 PCs of the data.

1. Is there a theoretical reason for performing a Diffusion Map on PCs rather than on the data itself ( other than speed)?
2. Is it reasonable to do the diffusion map on the normalized data itself rather than on Principal Components of the data, if that differentiates the states better ?

---

Techinically, I am using the R Bioconductor package destiny. The help page on destiny::DiffusionMap tells "Using e.g. 50 PCs results in more regular looking diffusion maps", but it's not clear compared to what will they be more regularly looking.

# Answer
### The empirical statistical point of view

The technique is uniform manifold approximation and projection (UMAP) and is comparable to t-SNE. UMAP is considered a slight improvement, personally I think they have similar behaviour. They are dimensionality reduction methods.

The reason UMAP or t-SNE is often not applied directly is it is weaker than PCA are capturing the multi-dimensional variance of the data. Both methods can be applied directly to the data, but will not necessarily deliver stronger clusters than PCA (which is the goal of the technique).

However, the weakness of PCA is that the variance can be spread across 50 principle components, making it hard to recover the representative clusters in certain situations. These situations are commonly found in RNA expression analytics, notably where the first 3 or 4 PCs only comprise ~30% of the variance. PCA should be trapping \>80% of the variance within the first 3 or 4 PCs. In expression data (RNAseq) the next 46 PCs might comprise a total of ~30% of the variance and omitting that is mistake, but manually comparing 50 PCs in bivariate plots is not practical. So at best a good unsupervised learning algorithm (UMAP/t-SNE) following PCA in RNA expression via RNAseq might trap 60% of the variance. It's not ideal but it's okay. Thus you can't ignore the 30% of the variance in the bottom 46 PCs and UMAP/t-SNE is recovering that. If 80% of the variance was in the top 4 PCs, then the remaining 20% of the variance in the lower 46 PCs can be ignored and t-SNE or UMAP is not needed.

**Resolution** The way you can capture the variance from all 50 PCs is to perform a clustering method over the top of PCA, this being t-SNE or UMAP.

Summary,

* t-SNE or UMAP without PCA might trap at best 20% of the variance.
* The first 4 PCs in a PCA might trap 30% of the variance (it will be a better result than t-SNE or UMAP).
* Combining PCA and then UMAP/t-SNE will capture 60% of the variance.

60% of the variance, its just about acceptable. Thus thats why this approach is done from an empirical point of view.

### Practical point of view

PCA and then UMAP (or t-SNE) will give tighter, neater clusters (cleaner, less messy groups), but thats not the reason why you are using it.

---

The publications in the link,

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC10312943/ and the PLoS Biology paper https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001532

... both use UMAP across PCA. Its a major approach to cluster resolution.

Diffusion maps are an attempt to extrapolate a variance within Eigenvalues via a random walk, i.e. to defuse within Euclidean space. It is not a primary calculation in my personal opinion because you are sort of introducing more variance into a variance dimensionality reduction calculation. Once you have very tight clusters using PCA/UMAP or PCA/tSNE then you could check to see whether Diffusion map clusters are similar.

This is my perspective only, not to replace PCA with Diffusion map nor to consider a "down-stream" calculation. Its separate calculation used to verify if PCA-UMAP is okay.

If you doing graph theory a Diffusion map could be really important as a null hypothesis and a really important way to obtain a different understanding of the graph (i.e. this is a technical term for a network). Using a diffusion map without any directionality (the papers you cited)\*, I struggle to see this as a primary analysis. In the PLOS Biology paper is part of the appendix.

\*, graph says A-\>B-\>C The publications treat A,B, C as unconnected. If you can hypotheses how A, B and C are connected then diffusion map could be cool. Thats my opinion from the perspective of graph theory.

**Summary** I would definitely choose PCA-UMAP over a diffusion map every time. This is what has happened in the two publications cited. The diffusion map is secondary analysis, so by all means do it ... but present PCA-UMAP first, and base the results around that.

---

From the comments there's some sort of time-series. So if you are performing a time-series - it's not datamining - this would be the justification of using a Diffusion map over PCA. This therefore answers points 1 and 2. However, it would imply you are pooling your data.
> 2 votes

---
Tags: dimensionality-reduction, diffusion-map
---


# Question
Title: Discrepancy in Read Counts Between FastQ and BAM Files in Adapter-Trimmed Pipeline
In a FastQ to BAM pipeline where only adapter trimming is performed, I've noticed a potential discrepancy in read counts between the initial FastQ files and their resulting BAM file. Specifically, I'm seeking clarification on whether the following statement holds true:

"Total number of reads in R1 and R2 FastQ files = mapped and unmapped reads in the corresponding BAM file."

I understand that adapter trimming is the sole preprocessing step in this pipeline. However, I'm uncertain about the impact of factors such as PCR duplicates, mapping efficiency, and any potential loss during processing. The BAM file does contain PCR duplicates, and I'm wondering if these factors contribute to the observed differences in read counts.

Can someone shed light on the expected relationship between the read counts in the FastQ files and the BAM file under these conditions? Any insights into how PCR duplicates may affect the discrepancy would be particularly helpful. Additionally, if there are common reasons for the observed differences in read counts between FastQ and BAM files, I'd appreciate guidance on identifying and addressing them.

**Additional information:**

The pipeline used:

pre-trimming: `fastqc R1.fastq.gz R2.fastq.gz -o pre_trimming`

adapter trimming: `fastp -i R1.fastq.gz -o trimmed_R1.fastq.gz -I R2.fastq.gz -O trimmed_R2.fastq.gz -h file.html -j file.json`

alignment: `bwa mem -M -t 4 -R “@RG\tID:FE_S211\tPL:ILLUMINA\tSM:A4” hg19.fa trimmed_R1.fastq.gz trimmed_R2.fastq.gz | samtools sort -o output.bam -O BAM -@ 4 -`

# Answer
> 2 votes
My main question is whether you refer to the fastq files after adapter trimming or before. It seems you mean the original fastq files. In that case, any read filtering done by fastp will remove low quality or totally-adapter reads.

Note from the fastp docs that length trimming is enabled by default, so trimmed reads will be removed if they are shorter than the length setting, even if they are otherwise high quality. And if they are low quality then they will be filtered out for sure. I refer you to the docs to see what all the various default settings are. Possibly you could turn them all off and still see the initial number of reads, but I very much doubt that they would be useful for you.

UPDATE: see also @Yossi Farjoun's excellent point about alignments vs. reads- not at all the same thing! though you can learn the exact number of reads from a BAM (using e.g. `samtools stats`), the BAM is in fact a record of alignments. So it's important to be looking at the right number when you're doing your comparison!

# Answer
> 2 votes
"Total number of reads in R1 and R2 FastQ files = mapped and unmapped reads in the corresponding BAM file."

there might be confusion regarding the word "read".

"record" refers to an entry in the BAM file (a row in SAM) "template" refers to the original sequence from the fastq/sequencer.

after mapping a single template may be the source of two different *records* due to secondary and supplementary alignments.

So the general answer to your question is: Not if you count "reads" as rows in the SAM file Yes, if you count carefully and make sure to not double-count records that have the same "query name" (the name from the fastq file)

For more information about SAM/BAM/CRAM the spec can provide hours of light reading... :-)

---
Tags: sam, fastq
---


title: tBLASTn vs NCBI TSA

text:
# Question
Title: tBLASTn vs NCBI TSA
Is there a way to blast my proteins against all the TSA (Transcriptome Shotgun Assembly Sequence Database) belonging to a specific taxonomy using command line?

Thanks,
Marco

# Answer
> 0 votes
Just go to Blast here

1. Click where is says Database go to Transcriptome Shotgun Assembly Sequence Database (TSA).
2. Click where is says "Limited by" click organism and thereon specify the group of interest.

Thats it.

---
Tags: proteins, taxonomy, blastn, tblastn, tsa
---