text:
I have really skewed RNA-seq data, what's the best way to normalise it? Preferably in R!
My data frame compares the RNA-seq reads from many genes in different tissues.
The reads look as follows.

I tried using log to make it better but still looks pretty skewed.

Are there any better transformations/functions for this in R?
This is part of an assignment I have, I have used in the past edge software and agree it will do a better job but in this work I'm not meant to use any of the pipelines.

I assume that the eventual goal is differential expression or something along those line, in which case you're much better off using one of the many existing Bioconductor packages.

Welcome to RNA-seq. This is not at all unusual, as RNA-seq raw data are not normally distributed. People typically approximate it via the Negative Binomial distribution. I suggest you check existing packages as Devon suggests in his comment for normalization and differential expression (or anything else). Typical choices are DESeq2, edgeR or limma, but there are others as well.
This guide might help you get started. Please also use google to explore existing threads both here at SE and over at biostars.org.
Most basic questions have already been discussed there and are informative to get started.

I like using DESeq2. There's a great document written by the developers about how to process gene count tables and look at differential expression:
https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The document also mentions other tools that have similar approaches:

Other Bioconductor packages with similar aims are edgeR, limma, DSS, EBSeq, and baySeq.

Another approach, that does not replace the methods already mentions, but could be useful if you don't have raw counts to start with, or if you want to screen hundreds of genes, and you want to focus on the most interesting.

either log transform or zscale your expression values
do a boxplot to see how your gene is expressed in the different tissues
run a non-parametric test to compare the distributions. I usually use a kruskal.test with multiple comparisons and a Mann-Whitney U test (wilcox.test) if I only have two groups.

This should give you an idea if the expression of a gene is different between tissue types.

Thanks for all your answers.
Although the pipeline packages would have been useful I was trying to normalise the data before analysing it later and so I used the bestnormalize package which looked like it did a pretty good job.

What is your goal with this analysis? Your RNA-Seq data looked relatively typical at a glance and I would really recommend not using a general normalization as you did. A lot of research has gone into normalizing RNA-Seq and you are at best hurting your analysis, at worst doing something outright wrong. Use - as recommended before - edgeR, limma or DESeq2 unless you have a good reason not to.

See e.g. Soneson & Delorenzi, 2013

This is part of an assignment I have, I have used in the past edge software and agree it will do a better job but in this work I'm not meant to use any of the pipelines

@Connor you should edit your question and include that information, so that people trying to help don't offer you a solution you can't use.


text:
less than 80% of your list has mapped to your chosen identifier type
I am trying to do GO analysis. However, when I run DAVID, I am getting this error:

You are either not sure which identifier type your list contains, or less than 80% of your list has mapped to your chosen identifier type. Please use the Gene Conversion Tool to determine the identifier type.
In addition, when I am trying to add a column with gene names next to esemble ids in R in the file with differential expression genes, many genes say NA under gene name column.

I have used STAR and Salmon for the alignment,and DESEQ2 for differential expression analysis.
I have used premade indexes for both tools.
Thank you in advance!
Some of the ENSEMBL codes I used:

ENSG00000249200 ENSG00000239920 ENSG00000242295 ENSG00000235763 ENSG00000226933 ENSG00000283632 ENSG00000235838 ENSG00000232429 ENSG00000281076 ENSG00000271042 ENSG00000229119 ENSG00000253384 ENSG00000285645 ENSG00000240627 ENSG00000213590 ENSG00000270739 ENSG00000240201 ENSG00000203531 ENSG00000271626 ENSG00000250316 ENSG00000218713 ENSG00000231548 ENSG00000231049 ENSG00000286185 ENSG00000214535 ENSG00000212961 ENSG00000219747 ENSG00000259100 ENSG00000231767 ENSG00000243175 ENSG00000236804 ENSG00000240634 ENSG00000270773 ENSG00000230267

Give some examples of identifiers that didn't map. Please note that many genes will lack GO annotations, though 80% feels a bit high there.

ENSG00000249200ENSG00000239920
ENSG00000242295
ENSG00000235763
ENSG00000226933
ENSG00000283632
ENSG00000235838
ENSG00000232429
ENSG00000281076
ENSG00000271042
ENSG00000229119
ENSG00000253384
ENSG00000285645
ENSG00000240627 ENSG00000213590
ENSG00000270739
ENSG00000240201
ENSG00000203531
ENSG00000271626
ENSG00000250316
ENSG00000218713
ENSG00000231548
ENSG00000231049
ENSG00000286185
ENSG00000214535
ENSG00000212961
ENSG00000219747
ENSG00000259100
ENSG00000231767
ENSG00000243175
ENSG00000236804
ENSG00000240634
ENSG00000270773
ENSG00000230267

DAVID last received an update in May 2016 (https://david.ncifcrf.gov/content.jsp?file=update.html), so any additions to the Gene Ontology database since then won't be present. I checked the examples you posted in Biomart and found that only 2 have GO annotations, so you're likely not missing an immense number of genes. Having said that, you'd be best served to use a more frequently updated tool than DAVID. Around 2015 there was an outcry because it hadn't been updated in 6 years and the database was so out of date that people were getting questionable results. Have a look at one of the many Bioconductor packages for convenient free solutions that are better maintained.

Thank you so much for your reply!! So, there is not a problem with my index, right? I have tried gprofiler as well, but I wanted to check more tools.

Your index is fine.

Thank you a lot!

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
Template files for bioattributes and meta-data
I am in process of submiting my FASTQ files to the SRA database. I am a little confused about how to fill the biosample attributes and SRA metadata.
I would really appreciate if you anyone could share example templates or help with the given database fields.
Thank you

Hi @Riya this is certainly a good question for the site, but could you explain where you have got to so far? Screenshots would be beneficial.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
gmx mdrun error: "Error in user input: Invalid input values In option s"
I'm new to bioinformatics and gromacs I tried to run the code and got this error
Program:     gmx mdrun, version 2020
Source file: src/gromacs/options/options.cpp (line 179)
Function:    void gmx::internal::OptionSectionImpl::finish()

Error in user input:
Invalid input values
In option s
Required option was not provided, and the default file 'topol' does not
exist or is not accessible.
The following extensions were tried to complete the file name:
  .tpr

I already have the file topol. I searched and found that I need to run mdrun-h to ensure that the command work but when I check
mdrun -h

I got
mdrun: command not found

try 'which mdrun', it might not simply be installed

thanks for replying, I tried it and got nothing

Is the tool installed, is it in PATH?

@ATpoint has a good point! The problem is your environment is wrongly configured

If the command line says "command not found", it means that the binary is not executable or not found in any of the paths mentioned in the $PATH environmental variable. You can check for example this tutorial for details on how it works and how to resolve problems related to this variable.
Which probably means, that you need to go a step back, to the installation instructions of gromacs and try to install the mdrun binary. However, if the mdrun was not installed successfully, it's very likely that more stuff of the package will be missing, so it might be worth reinstalling it all.
P.S. It's not obvious to me if getting mdrun it will resolve the first posted problem (I have never worked with gromacs). If you would share the link to the research you have done, that would help others to help you.

This was a great answer Kamil! We have a few unanswered GROMACS questions here, are you able to help us out at all with those? I'm trying to help take care of the unanswered queue.

I think there was a misunderstanding here. Two errors are in user5520049's post.

A .tpr is not present in the working
directory so gromacs complains, if your .tpr is topol then you have to give this file
to gmx mdrun:  gmx mdrun -s topol.
mdrun is not a command, you should have called instead: gmx mdrun -h. That's why you had command not found and not because gmx is not in your $PATH. You couldn't have run gmx mdrun in the first place otherwise.

Side tip: run gmx mdrun -h -hidden to see hidden options.

thanks for helping . the path is also was a problem but your answer is important for me . i tried to execute those commands gmx grompp -f run_$count.mdp -p topol.top -c em.gro -o run_$count.tpr -nb gpu;
gmx mdrun -deffnm run_$count -nb gpu and got the same error although topol.tpr is exists in the directory

I can assure you that gmx is in your $PATH, you could not launch any gmx commands otherwise (gmx grompp or mdrun etc.).
The problem is that you use -deffnm (which you didn't mention in your orginal post btw) with gmx mdrun which is the option to setup default input and output names. So with -deffnm un_count gromacs is looking for a tpr named: un_count.tpr and not topol.tpr. Either:

use -deffnm un_count and name your tpr un_count.tpr
do not use -deffnm and not use -s and have default topol.tpr in the working directory.
do not use -deffnm & use -s yourcustomtpr.tpr

Here is a link to good gromacs tutorials: http://www.mdtutorials.com/gmx/ , maybe you could take a look to be familiarized with this MD engine ;)

This was a great answer Jeremy! We have a few unanswered GROMACS questions here, are you able to help us out at all with those? I'm trying to help take care of the unanswered queue.


text:
How do I download the mitochondrial haplogroup datasets for human genetics online?
I seem to have landed on the mitomap.org site, but I don't know what to make of it or what do with it / how to get the genomes onto my computer. It sounds like the genomes are stored in GenBank, but that mitomap simply lists which of the GenBank genomes are human mitochondrial haplogroup genomes.
First of all, I can't find a simple list of links or list of names of the haplogroups the mitomap has... Can you show me where I can get or download that? They say:

This brings our total number of FL sequences to 51,673, and the number of CR sequences to 74,660. Our SNVs now total 19,227.

I don't know what FL or CR (Control Region sequences?) means, but I am expecting to see lists of genomes or something. What am I to do with this information? Where do I find this data? Is this it? If I click on one of the results, I see this. It has path nuccore/MT742594.1, can I find this in some NCBI FTP server somewhere? I don't see anything of that folder structure in the GenBank FTP server...
Basically my question is:

Where is the complete list of human mitochondrial haplogroup genomes?
Where can I download all of those human mitochondrial haplogroup genomes so I can get them onto my computer?

There is no such thing as "list of human mitochondrial haplogroup genomes". A haplogroup is not a genome but a list of variants common for many genomes.

Data is located here: https://www.mitomap.org/foswiki/bin/view/MITOMAP/Mitobank. Under the Sequences available column, you can download a text file containing genbank_id. You need to write a script to download all these files. My back of the envelope calculation suggest you will download about 1GB data. I suggest test few entry first.

Just realized this entry is very old :)

Hmmm, I am not experienced enough in biology or bioinformatics to truly understand the question, but I am guessing you want to download the FASTA file or genebank file containing the sequence of the accession number you listed above.
There are 2 ways to do this -

Using your browser, if you want to download a small sample of files. On the right hand side of your screen you can see the "send to" button on this link, click on it and choose the what's appropriate for you and then click on create file.

https://www.ncbi.nlm.nih.gov/nuccore/MT742594.1

If you familiar with Biopython and python programming you can download multiple files at once in whatever format you want, there is an easy script to do this and has been taken from Biopython's documentation on biopython.org
import os
from Bio import SeqIO
from Bio import Entrez
Entrez.email = "A.N.Other@example.com"  # Always tell NCBI who you are
filename = "MT742594.1.gb"
if not os.path.isfile(filename):
 # Downloading...

 net_handle = Entrez.efetch(

     db="nucleotide", id="MG762674", rettype="gb", retmode="text"

 )
 out_handle = open(filename, "w")

 out_handle.write(net_handle.read())

 out_handle.close()

 net_handle.close()

 print("Saved")

print("Parsing...")
record = SeqIO.read(filename, "gb")
print(record)

You can make a list of the accession numbers of the sequences you want to download and pass the list in the id, if you are curious as to how the code works visit biopython.org .

I still cant figure out this code block thingy!!!

The easiest thing to do is just wrap 3 ``` symbols before and after the code.


text:
Get protein names corresponding to PDB ID
I have a list of about 4000 PDB IDs and would like to get the actual names of the proteins (e.g. lactate dehydrogenase, cytochrome c). I tried the batch header section at the Protein Databank Download page but it refused to accept my PDB IDs in formats (xxxx or xxxx.pdb, individually or space-separated) that worked in an interactive search for the protein structure.
Any suggestions?

Have you tried writing a PyMol script?

Perhaps you could use ftp://ftp.wwpdb.org/pub/pdb/derived_data/index/entries.idx which contains HEADER and COMPOUND fields extracted from all PDB entries.

@Michael — Could you clarify please. A script tTo do what on what? And why in PyMol?

@marcin — Thanks, but the problem with the HEADER and COMPOUND fields is that they are such utter unparsable anarchy. I could write a script to get these fields from all my pdb files, but take for example 102L — I would get as HEADER "HYDROLASE(O-GLYCOSYL)" (useless) and as COMPOUND "HOW AMINO-ACID INSERTIONS ARE ALLOWED IN AN ALPHA-HELIX OF T4 LYSOZYME ENTEROBACTERIA PHAGE T4" where the name of the protein — Lysozyme — is buried. The prehistoric format of pdb files is the problem, but surely someone, somewhere has done something with them.

In which case, going from PDB to Uniprot using SIFTS https://www.ebi.ac.uk/pdbe/docs/sifts/quick.html is the best bet as the Uniprot entries will have all the synonyms you could have for.

@MatteoFerla Thanks. Will check that out.

You can use one of the UniProt Protein APIs.
As you said you have your pdb entries in a text file line by line you can, like this example.txt containing:
1brr
4lzm
2dyi

Using the commandline, you can use a little script like this to download the name, if it is available for the given pdb entry.
while read line;
do
 curl -X GET --header 'Accept:application/json' "https://www.ebi.ac.uk/proteins/api/proteins/pdb:$line" |
 jq -r '.[].protein.recommendedName.fullName.value'  |
 sed "s/^/$line\t/" >> pdb_names.txt;
done < example.txt;

You need to have curl, sed and jq installed on your system.
This gives you following output in pdb_names.txt
1brr    Bacteriorhodopsin
4lzm    Endolysin
2dyi    Ribosome maturation factor RimM

Update
if you want to speed it up, you can run it with parallel.
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" |   jq -r ".[]. .protein.recommendedName.fullName.value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

With the -j option you call how many jobs should run in parallel. The limit of the UniProt API is 200 request per second per user.

Update 7. Nov 2020
To get another info beside the protein name, you need to know how the JSON response from UniProt looks like.
To get also the scientific name, you can run following command:
parallel -j 4 'curl -X GET --header "Accept:application/json" "https://www.ebi.ac.uk/proteins/api/proteins/pdb:{}" | jq -r ".[] | .protein.recommendedName.fullName.value + \" - \" + .organism.names[0].value"  |  sed "s/^/{}\t/" >> pdb_names_parallel.txt' :::: example.txt

As result you get this:
1brr    Bacteriorhodopsin - Halobacterium salinarum (strain ATCC 700922 / JCM 11081 / NRC-1)
4lzm    Endolysin - Enterobacteria phage T4
2dyi    Ribosome maturation factor RimM - Thermus thermophilus (strain HB8 / ATCC 27634 / DSM 579)

sed should also be available for you. You can install jq and parallel with brew

Many thanks. Got this running. Only problem is that on OS X the \t tab escape does not work (it is apparently non-standard) and none of the answers here worked for me. However it is easy enough to insert a tab after the four-character PDB ID programmatically.

Check this answer https://stackoverflow.com/a/6392271 for your tab problem.

Thanks again. I know it is bad form to ask supplementary questions, but could you possibly provide a link in your answer to the particular page in the API which describes the need to specify all of the four "dotted" terms. I went to https://www.ebi.ac.uk/proteins/api to try to find out how to construct a query to get the species (so I didn't have to ask) but must be looking in the wrong place. I can find an example XML for a protein under proteins/proteins that includes some of these terms, and various unpromising json thingies. NCBI provides intelligible APIs (or used to) — whereas EBI…

I've updated my answer. The dotted format is not a problem of the API. It is just how jq is getting the keys from the JSON formatted result. https://stedolan.github.io/jq/manual/#Basicfilters gives you a basic overview of how to use jq.

Assuming you can use R, have you tried with biomaRt? For example, using 2bhl (my PhD lover :D)
library(biomaRt)
ensembl <- useMart("ensembl",dataset="hsapiens_gene_ensembl")

# get list of all available info
filters <- listFilters(ensembl)
attributes <- listAttributes(ensembl)

getBM(attributes=c('hgnc_symbol','ensembl_gene_id','entrezgene_id',
                'protein_id','description',"superfamily"),
      filters = 'pdb',
      values = "2bhl",
      mart = ensembl)

Returns
hgnc_symbol ensembl_gene_id entrezgene_id protein_id description
1         G6PD ENSG00000160211          2539   ADO22353 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
2         G6PD ENSG00000160211          2539   CAA27309 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
3         G6PD ENSG00000160211          2539   AAA63175 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
4         G6PD ENSG00000160211          2539   AAA52500 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
5         G6PD ENSG00000160211          2539   AAA52501 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
6         G6PD ENSG00000160211          2539   CAA39089 glucose-6-phosphate dehydrogenase [Source:HGNC Symbol;Acc:HGNC:4057]
7 ...

I am sure if you look at all the available filters and attributes, you can pinpoint the ID you are looking for.

Thanks. I can reproduce that, but only paint by numbers with R. How would I modify that to loop through a list of 4000 proteins and ensure the pdb ID was included in the output?

You can add 'pdb' to the list of attributes. That will add a column with the pdb id. c('pdb','hgnc_symbol',....),

Ok. Get it. Thanks.

I am not sure how your input file is (csv, vector, data.frame), but to run on multiple pdbs, you need to change yourvalues (values= c("2bhl","pdbID2","pdbID3",...))

My input is just a text file of filenames, one per line, generated from a unix ls command. I can convert it into a string of the type you suggest, although I wonder if there is any limit on the length of such strings. Referencing a text file would be preferable if it is possible.

I think the limit is 500 values x query, bu the newest versions of biomaRt should automatically split for you. I am not 100% sure, so you may want to double check. You can ls *.pdb >> all_pdbs.txt , use all_pdb_list <- read.delim("./all_pdbs.txt") to read all_pdbs.txt , and gsub to remove the extension from your entries (gsub(".pdb","", all_pdb_list)) .


text:
DESEQ2:Error in rownames, what is the problem?
I am using DESEQ2 for RNAseq analysis but i do not understand why I get this error.
> dds <- DESeqDataSetFromMatrix(countData = data,
+                               colData = labeldata,
+                               design =~ subtype, tidy = TRUE)

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

I really do not know what is the issue. Could you please help?
Thanks

Hi @athanasia and thansk for your post. Members here do experience of DESEQ2 for transformation (I don't) and whilst the information could be enough, it is etiquette to provide as much detail as you can. For example, your data, size dimension, runs etc. ...

Please edit your question to include the output of dput(head(data)). thanks

The hint is here:
duplicate 'row.names' are not allowed
rownames(data) should return a vector of only unique values, since it doesn't make sense to have the same gene/transcript in two different rows of your count matrix. How did you generate it?
What does which(duplicated(rownames(data))) return? Maybe that'll give you a hint on what went wrong. In general, if you use salmon or kallisto for counting, you should use the excellent tximport package to get the counts into R, then create a DESeq dataset with DESeqDataSetFromTximport().
EDITS:
The replies to my answer prompted me to re-read your post and I saw that you had set tidy=TRUE, which means the duplicates are actually in your first column of your data.

I thought along the same lines, but data cannot have duplicate row.names - no object can have duplicate row.names, so this is some generated object that is being assigned duplicate row names owing to something wrong with data or labelData that has nothing to do with the row names of either object.

Of course, they can. Here's a oneliner for you: m <- matrix(0, nrow=3, ncol=3, dimnames=list(c('a', 'a', 'a'), c('b', 'b', 'b'))). Anyway, I tried and DESeq2 handles this case semi-gracefully by appending numbers to the dupes. However, @Athanasia had set tidy=TRUE, which prompts DESeq2 to look for row names in the first column of the input, where duplicates trigger an error. I can post an MRE for this exact behavior if you'd like.

Thank you, TIL.


text:
Meta-analysis in R
I am completely new to meta-analysis and am a bit lost on how to proceed.
I have two biomarker datasets; each one is using a different microarray platform. I have the protein names, the concentration in each subject (one dataset has 15 subjects while another has 22) as well as fold change, p-value, and AUC values.
Is there a simple way to do a meta-analysis with this information? For example, is there any way to join the datasets and calculate a combined fold change and AUC? I was thinking of using sumlog() from the metap package in R to do a combined p-value, but I am lost on how to combine the other two values, or if this would even be appropriate. I've attached a screenshot of one of the datasets

Any help or guidance would be greatly appreciated!

You seem new to analyze microarrays. Look at the wonderful limma vignette about how to analyze one dataset. Then I suggest you to analyze both  datasets and see how much do they agree. However what you need to do depends on what you want to do (besides a meta-analysis). What is your biological question (or reason of this analysis)?

Thats a lot of questions rolled into one. First consider data mining to get an overview


text:
Need help with Pairwise Alignment module in iterating over the alignment
from Bio import SeqIO
from Bio import Align

ref_seq_1 = SeqIO.read('C:/Users/King_Arthur/Desktop/ref_seq/segment 1/ref_seq_8.fasta','fasta')

seq1 = SeqIO.read('C:/Users/King_Arthur/Desktop/file/segment 1/Myfile_1 (1).fasta','fasta')

aligner = Align.PairwiseAligner()
aligner.mode = 'global'
aligner.match_score = 1
aligner.mismatch_score = -2
aligner.gap_score = -2

alignments = aligner.score(ref_seq_1.seq , seq1.seq)

print(alignments)
for alignment in sorted(alignments):
    print(alignment)

So this is my code and as you can see in the last section i am trying to iterate over my alignment but I am getting this error
TypeError: 'float' object is not iterable

I have tried various things like using str() but it gives some strange values and I also tried to read the source code by using the inspect module but I can't figure out the problem.
Any help would be really appreciated.
My final objective is to find out how many matches, mismatches and gaps are present in the final alignment using biopython.
if there is any other better way to do it in python please feel free to suggest.

aligner.score(ref_seq_1.seq , seq1.seq) doesn't produce an alignment. It just gives you the score of the alignment as float value.
Instead you need to run:
alignments = aligner.align(ref_seq_1.seq , seq1.seq)
Now alignments is an Bio.Align.PairwiseAlignments you can iterate over alignments or get the score with alignments.score which gives you the same as aligner.score(ref_seq_1.seq , seq1.seq)`
If you need more information, you could check the documentation.

Thank you so much.
That helped a lot.
Do you know anything to solve this- My final objective is to find out how many matches, mismatches and gaps are present in the final alignment using biopython.

For your final objective check http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec100


text:
Are sequencing reads written differently when inserts and indexes are sequenced seperately?
I was going over different sequencing library protocols and noticed that the adapter sequence can vary between protocols. In some adapter sequences the index sits between the primer binding site and the insert. In this case during sequencing the index will be sequenced and then the insert that comes right after it. In other cases, the adapter sequence is built so that the index is outside the primer binding site, i.e. index - primer binding site - insert (see two leftmost parts of the attached figure). Here the insert is sequenced using the Read1 primer and afterwards the index is sequenced in a separate reaction using the 17 index primer.

In this case are the index reads just added on to the ends of the insert reads in the FastQ files? Is there a way to tell from the FastQ files where the index was positioned within the adapter?

Software FastQC reports overrepresented kmers per position in read, perhaps that could help you figure out the order? Also, note that some sequencing providers are doing some of the trimming before they hand users the data.

I do think this is a bioinformatics question.  Ideally you could ignore this as an implementation detail that doesn't concern bioinformatics, but in reality I've had to stare at diagrams exactly like this when trying to troubleshoot strange illumina runs and interpret changes in base quality, for example.  I'd like to see a clear answer to this myself.

The order of sequencing in a dual index flow cell is:

Read 1
Index 1 (I7)
Index 2 (I5)
Read 2

Between each step there is a switch of the primer used (Read1-primer -> Index1-primer (I7) -> Index2-primer (I5) -> Read2-primer)
There are differences in how the I5 is dealt with between dual- and single-read flow cells as explained in the Illumina Indexed Sequencing Overview. In this document there is a lot more detailed explanation of the sequencing workflow that may be helpful to you.
But as explained by the other answers, none of this will matter in most cases. Even if you decide to output the Index sequences they will be exported as separate FASTQ files.

This is definitely a bioinformatics question! Understanding how sequencing libraries are created in the wet-lab and how they are sequenced is fundamental to bioinformatics.
To answer the question, the index sequences are not included in the fastq file output from Illumina, or other machines. Illumina sequencers just use the index sequence to demultiplex the complex mixture of molecules in the sequencing pool. If a small "dot" on the glass sequencing slide has some sequence, but the indices don't match what the technician specified for a library, then that read is not included in the fastq output of that library since it could be contamination of some sort.
The sequences in the fasta file theoretically all should be the insert from the library that you prepared, but sometimes the molecule being sequenced is short, and there is read-through of the adapter sequence. This is the only circumstance under which the i7 or i5 index sequences might be in the normal fastq output sequence line. This is not ideal, and this is why people perform adapter trimming.

The index sequences can be returned as separate fastqs, but most of the time, you don't need them.  And cluster whose index doesn't match what's on the sample sheet does get outputted, in the "Undetermined" fastq.


text:
Nonbonded interactions on the GPU
i'm using this command
gmx mdrun -deffnm nvt0 -nb gpu

but got
Nonbonded interactions on the GPU were requested with -nb gpu, but the GROMACS
binary has been built without GPU support. Either run without selecting GPU
options or recompile GROMACS with GPU support enabled

I've Nvidia rtx2060 with cuda 11.1 on ubnutu 18.4.5

I don't know gromacs but that would be normal for GPU in ML, CPU and GPU uses two different compilations.

Someone (with the same GPU) asked the same question during the weekend that they have now deleted, where I said about CUDA and cmake with -DGMX_GPU=ON. If this was you, please do not do this as you are wasting volunteers' time —and I believe may get you blocked (maybe?). A lesser evil is to edit the Q to bump it.

I've put a post on the Meta forum https://bioinformatics.meta.stackexchange.com/questions/251/duplicate-ops-and-op-deletion-of-response-postings

Try here: https://mattermodeling.stackexchange.com/. There's an entire GROMACS tag there with way more questions.

As the error tells you: the GROMACS binary has been built without GPU support.
Did you install gromacs with apt-get? Usually these kind of programs are shipped with a basic built, with a lot of features disabled.
If you did built gromacs from source, then cuda was not found during configuration.
I would try building gromacs from source and:

be sure cuda is in your path and check that nvcc and the drivers have a matching version
run cmake and check the output, to see if the video card is recognised
run cmake with -DGMX_GPU=ON

Note that it seems that gromacs doesn't support multiple gpu accelerations( e.g. it is either openCL or cuda, you can't have both)

thanks for helping . i did it already used this command -DGMX_GPU=ON but still got error . i'm now removing gromacs and will install it from scratch taking your answer in my consideration and will get back .. thanks again

@user5520049 When you have time, please note whether you were able to finish compiler with -DGMX_GPU=ON, noting that you would need to still copy the resulting recompiled binary to somewhere in your $PATH.

This was a great answer fra! We have a few unanswered GROMACS questions here, are you able to help us out at all with those? I'm trying to help take care of the unanswered queue.


text:
Excess average estimated identical-by-descent in genotype data
(Cross-post with Biostars: https://www.biostars.org/p/470788/)
We have some genotype data that we are putting through quality control in PLINK 1.9. As part of this QC, we have limited the data to subjects of self-reported and genomic European ancestries, and have run basic QC steps, removing call rates below 95% for individuals and variants, and removing variants excessively out of Hardy-Weinberg equilibrium (< 10^-10).
Pruning for LD and limiting to the autosomes, we have calculated estimated IBD between our participants and the estimated the average pairwise IBD for each participant. We have found that an unusual number of our participants (233/13049) show a high mean pairwise IBD with others (>3SD from the cohort mean, > 0.028 in this instance, maximum 0.075). Examining examples of the individuals in question shows that this is driven by low-level relatedness (~0.03-0.1) with many other individuals, rather than a small subset of relatives in the data. In addition, the relatedness in each case stems entirely from IBD1 (one shared variant) and not at all from IBD2 (two shared variants). Individual average IBD correlates negatively (-0.69) with the genome-wide inbreeding coefficient.
I have one idea about what could be causing this (partial cross-contamination of the samples), but I am interested to know if others have other ideas what would be the cause?

What population are you conducting this in? Some populations e.g. Iceland have high levels of IBD due to small effective population size. Also, what length (in cM) are these 233/13049 individuals sharing?

UK population (although no guarantee of British ancestry, so country of birth may be an interesting thing to look at). Shared segments is also a good point - haven't looked at that but will.

A clear negative association within an inbreeding coefficient suggests convergent evolution driven by shared selective pressure. There is an outside chance of purifying selection, or else balancing selection, but selection is the key word. The phenotype of the locus needs investigating, i.e. outside nearly neutral theory what is constitutes positive selection?
Allopatry caused by genetic isolation and genetic drift is your null hypothesis. However Darwinian positive selection is the alternative hypothesis.
I would be wary against suggesting contamination unless you've a clear suspicion this is occuring. I do accept that in human genetics its often difficult to avoid, particularly in ancient DNA type work. In this case you would need to try and control for it. The argument I would put against contamination is that 'genetic relatedness' is localised within the genome, whereas contamination would be a genome-wide phenomena.
Convergence is common within for example HLA genotypes (T-cell) loci between global populations.
You could perform formal tests for postive selection either within pop gen (the locus will associate with significant homozygosity) or phylogeny (point mutation).


text:
How to calculate GenePy score for vcf files in all subdirectories?
In a parent directory, I have several sub-directories, each of them contain one .vcf.gz file on which analysis should be performed. I want to run following commands on all .vcf.gz files in all subdirectories. What is the best way of doing it?? I am the beginner with shell scripting.
I will add short explanations. The input is .vcf file from every subdirectory. at the end, GenePy score is calculated.
To improve score consistency, filtering step:
vcftools --gzvcf GENOTYPED_ALL.vcf.gz --minGQ 20 --recode --out GENO_FILT_GQ20 # Set GQ<20 as missing ( retain high confidence calls)
vcftools --vcf GENO_FILT_GQ20.vcf --max-missing 0.7 --out filtered # Remove SNVs with missing rate >30%
vcftools --vcf filtered.recode.vcf --min-alleles 2 --max-alleles 2 # Keep only Biallelic SNVs

Make annovar-ready file:
./annovar/convert2annovar.pl \
    -format vcf4 FINAL_GQ20M30_BIALL.recode.vcf.gz \
    -outfile ALL_genepy.input \
    -allsample \
    -withfreq \
    -include 2>annovar.log

Annotate against refGene,gnomad_exome,cadd13,eigen,revel,gwava,dann using Annovar
./annovar/table_annovar.pl \
        ALL_genepy.input \
        ./annovar/humandb/ \
        -buildver hg19 \
        -out ALL_genepy \
        -remove \
        -protocol refGene,gnomad_exome,cadd13,eigen,revel,gwava,dann \
        -operation g,f,f,f,f,f,f \
        --thread 40 \
        --maxgenethread 40 \
        -nastring . >>annovar.log

Combine Genotypes and annotations
cut -f 18- ALL_genepy.input > a1
zgrep '^#CHR' FINAL_GQ20M30_BIALL.recode.vcf.gz | cut -f 10- > b1
cat b1 a1 > ALL_temp
paste ALL_genepy.hg19_multianno.txt ALL_temp > ALL_genepy.meta
rm a1 b1 ALL_temp

Prepare output folders/dependecies, making new folders in a current directory to store raw GenePy score files
mkdir CADD13_RawScore Eigen GWAVA_region_score GWAVA_tss_score dann REVEL

Take the header from the ALL_genepy.meta file and stores it in a newly created header file
grep "^Chr" ALL_genepy.meta> header

Compute GenePy scores
Once the ALL_genepy.meta file is created, GenePy_1.2.sh can be run by simply iterating through the list of deisred genes. The make_scores_mat_5.py file is in the same directory of GenePy_1.2.sh (every subdirectory)
chmod +x GENEPY_1.2.sh
while read gene; do
  sh GENEPY_1.2.sh $gene
done < gene.list

This is code that I have run for .vcf files when I had only one .vcf file. Now, I have over 200 subdirectories, in each there is one .vcf file for analysis.
How to wrap up this code to be executable for all subdirectories at once?
I tried with the folowing code. But the result file is not saved in the corresponding subdirectory. It is saved in main directory, and every time it got overwritten with next one.
    find . -type f -name '*.vcf'      \
   -exec vcftools --gzvcf {} --minGQ 20 --recode --out GENO_FILT_GQ20 ';'   \
   -exec vcftools --vcf GENO_FILT_GQ20.recode.vcf --max-missing 0.7 --out filtered  ';'\
   -exec vcftools --vcf GENO_FILT_GQ20.recode.vcf --min-alleles 2 --max-alleles 2  ';'\
   -exec /home/bioinformatics/Dokumente/pipeline_test/annovar/convert2annovar.pl \
    -format vcf4 GENO_FILT_GQ20.recode.vcf \
    -outfile ALL_genepy.input \
    -allsample \
    -withfreq \
    -include 2>annovar.log ';'\
   -exec /home/bioinformatics/Dokumente/pipeline_test/annovar/table_annovar.pl \
        ALL_genepy.input \
        /home/bioinformatics/Dokumente/pipeline_test/annovar/humandb/ \
        -buildver hg19 \
        -out ALL_genepy \
        -remove \
        -protocol refGene,gnomad_exome,cadd13,eigen,revel,gwava,dann \
        -operation g,f,f,f,f,f,f \
        --thread 40 \
        --maxgenethread 40 \
        -nastring . >>annovar.log';'\
    -exec cut -f 18- ALL_genepy.input > a1';'\
    -exec zgrep '^#CHR' GENO_FILT_GQ20.recode.vcf | cut -f 10- > b1';'\
    -exec cat b1 a1 > ALL_temp';'\
    -exec paste ALL_genepy.hg19_multianno.txt ALL_temp > ALL_genepy.meta';'\
    -exec rm a1 b1 ALL_temp #remove';'\
    -exec mkdir CADD13_RawScore Eigen GWAVA_region_score GWAVA_tss_score dann REVEL ';'\
    -exec grep "^Chr" ALL_genepy.meta> header ';'\
    -exec chmod +x GENEPY_1.2.sh ';'\
    -exec while read gene; do
  sh GENEPY_1.2.sh $gene
done < gene.list ';'\
done

Thank you for all suggestion, I appreciate it.

btw, I would reccomend not using vcftools and switch to using bcftools, since I believe that vcftools isn't under active development any more and may therefore contain bugs of some kind.


text:
Interpreting contig alignments to a reference genome
I have applied two de novo genome assembly tools to data from the run SRR12707453, corresponding to a phage (I downloaded the data and have no relation with the authors of the study).
Using rnaSPAdes I get many long contigs (>50.000 bp). When analyzing the contigs using Quast against a reference genome (MN813694.1) most long contigs are discarded and the seventh-longest one is the only one selected. If I blast the other contigs I find correspondence with other phages, which makes me conclude that many phages were present in the same probe.
This seventh contig has 47686 bp length, while the reference has 47713 bp length. But they are not equivalently arranged. Apparently, the genome was circularized and then it was broken around the position 32400 in the reference. Also, the sequence does not match between positions 45000 and 45951 (in the reference), and a little misassembly seems to occur around position 46300.
See here some screen captures for the Icarus visualization obtained with Quast:

In contrast, using Quast with the contigs obtained from the other assembly tool, I get some differences. There is also a contig with almost the same size as the reference genome (47573 bp) that has exactly the same structure than the contig obtained with rnaSPAdes: circularized, broken around position 32400, does not match between position 45000 and 45951, and small misassemblie aronund position 46300. But there are also some contigs in which there is no break around position 32400.

May I interpret that these contigs indicate that the phage genome is sometimes circularized and sometimes broken at position ~32400 and that's why some contigs don't find a sequence break at that position?
Am I interpreting anything wrong? Any other conclusion can be obtained from the data?
You can (temporaly) download here the Quast results.

some additional information to your request.
https://ftp.ncbi.nlm.nih.gov/genomes/Viruses/FamilyPhylogeneticTree/Graph/SRR12707453.MN813694.1.gfa
it contain assembled by reference graph
using ncbi tool saute(skesa) which will be public later.
install bandage app and load this graph
https://rrwick.github.io/Bandage/
picture looks intriguing,I would say , run contain many errors located in one position.
Reason? I don’t really know.

Yep... I don't know how circularization works. I guess it could be that binding the sequence generates every time some garbage in the union and so every RNA chain has a slightly different sequence at this point.


text:
Find CRISPR PAM-sites with python
I am a college student and just started bioinformatics. I am trying to write a script (not for school) that finds all potential NGG and TTTN protospacer adjacent motif sequences in a genome string. I wrote the following script:
TargetFile = open('file.fasta')
readTF = TargetFile.read()
useTF = "".join(readTF[66:].split())

#no_of_basepairs = len(useTF)

sites = {"ngg_strand" : [],
         "ngg_antistrand" : [],     #location counted from reverse direction
         "tttn_strand" : [],
         "tttn_antistrand" : []}    #location counted from reverse direction

for i in range(0, len(useTF)-2):
    gg_mer = useTF[i:i + 2]
    ttt_mer = useTF[i:i + 3]

    if gg_mer == "GG":
        sites["ngg_strand"].append(i)
    elif gg_mer == "CC":
        #look for reverse complement, record location on antistrand from reverse direction
        sites["ngg_antistrand"].append(len(useTF)-i-2)
    elif ttt_mer == "TTT":
        sites["tttn_strand"].append(i+1)
    elif ttt_mer == "AAA":
        # look for reverse complement, record location on antistrand from reverse direction
        sites["tttn_antistrand"].append(len(useTF)-i-2)

print(len(sites["ngg_strand"]))
print(useTF.count("GG"))

At the bottom I used the .count() method to compare results. I noticed that the method returns a smaller number than my code. I believe this is due to the .count() method not counting overlap, whereas my code does. I.e. in the string "NNNGGGGNNN" my code would count 3 occurences of the "GG" substring, whereas the .count() method returns 2 occurences.
My Questions:

If I want to find all possible PAM sites, which way of counting is correct? I believe my way would be correct, since as far as I know the RNA molecule doesn't advance in "triplet-steps". However, I am very much new to this and might be missing something.

If the reverse complement is found, is the location usually recorded from the same direction as on the strand, or in the reverse direction (as I implemented in the code)?

Like I said, I am new to this so if you see any other glaring problems I would appreciate any and all help. :-)

Thank you!

Hi @heibai thanks for your question. If you have  a moment could you please clarify  "NGG and TTTN protospacer"? Its just question 2 is difficult to understand and I suspect relates to the definitinon of  NGG and TTTN protospacer motifs. Do you mean in question 2 that AAA and TTT  are reverse complements?

Hi @Michael! So I am trying to find all possible editing sites for Cas9 and Cas12, which in this case have PAM sequences of NGG and TTTN. As far as I understand the protein complex reads the strands in opposite directions, meaning that if the .fasta file contains a substring "NAAA", the antistrand at that position would be "NTTT". Read in the reverse direction it becomes "TTTN", i.e. a potential editing site. Since the "N" can be any nucleotide, I am not including it in the search. I hope this clarifies things and doesn't just further confuse what I mean. :-D

@Michael If you have time could you expand on that? What kind of situations would result in false hits? What do you mean by defining forward and reverse direction, and how would I do that?

If this is not an assignment but for productive use then use specialized software such as ChopChop of CCTop (or many others) as you also have to check for genome-wide similarity regions (off-targets) plus there are many factors that can influence efficacy of the gRNA to actually recruit the Cas protein to the target site which is implemented in these tools. It is not as simple as a mattern match. For example: https://cctop.cos.uni-heidelberg.de:8043/ It has a stand-alone version.

Okay @Michael, thank you for your help.

@ATpoint it is an informal assignment given to me by one of the postdocs at the department. I'll have a look at the specialized software out of curiosity anyway, thank you!

The easiest solution is finditer
use re

# Compile you search term where [A|T|G|C] means A or G or C or T
p = re.compile("[A|T|G|C]GG")

iterator = p.finditer(useTF)
for match in iterator:
    print(match.span())

** Example Output**
(0, 2)
(22, 24)
(29, 31)

Easy does it. Its worth looking at the documentation for re here, https://docs.python.org/3/howto/regex.html for alternative approaches using match or search below. finditer is definitely the solution in context.

I support your coding approach because it provides lots of flexibility. For example, you might want to target GG or TTT around a given locus and you can easily modify the code to accomodate this. Someone elses code would be less flexible.
In summary, I kinda like your substrings and do like the count method. I'm not a fan of your == in this instance and your C-style matching should be substituted for reg-ex which in Python is implemented via re.
The basic issues are,

Biopython would be used via AlignIO or SeqIO. This will allow any given format to be analysed. Yours approach is okay but prone to bugs (below). [probably a minor issue]
Your sequence is in some sort of weird format, like extended phylip format and it is not fasta (that wouldn't work). It needs verifying because removing the first 65 characters might remove bits of sequence [probably a minor issue]
Chopping of -2 from the length ... I dunno I'd check out -1 [minor]
C-style matching. I don't like this. I would definitely use re. One solution uses match or search, but the solution above uses finditer, Reg-ex will allow to search for NGG and TTTN, which you are not doing
I think you are looping a locus of 3 and 4 bp and I'm not sure how the C-style == would match this because thats 2 and 3 bp respectively. Just thoughts. You will have 'edge effects' at the beginning and end of sequences, e.g. TTT 3' will match but that isn't TTTN.
I like your count method, but I think finditer is better because it gives you location information
The dictionary is sort of fine. Personally, I would dump the dictionary into a pandas dataframe This way the results are output as e.g. CSV, spreadsheet flat-files for your colleagues to view.
the if and elif would be placed in a subroutine, this is good practice to convert to an object orientated approach.

Great, thank you @Michael! This already gives me lots of stuff to research.

Upvotes and 'accept' preferred.


text:
How to filter allele frequency VEP with bcftools?
I have a VEP annotated VCF file, and I wish to use bcftools to filter it. bcftools work very well on standard VCF fields, but I'm not able to make it work on VEP fields.
For example, the following command would report "not defined in the vcf header" (it's not a standard VCF field but part of the VEP fields in the VCF file). split-vep is a VEP plugin offered by bcftools.
bcftools +split-vep MyFile.vcf -i 'gnomAD_exomes_AF<0.01'

Thanks for the edit! Now, please also give us a file we can test on. Try to make it as minimal as possible while reproducing the error. Note that if the field isn't actually included in the headers, then the VCF is not valid, and bcftools won't know how to deal with it. You might, however, be able to script something, but for that we will need to see at least some part of the file.


text:
Best way to align to find inserted sequence
We have some RNA from knock-in mice, there are two different sequences we're looking for. We have aligned to the mouse genome using STAR but the sequence isn't there which isn't too surprising What is the best way to check for the presence of the sequence? Should I just turn them into a sequence and use them as the reference genome for STAR to align to? Is there a way I can easily incorporate it into the existing genome so I could get read counts for everyhting including those reads? Should I just BLAST using something like magicblast?

Please search biostars.org, this has been asked before and discussed there. I think the naive approach is to include the sequences as extra chromosomes, then extract reads mapping to it and then check the soft-clipped part of the reads (that actually come from the genomic neighborhood) and see whether you can track it back to get a proxy for the insertion site.

Depends on the throughput. High throughput download the mouse genome and use Blat (a bit out of date now)? Megablast just make sure you select the mouse genome to speed the blasting.

@ATpoint could you give me some suggestions on what to search by any change? I couldn't really find anything useful partly because it kept talking about Indels.

@Michael I'll look into that. Thanks.

@Sethzard
https://www.biostars.org/p/374583/


text:
How to find all WGS assemblies accessions of a species
Some background
Similar to the OP of https://www.biostars.org/p/377840/, I would like to programmatically BLAST a sequence to a local database of all WGS assemblies.
Since this isn't feasible for the average biology lab server (correct me if I am wrong), I plan to use ncbi-acc-download to download all WGS assemblies of the species of interest (not a popular species like E. coli, so it should be feasible). Then, I will create a BLAST database for the downloaded assemblies and BLAST the sequence to it.
My question
How can I find all WGS assemblies accessions of a species?

My current plan is to search the NCBI Assembly database using Entrez and a search term such as "wgs"[Properties] AND txid1337[orgn:exp].
EDIT: IIUC, this approach might miss some WGS assemblies. See my answer.
I am worried (and thus ask for your help) this isn't the right approach because there seem to be at least 3 other places in which assemblies can be found:

The assembly_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt.
https://www.ncbi.nlm.nih.gov/Traces/wgs/. This is the one I am most worried about. Both because it seems to have much more entries, and also because https://www.ncbi.nlm.nih.gov/genbank/wgs/ says: "The public WGS projects are at the list of WGS projects ."
taxid2wgs.pl. (As mentioned in https://www.biostars.org/p/377840/#471535, this one seems to miss some assemblies.)

There seem to be WGS assemblies that can't be found in NCBI Assembly database, e.g.: https://www.ncbi.nlm.nih.gov/nuccore/1779902990. I guess that such assemblies also cannot be found in The assembly_summary.txt files that are described in https://ftp.ncbi.nlm.nih.gov/genomes/README_assembly_summary.txt.
My current best guess is that each WGS assembly has a "WGS master record" in NCBI Nuccore database. To find all WGS assemblies of a taxon whose uid in NCBI Taxonomy database is 1337, search NCBI Nuccore database using Entrez and the search term "wgs master"[Properties] AND txid1337[orgn:exp].


text:
Get raw genome sequence
I'm currently working on testing some matching algorithms for strings. I would like to do some tests on raw genomic data as I expect different results from random strings given the lower entropy. How can I get a whole human genome sequence? (just the raw basis)
I have tried websites like https://www.ncbi.nlm.nih.gov/ it seems very confusing and full of unnecessary information given that I just want to have the raw basis without annotations.

Its potentially an interesting project. I'm aware of comparable projects, but ... your question and investigation is ill-defined. Simply if you download 'fasta' format for any given sequence you will get the description on one line beginning with '>' and the sequence on the line below. Genbank is a good place to go

@Michael I'd be happy to have references to the comparable projects, thanks!

Just a small comment, "raw" usually in genomics means unprocessed sequencing reads: https://en.wikipedia.org/wiki/Read_(biology), the sequence you seem to be after is "genome assembly".

You want FASTA files from here.
If you want the entirety of the genome, get the dna-toplevel file. Files with rm and sm have repeat sequences masked, the former by replacing the repeats with Ns (hard mask) the latter by designating repeats as lower case letters (soft masked).
See the README for more info.


text:
Gene expression Table to Expression Matrix converstion
I have an RNAseq gene expression file (Count data) as follows, I need to conduct a Differential gene expression analysis between Patients, for that, I need to have this file as a Matrix of Rows as Genes and columns as counts (samples), how best can I get this as matrix file?
Thanks a lot for your help, Stay safe

Cross posted here:
https://www.biostars.org/p/471641/

Your question stated differently is as follows:

I have data in long format and need to change it to wide format. For that, I want to keep the gene_raw, use patient_id for the column names and gene_code for the row names. How can I do this easily?

Assuming you know how to read that into R as a data.frame that I'll call d:
library(tidyr)
library(dplyr)
reformatted = d %>% select(patient_id, gene_code, gene_raw) %>% spread(patient_id, gene_raw) %>% column_to_rownames(gene_code)

Or something along those lines.


text:
BinaryCIF vs MMTF formats, which one to choose?
Both file formats are advertised as more memory/parsing efficent than old PDB format and its successor mmCIF format, but I can't understand the context (I'm a developer, but pretty new to the field of bioinformatics).
Is BinaryCIF a later iteration of MMTF concepts/ideas (so it's a new competing format) or is the purpose of the two formats different ?
Refs :

BinaryCIF: https://github.com/dsehnal/BinaryCIF
MMTF: http://mmtf.rcsb.org/index.html

So binaryCIF is newer and is heavily inspired by MMTF.
NGL, the most "loved" JS library (IMO) to show protein fetches proteins as MMTF, while Mol* uses bCIF —Mol* has the official repo of the bCIF format. Alex Rose, developed the former at the RCSB PDB, and the latter in collaboration with the PDBe. He also worked with Antony Bradley, the main author on the MMTF format, who left the RCSB PDB, some 4+ years ago. So there is definitely no competition there.

The MMTF download of the RCSB works well: https://mmtf.rcsb.org/download.html
The MMTF format did not catch on —okay, mmCIF is only slowly replacing PDB, which is still the only accepted format for analyses—, hence the need for a bCIF format
The libraries to do the conversion aren't without glitches
In the distant future bCIF will be the sole binary format
mmCIF <--> bCIF conversion is easier

However, I really ought to stress that mmCIF is a better choice to bCIF or mmTF. And even then many tools still hate mmCIFs and use PDB format.
If memory footprint is not really a problem, then stick with mmCIF. If you are running a webserver, you may be sending data as gzipped making this a moot point.
Last year for my site, michelanglo, I was faced with a dilemma: Mol* has dotted lines for gaps and does not crash the browser with a virus or  someone wanting to show hyperball representation on all residues. However, it was not documented (and then glitchy) so I firmly stuck with NGL, which is excellent. For file formats, I stick with PDB and a metadata dictionary as I know the format well and is despite having the jankiest syntax and no real metadata support (REMARK... ha!), does not break isopetides, disulfides or for no apparent reason.

Good post. As stated, pdb (from the age of Fortran) remains a strong favourite

Thanks a lot for this useful and detailed response : a lot of context, exactly what I needed.
We'll use MMTF essentially to send files to NGL (as you have guessed ;) ), internally our tools will work with PDB and mmCIF.
I'll keep an eye on bCIF, as Mol* has impressive tools (but not needed at our current stage).
Have a nice day and thanks again !

Consider mmJSON.
Each of the three PDB sites had the same problem: using mmCIF files in web-based viewers is inefficient. To solve this problem three new file formats were introduced:

MMTF at RCSB,
BinaryCIF at PDBe,
mmJSON at PDBj.

MMTF has the smallest files – it's aggressively optimized for size.
It is based on MessagePack, but on top of it uses a number of custom compression schemes, which make the format rather complex.
Usually it's not a problem – MMTF developers made an effort to provide libraries for reading MMTF files in the most popular programming languages. But you need to install these libraries.
To make the file smaller, MMTF doesn't contain all the metadata from mmCIF, which may be a blocker in some cases.
What is nice is that the MMTF developers reached out to other developers in this field, organized some gathering and they were taking comments and suggestions from the outside.
As the result, a few popular programs such as PyMOL and JMol can also read MMTF.
BinaryCIF is clearly based on (or at least inspired by) MMTF.
So it's also built on top of MessagePack with a number of custom compression schemes.
It includes all the mmCIF metadata, which makes the files slightly bigger, but I think it's a good trade-off.
Initially there was only TypeScript code for reading it, now I see that there are also implementations in Java in Python.
mmJSON is 1:1 equivalent of mmCIF, but uses JSON instead of CIF.
There is no extra compression added, but having the data from mmCIF re-ordered column-wise makes gzip compression much more efficient (this was also noted by MMTF developers).
If the file is stored as json.gz or downloaded in web browser (with the usual compression between the server and the client) its size is 30-40% of gzipped mmCIF. It's not as small as MMTF and BinaryCIF (which are several times smaller than mmCIF) and this is the only downside.
The advantage is that you can read in it any language like all other JSON files.
No need to install/update MessagePack + MMTF/BinaryCIF libraries.
A few years ago I compared reading MMTF and mmJSON in Python – reading mmJSON was order of magnitude faster. It may be very different with libraries in other languages, but generally JSON parsers are fast.
The fastest ones parse GBs per second. Reading binary formats can be faster, but decompression schemes can also make it slower.
Anyway, all the three formats are fast enough to read.
Sadly, I don't think that mmJSON is as popular as it should be.
Perhaps I'm the only person using it outside of PDBj.
That's why I'm writing this answer.
The mmJSON files are available at:
https://pdbj.org/help/mmjson

Thanks ! Really useful added context and experience.
Your POV on mmJSON is really interesting if one want to read/parse it in a JS or Python based stack.
Do you know if pdbj.org is in good sync with rcsb.org (on which we sync regularly)?
Is the transformation mmCIF => mmJSON reliable (Matteo said above that some conversions are glitchy)?

@ashas all the wwPDB sites are updated on Wednesdays. I suppose the conversion is reliable, I don't know about any problems.


text:
Pymol: color spheres with different colors
Hello I am a Pymol user
I want to color all alpha carbons (ca) and beta carbons (cb) from a certain ressidue in different color. How can I do this?
I selected the spheres with:
show spheres, C/1-9/ca
show spheres, C/1-9/cb

and if I use
set sphere_color, theColor

I color all the spheres what I don't want.
according to: color    gold,    resi    200-‐221
I tried
gold,    spheres, C/1-9/ca
but this didn't work

Hi @Elias whilst your question is certainly answerable it is etiquette to provide a few more details about your investigation. Past posts do answer this specific question

Spheres
So how PyMOL stores how it shows atoms is that each atom has a colour and a represention, which is a binary number telling what to show. Therefore the colour for spheres is the same as the colour for sticks. The only way round this is to copy the object (create copycat, original and you set copycat differently)
To set the atom you do
select iliketheschaps, resi 200‐221 and name CA+CB
color tv_green, iliketheschaps
show spheres, iliketheschaps

The object iliketheschaps is a selection. Do make sure the range is a hyphen-minus (character next to 0) as opposed to an en-dash (–) which Word subs and on a mac you get with cmd+-
non-PyMOL colours
Now you asked for "gold". This is not a colour in the menu —famous for the colour limon. Gold is a CSS colour name, but it won't work.
You will have to colour by hex code. Say you want "gold" you need to google the colour name and hex. #FFD700.
color 0xFFD700, iliketheschaps

The Fashion Police have been alerted

This paragraph is in jest.

However, as you may know bioinformaticians are very fashion savvy for colours. Mustard was the in-colour from last year, so you really ought to use mustard (0xffdb58). This year is all about nectarine (0xffa168) (which is not salmon —this ain't 2015!). Coral (0xff7f50) and teal (0x008080) were from previous years and are still are going strong!

show sphere, A/1/C
set sphere_color, red, A/1/C

Please add some explanation to your answer, and format code lines better. As it stands, this post is low quality and might be deleted.


text:
mRNA cap identification by sequencing?
If you sequence a positive-sense  RNA virus genome using NGS via Nanopore, Illumina or even the old Sanger method, is it possible to tell if the virus has 5' (five prime) cap at the 5' end of the RNA genome?

Hi @LRJoshi, please do consider gringer's answer below. If this addresses your question please do "accept" it.

It's possible to do 5' cap enrichment sequencing (e.g. see here), and presumably you could look at differences between an enriched vs non-enriched sequencing run to determine whether or not there's a substantial number of sequences with the cap.
Unfortunately, because the cap is right at the end of the sequence, it zips through a nanopore very quickly, so may be quite difficult to distinguish from any other base at the end. There are also no official models that have been trained to specifically detect the cap, so you'd be initially restricted to looking for signal deviations. If you could find a way to extend the sequence for 20 or so based beyond the cap site, it might be possible to detect it at a single-molecule level.
Despite all that, this is something that UCSC is working on; see their repository here:
https://github.com/mitenjain/dRNA_capping_analysis

In short the answer is no, however the 'car-jumping method' described here it is possible, but it is about the enzymic MMLV reverse transcription in response to the capping structure. I assume that an alternative reverse transcriptase doesn't do 'car jumping' therefore any sequencing method will pick this up, but the difference between MMLV verse the heat-tolerant style reverse transcriptases. I would assume the secondary structure associated with viral capping prevents a low temperature reverse transcriptase, i.e. MMLV from correctly reproducing the vRNA into cDNA.
It is important to note that a cap, isn't just a virus thing, its mRNA eukaryotes per se. The cap is a 7-methylguanosine (m7G) linked to the 5′ end of the transcript by a 5′-5′ triphosphate bridge [m7G(5′)ppp(5′)N- or m7G(5′)ppp(5′)Nm in which Nm is a 2′-O methylated nucleoside].


text:
Why does Coronavirus genome contain poly A region at the 3' end in GenBank?
If you look at the Coronavirus genome below, it contains poly-A region at the end. I know that the viral RNA contains a poly-A tail. But, the genome in GenBank is stored in cDNA format. That means the end should have poly-T instead of A. Can anybody help me to understand this, please ?
>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
--------------------------------------------------------------------
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

I figured out that it has to do with second strand of cDNA. Usually genbank sequences are equivalent to 5’ to 3’ strand of cDNA. This strand is generated using second strand synthesis enzyme.

No its definitely a poly-A tail for both the in situ virus and the cDNA sequence. It couldn't be a poly-T in situ because,

the Ts are uracils (Us) because its an RNA virus. Thus if you substituted all the 'T's for 'U's you'd recover the in situ genome
specifically its a positive single-stranded RNA virus and the cDNA isn't a reverse complement but the positive strand.

If you look at the poly-A together with the 5' capping then you get the idea its to assist translocation and transcription to the endoplasmic reticulum. The genome pulls a few tricks which make it a bit more complicated than that, but thats the gist.


text:
Drug classification database information
Is there any database that give information about the drug class. Right now I have all the MeSH id of the drugs from different publications. So right now I'm looking one by one drug looking at the drug description,usage as well as I'm looking at this MeSh tree structure. There are many drug which overlaps mostly due to the core structure . For example Doxorubicin and Aclarubicin both have the same origin which is Anthracyclines. I would like to group them together based on this first.
So Is there any database where I can make query in batches to get my result.

drug description
usage
MeSH tree structure or similar information which gives output about the drugs hierarchical classification.


text:
co-expression analysis for two different tissues of same sample
I have gene expression datasets A and B that contain as many rows as genes and as many columns as samples. The rows in A and B represent a common set of genes measured in different tissues of the same 5 patients. My objective is to do a correlation between the set of genes in A and B to infer whether genes in both tissues are co-expressed based on the given condition or not?
Here are sample and patients IDs information with their Expression profile
Patient_ID <- as.character(quote(c(p_1, p_2,p_3, p_4, p_5)))[-1]
SampleA <- as.character(quote(c(A_1, A_2, A_3, A_4, A_5)))[-1]
SampleB <- as.character(quote(c(B_1, B_2, B_3, B_4, B_5)))[-1]

Expression profile of the two tissues
set.seed(1)
A <- matrix(rnorm(250), nrow = 50)
set.seed(2)
B <- matrix(rnorm(250), nrow = 50)

I have tried two approaches
first approach
library(tidyverse)
rownames(A) <- paste0("g_",1:50) ## genes
colnames(A) <- Patient_ID # patient IDs

prefix "b." is added to every gene in Tissue B to identify that they are from tissue B.
rownames(B) <- paste0("b.g_",1:50) ## genes
colnames(B) <- Patient_ID  ## patient IDs

Since both have the same column IDs, I used rbind function to join the two datasets.
Then I did correlation as a transposed matrix and extracted the top right bock of the matrix where I can get correlation values between Tissue A and Tissue B
bind <- Hmisc::rcorr(t(rbind(A,B)))
n = nrow(A)
final <- bind$r[1:n, (n+1):(2*n)]

Second approach
rownames(A) <- paste0("g_",1:50) ## genes
colnames(A) <- SampleA ## Sample IDs in A

rownames(B) <- paste0("g_",1:50) ## genes in B tissue
colnames(B) <- SampleB  ## Sample IDs in B
bind.col <- Hmisc::rcorr(t(cbind(A,B)))
bind.col <- bind.col$r

Therefore, I would like to ask which approach is correct to determine if genes of Tissue A and Tissue B are co-expressed based on the given condition?
If both are incorrect I wish to get feedback on how to proceed?
Thank you!


text:
I want to write python program to cut a DNA sequence at an EcoRI restriction site and print the two fragments after cutting
Seq =
AAAAAAAAAAATTTTTTTTTTTGAATTCCCCCCCCCCCGGGGGGGGGGG
Restriction site = GAATTC
I want to cut the sequence at GA/ATTC
I tried split method in python but It does not cut at GA/ATTC.
Please help

Have you looked up the documentation on working with restriction enzymes in biopython? http://biopython.org/DIST/docs/cookbook/Restriction.html#1.3

Nice documentation It really helped me. Thanks.

If you're able to do this after actually reading the documentation then please post a working example as an answer for others.

The obvious bare native Python 'AAAAAAAAAAATTTTTTTTTTTGAATTCCCCCCCCCCCGGGGGGGGGGG'.split('GAATTC') produces what I expect.

@tripleee, your solution would not reflect the "true biology" The restriction enzyme cuts within the recognized bases and the two resulting fragments would contain a part (different parts obviously) of the restriction site.

import Bio
from Bio.Seq import Seq
my_seq = Seq("aaaaaaaaaaatttttttttttgaattccccccccccggggggggggg")
from Bio.Restriction import *
print(EcoRI.search(my_seq))
print(EcoRI.catalyse(my_seq))

Actully the EcoRI cuts the DNA after G in restriction site GAATTC.

Okay my apologies, optical analysis is one of the reasons why bioinformatics is so important. You can mark your post as 'accepted'.


text:
Metapath2vec on Drug-ADR Heterogeneous Graph
I am new in this field and I am having some problems regarding a new project.
I built a graph using Drugbank Data connected to SIDER Adverse Reactions.
I used Organ- level Terms to classify the ADRs in 25 different groups.
I aim to use the algorithm Metapath2vec in order to cluster drugs based on the group of ADRs they lead to.
I am currently using the python package for Machine Learning on graphs Stellargraph and the related implementation of Metapath2vec.
I chose Clustering as downstream task with the resulting node embeddings, e.g. DBSCAN but the results are not promising.
I would like to have some clusters of drugs related to different adverse drug reaction.
Since I am new in this field everything I am trying is based on scientific literature but I don't know if this is the right approach for my objective.
This is the code related to the metapath2vec algorithm:
walk_length = 100  # maximum length of a random walk to use throughout this notebook

# specify the metapath schemas as a list of lists of node types.
metapaths = [
    ["drug", "adr", "drug"],
    ["drug", "adr", "drug", "drug"],
    ["drug", "drug"],
    ["drug", "adr", "group_adr", "adr", "drug"],
    ]

# Create the random walker
rw = UniformRandomMetaPathWalk(graph)

walks = rw.run(
    nodes=list(graph.nodes()),  # root nodes
    length=walk_length,  # maximum length of a random walk
    n=1,  # number of random walks per root node
    metapaths=metapaths,  # the metapaths
)

from gensim.models import Word2Vec

model = Word2Vec(walks, size=128, window=5, min_count=0, sg=1, workers=2, iter=1)

Should I change approach?
what other algorithm for representation learning can I use in order to reach my goal?
what can I improve in the presented approach to have better node embeddings?

Interesting, will be a while before I answer. I'm not 100% sure, intuititively the algorithm you propose appears to be a vectorisation. However, you are using random walks ... I've seen random walks in Monte Carlo yes .... but classification no. Could you explain the purpose of the random walk? Are you doing a simulation to assess the classification then applying to real data?

The meta-path-based random walk strategy ensures that the semantic relationships between different types of nodes can be properly incorporated into skip-gram.
A better explanation is taken from this paper: link.
I didn't take into account any simulation to assess the classification, do you have any suggestion? @Michael

Thanks,  I get the idea its for 1:1 relationships. This appears to be for natural language processing rather than drugs and  will affect vectorisation. I'll respond towards the end of the week. I need to get my head around the biology of your stuff

@MIchael It would be great, thanks .
I provide you some papers I have been referring to for this project: link link
I hope it's quite clear what's my goal.

Clearly the additional references do put Metapath2vec into good use. Its a many to many ML relationship. Most catagorization/classification particularly in deep learning is many  to one. So thats the good thing about the technique you  are using, you've 25 catagories, which is quite large for machine learning(ML) and large for deep learning. The technique appears to be half way between supervised and unsupervised learning. The random walk establishes a zone (circle) of catagorisation.
NLP Background
The situation where it would work in NLP (natural language processing) is for example a spell checker, where you have one correct word and loads of possible typos, some of which you not observed in training. So what will happen is the random walk will pick up the local vacinity around the 'correct spelling' by the vectorization (converting a word into numbers) of the incorrect spelling will land it close to 'true spelling' in a 2D vectorization plot. The training following the random walk says  'yep you belong with me'.
The papers given talk about being able to identify the context of a word within a sentence which is a very NLP thing, one vectorisation being 'bags of words'. So here the sentence represents one of many categories of information, i.e. each is the central node within the graph, the vectorisation of the sentence places of meaning closer to one 'hub' than another and again the local vacinity is identified by the random walk. In training this 'random walk vacinity' is identified, so it will algorithmically recreate the 'random walk locality' so that sentences that have not been trained can have their meaning identified.
The central theory of deep learning, which not quite true of other ML, is any algorithm can be recreated through training. So like it will sort of recreate a 'random walk zone'
Difference with drug discovery
The papers you have map drugs to genes, but it means you have to have a "true drug to gene" relationship. The problem I see here is what is 'true' for a drug-protein interaction? With a spell-checker there is one 'true spelling', with a sentence there is pretty much one 'type of meaning', like a movie was great, average, terrible based on a ratings score and stuff like gender classification. If its not 'true' the random walk will get it wrong. So its a big assumption. If you have a circle and pick a point on the edge of the circle then try and draw a circle around the original circle, you can't unless the circle is so large it risks including stuff that is not part of the original circle. The problem the papers is how do they know what is the 'centre' of the circle? Its an a priori assumption and you need really good biological rationale.
With your stuff I agree its easier because adverse drug reaction, some drugs could be fatal or toxic (traditional cancer drugs or parasite drugs), so you'd pick the worst case adverse reaction. That then centres the definition around which everything else is relative. If however your 'true catagory' turned out to be an outlier its not going to work, because the random walk only works if its starts bang in the centre of the 'catagory'. Random walks form spheres, in case you are not aware. I can see this working as an advanced 'unsupervised learning' technique.
One very real advantage is if a drug doesn't have an adverse reaction it MIGHT (I'm guessing) not give a result, in which case you are good to go - that a good drug. If that is what happened that would be amazing because its not how ML, or particularly DL works. Traditional ML/DL will force a classification even if it doesn't fit. There are methods to accept some sort of "half-classification" but I don't know how to do it. This is one of the reaons you need lots of things to train.
I am assuming that vectorisation of your drug is fine because there is a huge history of this in drug-discovery going right back to the origins of the field.
Summary
So yeah I get the idea, I see the flaws and the advantages. The final issue is that most ML requires pretty large data sets and DL require huge data sets and when you've got 25 catagories. For DL approaches the data set size might need to be huge unless you can do 'transfer learning' (backing your training onto a pre-trained model). If your method resulted in a zero classification i.e. "doesn't fit anything I've got" - that would be cool and that is the key advantage I can see in the 'random walk' approach which is not available using traditional methods.
As long as you are aware, you do not have to use this algorithm to get a classification. This final point what appears missing in the papers is using a 'traditional' method by comparison. The ML before deep learning, stuff like random forests, routinely takes loads of catagories. If you are using a non-deep learning ML you would traditional use several different approaches to justify your result.

My answer is a bit rambly but I think you'll get the gist. Peronally I would complement with a standard ML, but all methods have their flaws. Deep learning (DL) has a big history in drug discovery (alot of work was done via Tensorflow), but I'm not sure you've data set size for this.

It is all clear enough and I thank you for your precious answer.
So if I would like to classify drugs based on their ADRs I should give up on this approach and try out a standard ML classification algorithm, right?
I am quite confused about it, it is the first time working in this field and I don't really have much experience. If you could give me a more detailed explanation of how to achieve my goal, I would be grateful. @Michael

I would do the two side by side. I think this approaches has merits. All I am pointing out is there are weaknesses in the approach, but you need to be aware that ML approaches, outside deep learning, give varying answers and there's no real way to say one is better than another. This method also might have advantages that 'traditional ML' does not have (possibly a 'zero classification'). The thing about this one is at least there is a known basis where it might not give the right answer.


text:
Visualizing FindMarkers result in Seurat using Heatmap
I used FindMarkers to find the DEGs between two clusters in my dataset using Seurat. I am trying to visualize the outcome using a heatmap but I failed to write the command in R.
dput(cluster4_3.markers)
              p_val  avg_logFC pct.1 pct.2    p_val_adj
Vim    1.803527e-44 -1.1689785 0.996 0.999 3.458985e-40
Lgals7 1.479278e-38  1.7615232 0.993 0.552 2.837108e-34
Anxa2  1.797016e-37 -1.3971764 0.629 0.991 3.446497e-33
Tmsb10 6.936830e-31  0.7937881 1.000 1.000 1.330415e-26
Ptma   7.518439e-27  0.5622808 1.000 1.000 1.441961e-22
Jund   2.044576e-25 -0.9881766 0.993 0.999 3.921293e-21

When I tried:
DoHeatmap(cluster4_3.markers, features = NULL) +
scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))
Error in UseMethod(generic = "DefaultAssay", object = object) :
  no applicable method for 'DefaultAssay' applied to an object of class "data.frame"

It didn't work. I was trying to replicate this command:
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend() +
  scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))

I tried to extract some info from the csv file but I couldn't determine which parameter is important (p_val avg_logFC   pct.1   pct.2   p_val_adj). Also, I didn't know which parameter to use to use the group_by function. That's one of my attempts:
clus4vs3 <- cluster4_3.markers$p_val %>%
    group_by(pct.1) %>%
    top_n(n = 10, wt = avg_logFC)
Error in UseMethod("group_by_") :
  no applicable method for 'group_by_' applied to an object of class "c('double', 'numeric')"

Now that I understand 'group_by' only accepts a categorical variable, how can I write a command that displays the top20 differentially expressed genes between cluster 4 and 3 in a heatmap ?

Please a) include a reproducible example of your data, (i.e. by using dput(cluster4_3.markers) b) tell us what didn't work because it's not 'obvious' to us since we can't see your data. Include details of all error messages. We can't help you otherwise.

@user438383 Sorry for the ambiguity. I just attached the result of dput(cluster4_3.markers).

In your DoHeatmap() call, you do not provide features so the function does not know which genes/features to use for the heatmap.
In your last function call, you are trying to group based on a continuous variable pct.1 whereas group_by expects a categorical variable.

I understand a little bit more now. I just want to show the DEGs expressed in both clusters in 1 heatmap. But I failed (still failing) to write the command logically.

I am having a hard time understanding "DEGs expressed in both clusters", by definition a DEG is expressed in one of the clusters, not both. You can feed rownames(cluster4_3.markers) to the features variable of DoHeatmap().

What I mean is to compare the DEGs expressed between those two clusters. Does that make sense?

> DoHeatmap(cluster4_3.markers, features = rownames(cluster4_3.markers)) +
+   scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))
Error in UseMethod(generic = "DefaultAssay", object = object) :
  no applicable method for 'DefaultAssay' applied to an object of class "data.frame"

DoHeatmap() takes a Seurat object as a first argument. Please go through the documentation and vignettes.

Thank you for bearing with me.
I edited the command using the only seurat object I've created with the main set of data:
> DoHeatmap(tg14, features = rownames(cluster4_3.markers)) +
+   scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))
Error: cannot allocate vector of size 580.8 Mb

I am currently looking through the vignettes but if you can still give me more clues, that'll be great!

Thanks a lot for helping. This fixed my problem:
top10 <- cluster4_3.markers %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(tg14, features = top10$gene) + NoLegend() +
  scale_fill_gradientn(colors = rev(RColorBrewer::brewer.pal(n = 7, name = "RdYlBu")))


text:
gprofiler output ENSG to gene name
I would like to ask if there is any tool or R script that translates Ensembl IDs to gene names in a gprofiler GO output. An example of this file is the following one:

This is usually done using AnnotationDbi in combination with the annotation database for human org.Hs.eg.db.
library(org.Hs.eg.db)
genes <- c("ENSG00000102055", "ENSG00000232605", "ENSG00000136828",
           "ENSG00000159197", "ENSG00000201394", "ENSG00000278580",
           "ENSG00000213440", "ENSG00000166816", "ENSG00000275895",
           "ENSG00000105723")
select(org.Hs.eg.db, genes, c("ENSEMBL", "ALIAS"), "ENSEMBL")

This returns:
           ENSEMBL     ALIAS
1  ENSG00000102055       I-4
2  ENSG00000102055  PPP1R2P9
3  ENSG00000102055   PPP1R2C
4  ENSG00000232605      <NA>
5  ENSG00000136828   RALGEF2
6  ENSG00000136828  RALGPS1A
7  ENSG00000136828   RALGPS1
8  ENSG00000159197     ATFB4
9  ENSG00000159197      LQT5
10 ENSG00000159197      LQT6
11 ENSG00000159197     MIRP1
12 ENSG00000159197     KCNE2
13 ENSG00000201394   RN5S140
14 ENSG00000201394 RNA5SP140
15 ENSG00000278580      <NA>
16 ENSG00000213440      <NA>
17 ENSG00000166816     DLACD
18 ENSG00000166816       DLD
19 ENSG00000166816      LDHD
20 ENSG00000275895   U2AF1L5
21 ENSG00000105723     GSK3A

If you want to apply a 1:1 mapping of gene IDs to common names, you can use the "SYMBOL" column. In a data.frame, you can do something like this:
ensembl_to_symbol <- function(e) {
  hits <- select(org.Hs.eg.db, e, "SYMBOL", "ENSEMBL")$SYMBOL
  # Gene IDs without symbols should just stay as IDs
  no_hit <- is.na(hits)
  hits[no_hit] <- e[no_hit]
  hits
}

recode_symbols <- function(s) {
  e <- ensembl_to_symbol(unlist(strsplit(s, split = ",")))
  paste(e, collapse = ",")
}

df <- data.frame(
  genes = c(paste(genes[1:5], collapse = ","),
            paste(genes[6:10], collapse = ",")
  )
)

df$genes <- sapply(df$genes, recode_symbols)

This isn't the most vectorized (fast) solution, but it's probably easier to understand this way.

Yes I know the package. My problem is to use it in the format of gprofiler file. How could I  isolate each ensemble ID and replace it with a gene name, when ensembl ids are in one cell of a column of dataframe?THANKS!

I added an example of how you can use this in a data.frame.

Thnak you so much!!!


text:
How can atomic RMS differences / atomic RMS distributions help us understand the structure of a protein?
For example, I have a paper that says: "the atomic RMS distribution about the mean coordinate positions is 0.41 +/- 0.04 A for the backbone atoms and 0.82 +/- 0.04 A for all atoms"
And I have no idea what any of this means, I found several papers with this kind of information but none of them explain why this is relevant in order to understand the structure of a protein.

Hi @Ig009 there is nothing wrong with the basis of your question. The resolution in question is excellent if correctly reported. Could you elaborate on the basis of your question/investigation please?

The "RMS distribution about the mean coordinate positions" could be restated as "certainty of the structure". In other words, the main backbone that defines the large-scale shape in the structure is expected to be accurate to within about 0.41 angstrom, while the atoms more on the periphery (or otherwise not forming the backbone of the structure) are less accurately solved. How important this is will depend on what you want to do with the structure. If you want to get an idea for what it looks like then the accuracy is less important, if you want to do docking then the more accurate the atoms around any sort of pocket/interaction surface are the better.


text:
Why is there a minus in the 2ˆ(–delta delta CT) method (qPCR)
Question: Why is there a minus in the 2ˆ(–delta delta CT) formula?
My line of thoughs:
Consider Ct values of some pPCR experiment (some gene of interest X and a reference gene REF is measured under some experiment condition and a control condition).

The delta CT (dCt) would correspond to normalized expression Ct of X against REF of both conditions. (dCt = X-REF)

The delta delta CT (ddCt) would then correspond the the difference of the dCt of the experiment condition minus the dCt of the control condition.

=> If ddCT > 0 then the gene of interest is upregulated downregulated (Ct is antiproportional to gene expression). (the difference of the dCt in the experiment condition is larger than in the control condition)
But now: the ratio R = 2ˆ(–delta delta CT) would then correspond to a value smaller than 1 (because of the minus).
I read sentences like

"We see a 10% decrease of expression (R=0.9)..."

But actually the raw data clearly shows a upregulation... This is confusing since only the ratio R decreases by 10% but the this would correspond to an increase of raw expression...
I could interprete R as the ratio of the control against the gene of interest. But why not the other way around?
Edit: Ct is antiproportional to gene expression => dCt does not directly reflect the gene expression.

Sources:
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-020-03696-y
"We can predict that a plant with a hormone concentration of
8.85 pg/mL would have an expression of Gene G that is 27% (̂ =0.73) lower than that of plants with average hormone concentration" (just
below equation 8)
But clearly in Fig. 3a you can see that the expression at 8.85 is
obviously higher than with average hormone concentration (x0 = 9.85)

Perhaps this pdf of the ddCt package helps: https://bioconductor.org/packages/release/bioc/vignettes/ddCt/inst/doc/rtPCR.pdf . See equations 6, 7 and 8.

ok I found the reason: Ct is antiproportional to the gene expression. Therefore the statements from the paper are correct. But thank you for your link, very nice derivation


text:
is it ok to use gene sets from an organism to enrich genes from a different organism?
Lets say I performed a single cell RNA sequencing analysis on mice blood samples. Eventually I am able to cluster cells together and get a set of genes of interest for each cell type since I want to study the gene expression of specific cells after a treatment.
At this point I want to perform an enrichment analysis to get a list of those GO biological processes and/or molecular functions that might be affected by the novel genes that are up or down regulated after the treatment. I can use several tools such as ClusterProfiler, Panther, EnrichR but, and this is the point, should I use gene sets that are compiled for mouse only or can I enrich my genes using gene sets compiled from other organisms too such as, for instance, human gene sets? Is it ok to mix genes and gene sets from different organisms when performing an enrichment?

I am not very familiar with those programs but maybe you are searching for the term ortholog? Orthologs are "same" genes/proteins of different species (roughly speaking). So you could convert your genes into the respective ortholog gene of another species (like human) and then perform your analysis on human. You will loose some genes that do not find a counterpart and you have to resolve some ambiguities if your gene is mapped to multiple candidates (choose all/ any/ best).

@Paul, ok this could be done, of course but is it common practice? also, I need to study the genes in that specific organism not in its ortholog version.