text:
Best resource to link a gene/protein to a druggability score?
Could I ask, is anyone familiar with the Drugebility score provided by Chembl?
I wanted to add a druggability score to a set of human targets, to understand the likelihood that they are druggable.
I found Chembl's drugebility score here, but I don't fully understand which file I'm meant to use that adds a drugebility score to each human protein? In particular, should I be looking at the drugebility.dmp.gz file (that was created in 2019) or the site_drugebility.txt file that was created in 2013 (seems very long ago)?
I'm unclear in drugebility.dmp how do understand what ensembl/NCBI/uniprot ID each score is assigned to?
Or does someone have a better suggestion for how to link ideally Ensembl/UniProt IDs to a score indicating the likelihood that the protein is druggable?

The DrugEBIlity score is multiple scores (tractable, druggable & confidence etc.) and actually uses an even older pocket finding algorithm IIRC (probably written in here), but the major issue that it is solved structures only —not even close threaded models. Now with AlphaFold2 there's a whole new dataset to play with. I would bet they are going to or are currently rescoring... May be worth checking with them by email/Twitter/etc.

@MatteoFerla while AlphaFold2 has generated a lot of media attention, as far as I can tell
it has actually generated very little new data that wasn't already available (and likely more reliable) in SwissModel.

@terdon. The oligomerisation in SwissModel is a real winner and being threaded the other chains and ligands from the template can be stolen. But bar for oligomers and the spaghetti loops of AF2, AF2 is more accurate (for starters it does not blow-up in minimisation) and there are a lot of horrors in the Felsenstein zone of SwissModel (<30% id). Also, there are definitely plenty of protein which we have no idea what they do that have enzyme-like cavities solved in AF2, but not SM, e.g. Q9NX78, simply because it's nobody's favourite protein.


text:
Transfer learning for a Convolutional NN for Recursion dataset
As a practice data analysis I am trying to train a convolutional neural network (CNN) on some cellular images made publicly available on Kaggle in 2019 by the company Recursion via supervised deep learning.
I want to increase the predictive power using transfer learning at the lower layers of my CNN preferably from a public imaging model.
Are there any appropriate models available?

Good question, the HUGE advantage of CNN is transfer learning. Google do a huge photo database for transfer-learning and if that dataset had cell images then you've made it. If not - you could start the trend and train up a database.

Hmm very cool thanks for the suggestion

Would you consider andwering so I can lay this questuion to rest?

learn = cnn_learner(data,
                    models.resnet50,
                    metrics=error_rate)

resnet50 is one of many pre-trained 2D CNN computer vision models which use transfer learning. If you've classifying a common image such as a given mammal you're sorted because they'll be loads of pictures of that mammal all over resnet50. However, 'a cell' I dunno.
When the computer vision challenges were in full swing between the tech giants each put their trained models onto publication release and are easily accessible e.g. via keras, fastai etc... I'd have to look them all up and each has different architecture, they are serious architecture like 50 layer NN. The competition was eventually stopped because computer vision had a lower error rate than human vision, at which point its pointless judging between DL algorithms because we can't say what it is how can we judge an algorithm?
Transfer learning places CNN at significant advantage over RNN.
The other possibility is google images
pip install google_images_download

I don't think its on conda. Then see if you get lucky with the following
from google_images_download import google_images_download
response = google_images_download.googleimagesdownload()
arguments = {"keywords":"cells,organelles, membranes","limit":60,"print_urls":False}
paths = response.download(arguments)
print(paths)

Please note I don't do 2D CNN,i.e. images.

source: Stack Exchange

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

sort: votes


text:
Difference between bioinformatics and cheminformatics
Anyone can answer, and should be in simple way

A reminder to not downvote newbies, but inform them why their Q is problematic. In this case it's a LMGTFY situation w/ no info of what may be the source of confusion and the body is not informative.

There are four terms to consider, all are computational

cheminformatics — Big data to explore relationship between small molecules (e.g. enumeration and machine learning stuff)
computational chemistry ("compchem") - thermodynamics-based methods to explore small molecules (e.g. docking)
bioinformatics — big data primarily for exploring -omics data (genomics, proteomics, transcriptomics, metabolomics)
computational biochemistry - thermodynamics-based methods to study structure of protein

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


text:
Multiple Sequence Alignment
I applied three multiple sequence alignment methods of muscle, mafft and clistal omega on 30 sequences and 400 sequences . In the case of 30 seq, the output of alignment from the 3 methods was in the same length but in the case of larger data with different length of sequences (input data) the output from alignment methods was different length clustal length is different than the two others. My question is how can i assure that the length from MSA should be equal or different?

The length of the alignment will depend on the model. For example, high gap costs will typically result in shorter alignments, whereas low gap-costs will do the opposite.
Consider the sequences below:
AAAAA
TTTTT
GGGGG

with a very large gap cost, these sequences will result in the following MSA
AAAAA
|||||
TTTTT
|||||
GGGGG

but with lower gap costs, we'd observe
AAAAA----------
-----TTTTT-----
----------GGGGG

Of course, this is an extreme case, but it goes to show that the alignment model will vary how long the resulting alignment is.

I noticed also one thing when the input (length of sequences)  before the MSA method are on the same length, the output after alignment methods are all on equal length, however, when the input differs, the output differs as well. Is it right?

For a single run of MSA on a set of sequences, the output lengths should all be the same, regardless of input length.

It is likely the sequence data in this instance, the more sequences the more chance an issue with the alignment where variations in the gap penalty are observed.
MAFFT outperforms Clustal Omega to some extent over significantly larger alignments >50 000 sequences from published benchmark studies (FAMSA paper), but this is not relevant here.


text:
blastn returning an inferior alignment
In my work I've found that in circumstances where there is a mismatch near the edges of the query sequence, blastn prefers to return a shorter contiguous alignment, rather than allowing for a mismatch near the 5' or 3' end of the query sequence.
For example, consider the subject sequence (bold indicates alignment with mismatch in brackets):
GACCAACCTGCTGGATCTGCTGATCGGCG[G]TTTGCGAAATA
and query sequence:
CTGGATCTGCTGATCGGCG[A]T
Blasting this will return the 19bp contiguous alignment upstream of the mismatch, rather than returning a 21bp alignment with a single mismatch at the G/A base at position 20 in the query.
I'm using the command line version of blastn, and I've tried messing with the penalties as well as setting -task blastn-short, but nothing seems to fix this issue.
I speculate that the various gap and mismatch penalties are working out such that a shorter contiguous alignment is preferred over a slightly longer alignment with a mismatch, but I haven't been able to figure out how to modify the penalties in such a way as to change this preference.

Nope, altering the weights is the solution as you said. So please do answer your own question —this a totally encouraged behaviour.

Here's my working solution,
So I found a choice of parameters that resolves this, -reward 3 -gapopen 5 -gapextend 5.


text:
tmhmm installation and running
I am trying to predict transmembrane protein using tmhmm-2.0c. According to Package instruction, I have changed the path of perl inside the program like as following:
First i checked the location of perl
which perl
/usr/bin/perl

Then I checked the version also to make sure it is above 5.
perl -v
This is perl 5, version 16, subversion 3 (v5.16.3) built for x86_64-linux-thread-multi (with 44 registered patches, see perl -V for more detail)

Copyright 1987-2012, Larry Wall

I have pasted the path (/usr/bin/perl) in the tmhmm and tmhmmformat.pl
Run the cmd:
./tmhmm --short < /mnt/genome3/Lab_Users/Kishor/DISK_3/mnspt1/blood/rsem_outdir/blood.Trinity.RSEM.retained.clustered.fasta.transdecoder.pep > tmhmm.out

But I only received a blank folder. It will be great for me if I get any suggestions to solve the problem.

You're right about the Perl path, both in the tmhmm and the tmhmmformat.pl file, but if you select the short option you will not get the output to create the output plot in the folder, that is why is empty.

Thanks Zora. But also tried without setting parameter '--short' and still didnt find output except a blank folder. I actually need following typr of output to fed into Trinotate:
TRINITY_DN100_c0_g1_i1.p1 len=115 ExpAA=0.00 First60=0.00 PredHel=0 Topology=i
TRINITY_DN102_c0_g2_i1.p1 len=256 ExpAA=0.16 First60=0.07 PredHel=0 Topology=o
TRINITY_DN103_c0_g1_i1.p1 len=107 ExpAA=16.73 First60=0.00 PredHel=0 Topology=o
Can you give me any other suggestion?

can you share minimal input to reproduce your error?

@zorbax
Thanks for your concern. I have solved the issue. Needed to do some change inside the  tmhmm program.

Hi Kishor Kumar Sarker What changes you had done inside the tmhmm program?????? I also got the empty folder.
Dimple

Hi @Dimple,
In following link you will find the detail about the problem and as well as the solution at the end.
https://github.com/Trinotate/Trinotate/issues/77


text:
Snakemake: Accept any of a possible set of inputs (any one being sufficient)
I have a little pipeline that I want to execute using Snakemake.
One of my rules (let's call it create_somefile) creates a file somefile_{n_dims}.json, and when running the pipeline can specify that I want to run the pipeline for a pre-defined set of values for n_dims, which is why I specified this as a wildcard.
Now I want to have another rule, that takes as input any of these files. It can take either somefile_1.json or somefile_2.json or somefile_32402348.json, it doesn't matter which one as long as any of these exists. The rule then creates a file that does not depend on the number of dimensions, and thus doesn't want and need the n_dims in its output.
When I try to have the wildcard only in the input, like this:
rule do_stuff:
    input:
        "somefile_{n_dims}.json"
    output:
        "general_file.json"
    shell:
        "echo {input} &&"
        "touch {output}"

Snakemake gives me the error Wildcards in input files cannot be determined from output files: n_dims.
I don't want the output file to contain n_dims whatsoever, and unfortunately it's also impossible in snakemake to just create an additional throwaway temporary(..)-output that does depend on the n_dims just so snakemake is happy (error message then: ot all output, log and benchmark files of rule do_stuff contain the same wildcards. This is crucial though, in order to avoid that two or more jobs write to the same file..
Is there any way in snakemake to specify a rule's inputs as "take any of these", without wildcards?
Looking at the documentation, the only thing I found is the concept of checkpoints (https://snakemake.readthedocs.io/en/stable/snakefiles/rules.html#data-dependent-conditional-execution), however that doesn't seem to work as well - I tried making the rule create_somefile a checkpoint and changing the do_stuff rule to
rule do_stuff:
    input:
        lambda wildcards: checkpoints.create_somefile.get(**wildcards).output
    output:
        "general_file.json",
    shell:
        "echo {input} &&"
        "touch {output}"

however then I get the error WorkflowError: Missing wildcard values for n_dims.
Is there any concept in snakemake to allow me to do what I want? To have "wildcards" in inputs, that just stand for "execute this rule as soon as the input for any value of the wildcard is fulfilled, execute only once and be done with it"?

You might be able to make this work with some combination of the checkpoints/directories as outputs/functions as inputs features, but it depends on the details of the behavior you want.  What should happen if you request general_file.json with none of the possible input files already there?  Or, what if more than one of them already exist?

Well, it shouldn't be scheduled with none of the input files already there, and if it's more than one that doesn't matter, it will take just one of them. Could be specified in the shell/run section which one, but doesn't even need to be, the code itself could also just take a random one

I don't think Snakemake can follow a "use any of these" pattern with a unified dependency graph and without additional assumptions, but you can do it if you break the workflow into two pieces.
This will take whatever somefile_*.json shows up first, or throw an exception if there aren't any:
from pathlib import Path

def input_for_do_stuff(wildcards):
    try:
        return next(Path(".").glob("somefile_*.json"))
    except StopIteration as err:
        raise Exception("No n_dims JSON files found") from err

rule do_stuff:
    input: input_for_do_stuff
    output:
        "general_file.json",
    shell:
        "echo {input} &&"
        "touch {output}"

rule create_somefile:
    output: "somefile_{n_dims}.json"
    shell: "touch {output}"

For example if you have somefile_2.json already there:
$ snakemake -nrp
Building DAG of jobs...
Job stats:
job         count    min threads    max threads
--------  -------  -------------  -------------
do_stuff        1              1              1
total           1              1              1

[Sun Oct 10 12:22:14 2021]
rule do_stuff:
    input: somefile_2.json
    output: general_file.json
    jobid: 0
    reason: Missing output files: general_file.json
    resources: tmpdir=/tmp

echo somefile_2.json &&touch general_file.json
Job stats:
job         count    min threads    max threads
--------  -------  -------------  -------------
do_stuff        1              1              1
total           1              1              1

This was a dry-run (flag -n). The order of jobs does not reflect the order of execution.

With that, you'll never be able to run the workflow by just requesting the final output, which kind of breaks Snakemake's way of thinking about things.  But at least it'll work like you want if you run it in two pieces.  If having a default makes sense for your case, you could also do something like this:
def input_for_do_stuff(wildcards):
    try:
        return next(Path(".").glob("somefile_*.json"))
    except StopIteration as err:
        return "somefile_1.json"

That way, even if no files exist to start with, Snakemake can come up with a way to make the output you ask for.


text:
DEseq2 contrast in two time points
I need help with some clarifications please
I have control versus treatment in two time points like
> y
                             X condition time
1       CTRL_24_hrs_replicate1   control   24
2       CTRL_24_hrs_replicate2   control   24
3       CTRL_24_hrs_replicate3   control   24
4  treatment_24_hrs_replicate1        t   24
5  treatment_24_hrs_replicate2        t   24
6  treatment_24_hrs_replicate3        t   24
7       CTRL_48_hrs_replicate1   control   48
8       CTRL_48_hrs_replicate2   control   48
9       CTRL_48_hrs_replicate3   control   48
10 treatment_48_hrs_replicate1        t   48
11 treatment_48_hrs_replicate2        t   48
12 treatment_48_hrs_replicate3        t   48
>

I want to see the difference of treatment versus control from 24 hours to 48 hours, so this is not just the difference of treatment versus control or difference between 24 versus 48 hours rather the treatment versus control from 24 hours to 48 hours
This is what I have done so far
>  dds <- DESeqDataSetFromMatrix(countData=count,
+                                 colData=design,
+                                 design=~condition, tidy = FALSE)
Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors
>   dds$group <- factor(paste0(dds$time, dds$condition))
>   design(dds) <- ~ group
>   dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing
>   resultsNames(dds)
"Intercept"                    "group_24IT_vs_24control"         "group_48control_vs_24control" "group_48IT_vs_24control"
>

How to find the difference at time point 24 vs 48 in condition ?
I searched and I found a similar thing but I don't know how to adopt that to my data
results(dds, contrast=list(c("time_4_vs_0", "conditionLD.time4")), test="Wald")

I tried this and I got error
> results(dds, contrast=list(c("time_48_vs_24", "conditionIT.time24")), test="Wald")
Error in checkContrast(contrast, resNames) :
  all elements of the contrast as a list of length 2 should be elements of 'resultsNames(object)'
>

Can you help me in adopting this to my data please?

You can't figure out how to adapt this to your data?
http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts
Ignore resultsNames.  You don't need them for simple comparisons like you want to do.

Sorry even by ignoring resultsNames, I don't know what to put into this formula. For instance what I should put instead of time_4_vs_0 ?   This gave me error  > results(dds, contrast=list(c("time_48_vs_24", "conditionIT.time24")), test="Wald")
Error in checkContrast(contrast, resNames) :
all elements of the contrast as a list of length 2 should be elements of 'resultsNames(object)'

Where have you looked?  Did you read the link I posted?  Did you look at the entry for deseq2 results?  https://rdrr.io/bioc/DESeq2/man/results.html   Can you not find a single tutorial anywhere which shows you how to use specify the contrast?

you can do:
data$condition_new <- as.factor(paste0(data$condition, " - ", data$time))
dds <- DESeqDataSetFromMatrix(countData=count,
                             colData=design,
                             design=~condition_new
                             , tidy = FALSE)
enter code here
results(dds, contrast=c("condition_new", "control - 24", "condition - 4"))

I hope it helps.

You can find here: http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#contrasts
the tip as pointed from @swbarnes2


text:
It is possible to have a UniProt protein entry that does not have any genomic coordinates, because lacks of an Ensembl link?
I have a list of UniProt codes and I am trying to fetch the genomic coordinates through the API from programmatic access to UniProtKB (https://www.ebi.ac.uk/proteins/api/doc/).
Moreover, I read that UniProt in collaboration with Ensembl have mapped the protein sequences on the GRCh38 genome (doi: 10.1093/nar/gkw1099)
However, it is true that the number of protein codes in UniProt is updated more frequently than the GRCh38 assembly.
Considering only the protein genomic coordinates data:
It is possible that a new submission in UniProt (in TrEMBL for instance) does not have any correspondence with Ensembl, due to the possibility that the genomic coordinates of a new protein do not match any gene model stored in GRCh38?
If not:
It is possible to say that all the protein genomic coordinates in UniProt have been mapped on the GRCh38 (in connection with Ensembl) and it is impossible that a new submission in UniProt lacks of this data?
Thanks

If you look at this query
https://www.uniprot.org/uniprot/?query=NOT%20database:(type:ensembl)&fil=organism:%22Homo%20sapiens%20(Human)%20[9606]%22&sort=score
and restrict it to the reference proteome (by using the link on the left hand side):
https://www.uniprot.org/uniprot/?query=NOT%20database:(type:ensembl)&fil=proteome%3AUP000005640+AND+organism%3A%22Homo+sapiens+%28Human%29+%5B9606%5D%22&sort=score
you will see that there are "only" 1100 entries, and that they are all in UniProtKB/Swiss-Prot, but none in UniProtKB/TrEMBL, the unreviewed section.
The manual review of sequences for human and other proteomes is one of
our most important and crucial tasks but as you can imagine this is an
ongoing process. First because it represents an immense work, second
because knowledge is still evolving, the human reference genome being
for instance constantly improved by the Genome Reference Consortium. For
the human proteome, there are still a small but non negligible number of
UniProtKB/Swiss-Prot entries for which the sequence is not identical or
is not represented at all by Ensembl sequences. We are constantly
collaborating with Ensembl and other genome annotation resources trying
to agree on identical sequences to fill the gaps. UniProtKB human
entries are automatically mapped to the latest Ensembl release but you
have to take into account that UniProtKB and Ensembl releases are not
synchronized.
Starting with release 2022_01, we will have cross-references from relevant UniProtKB/Swiss-Prot entries to MANE-Select (Matched Annotation from NCBI and EMBL-EBI, see https://www.ensembl.org/info/genome/genebuild/mane.html and
https://www.ncbi.nlm.nih.gov/refseq/MANE/)


text:
Removing columns from meta.data in Seurat
I'm trying to remove columns from the meta.data from a Seurat object. Is there a way to do that or would I have to re-create the Seurat object without that column?

Set the column you don't want to NULL
x$remove.column <- NULL


text:
NMF - Subscript out of bounds
I've calculated the cophenetic coefficients for the NMF of gene expression data, but is giving error when performing the clustering information step as given below:
clustered_data <- NMF_estimate[[num.opt - 1]][[2]]

Error in NMF_estimate[[num.opt - 1]] : subscript out of bounds

Here num.opt is num.opt <- coph[which(coph[,2]==max(coph[-1,2])),1] which gives 2, 5 as outputs, and NMF Estimate is a list of 6 elements.
Please suggest a solution as soon as possible.


text:
Why does it require 4 vectors to calculate a torsion angle?
The BioPython documentation shows the following:

11.6.3 Measuring torsion angles Use the vector representation of the atomic coordinates, and the calc dihedral function from the Vector module:
>>> vector1 = atom1.get_vector()
>>> vector2 = atom2.get_vector()
>>> vector3 = atom3.get_vector()
>>> vector4 = atom4.get_vector()
>>> angle = calc_dihedral(vector1, vector2, vector3, vector4)

We know that a torsion angle between N-C-alpha and C-alpha-C planes.
So, the minimum number of points required to form a plane is 3.
Then, Why does BioPython require 4 vectors?

You say you need two planes and that you can describe each plane with three points. That's 6 points total, which in three-dimensional space require 18 coordinates. Four vectors require only 12 coordinates. What is there to complain about?

If you want to check out the actual source code to clarify this further see https://docs.python.org/3/library/inspect.html#inspect.getsource —a handy function in a great package.

Vector spaces include the origin, by definition. The origin supplies two of the six points need to define the two planes.

From wikipedia on torsion angles:

A dihedral angle is the angle between two intersecting planes or half-planes. In chemistry, it is the clockwise angle between half-planes through two sets of three atoms, having two atoms in common.

Two sets of three with two in common makes four total atoms.
You are defining not one plane but two half-planes or planes.
Wikipedia includes the following as an illustrative image:

"Angle between two half-planes (α, β, green) in a third plane (pink) which cuts the line of intersection at right angles"
In this image I interpret the atoms to be at the vertices of the planes, e.g. the two atoms that are in common between the two planes are those at the points of intersection of the planes, and the two other atoms extend the planes $\alpha$ and $\beta$.

Could you be kind enough to add a picture/diagram?

@user366312 I have added the image from the wikipedia page.

To compute dihedral angles you require 4 atoms and hence for vectors to determine two planes (see this figure here) whose angle you want to measure - dihedral angle.
You can also check Fig 55 from the GROMACS documentation
https://manual.gromacs.org/documentation/current/reference-manual/analysis/bond-angle-dihedral.html
Definition of (virtual) Cα backbone bond (κ) and torsion (τ) angles.
from https://journals.plos.org/plosone/article/figure?id=10.1371/journal.pone.0215141.g002

this image is a much better illustration than the one I got from wikipedia


text:
How to convert the chromosome position for insertions/deletions to rsIDs in a GWAS summary statistics?
I am trying to convert chromosome positions to rsIDs of a GWAS summary statistics file. I used bedtools intersect to merge the reference genome GRCh37 and the summary statistics based on the same chromosome positions.
But, INDELs regions resulted in several rsIDs. Some regions counting +1 base pairs regions resulted in duplicate rsIDs. In some cases, a single SNP related to many positions were not included in the summary statistics file.
I used this command intersectBed -sorted -a sumstats.bed -b 00-All.vcf -wb but it did not work. Please let me know if anybody who had a similar experience and solved it before.


text:
Colouring amino acid substitutions by scoring matrix
I have a list of site-specific amino acid substitutions from a mutagenesis experiment which I am visualising by colouring the structure of the protein in PyMOL. I would like a way to colour the substitutions by some approximation of the biochemical (dis)similarity between them.
So far, I have tried to create a colour gradient of scores from an alignment scoring matrix (e.g. BLOSUM62, therefore a 16-shade white to red gradient; a substitution of A to R for example gets a -1). However, I am fairly certain this is not the best interpretation of of the log-odds scores in the BLOSUM matrix, and I wondered if anyone had any thoughts on how I might better accomplish this?
Note that other colour schemes such as hydrophobicity, Tailor etc. aren't suitable for my purpose.

If you wanted to plot the diversity of an MSA Cosurf (https://consurf.tau.ac.il/) would the tool to go for (consurfDB is the precalculated version). For a set of point mutations and using PyMOL only (ouch!), I would suggest checking this PyMOL snippet, which gives a rough ∆∆G, which could be set as a b factor (alter resi xx, b=1.23) as played with that way.

*) I mention ∆∆G because that would be a metric of change in stability. Whereas a substitution matrix is a metric of how uncommon such a change is from a selection point of view. If your experiment is a Lenski style experiment (long term evolution experiment), that makes sense, but not if it's directed evolution following in vitro mutagenesis (mutazyme, magnanese, wobbly analogues etc) even if with in vivo selection (e.g. phage display, biosensor etc etc.) as the mutational spectrum is a product of the in vitro step not biological selection of a population.

Hi @MatteoFerla many thanks for your thoughts. I think the Gibbs energy approach is maybe not quite what I am looking for because we know that both the start and the end point fold and are active - we are particularly interested in IDing important substitutions for changes in activity between the two isoforms (this being a somewhat unguided experiment, based on ancestral sequence reconstruction)

Ah. That makes sense if it's an ASR —plus you probably have some rad epistatic mutations, which the ∆∆G of each in isolation would have been dismal. A substitution matrix is poor metric for dissimilarity because it would just light up the surface while in reality the slow evolving core mutations are the interesting ones. Consurf would reflect this, but is not specific for your changes. Is surface/core the reason why hydrophobicity, volume etc etc. are not suitable?

Indeed. So ideally we would like some way to winnow this (obviously other than exhaustive combinatorial mutagenesis, which we have in the pipeline).


text:
Different contrasts in DESeq2
I have a treatment and control in two time points like this
 > design
                                 X condition time
    1       CTRL_24_hrs_replicate1   control   24
    2       CTRL_24_hrs_replicate2   control   24
    3       CTRL_24_hrs_replicate3   control   24
    4  treatment_24_hrs_replicate1        t   24
    5  treatment_24_hrs_replicate2        t   24
    6  treatment_24_hrs_replicate3        t   24
    7       CTRL_48_hrs_replicate1   control   48
    8       CTRL_48_hrs_replicate2   control   48
    9       CTRL_48_hrs_replicate3   control   48
    10 treatment_48_hrs_replicate1        t   48
    11 treatment_48_hrs_replicate2        t   48
    12 treatment_48_hrs_replicate3        t   48
    >

I want to test between between treatment and control considering time point 2 hours to 4 hours
I have done like this
dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)

ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)

> resultsNames(ddsTC)
[1] "Intercept"                "condition_t_vs_control"
[3] "time"                     "conditiont.time"

What is the difference of the results of these codes
results(ddsTC, name="condition_t_vs_control", test="Wald")

versus
results(ddsTC, name="conditiont.time", test="Wald")

Actually I want to know what condition_t_vs_control gives and what conditiont.time gives
Thank you so much for any intuition

condition_t_vs_control gives the effect of condition, conditiont.time is the interaction of condition and time. If you wanted to test the effect of time, use name="time".

Thank you, using which contrast I can get the difference of IT versus control but considering time goes from 24 hours to 48 hours?

Just condition_IT_vs_control, since time is compensated for by fitting all coefficients at the same time (let me know if that phrasing is confusing).

Sorry @Devon Ryan I get  error for ddsTC <- DESeq(dds, test="LRT", reduced = ~ time + time:condition)                                                                                error in nbinomLRT(object, full = full, reduced = reduced, quiet = quiet,  :
less than one degree of freedom, perhaps full and reduced models are not in the correct order  what should I do please ?

Make sure to specify the full model, I don't know what that defaults to.

Thank you so much for your answer, I think this is my full model                                              dds <- DESeqDataSetFromMatrix(countData=a,colData=design, design=~time + condition + time:condition)   am I right? design=~time + condition + time:condition  When I try to reduce that by removing reduced = ~ time + time:condition I get error

Then convert time to a factor. Perhaps that will solve the problem. You have enough samples that you shouldn’t run into this issue. BTW you can use the default Wald test instead and avoid the error.

Thank you so much says that > ddsTC <- DESeq(dds,reduced = ~ time + time:condition)
Error in DESeq(dds, reduced = ~time + time:condition) :
'reduced' ignored when test='Wald'

Correct, there is no reduced formula with a Wald test, you just extract the appropriate coefficients with lfcShrink()


text:
I can't launch FastqCleaner I always get a warning message and the application never starts
I tried to install all the needed and related packages but I still did not know what the problem is, Can anyone please help if anything else I can do??
I always get this over and over:
Warning: Error in : Navigation containers expect a collection of `bslib::nav()`/`shiny::tabPanel()`s and/or `bslib::nav_menu()`/`shiny::navbarMenu()`s. Consider using `header` or `footer` if you wish to place content above (or below) every panel's contents.
  85: stop
  84: FUN
  83: lapply
  82: buildTabset
  81: tabsetPanel_
  80: func
  78: tabsetPanel

This error message looks incomplete. It would help for you to post your exact workflow and the full output. Links to the tool/documentation that you used are also always helpful.


text:
UPGMA computation
The unweighted pair group method with arithmetic mean (UPGMA) is a hierarchical clustering method, for example used in phylogeny. In a phylogenetic tree it would  result in a global molecular clock for taxa under consideration.
Two clusters $A$ and $B$ with the minimum distance are selected and merged to form a new cluster $C$, where the distance between $C$ = $A \cup B$. Assuming another $D$, we can compute:
$$
d(C, D) = \frac{1}{|C| \cdot |D|} \sum_{x \in A \cup B,\, y \in D} d(x, y).
$$
I think $d(x,y)$ is the distance between objects $x$ and $y$ in the input matrix.
Given the information, is there a way to prove/disprove the following statement?
$$
d(C, D) = \frac{1}{|A| + |B|} \Big( |A| d(A, D) + |B| d(B, D) \Big).
$$

I’m voting to close (=migrate) this question because this is on-topic for CrossValidated rather than Bioinfo.

UPGMA has been widely used in biological systems involving multivariate statistics and phylogenetics. It has been widely used in the early phases of genomics and might still be used in commercial packages. Understanding it and why it has limited application is kinda important.

The key words in UPGMA are 'pairwise' and 'mean'.
In phylogeny the format is a bit unusual because $C$ here is the MRCA (most recent common ancestor) for $A$ and $B$. In other words this represents a three taxa tree. I agree in UPGMA this works, because of its use of the mean (which is also its criticism - below).
Under UPGMA - in your equation $d(C) = \frac{d(A,B)}{N}$ ...where $N=2$, i.e. mean.
Thus under UPGMA
$d(C,D) = d(A,D) - \frac{d(A,B)}{2}$
There is no need to include $d(B,D)$ because $d(B,D) = d(A,D)$, its the mean that was resolved in establishing $C$. If you want to assign $E = D U A$ that easy too.
Your equation doesn't work because |A| and |B| are both 0 under UPGMA, i.e. is pairwise so d(A,A) is always zero. Thus under your equation $d(C,D) = 0$ which simply means $D$, $A$, $B$ are identical and results in a matrix of 0.

UPGMA isn't really used in general in biological systems, because the use of the mean rarely makes it a good assumption. I 100% agree has been used in genomics in the modern era, absolutely true, its use is not encouraged because in results in artifactual clusters due its inability to accommodate unequal variation in the underlying genetic process. It could be used in for example in agricultural systems because these conform to uniform distributions, e.g. the distances between apple trees in an orchard. In the natural environment species clump and don't show uniform distributions. UPGMA was used in phylogenetic but the assumption of a global clock very rarely held, so it was discontinued.


text:
Calculating percent mouse cells per sample
I want to see the correlation between % mouse reads per sample vs % mouse cells per sample. I've already calculated the % mouse reads per sample but I'm stuck on calculating the % of mouse cells per sample and creating a scatterplot.
I calculated the % mouse reads per sample by doing the following loop:
sampleid <- levels(df_combined$orig.ident)
df_fractions <- as.data.frame(sampleid)
for (i in 1:nrow(df_fractions)) {
sample <- df_fractions[i,1]
sample_tmp <- subset(df_combined,orig.ident == sample)
tmp_mouse_total_UMI = sum(sample_tmp$mouse_gene_UMI)
tmp_human_total_UMI = sum(sample_tmp$human_gene_UMI)
tmp_Total_UMI <- tmp_mouse_total_UMI+tmp_human_total_UMI
df_fractions$fraction_mouse[i] = tmp_mouse_total_UMI/tmp_Total_UMI
df_fractions$fraction_human[i] = tmp_human_total_UMI/tmp_Total_UMI
}

Any help is greatly appreciated it.


text:
Error to running "CIRIquant" function in CIRIquant package at python Language
"CIRIquant" is a Python package for detecting Circular RNA. I installed "CIRIquant" through the below command in terminal:

pip install CIRIquant

My Linux distribution is Ubuntu 16.04 and my python version is Python 2.7.12. After installing, I got the below message:
“Installing collected packages: numpy, scipy, scikit-learn, pysam, PyYAML, argparse, numexpr, CIRIquant
Successfully installed CIRIquant-1.1.2 PyYAML-5.1.1 argparse-1.2.1 numexpr-2.6.9 numpy-1.16.4 pysam-0.15.2 scikit-learn-0.20.3 scipy-1.2.2”
So, When I call CIRIquant for running its function in Terminal, I got the below message in Terminal:
CIRIquant: command not found

Now, I don't know what is my problem?
I appreciate it if anybody shares his/her comment with me.
Best,

did you activate your virtualenv environment? source ./venv/bin/activate

Dear @zorbax, No. I didn't activate virtualenv. is it necessary? because I didn't install CIRIquant from the source code. could you please guide me?

When I try this as a regular user on our Ubuntu 16 LTS system, I end up with ~/.local/bin/CIRIquant from the pip install, so maybe just add that directory to your path:
$ export PATH="$HOME/.local/bin:$PATH"
$ CIRIquant -h
usage: CIRIquant [-h] [--config FILE] [-1 MATE1] [-2 MATE2] [-o DIR]
                 [-p PREFIX] [-t INT] [-a INT] [-l INT] [-v] [--version]
                 [-e LOG] [--bed FILE] [--circ FILE] [--tool TOOL]
                 [--RNaseR FILE] [--bam BAM] [--no-gene]

Dear @Jesse I set export PATH="$HOME/.local/bin:$PATH" and I could run CIRIquant -h  but now when I run e.g: ls  then the system give me below comment: "ls: command not found" and error related to PATH environment variable

Make sure you kept your old $PATH as part of the new string, rather than having only the new path (which is what that sounds like to me).  You can echo $PATH before and after the change to check it.  You should still see things like /usr/bin in there afterward.


text:
Is there a sam flag for all&none of the reads?
a. Is there a SAM flag that specifies all reads?
b. Is there a SAM flag that specifies none of the reads?
So that if I run samtools view -f (a) -F (b)
the result will be all reads of the file (as if omitting the specification)  ?

Instead of doing this kind of tweaking I would write my script in a way that these options are passed as strings that can be empty. Something like:
#/ dummy syntax just for illustration:
if(you_want_filtering){
  filter_f='-f something'
  filter_F='-F something_else'
} else {
  filter_f=filter_F=''
}

samtools view (...) ${filter_f} ${filter_F}

So if you want all reads then make these strings empty. That is far more generic than your idea because it can well be that the flag values you need do not exist.

This is the proper approach and I'll add that you can't specify all reads using -f unless all of the reads happen to have a given flag set (or all not have a particular flag set). The best you could do was exclude a flag of 4 to ensure only mapped reads.


text:
What is the 4th "level" of version number for the S. cerevisiae reference genome?
yeastgenome.org has releases called R64-1-1, R64-2-1, etc, with 3 "levels" of versioning but in a few places I have seen a "fourth level". For example, Ensemble calls their hosted version of R64-1-1 "version 51" (http://ftp.ensemblgenomes.org/pub/fungi/release-51/fasta/saccharomyces_cerevisiae/dna/). And the pre-created databases available with snpEff 4.3 include "R64-1-1.75" and "R64-1-1.99". I assume the "75" and "99" here correspond to Ensembl's "51". What is this fourth number and what does it mean?


text:
Understanding the different designs in DESeq2
I am using DEseq2 and trying to understand the results obtained using different models.
I have a data design with 2 genotypes and 2 time points.
    sample genotype time
1   WT_S1       WT   T1
2   WT_S2       WT   T1
3   WT_S3       WT   T1
4   WT_S4       WT   T2
5   WT_S5       WT   T2
6   WT_S6       WT   T2
7   KO_S1       KO   T1
8   KO_S2       KO   T1
9   KO_S3       KO   T1
10  KO_S4       KO   T2
11  KO_S5       KO   T2
12  KO_S6       KO   T2

I want to know the differences in results obtained from genotype coefficient and time coefficient when using different models.
Model1) ~ genotype
Since it is comparing the differences in genotype regardless of time , samples 1-6 vs 7-12 are being compared.
Model2) ~ time
This model compares time T2 vs T1 regardless of genotypes. So it is comparing samples (1-3 + 7-9) vs (4-6 +10-12), is this correct?
Model3) ~ genotype + time
My understanding is that this model assumes the genotype effect is the same at both time points and so it adds a time effect to both genotypes.
Does that mean results(obj3, name="genotype_KO_vs_WT") give differences in genotypes by comparing samples 4-6 vs 10-12 ?
What samples are being compared in the results obtained from results(obj3, name="time_T2_vs_T1")? How is it different from model 2?
Model 4) ~ genotype + time + genotype:time
Here I understand results(obj4, name="genotype_KO_vs_WT") gives the differences in genotypes at reference levels ie, samples 1-3 vs 7-9
results(obj4, name="time_T2_vs_T1"). What samples are compared here? How is it different from ‘time_T2_vs_T1’ results in model 2 or model3?
The interaction term as I understand is giving the specific effect due to KO at time T2 controlling for the baseline differences in genotypes.
results(obj5, name="genotypeKO.timeT2")
What samples are being compared from the model matrix?

I strongly suggest that you take your normalized counts, get the averages of different groups, and compare those ratios to the different designs.
That said, you might not be able to replicate the fold changes from ~ genotype + time well in Excel.  What that design does is basically say "Lots of the variance within each genotype/time point is due to the time point/genotype.  Model that so that those differences are taken into account."  Whereas a design of just ~ time, the software will go "Oh, there's a lot of variance within each time point".  You won't be able to replicate that well in Excel...but that's why you are using DESeq.
That design is probably better than just time or genotype alone, if you want to look at all the KO versus all the WT

Does that mean results(obj3, name="genotype_KO_vs_WT") give
differences in genotypes by comparing samples 4-6 vs 10-12 ?

No.  The simplest way to do that comparison, and comparisons like it is to make a third data column that's genotype_time, and use contrasts to specify what subgroup to compare to what subgroup.

Here I understand results(obj4, name="genotype_KO_vs_WT") gives the
differences in genotypes at reference levels ie, samples 1-3 vs 7-9

Yes.  But since that's super confusing, I'd recommend against it.  The way I outlined above is far more readable, and should work the same.
The interaction design is best used when you want to see of the changes caused by time are different in the KO than the WT.  In theory, you could just generate the two timepoint fold changes, and subtract one from the other, but you'd have no p-value.  Using a design with interactions will get you that fold change difference with a p-value.

Thank you very much for the detailed reply. I want to make sure if I am understanding it correctly.  So in the (obj3, name="genotype_KO_vs_WT") from ~genotype+time model, time differences are included  as well, is that correct? Also, how is (obj3, name="timeT2vs_T1) in this model is different from just ~time model?

In the same way.  Remember that the order of the elements in the design does not matter at all if you specify what you want in the results command.


text:
Shotgun analysis: Where and how to start?
I am a complete rookie in shotgun analysis, having only performed marker gene analysis before. I have a pipeline that I want to try but I can't wrap my head around how it all works.
Please correct me if I'm wrong:

I start with a paired-end Illumina short reads fastq files (.R1,.R2)

I will use FASTQC to trim according to quality. I assume that since data has been demultiplexed, the primers and and adaptors are removed as well.

? First question, a lot of people advise to merge all the samples together and then do the co-assembly since apparently it helps with the assembly. Would this be necessary if I also have nanopore reads and only 9 samples? If if I did this, how do I then separate the samples from each other?

Because I also have nanopore, I will use a hybrid assembler to get contigs. The output will be a file with contigs.

I then plan to use a binning tool (Metawrap) to create bins.

? I believe the output of these will be a file containing bins per sample. How do I estimate the abundance of each bin in the sample? And how do I factor in the depth of sequencing to normalise them? This is my biggest problem, I can't wrap my head around it.
?Can some explain VERY plainly what is sequencing depth and how do I factor it into my analysis? I've seen some people using it for normalisation of the bin/contigs abundance but again, this is idea is not intuitive to me at all.

From then on, taxonomy and functional annotation is pretty straightforward.

#Update
So a little background on the project. I have samples sequenced from a bioreactor containing a small community of bacteria (no host removal required). The samples derive from 6 different reactors = 6 samples as well as one reactor will be samples three times (time series) so I have 9 samples in total. The shotgun Illumina will be combined with nanopore. The idea is to get the composition of the community (what kind of bacteria are there?), abundance (what is their relative quantity?) but also functions (what do they do?).
I'm sorry for a very generic question, I would really appreciate somebody answering very plainly without too much jargon. Hopefully it will help me and other people who are just starting out with the shotgun.
Thank you!

Welcome to the forum.  I think that in order to get better responses it would help for you to define what your goal is for this data and analysis.

Thank you Bioathlete! I added it to the original post.

Question 1
For co-assembly vs individual assembly, it is sometimes advisable to do them separately, perhaps if you expect samples to be very different or if you have huge amounts of data. For <10 species I imagine that you will have quite high coverage, especially with nanopore.
Some metagenomic binners (e.g. MetaBAT2, a standard from a few years ago) can even incorporate information across samples, which might argue for maintaining at least some separation.
"How do I separate the samples from each other" I guess the point is that you don't, for assembly at least. You can then go in post hoc with your sequencing data from each sample and try to figure out with abundance analysis which of your genomes are in each sample.
Question 2
For abundance analysis/depth, there are strong statistical parallels in RNA-seq analysis, as Lior Pachter is fond of saying, if that is helpful. Any number of metagenomic shotgun abundance computations can be used to estimate abundance of genomes, including Pachter's metakallisto (linked therein). I found metakallisto fairly straightforward to use on contigs resulting from assembly.
Basically, you have some set of DNAs in your sample with some stoichiometry (relative abundance with respect to each other). The idea is to use DNA sequencing to estimate that stoichiometry (absolute abundance is much harder so most people don't bother).
Here is a broad view of the workflow:

You come up with some estimate of the DNAs present in the community. This can be a shotgun assembly, or some set of reference genomes, or whatever.

You count how many sequencing reads derive from each of those DNAs, usually via alignment or pseudoalignment.

YOu do some numerical corrections for systematic covariates such as DNA length. (These can look like RPKMs or TPMs from RNA-seq; metakallisto literally uses TPM calculations)

The resulting counts should be abundances, more or less.

You don't have to figure this stuff out yourself. There are various end to end pipelines that do this.
You will almost certainly need to compute abundances of contigs if you want to do metagenomic binning.


text:
Different structures of same protein from two softwares
I used AlphaFold and Swiss-Model to get the structure of a protein (KiSS1R) and then compare them. The structures don't match and  I was wondering why the softwares produce models that are different in some regions?
Thankful for all help!

Superposition

https://michelanglo.sgc.ox.ac.uk/data/d6e427bb-a0cb-4a65-b988-b08237f5054a

SwissModel

https://swissmodel.expasy.org/repository/uniprot/Q969F8

SwissModel models are threaded against a template. In this case the homologue has 31% identity, which isn't great. The template is 6TP3 which has a rather high R_free of 0.32 and it's solved at 3 Å, but it's a membrane structure so I personally would forgive it —others may not. SwissModel has switched primary metric to its QMEANDisCo, which is present as global (0.6 in this case) and per residue (saved in the PDB file as the temperature factor column), which is roughly similar conceptually to AlphaFold2's pLDDT (vide infra), i.e. it tells you what resides to trust. 0.6 not great, not terrible. A nice thing in SwissModel is that, by virtue of being threaded, it's easy to steal ligands and other protein and it shows models of oligomers.
AlphaFold2

https://alphafold.ebi.ac.uk/entry/Q969F8

AlphaFold2 is by far the best modelling algorithm as see in the latest CASP competition (#14). However it natively doesn't do oligomers —ColabFold notebook can.
AlphaFold2 gives a per residue score (plDDT). AlphaFold2 models the whole sequence, which results in some nasty "spaghetti loops": these have really bad scores, so need to be ignored. These are not trash, but are likely unstructured in isolation, but bind other protein etc. (more info in this blog post). This structure has N- and C-terminus loops that need to be ignored. The N in particular is the signal that would be cleaved off —to get the membrane embedding dots, one can use the OPM server.
Choice
The difference in the two models only involves low certainly regions and the support for the SwissModel is problematic. So the AlphaFold2 model is better, but requires removal of the C- and N-terminal loops.


text:
Convert Entrez gene ids and HUGO symbols, whose genome assembly is unknown, to GRCh38 Ensembl gene ids
I want to use a file available on a database called DepMap (https://depmap.org/portal/download/). This file called CRISPR_gene_effect.csv contains gene ids in the format HUGO (Entrez) eg. TP53 (7157), without having any complementary information regarding which genome assembly (GRCh37/38) the Entrez ids (and also HUGO symbols) belong to.
I want to convert these gene ids in the HUGO (Entrez) format to  GRCh38 Ensembl ids (specifically to Ensembl release 100), so that I can map the data to Ensembl's annotations. I would like to get suggestions on this.
I feel that figuring out the genome assembly version (GRCh37/38) of the Entrez ids would be critical in this conversion. Although, I am not aware if Entrez ids change from one genome assembly to other.

Hugo names don't really have an associated assembly. They are just the official name of the gene in human. You will find them in all assemblies (except if they hadn't been annotated/identified in one of the older ones). Entrez entries are also not tied to a specific genome assembly. They usually have mappings to multiple ones. For example, Hugo gene name A1BG, with Entrez ID 1 has:

Now, if you want to get Ensembl gene IDs for a list of Hugo gene names, you can use BioMart to convert. First, get the list of target IDs. You can use either Hugo or Entrez. I found the Entrez simpler to parse out, so I went for those (note that this assumes GNU grep):
grep -oP '(?<=\()\d+(?=\))' CRISPR_gene_effect.csv > entrez.ids

I then went to BioMart and selected "Ensembl genes 104" (note that I could only get Ensembl version 104, I assume you're not adamant about getting version 100 given that it doesn't seem to be current) and "Human genes GRCh38.p13" :

Next, I set the filters to "NCBI gene (formerly Entrezgene) ID(s)" and uploaded my file:

And, finally, selected the Attributes I want to see, namely Ensembl gene ID, HGNC (Hugo) symbol and Entrez ID so I can map everything to each other:

I then clicked on "Results" and got:
[
You can now download this to a file and have your mapping. You might even be able to use the one I created: https://www.ensembl.org/biomart/martview/352aae8e0b5c774cbbee73d67634d577.

Your answer is very helpful. Many thanks!


text:
Error in running CIRIquant package for 2 RNA-seq in Ubuntu 16.04
"CIRIquant" is a Python package for detecting Circular RNA. I ran "CIRIquant" individually and it was run. consequently, I ran CIRIquant through the below command and parameters in the terminal:
 ~/.local/bin/CIRIquant -t  4

     -1 ./SRR1176949_1.fastq.gz
     -2 ./SRR1176949_2.fastq.gz
     --config ./CIRIquant-Config.yml
     -o ./SRR1176949
     -p  SRR1176949

But I got the below message:
Traceback (most recent call last):
  File "/home/genetic/.local/bin/CIRIquant", line 11, in <module>
    sys.exit(main())
  File "/home/genetic/.local/lib/python2.7/site-packages/CIRIquant/main.py", line 89, in main
    config = check_config(check_file(args.config_file))
  File "/home/genetic/.local/lib/python2.7/site-packages/CIRIquant/utils.py", line 82, in check_config
    globals()[i.upper()] = check_file(config['tools'][i])
  File "/home/genetic/.local/lib/python2.7/site-packages/CIRIquant/utils.py", line 49, in check_file
    raise ConfigError('File: {}, not found'.format(file_name))
CIRIquant.utils.ConfigError: File: /home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI/bwa, not found

###########
My "CIRIquant-Config.yml" content is as below:
tools:
  bwa: /home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI/bwa
  hisat2: /usr/bin/hisat2-2.2.1
  stringtie: /usr/bin/stringtie
  samtools: /usr/bin/samtools-1.9

reference:
  fasta: /home/genetic/ncbi/dbGaP-26368/fasta_file/hg19.fa
  gtf: /home/genetic/ncbi/dbGaP-26368/gtf_file/gencode.v19.annotation.gtf
  bwa_index: /home/genetic/ncbi/dbGaP-26368/bwa-sam-CIRI/bwa-0.7.17/hg19
  hisat_index: /usr/bin/hisat2-2.2.1/hg19

I apprecite it if anyone share his/her comments with me.

Do you see bwa under /home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI?

Dear @haci , I installed bwa aligner in bwa folder. In fact, in "/home/genetic/ncbi/dbGaP_26368/bwa_sam_CIRI/bwa" bwa is a folder and the executive file of bwa is in  bwa folder.  should I Change this path structure?

I believe adding /bwa right after bwa in your tools: section of the yaml file is worth a try.

Dear @haci, I did your recommendation in tools :  section and was passed but Now I have the same error for refrence :  section for fasta file as below:        CIRIquant.utils.ConfigError: File:/home/genetic/ncbi/dbGaP-26368/fasta_file/hg19.fa , not found . I appreciate it if you share your comments with me

@haci I downloaded hg19.fa from : [link] (https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/)     through this link:   [link] (https://hgdownload.soe.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz)

Once again, do you see hg19.fa under /home/genetic/ncbi/dbGaP-26368/fasta_file/? h19.fa.gz is not the same as hg19.fa, the former is compressed. Most aligners if not all are capable of working on compressed fasta files so you can give a try by altering reference section in your yaml file to reflect that.


text:
Combine all alternate alleles
I have a vcf with multi-allelic SNPs, e.g. REF: A and ALT: G, C.
I would like to combine all reference alleles to the most frequent reference alleles. This should reflect both in the REF field and in the genotypes themselves.
Can this be done with an existing tool?

Please [edit] your question and add a few lines from a VCF file that we can use to reproduce your issue. Make sure to also show us the output you expect from those lines.

I would suggest computing the population frequencies of each allele, using suggestions here: https://www.biostars.org/p/366694/. You will then have the frequencies of each allele, you can make up your own mind which alleles to combine.


text:
Bowtie2 gets stuck on alignment
I am aligning a fastq file as follows:
cutadapt -j 8 --quality-base 33 -l 20 \
     -o - /home/user/Documents/analyses/CRISPR-screens/screen3/D16.fq.gz |
       ~/bowtie2-2.4.3-linux-x86_64/bowtie2 --no-hd -p 8 \
         -t -N 0 -x ~/moffat_tko3 - |
           sed '/XS:/d' | cut -f3 | sort | uniq -c > ~/D16.guidecounts.txt

The expected output is similar to this:
head D16.guidecounts.txt
329 A1BG_sgA1BG_1
669 A1BG_sgA1BG_2
65 A1BG_sgA1BG_3
397 A1BG_sgA1BG_4
100 A1CF_sgA1CF_1
731 A1CF_sgA1CF_2
396 A1CF_sgA1CF_3
536 A1CF_sgA1CF_4
427 A2ML1_sgA2ML1_1
978 A2ML1_sgA2ML1_2

However, there is no output from this command as bowtie2 gets stuck. cutadapt works fine, but when that step is done, nothing happens for hours and it should only take a few minutes.
I have tried to take a subset of reads using head and tail as follows:
zcat D16.fq.gz | head -10000 | gzip -c > D16head.fq.gz
zcat D16.fq.gz | tail -10000 | gzip -c > D16tail.fq.gz

What confuses me is that these files don't freeze bowtie2 and align well (> 90% alignment rates). Other files from the same sequencing run also align fine (NovaSeq SP flow cell 50 bp paired-end).
As far as I can tell the original fastq file looks fine:
zcat D16.fq.gz | head -12
@A00489:1088:HKN2TDRXY:2:2101:1334:1172 1:N:0:GCCAAT
GCCCAGTGACCCTTTCACGGGTTTTAGAGCTAGAAATAGCAAGTTAAAAT
+
FFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,,FFFF
@A00489:1088:HKN2TDRXY:2:2101:1606:1172 1:N:0:GCCAAT
AAAGTTCTATCTGTTGGGCGGTTTTAGAGCTAGAAATAGCAAGTTCAAGT
+
FFFFFFFFFFFFFFFFFFFFF:FFFFFFF:FFFFF:FFFFFFFFF,FFFF
@A00489:1088:HKN2TDRXY:2:2101:1624:1172 1:N:0:GCCAAT
GCGACGGTTCATCTCCGCAAGTTTTAGAGCTAGAAATAGCAAGTTAAAAT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFF:FFFFFFFF,FFFF

zcat D16.fq.gz | tail -12
@A00489:1088:HKN2TDRXY:2:2278:30644:37059 1:N:0:GCCAAT
AGGTCATCCACTCTAGACCTGTTTTAGAGCTAGAAATAGCAAGTTAAAAT
+
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
@A00489:1088:HKN2TDRXY:2:2278:30915:37059 1:N:0:GCCAAT
ATGAGGAGGGTGATGTTCCGGAGTTTTAGAGCTAGAAATAGCAAGTTAAA
+
FF:FFFFFFFFFFFFFFFFFFFFFFF,FFFFFF:FFFFFFFFFFFFFFFF
@A00489:1088:HKN2TDRXY:2:2278:31638:37059 1:N:0:GCCAAT
GAGTGTGGTTGAGATCCCAAGTTTTAGGCTAGAAATAGCAAGTTAAAATA
+
FF:FF:F:FFFFFFFF,FFFFFF:FF::FFFFFFFFFFFFFFFFFFFFFF

Does anyone know what could be wrong?

I would suggest breaking up the piped commands and see if you get some errors when the commands are run separately


text:
Is there any genetic sequence, such that, for all the kinds of mutation, it still has a 100% probability to be passed on?
This question is more on the theoretical-encoding side. Given the known kinds of genetic mutation, transposition, deletion, and such, what can be said about genes that would be left untouched by the mutation, left invariant by the permutation, if you will? Are there any such sequences? Would they be more common in a population, if not detrimental?

By definition, mutations are genetic variation, so a gene cannot be “left invariant by permutation”. Are you perhaps referring to variation in the proteins encoded by protein-coding genes or the downstream phenotypes?

You are correct, I mislabeled what I meant. Technically not mutations, but defects than don't affect the sequence, i.e.: $(ACT)_1(ACT)_2 -> (ACT)_2(ACT)_1$ the resulting sequence is the same, but there was a change of the 'physical' bases. It occurred to me that, given the possible mechanisms for error production, some sequences would be more stable, and therefore, less prone to mutation.

@ElkinMontoya - you can't distinguish the examples you gave, they are literally chemically the same thing. I suggest you look up a codon table and see the redundancy in the codons and think about how that affects your question.

@bob1 That is precisely the point, the point I'm trying to get across is that some sequences, like $CGA$ can, due to imperfect replication, translate into $AGC$ (the first base is transposed with the third), but others, like $CAC$, will still produce the same sequence even if the same physico-chemical events that made the other sequence change where to happen. And so, my reasoning goes, proteins that are encoded by sequences like $CAC$ will not mutate at the same rate as others Because even if something were to happen, it will still be the same sequence

@ElkinMontoya - I'm not sure that you can study mutations where the "error" doesn't actually produce an error, because it then isn't a mutation at all. You can examine the mutation rate for individual bases in a protein coding sequence though.

@bob1 Thanks for advice. I'll try that.

Have you looked up the chemical nature of mutations? Mixups like "swap the first base for the third" are not the most common form of mutation, not by a long shot.

Not in the way you suggest: a sequence of letters with the property of being immutable by any type of permutation on its constituent bases. I would be interested if anything of the likes exists in mathematics.
Having said that, there are examples of sequences that seem to have such property but only in the light of evolution and the biological processes that come with it:
Ultra-conserved regions have been described as sequences of about 200 base pairs that appear to be exactly the same from humans to mice. They have not experience change in about ~300 million years of evolution (!) and exposure to the various forces of mutation.
The reason for this is not an underlying 'invariant to permutation' property, but rather strong negative selection that prevents individuals that acquire changes on the sequence from being able to pass the changes over to their progeny, thus effectively maintaining only the original version in the population.
The exact mechanisms, functions and properties of these regions are under research. You might want to see if the invariant properties you have in mind have anything to do with these interesting regions.

Enzymes' catalytic site(s) are also very often left unchanged for the same reason. If there happen to be a mutation there, the protein's function might be compromised. But sometimes it is  instead enhanced....

Very interesting. Based on the criteria that their existence is due to very strong selective pressure, and not due to external mechanisms, I figure it is possible that, given they are relative few in number, the UCR sequences that are not completely essential could be something of this kind, were some to be actually inessential, but it is also possible that the UCR change with species, so it may actually take an enormous to prove/disprove.

Yes, I guess UCRs are specific to particular groups of species and/or have certain upper limits to the evolutionary time within which they remain conserved. The UCR shared by humans and mice for instance were present in their common ancestor, so therefore all species derived from this common ancestor have them too. As usual, proving a certain property in this or other UCR example might well lead to some degree of generalization,


text:
Selection analysis based on each leaves as foreground with FastCodeML
I am using positive selection analysis for each leaves in a phylogenetic tree.
In the routine calculation via codeml, I will be running model zero for all leaves and then model 2 for each leaves as foreground (marking leaves as #1 in phylogenetic tree), then calculate selection using a likelihood ratio test.
Is there a method for using FastCodeML for each leaves as foreground instead of the branch?

First and foremost, the fastcodeml guys didn't replace codeml via refactorising all the functionality of codeml, they refactorised one particular calculation and hence the confusion. fastcodeml isn't a particularly name helpful in truth, because it is not a faster version of all the functionality codeml, its just one of the calculations.
In the authors defence they (notably) improved the runtime of the most intensive calculation in codeml and thats their point - its proof of concept.
Thus, basically no fastcodeml can not be used in this context. fastcodeml runs the branch-site model of codeml only.
What the OP describes isn't how the branch-sites model works, they are describing the branch model.
The branch-site null model sets dN/dS to 1, and the LRT compares dN/dS likelihood between the designated branch(es) between dN/dS fixed to neutral against allowing dN/dS to be any value. Thus it's measuring whether the branch (branches) deviate from neutrality. What the OP describes via, model 0 one rate for the tree versus model 2 meaning two or more rates for the tree (basically two rates) and that is a different analysis and it is called a branch model.

If fastcodeml can be shoehorned to perform sites and branch-models that would be good to know.

There is an issue about 'tip selection' but the key answer concerns the basic difference between codeml and fastcodeml.


text:
How to convert Sars-CoV-2 ORF1ab codon positions to ORF1b positions?
Some websites mention amino acid mutations as ORF1ab, most however use ORF1a/b instead.
So how can I translate from one to the other? What do I need to subtract from say ORF1ab:4588I to get the corresponding N in ORF1b:<N>I?

ORF1ab:4588I corresponds to ORF1b:187I so one needs to subtract 4401 to go from 1ab to 1b and add 4401 to go the other way.

Just for clarification: 1a or 1b? You reference 1b in the question but answer with 1a.

You're right, I fixed it, thanks!

There is now an Open Source online tool that does the conversion, Codon2Nucleotide by Theo Sanderson.


text:
pysam "Exec format" error
I am a beginner and trying to read a bam file in Python.
The lines below throw the error
OSError: [Errno 8] could not open alignment file `mybam.bam`: Exec format error

What could be the possible solution?
user@server:~/Documents$ python3
Python 3.8.10 (default, Jun  2 2021, 10:49:15)
[GCC 9.4.0] on linux
Type "help", "copyright", "credits" or "license" for more information.

>>> import pysam
>>> samfile = pysam.AlignmentFile("mybam.bam", "rb")

[E::hts_hopen] Failed to open file mybam.bam
[E::hts_open_format] Failed to open file "mybam.bam" : Exec format error
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "pysam/libcalignmentfile.pyx", line 742, in pysam.libcalignmentfile.AlignmentFile.__cinit__
  File "pysam/libcalignmentfile.pyx", line 941, in pysam.libcalignmentfile.AlignmentFile._open
OSError: [Errno 8] could not open alignment file `mybam.bam`: Exec format error

this is something related to the file not with python. Check the file with the linux command file mybam.bam just to know that you have a bam instead of a sam, text format. If you have a compressed file in bam format, check the head with samtools view mybam.bam | head -100, if this fails check the mapping software again where you got the alignment from.

https://www.biostars.org/p/284417/

Thank you for the help. samtools view mybam.bam | head -100
outputed a hundred lines, which have no index. As there is no index, I was able to retrieve the content with:
for read in samfile.fetch(until_eof=True):
print(read)


text:
Tool that can sort NCBI taxonomy ids based on phylogeny
I have a set of NCBI txids (of varying ranks) that I want to sort in some way (e.g. a table or a tree, the way it is sorted isn't super important), based on the current NCBI data for phylogenetic relatedness.
I don't have so many txids that I can map phylogenetic relationships at a really high resolution, but if, say, the txids for plants, mammals, chordates, humans could be sorted so that humans are nested within mammals which is nested within chordates, which is separated from plants.
Is there a tool that can do this?

not sure that i understand. what do you mean by "sort"? if you just want a tree of your taxids, what is wrong with tools like this: https://www.ncbi.nlm.nih.gov/guide/howto/gen-com-tree/


text:
Why does it require four vectors to calculate a torsion angle?
The BioPython documentation shows the following:

11.6.3 Measuring torsion angles Use the vector representation of the atomic coordinates and the calc dihedral function from the Vector module:
>>> vector1 = atom1.get_vector()
>>> vector2 = atom2.get_vector()
>>> vector3 = atom3.get_vector()
>>> vector4 = atom4.get_vector()
>>> angle = calc_dihedral(vector1, vector2, vector3, vector4)

We know that a torsion angle between N-Cα and Cα-C planes.
So, the minimum number of points required to form a plane is 3.
Then, Why does BioPython require 4 vectors?


text:
Help with performing PCA analysis on data from 3 different datasets of normalized counts data for RNA-seq experiment
So, I have limited knowledge of R but I need to do a PCA analysis of 3 different datasets of gene expression as a result of combined growth or mono-culture growth.
The 3 different datasets I performed DESeq2 analysis on are as follows:
Dataset 1: This was sequenced by Ion Torrent and it was single-end reads. Samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown in direct contact with yeast species B. This was all done for 1 time point @ 24h.
Dataset 2: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.
Samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown in direct contact with yeast species B. This dataset contains 2 time points @ 2h and 24h.
For the other part of this, samples that were compared were biological quadruplicates of yeast species A grown as pure culture versus species A grown via indirect (metabolic exchange possible) contact with yeast species B. This dataset contains 2 time points @ 2h and 24h.
Dataset 3: This was sequenced by Illumina MiSeq and paired-end reads were generated, but this dataset had two parts to it.
Here, two biological repeats of samples of yeast species A grown as pure culture versus two biological repeats of species A grown in direct contact with yeast species B, under AEROBIC conditions were performed. This dataset contains 1 time point @ 48h.
Then one biological repeat of sample of yeast species A grown as pure culture versus two biological repeats of species A grown in direct contact with yeast species B, under ANAEROBIC conditions were performed. This dataset contains 1 time point @ 48h.
I wrote the results out of the normalized counts per sample combined with the DESeq2 results as follows:
res <- results(dds)
table(res$padj<0.05)
res <- res[order(res$padj), ]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)), by="row.names", sort=FALSE)
names(resdata)[1] <- "Gene"
head(resdata)
write.csv(resdata, file="diffexpr-results.csv")

So the issue I am having is, because these are different datasets and conditions differ across them there are different numbers of genes expressed in each dataset that are not present across all datasets. So I need to combine all the datasets to only compare the genes found in all the datasets mentioned and then perform a PCA analysis on those genes only.
An example of the setup of the excel files I have, without the other info that you get with DESeq2 analysis e.g. columns baseMean, log2FoldChange, lfcSE, stat, pvalue and padj, is as follows:

Gene_ID
Sample1_PC
Sample2_PC
Sample3_PC
Sample4_PC
Sample5_Mixed
Sample6_Mixed
Sample7_Mixed

YL09980
1240,016
759,683
985,322
1025,855
5002,654
4200,963
4002,855

YLW9980
589,456
452,236
630,228
523,821
2325,556
2566,987
1800,633

So these tables are just the normalized counts of the genes per sample.
I was told that the data in the table has to be transposed basically, so that the genes are the variables and the samples are the rows. I am just not sure how to do this in R.
I was also told there is a good chance that the noise resulting from the use of different platforms would likely cause these datasets to separate from each other based on that, making this analysis meaningless or not allowing for comparison of gene expression that is conserved across the datasets as a result of mixing the two species during growth.
However, I just need to do this analysis and find that out as soon as possible, as this is a bottleneck to my progress at this point. Please any assistance would be greatly appreciated. I am currently teaching myself more about "tidyverse" in an effort to try and do this, but I am on a deadline with this.
Edit: So I did the PCA plot using the following function:
data_total.pca <- prcomp(data_total[,c(2:32)], center = TRUE, scale. = TRUE)

and this was the output I got:
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7     PC8     PC9    PC10    PC11    PC12    PC13    PC14
Standard deviation     4.2180 2.5400 1.9188 1.2866 0.65393 0.47799 0.40826 0.38879 0.35585 0.34915 0.27036 0.18549 0.14335 0.11830
Proportion of Variance 0.5739 0.2081 0.1188 0.0534 0.01379 0.00737 0.00538 0.00488 0.00408 0.00393 0.00236 0.00111 0.00066 0.00045
Cumulative Proportion  0.5739 0.7820 0.9008 0.9542 0.96800 0.97537 0.98074 0.98562 0.98971 0.99364 0.99600 0.99711 0.99777 0.99822
                          PC15    PC16    PC17    PC18    PC19    PC20    PC21    PC22    PC23    PC24    PC25    PC26    PC27    PC28
Standard deviation     0.10153 0.09239 0.08407 0.07765 0.07536 0.05978 0.04883 0.04556 0.04406 0.03960 0.03764 0.03378 0.03167 0.02897
Proportion of Variance 0.00033 0.00028 0.00023 0.00019 0.00018 0.00012 0.00008 0.00007 0.00006 0.00005 0.00005 0.00004 0.00003 0.00003
Cumulative Proportion  0.99855 0.99883 0.99906 0.99925 0.99943 0.99955 0.99963 0.99969 0.99976 0.99981 0.99985 0.99989 0.99992 0.99995
                          PC29    PC30    PC31
Standard deviation     0.02695 0.02126 0.02085
Proportion of Variance 0.00002 0.00001 0.00001
Cumulative Proportion  0.99997 0.99999 1.00000

So according to the tutorial I am following a general rule of thumb is that the first 3 PCs should explain at least 70% of the variance. The issue is my PCA plot looks weird all of the points are squished up at the x = 0; y = 0 axes of my PCA plot. So I think I need to transpose the data so that the genes are the columns (variables) and the samples are the rows (objects). So I tried this:
data_total_t <- transpose(data_total)
rownames(data_total_t) <- colnames(data_total)
colnames(data_total_t) <- rownames(data_total)

And this worked, but now I need to exclude row 1 from the function I used to plot the PCA as the PCA cannot have character (gene name) info in it.
Any help on how I can do this please?

I think you are making this harder than you have to.  Why can't you follow the code used in the vignette for making PCA plots?  It doesn't even involve getting all the way to the DE genes.

Hi there @swbarnes2, I found a tutorial online to follow but now I am uncertain of the plot I am getting in terms of if it is correct or not. All of my data is now squished into the x = 0; y = 0 point of the PCA axes. Is this because I should rather be plotting the genes as columns and the sample names should make up the rows?

IF you are already using DESeq, why not just do what they do in the vignette?

I will give it a look thanks @swbarnes2, I have also generated PCA plots already now for the data but I am not sure if they are very meaningful or if the data can actually be compared in this way.

So I need to combine all the datasets to only compare the genes found in all the datasets mentioned and then perform a PCA analysis on those genes only.

You can merge all your datasets for common genes in all datasets and perform PCA.

I was told that the data in the table has to be transposed basically, so that the genes are the variables and the samples are the rows. I am just not sure how to do this in R.

Your choice of transposing data will depend on whether you want to show variability between genes or samples. However you can get results for both individuals (genes) as well as variables (samples), using the prcomp(), factoMineR and factoexta packages in 'R'. Usually I like to use $rotation values to plot the variability between samples using ggplot2 package. You can use this PCA Tutorial to learn more.

I was also told there is a good chance that the noise resulting from the use of different platforms would likely cause these datasets to separate from each other based on that, making this analysis meaningless or not allowing for comparison of gene expression that is conserved across the datasets as a result of mixing the two species during growth.

This is true as the datasets are from different platforms, there may be some noise or batch effects. However, you can plot the 'PCA' and decide for yourself whether it is good or not.


text:
R help in Jupyter notebook?
I'm using this Jupyter notebook to work on R... How do I do the equivalent of ?someFun?
When I run that in the notebook, I get:
Warning message in seq_len(head.end.idx):
“first element used of 'length.out' argument”
ERROR while rich displaying an object: Error in seq_len(head.end.idx): argument must be coercible to non-negative integer

Traceback:
1. FUN(X[[i]], ...)
...

Update
When I execute docker run jupyter/r-notebook I notice that it does pop up the widget (after throwing the error). Not it isn't a popup window, but a panel at the bottom of the workbook. My notebook is 'Trusted', which I think is significant, but I don't know how to turn it off.
When I execute docker run -e JUPYTER_ENABLE_LAB=yes jupyter/r-notebook I get the error and the widget does not pop up. I'm not sure how to set or inspect the JS trust in that environment. I suspect it's related. I don't see any blocked popups in the browser.
Strange that even when it works (the widget opens at the bottom of the browser) I still get an error in the output.

If you launch like docker run -p 8888:8888 jupyter/r-notebook you will get the help page like a widget and a warning in your cell output. Check your browser, maybe is blocking pop-up elements.

Thanks for showing me what the widget looks like. I do see this now (i.e. working correctly), but only when not setting JUPYTER_ENABLE_LAB=yes

I met the same problem. R after v4.1.1 seem to have problems with repr_html.
I fixed the problem by following the quotes.

This tiny commit solves a long standing problem while trying to read help from inside jupyterlab.
In order not to wait for the PR to be merged, you may just install from @renkun-ken 's fork:
devtools::install_github("renkun-ken/repr", ref = "fix-html-filter")

Reference
https://github.com/IRkernel/repr/pull/144#event-5589475824
https://github.com/IRkernel/repr/issues/143

Thanks. I mentioned this here:
https://github.com/jupyter/docker-stacks/issues/1494

The PR is now merged and installing repr==1.1.4 resolved the issue for me.

I think error is not related to ?. This is what I see when I try ?seq:

Yes, help(seq) and ?seq are both a problem. I think it's related to trusted JS execution... I think...


text:
How to simulate replicates for DGE analysis?
I am prototyping with data visualization of DGE results, and I would like to work on the analysis pipeline before the real data is available. Currently, I only have 3 samples for  wild type and 1 sample for a mutant. I would like to simulate 2 extras maples of the mutant with the same variance as the wild type sample. How can I simulate data to make up for the absence of replicates in the comparison group?
My idea was to use pandas and SciPy to import the count tables, use these tools to simulate the replicates and then run DESeq2 on it. Still, I am not sure which functions I should look for once I am not so proficient with SciPy and pandas yet.


text:
Installing snakemake following official docs just won't work. Help needed
Following the docs just won't work for me.
conda install -n base -c conda-forge mamba
conda activate base
mamba create -c conda-forge -c bioconda -n snakemake snakemake

This is the Dockerfile I use:
FROM ubuntu:20.04
ENV PATH /opt/conda/bin:${PATH}
ENV LANG C.UTF-8
ENV SHELL /bin/bash
COPY aws/docker/Miniconda3-py38_4.10.3-Linux-x86_64.sh /tmp/
RUN bash /tmp/Miniconda3-py38_4.10.3-Linux-x86_64.sh -b -p /opt/conda
RUN conda install -y -c conda-forge mamba
RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps

The error is:
WARNING Extracted package cache '/opt/conda/pkgs/_openmp_mutex-4.5-1_gnu' has invalid size
WARNING Extracted package cache '/opt/conda/pkgs/_openmp_mutex-4.5-1_gnu' has invalid SHA-256 checksum
WARNING Package tarball '/opt/conda/pkgs/_openmp_mutex-4.5-1_gnu.tar.bz2' is invalid
  Package                              Version  Build                Channel                     Size
───────────────────────────────────────────────────────────────────────────────────────────────────────
  Install:
───────────────────────────────────────────────────────────────────────────────────────────────────────

  ...Package list omitted

  Ignored:
───────────────────────────────────────────────────────────────────────────────────────────────────────

  = snakemake                            6.9.1  hdfd78af_0           bioconda/noarch          Ignored
  = snakemake-minimal                    6.9.1  pyhdfd78af_0         bioconda/noarch          Ignored

  Summary:

  Install: 215 packages
  Ignored: 2 packages

  Total download: 256 MB

───────────────────────────────────────────────────────────────────────────────────────────────────────

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/gettext-0.19.8.1-h73d1719_1008.tar.bz2"
        Expected: 1bb53c99b4943d210c881aad9158fb0235b348498bad1a7076d1f2bef6671922
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/gtk2-2.24.33-h539f30e_1.tar.bz2"
        Expected: 2a0aef09a947e3768a426b0ea414ddd1e015c9804a46de619a5e67a689ea1dec
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/botocore-1.21.64-pyhd8ed1ab_0.tar.bz2"
        Expected: ee5501c74eefd0396f141593933ef2cc3419ad21624ed140bbc1e04944af6e31
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/libopenblas-0.3.18-pthreads_h8fe5266_0.tar.bz2"
        Expected: 6d67f0f7bda1503af83659ca45d8c8ba34d94091fa4cd5b1ecd85332c14db3cc
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

ERROR   File not valid: SHA256 sum doesn't match expectation "/opt/conda/pkgs/scipy-1.7.1-py39hee8e79c_0.tar.bz2"
        Expected: 536500b79a5ccc91f6059fcaa6d7e959bddb37a0ed7a930797375c1ed51a6af1
        Actual: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855

# >>>>>>>>>>>>>>>>>>>>>> ERROR REPORT <<<<<<<<<<<<<<<<<<<<<<

    Traceback (most recent call last):
      File "/opt/conda/lib/python3.8/site-packages/conda/exceptions.py", line 1079, in __call__
        return func(*args, **kwargs)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 926, in exception_converter
        raise e
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 920, in exception_converter
        exit_code = _wrapped_main(*args, **kwargs)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 878, in _wrapped_main
        result = do_call(args, p)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 757, in do_call
        exit_code = create(args, parser)
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 660, in create
        return install(args, parser, "create")
      File "/opt/conda/lib/python3.8/site-packages/mamba/mamba.py", line 609, in install
        downloaded = transaction.prompt(repos)
    RuntimeError: Found incorrect download: libopenblas. Aborting

`$ /opt/conda/bin/mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps`

  environment variables:
                 CIO_TEST=<not set>
               CONDA_ROOT=/opt/conda
           CURL_CA_BUNDLE=<not set>
                     PATH=/opt/conda/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin
                          :/bin
       REQUESTS_CA_BUNDLE=<not set>
            SSL_CERT_FILE=<not set>

     active environment : None
       user config file : /root/.condarc
 populated config files :
          conda version : 4.10.3
    conda-build version : not installed
         python version : 3.8.10.final.0
       virtual packages : __linux=5.10.16.3=0
                          __glibc=2.31=0
                          __unix=0=0
                          __archspec=1=x86_64
       base environment : /opt/conda  (writable)
      conda av data dir : /opt/conda/etc/conda
  conda av metadata url : None
           channel URLs : https://conda.anaconda.org/conda-forge/linux-64
                          https://conda.anaconda.org/conda-forge/noarch
                          https://conda.anaconda.org/bioconda/linux-64
                          https://conda.anaconda.org/bioconda/noarch
                          https://repo.anaconda.com/pkgs/main/linux-64
                          https://repo.anaconda.com/pkgs/main/noarch
                          https://repo.anaconda.com/pkgs/r/linux-64
                          https://repo.anaconda.com/pkgs/r/noarch
          package cache : /opt/conda/pkgs
                          /root/.conda/pkgs
       envs directories : /opt/conda/envs
                          /root/.conda/envs
               platform : linux-64
             user-agent : conda/4.10.3 requests/2.25.1 CPython/3.8.10 Linux/5.10.16.3-microsoft-standard-WSL2 ubuntu/20.04.3 glibc/2.31
                UID:GID : 0:0
             netrc file : None
           offline mode : False

An unexpected error has occurred. Conda has prepared the above report.

Curious that the checksums are wrong. Is there something wonky with the network that it might be changing the packages? What happens if you manually download those packages and look at their checksums?

Yes I have found similar comment on Github page I think. Somehow it is hard for me to blame the network. I thought that I am missing something. When I manually download the file. SHA256 matches the expected one.

What happens if you run mamba outside of a container?

Yes, on fresh Ubuntu 20.04 EC2 machine it installs just fine. On Docker packages are randomly failing SHA256 checksums. Interesting fact is that I am on Windows WSL (company policy) and that thing gave just trouble.

Weird, I wonder if something is messed up in the docker networking interface. I never tried to set docker up under wsl, I just start a new wsl OS if needed (I have similar company restrictions).

I've ran a build using our Gitlab runner. Everything works just fine.

I select the Miniconda installer directly from the official repository:
FROM ubuntu:20.04
ENV LANG C.UTF-8
SHELL ["/bin/bash", "-c"]
RUN apt update -qq && apt -qq -y install curl ca-certificates --no-install-recommends
RUN ( cd /tmp/; curl -LO "http://repo.anaconda.com/miniconda/Miniconda3-py38_4.10.3-Linux-x86_64.sh")
RUN bash /tmp/Miniconda3-py38_4.10.3-Linux-x86_64.sh -b -p /opt/conda
ENV PATH=/opt/conda/bin:${PATH}
RUN conda install -y -c conda-forge mamba
RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps
RUN echo "source activate snakemake" > ~/.bashrc
ENV PATH /opt/conda/envs/snakemake/bin:${PATH}

I agree with this AWS was being used as the source of miniconda rather than a direct download from conda. The alternative is simply to download the official conda docker

Try
FROM ubuntu:latest
MAINTAINER Conda Development Team <conda@continuum.io>

RUN apt-get -qq update && apt-get -qq -y install curl bzip2 \
    && curl -sSL https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -o /tmp/miniconda.sh \
    && bash /tmp/miniconda.sh -bfp /usr/local \
    && rm -rf /tmp/miniconda.sh \
    && conda install -y python=3 \
    && conda update conda \
    && apt-get -qq -y remove curl bzip2 \
    && apt-get -qq -y autoremove \
    && apt-get autoclean \
    && rm -rf /var/lib/apt/lists/* /var/log/dpkg.log \
    && conda clean --all --yes
RUN conda install -y -c conda-forge mamba
RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake snakemake-minimal --only-deps
ENV PATH /opt/conda/envs/snakemake/bin:${PATH}

The && approach is better as it will minimise the container size (which is going to be big in any case). Then look at an alpine version once its working.
To be honest, I'm not sure about not activating conda base and personally I'd look to correct that.
Alternatively ...
docker pull continuumio/miniconda3
docker run -i -t continuumio/miniconda3 /bin/bash

Just for a chuckle, a cool approach is
docker pull frolvlad/alpine-miniconda3

This is super small.
The docker file being
FROM frolvlad/alpine-glibc:alpine-3.14

ARG CONDA_VERSION="4.7.12.1"
ARG CONDA_MD5="81c773ff87af5cfac79ab862942ab6b3"
ARG CONDA_DIR="/opt/conda"

ENV PATH="$CONDA_DIR/bin:$PATH"
ENV PYTHONDONTWRITEBYTECODE=1

# Install conda
RUN echo "**** install dev packages ****" && \
    apk add --no-cache --virtual .build-dependencies bash ca-certificates wget && \
    \
    echo "**** get Miniconda ****" && \
    mkdir -p "$CONDA_DIR" && \
    wget "http://repo.continuum.io/miniconda/Miniconda3-${CONDA_VERSION}-Linux-x86_64.sh" -O miniconda.sh && \
    echo "$CONDA_MD5  miniconda.sh" | md5sum -c && \
    \
    echo "**** install Miniconda ****" && \
    bash miniconda.sh -f -b -p "$CONDA_DIR" && \
    echo "export PATH=$CONDA_DIR/bin:\\\$PATH" > /etc/profile.d/conda.sh && \
    \
    echo "**** setup Miniconda ****" && \
    conda update --all --yes && \
    conda config --set auto_update_conda False && \
    \
    echo "**** cleanup ****" && \
    apk del --purge .build-dependencies && \
    rm -f miniconda.sh && \
    conda clean --all --force-pkgs-dirs --yes && \
    find "$CONDA_DIR" -follow -type f \( -iname '*.a' -o -iname '*.pyc' -o -iname '*.js.map' \) -delete && \
    \
    echo "**** finalize ****" && \
    mkdir -p "$CONDA_DIR/locks" && \
    chmod 777 "$CONDA_DIR/locks"

Whoever frolvlad is they certainly know their stuff, its a highly optimised image. This will absolutely nail a light-weight conda file. How well conda will work with Snakemake  ... it is  well worth trying because the alpine image is 5meg (Ubuntu is 100+ Megs) and I forget how big conda is. I need to do this more for sure, because almost certainly what they have done is streamlined conda to the absolute essentials. Get that working your looping the loop.
Why is this potentially super cool? The problem with alpine is its so minimalist getting anything to work is not easy. frolvlad in theory solved this bit via conda, because conda is completely self-contained once its running. Thus its doesn't need any system stuff at all, hence alpine and conda in THEORY would be a fantastic solution.
My honest opinion is you are sort of trying to reinvent the wheel  (which sometimes that's a good thing). However, I think you'd find the other stuff just works because someone has sat down with yaml and tested everything endlessly.
The problem with your file is you've got to clean it without additional layers being dumped onto it and that's why loads of && because that enables a clean without a layer collapsing onto it blocking the clean. The docker file you create has loads of unnecessary stuff, which you can't delete. I guess your priority is to get it working.


text:
general cheminformatics references
Can anyone explain what is cheminformatics and how does it help in the future?
I need to write an article that should be understandable for 14-20 yr age group. Please help with references, I don't know anything.

Hi Harsha, is this a homework of some sort? What's the context of figuring out what chemoinformatics is?

Hi, this is about writing an article on cheminformatics for people 14-20 yr age group should understand, I know some basics but can't figure what should I include like subheading and context

The definition of cheminformatics has been covered in Difference between bioinformatics and cheminformatics

For a reference, simply citing Essentials of Computational Chemistry by C. Cramer will do it. I do not know if the newer Clayden (nitrile cover) deals with computational methods, but the old one (water+formaldehyde on cover) did not. In generally, a textbook is a very good way to get familiar with a field and Cramer's book is good.


text:
In CIRIquant output, which column shows the expression level of circular?
I am running CIRIquant on Ubuntu 16.04 LTS following the instruction which indicated on https://ciriquant-cookbook.readthedocs.io/en/latest/de.html. I got .GTF file that contains detailed information of BSJ and FSJ reads of circRNAs and annotation of circRNA back-spliced regions.
In the following, I added a record of Circular RNA detected in .GTF file:
Chrom: chr1

Type: circRNA

start: 1289228

End:  1289448

score:0.02

Strand Information: -

circ_id: chr1:1289228|1289448

circ_type   :exon

bsj: 1000

fsj: 515.000

junc_ratio: 0.004

gene_id: ENSG00000162576.12

gene_name: MXRA8

gene_type: protein_coding

I need to know which column of .GTF file indicates the expression count of circular detected?
I appreciate it if you share your comment with me.
Best Regards,