text:
GWAS phenotype data format and preprocessing
I have a set of different phenotypes which I want to use for a GWAS analysis (general linear model). I have a couple of questions and uncertainty about the phenotype data input.

I have control and samples and usually we normalize all sample to the control for other calculations.

Example:
Mean control - 60%
Mean sample1 - 20 %
Mean sample1 normalized to control - 20*100/60 = 33,3 %
Is this normalization step correct for a GWAS analysis?

Should the phenotype data scaled into 0...1 range? I have different scales for different phenotypes but the results should be compared and compatible.

Some of my phenotypes are not normality distributed. What is the best way to handle such data? I have read that the data should be normally distributed for GWAS analysis.

Here are the histograms of two phenotypes:

Please could you clarify if by "GWAS analysis" you are referring to general linear modelling (GLM)? This is regression analysis assessed using regression weights and a residual to understand the association between alleles and phenotypes. Is that correct?

Normalisation is standard procedure for GLM

Yes it's general linear modelling (written in the first sentence).

Okay you need to do this via a Q-Q plot, or at least it is one way to do this. The an alternative is to minimise the residual. When I get a chance I'll post an explanation. You don't necessarily need to transform, it depends on the 'diagnostics', e.g. Q-Q.


text:
a proper Design Matrix for several drug treatments with both control negative and control positive
I have a dataset of RNA-seq samples for testing different drugs on the presence of another drug. One of my samples is the normal cells with no drugs (control negative) and another is the cells with just the base drug (control positive). All of the other samples are cells with the base drug and another testing drug. So, my question is how to do design a proper Design Matrix for comparing these samples. I can separately compare with control negative and control positive. But, my intention is to somehow neutralize the effect of one of the controls in all the samples and compare to the other one. Is that possible? How should I design the matrix?

Literally one undrugged control?

Yes, Literally! Unfortunately I don't have biological replicates for any of these samples.

That's unfortunate.  I guess you can do some ratios in excel, but you can't do anything sophisticated without replicates.  And your conclusions will be really shaky.  Wet lab people are going to have to test your conclusions with better experiments.

I know that. But the whole experiment's purpose was to generate some hypotheses for further testing with replicates. So, I guess it can give me some insights!

Or it might give who a whole lot of false positives that waste time and money.  Everything  hinges on your one replicate being perfectly representative of what the control should look like.

I totally agree and discussed this issue a lot in the lab. But people still come to me with their no replicate datasets and I have to give them some results or insight! What can I do?

Give them a big 'told you so' 2 years down the line :-)

Are there any technical replicates? Could you share a table or list of the amount of samples and their conditions?

Condition is as I described. One sample is Mock, another sample is Drug1 and all the other samples are Drug1 + Drug2-7. Deug2-7 means 6 different drugs which are being combined with drug1 separately to see their simultanous effect. What do you mean by the amount of samples? Like number of Reads?

source: Stack Exchange

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

sort: votes


text:
Removing Batch Effect in Heatmaps after Differential Gene Expression Analysis
I'm working on a dataset in which the first replicate of each group is one batch and the second replicate is in a second batch. After checking the PCA plot and seeing the batch effect in PC1, I used removeBatchEffect function from limma to subtract the batch effect from my count data. Then, using PCA gives me a plot that doesn't seem to have any apparent batch effect left! However, it is recommended not to use batch effect correction for differential gene analysis but use the batch variable along with the group variable in constructing the model.matrix. So, I did that, ran limma/voom on the normalized counts, and extracted differentially expressed genes. However, when I'm trying to make a heatmap from the DEGs, I still see that the samples from different batches are clustered separately, instead of seeing the clustering of replicates of the same sample. So, my question is whether I should use removeBatchEffect on the count data from DEGs and then use the transformed dataset for heatmaps or there is another way to fix this?

Maybe you could try more state-of-the-art of methods like "comBat" from "sva" package to remove the batch effects.

I read a few posts in biostar that confirmed my choice of removeBatchEffect from limma over comBat from sva. The gist of it was that we should not use tools like comBat because they are biased towards removing all batch effects. Moreover, comBat and also removeBatchEffect generate negative sign results which can't be processed by voom or edgeR and trying to remove these negative results is another level of manipulation in the dataset. So, it is better to use removeBatchEffect for visualizations but just use batch in model.matrix for differential gene expression analysis

Read this for more information.

It is true that comBat introduces negative values. This is because comBat should be used on a normalized matrix. If all values in the matrix is positive, I am assuming it is a count matrix (which a newer version of sva takes care of) or a TPM matrix (which should be normalized first.) I think maybe you could run comBat after voom if you decide to use that. But of course using limma in a right way could work as well.

Could you please clarify your comment a bit. my data is already log2 scaled, filtered for lowly expressed genes and also TMM normalized, but I still get negative values with comBat and also removeBatchEffect. However, I'm using batch as a covariat in my model.matrix and do the DEG analysis. You meant I can use comBat on DEGs for visualizations?

If your data is log2 scaled, are there negative values before removing any batch effects?

I just checked and Yes there are a few negative values before removing batch effects.

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


text:
Find all the bases for given reference position
Hi im looking for the aligned bases in the reads for a given reference position. im using the following script from the pysam documentataion.
I adjusted it to find the specified position. in this case 24793.
import pysam

samfile = pysam.AlignmentFile("generated_alignment_sorted.bam", "rb" )
for pileupcolumn in samfile.pileup("chr05_modified.copy0", 10, 52000000):
    if pileupcolumn.pos == 24793:
        print ("\ncoverage at base %s = %s" %
               (pileupcolumn.pos, pileupcolumn.n))
        for pileupread in pileupcolumn.pileups:
                print ('\tbase in read %s = %s' %
                      (pileupread.alignment.query_name.split(';')[0],
                       pileupread.alignment.query_sequence[pileupread.query_position]))

samfile.close()

This outputs:
coverage at base 24793 = 6
        base in read m419941/6207/CCS Read=419941 = C

For all the position I tried it doesnt print out the same number of bases as it says that the coverage is.
shouldnt this print out 6 times: base in read m419941/6207/CCS Read=419941 = C for different reads?
Why is this?

pileupcolumn.n shows the total number of the reads that cover this basepair (in this case 24793). However, pileupcolumn.pileups iterates over the reads that cover this base pair AND this base pair has a minimum quality in those reads. So it has an additional filter for the reads to be considered. You can modify the quality threshold by setting min_base_quality parameter when you call pileup.
To have a similar filter applied to both coverage and reads, you can use piliecolumn.get_num_aligned() to get the number of segments and pileupcolumn.pileups to get the segments themselves.


text:
Modifying protein_scale in ProteinAnalysis in SeqUtil in Biopython?
Is there a way to give weight arguments to each of the residues when using protein_scale in the ProteinAnalysis module of Bio.SeqUtils.ProtParam? I thought this would be done with the window and edge arguments, but from the documentation and the question How to allow "half windows" in biopython? I understand that the weights are symmetric and calculated automatically given a window size.
I'd like to understand what the idea behind this was, just to check that I understand this correctly. For example, if one calculates eg. hydrophobicity of a protein, why should it be a value per residue that depends on a particular window?
I was looking for a functionality like where the value given for each residue in the ParamData dictionary would be multiplied by a specified weight, like so:
```Y = ProteinAnalysis('MSTAEFAQLL')
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10)
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,edge[0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.])```

Of course I can code this myself, but was wondering if that functionality already existed.

There is not a way to control the edge weights as written, but making a modified version is straightforward. Namely make a class that inherits the class Bio.SeqUtils.ProtParam.ProteinAnalysis and overwrites the method relevant methods.
Having a gander at the files one sees the following:
def protein_scale(self, param_dict, window, edge_list: list):
    """[...]
    """
    # generate the weights
    #   _weight_list returns only one tail. If the list should be
    #   [0.4,0.7,1.0,0.7,0.4] what you actually get from _weights_list
    #   is [0.4,0.7]. The correct calculation is done in the loop.
    weights = self._weight_list(window, edge)
    scores = []
    ...

So the weights are generated by the 'private' method ._weight_list. Underscore means it is for internal use and for production code it is bad to assume that the name of the private method many not change, but Biopython has not changed in years so you'd be safe. And this is for a quick hack anyway...
class JankyProteinAnalysis(Bio.SeqUtils.ProtParam.ProteinAnalysis):
    overriding_weights = None

    def _weight_list(self, window, edge):
          if self.overriding_weights is not None: # dodgy override!
               print("Weight method hack done!")
               return overriding_weights
          else: # behave normally
               return super()._weight_list(window, edge)

So you'd do:
Y = JankyProteinAnalysis('MSTAEFAQLL')
Y.overriding_weights = [0,0.2,0.3,0.1,0.,0.4,0.5,0.8,0.8,1.]
Y.protein_scale(Bio.SeqUtils.ProtParamData.kd,10,float('nan'))

Note that the edge argument is no longer doing anything.

Thanks so much! Super helpful.


text:
Extract residue sequence from pdb file (in biopython, but open to recommendations)
I'm new to Biopython and I'd like to extract the sequence of residues from a pdb file.
My two questions are:

What is the simplest way to do this? (Esp. when there is more than one sequence)
and
Should I be worried about "PDBConstructionWarning: WARNING: Chain B is discontinuous"?

So far, I've obtained the residue sequence via:
p = PDBParser()
structure = p.get_structure("1ppi", "1ppi.pdb")
ppb=PPBuilder()
for pp in ppb.build_peptides(structure):
     print(pp.get_sequence())

seq = pp.get_sequence().__str__()

This seems to work well for this molecule. However, is there an easier way, especially when there is more than one sequence?
For example, I've read that one can also do
res_list = Bio.PDB.Selection.unfold_entities(structure, 'R')

but res_list is not a sequence of residues in str, and I don't know how to convert the output from res_list to a str sequence.
In addition (or perhaps because I'm going through the PPBuilder), I've recently gotten a lot of warnings of the type:
/usr/local/lib/python3.8/site-packages/Bio/PDB/StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line..
(For example, with structure = p.get_structure("5owu", "5owu.pdb"))
I've seen a discussion about how to silence such warnings, but should I be worried about this? I've also noticed that when I get these warnings, pp builder seems to give me more sequences then are there.

One thing to consider is whether you want the full sequence of the protein (equivalent to the SEQRES record) or the sequence of the residues with coordinate data in the file (from parsing the ATOM records). They often differ, especially for proteins with disordered regions, and the answers below cover both methods.

There are other ways to do it (many and some might be easier), but if you want to do it with BioPython PDB module you can iterate the residues like this:
 # You can use a dict to convert three letter code to one letter code
 d3to1 = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
 'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
 'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
 'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

 # Just an example input pdb
 record = '1pa2.pdb'

 # run parser
 parser = PDBParser(QUIET=True)
 structure = parser.get_structure('struct', record)

 # iterate each model, chain, and residue
 # printing out the sequence for each chain

 for model in structure:
     for chain in model:
         seq = []
         for residue in chain:
             seq.append(d3to1[residue.resname])
         print('>some_header\n',''.join(seq))

Not really, is very common, you can turn it off like this (taken from https://lists.open-bio.org/pipermail/biopython/2014-July/015371.html):
 import warnings
 from Bio import PDBConstructionWarning
 #your code
 with warnings.catch_warnings():
 warnings.simplefilter('ignore', PDBConstructionWarning)
 #your code which might trigger the warning
 #rest of your code here

 However, because this is so common, you can just use the
 QUIET=True option on the PDBParser object:

 from Bio.PDB.PDBParser import PDBParser
 struct = PDBParser(QUIET=True).get_structure('tmp', pdbf)

Try help(PDBParser) for more, or see:
http://biopython.org/DIST/docs/api/Bio.PDB.PDBParser%27.PDBParser-class.html#__init__

The discontinuity is normal: it is a missing stretch due to lack of density. In PyMOL you see them as dotted lines.
I am sorry to say, the easiest/safest way to get the sequences is not using the PDB files...
import requests
data = requests.get(f'https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/{code}').json()[code.lower()]
print(data[0]['sequence'])

One gets back a lot of information: example. This will tell you when you have multiple copies of the same peptide. The only information you do not get is what the PDB numbering is relative to the canonical Uniprot for that chain, which has to be taken from SIFTS.
A side note... Biopython.PDB is one of the best BioPython packages, but I personally I much prefer PyMOL as a python module (pymol2 module) to Biopython.PDB as it is functionally more complete and straightforward even if everything operates in a stage. You might want to check it out.
with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1UBQ', 'prot')
    print(pymol.cmd.get_fastastr('prot'))

Thank you! I've accepted the other answer because of the point with how to iterate over the residues, but the info about the density/unitpro were very helpful, and I really appreciate the pointer to pymol!

You can also use the SeqIO module from Biopython https://biopython.org/docs/1.75/api/Bio.SeqIO.PdbIO.html
and do something like this:
from Bio import SeqIO

for record in SeqIO.parse(target, "pdb-atom"):
    print(record.seq)

where target is the path to your pdb file
the parser has multiple input/output functions, and pdb-atom is the one that will read out of the coordinates itself.

you can try PDBx Python Parser too:
http://mmcif.wwpdb.org/docs/sw-examples/python/html/index.html
it has an example showing exactly what you need
http://mmcif.wwpdb.org/docs/sw-examples/python/html/fasta.html


text:
How to represent random sequence elements in SBOL?
In our lab, we work with synthetic biology components with partially random sequences (similar to work in directed evolution). So, for example, we have a plasmid design with several components, including one protein coding sequence that will be created from a wild-type sequence via mutagenic PCR.
At the design stage, we don't know the exact DNA sequence for any of the resulting variants of that protein coding sequence. We just know the wild-type sequence and the average number of mutations we'd like to apply.
How can we represent that design in SBOL (via the Component class)?
Then, when we implement the design, we end up with a single sample that contains plasmids with 10^5-10^6 (or more) different protein coding sequences. We then measure the complete DNA sequence for every plasmid variant, but I don't want to end up with >10^5 different SBOL definitions to describe the sample.
How can we best represent that sample in SBOL (via the Implementation class)?
Finally, we also end up measuring the performance of every plasmid variant, and we identify the specific protein coding sequences for useful and/or interesting variants. We then design and build specific, clonal plasmids with those useful/interesting protein coding sequences.
Representing those clonal plasmids as SBOL Components seems pretty straight-forward, but is there a good way to capture the relationships with the original random-library plasmid design and the specific implementation of that design that we built and measured?

This sounds like a job for the CombinatorialDerivation class, which is designed for representing libraries and the relationship between derived designs and the originals.
About CombinatorialDerivation
In SBOL 3, a CombinatorialDerivation takes a Component in its template property to represent the system that you intend to vary:

Some of the Features of the template Component are designated as the variables of the library via VariableFeature child objects of the CombinatorialDerivation, which point at the variables with their variable property.
The values that can be filled in for each variable are then designated via the variant* properties of the VariableFeature.
Finally, the distribution of values is determined by the strategy property of the CombinatorialDerivation and the cardinality properties of each variable.

Representing Random Mutagenesis
In the case of random mutagenesis, if I'm understanding your design correctly, you would represent it with a CombinatorialDerivation as follows:

The template is Component of type DNA with the original sequence. Each site designated for mutagenesis can be designated with a SequenceFeature.
For each of the mutagenesis sites, create a VariableFeature in the CombinatorialDerivation with the corresponding SequenceFeature as its variable and a set of variant values for each of the mutant options set by your primers. The cardinality is probably one (if every site has to be hit by something in your protocol) or zeroOrOne (if sites might or might not be hit).
The strategy for the CombinatorialDerivation is sample, since only a few of the many possibilities will actually be realized.

Representing Samples and Characterized Constructs
Now, when you take it to the lab, you can represent your aliquot with an Implementation whose wasDerivedFrom points to the CombinatorialDerivation, representing the fact that you've got a randomly generated library whose contents are not precisely known.
Finally, once you've got your isolated products represented as Components, there is a recommended best practice to use wasDerivedFrom links to connect the derived Component and its varied Features back to the CombinatorialDerivation and its template. If you make a Collection of your results, the Collection should also link to the CombinatorialDerivation.
This should then capture the relationship between all of the products, their associated data, and the original mutagenesis design.

I have some code for doing a similar task that I'm happy to share if that's of interest; it happens to not have yet been posted publicly yet, but I can quickly clean it up a bit and do so if desired.

No need to rush posting of the code. I'm working toward the framework outlined in my question, but have a lot of other things to do first. I just wanted to make sure I was pointed in the right direction.

I do have one clarifying question/comment: From the "Representing Random Mutagenesis" part of your answer, if we are going to use mutagenic PCR to make the protein variants, then every base in the protein coding sequence would be a potential mutagenesis site. So, I would need to define ~1000 SequenceFeatures, one for each base in the protein coding sequence. Is that right?
This seems a bit cumbersome, but but not too hard to write code for. Does the idea of a SequenceFeature for each base in a protein coding sequence raise any red flags for you though?

@DavidRoss I think I misunderstood your protocol: I thought it was something like this protocol for targeted swaps with a random subset of a fixed set of variations.  If you're randomly varying every base in a large region, then making a SequenceFeature for every base in the sequence does raise a red flag to me.  I don't think we've worked out a best practice for this case yet, but my thought would be to mark the whole region as one SequenceFeature and tag it with a variant of all NNNNN.... bases

That way, you're saying "replace this region with indeterminate bases" and any variant you find will be an appropriate realization of that specification.


text:
Alignment to verify plasmid from forward and reverse Sanger sequencing reads
I am analyzing Sanger sequencing data to verify whether or not the correct protein coding sequence was inserted during cloning.
For each sample, I have both forward and reverse Sanger sequencing reads that cover the coding sequence.
Would it be better if I first align the forward and reverse reads to get a consensus, and then align the consensus to the reference?
Or should I just align each read directly to the reference?

Align all reads from a sample at the same time IMO. This will give you all of the information at once, and it should not be difficult. It should be very obvious whenever there's a disagreement between reads or between reads and reference.
I would also recommend inspecting chromatograms directly at all sites of interest (such disagreements, insertion sites, possible errors or mutations).


text:
How to obtain NCBI Nucleotide (Entrez) IDs for all E coli sequences?
I would like to gather a list of ALL E coli NCBI (Entrez/Nucleotide) IDs so that I can download fna and faa files, through a function like this:
# python code #
from Bio import Entrez
from glob import glob

def dl_genome(id, folder='genomes'):
    files=glob('%s/*.fna'%folder)
    out_file = '%s/%s.fna'%(folder, id)

    if out_file in files:
        print (out_file, 'already downloaded')
        return
    else:
        print ('downloading %s from NCBI'%id)

    Entrez.email = ""
    handle = Entrez.efetch(db="nucleotide", id=id, rettype="fna", retmode="text")
    with open(out_file, 'w') as fout:
        fout.write(handle.read())

for strain in targetStrainIDs:
    dl_genome(strain, folder='genomes')

My problem is that I don't have a good way to get the E coli strain names and their corresponding NCBI IDs for Entrez. Is there a way to get these NCBI IDs without manually searching for them?
Context
I have been trying to develop an E coli pan-genome for one of my projects. I initially just downloaded the relevant genomes (fna, faa, & features.tab files) from the PATRIC database. However, I've come to realize that the reference sequence on PATRIC has issues. Specifically, the reference sequence on PATRIC's fna/faa files for E coli K12 MG1655 notably do NOT list genes like b2031 or b2033, despite both being present in this strain.
The NCBI reference genome DOES include these genes (with proper annotations), so I'd like to just download what I need from NCBI instead. However, I'd like a systematic way to download all the fna/faa files since my pan-genome is eventually going to be built from hundreds of of these fna/faa files.
Any help would be appreciated.

You should check out NCBI Datasets for this sort of thing. Since this is bacteria, you won't be able to use the web interface but there's an API and a command-line tool that you can use.
For example, you can use the datasets command line tool to download the E. coli RefSeq genome assemblies as follows:
## see how many genomes are in scope
$ datasets summary genome taxon 562 --limit none --refseq
{"total_count":22287}

## download data excluding gff3, rna and genome seq files
## as an example but drop --exclude-seq to download genome seq as well
$ datasets download genome taxon 562 --refseq --exclude-gff3 --exclude-rna --exclude-seq --filename ecoli_genomes.zip

Do you mean the NCBI taxonomy IDs? If so, I've personally always preferred the way UniProt displays NCBI taxonomy, and it's easier to download IDs in mass.
For example, if you find E.coli's entry: https://www.uniprot.org/taxonomy/562
You can then click on "all lower taxonomy nodes", which will show all E.coli strains. If you click 'Preview first 10' before downloading, you'll see that the first column is the NCBI taxonomy ID of each strain of E.coli.
For example, E.coli strain 0.1197 is reported to have the ID 1005443 from this preview page. Indeed, if we search this ID on NCBI's taxonomy page, we do get the correct entry: https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=1005443
If this is what you want, click 'download all' to get a nice tab-separated file of all the strains of E.coli and their NCBI taxonomy IDs. This file could then be read and processed before you send API requests to Entrez.
Hope this helps!


text:
Why do molecular generation models maximize “penalized logP” as a measure of drug-likeliness?
I found that Lipinski's rule of five states that Log P (octanol-water partition coefficient, lipophilicity measure) usually should not exceed 5.
Many papers about drug discovery machine learning models tell about maximization of "penalized logP", but the paper they refer in the end does not contain any information on it. Let me show that.
Let's check 1, here we read about penalized log P:

The second task is to produce novel molecules with desired properties. Following (Kusner et al., 2017), our target chemical
property y(·)is octanol-water partition coefficients (logP) penalized
by the synthetic accessibility (SA)score and number of long cycles
y(m) = logP(m)−SA(m)−cycle(m)

where cycle(m)counts the number of rings that have more than six atoms.

Then we go into Kusner et al., 2017:

For the second optimization problem, we follow (Gómez-Bombarelli et al., 2016b) and optimize the drug properties of molecules. Our goal
is to maximize the water-octanol partition coefficient (logP), an
important metric in drug design that characterizes the drug-likeness
of a molecule.As in Gómez-Bombarelli et al. (2016b) we consider a
penalized logP score that takes into account other molecular
properties such as ring size and synthetic accessibility (Ertl &
Schuffenhauer, 2009).

Here is Gómez-Bombarelli, 2016 and it seems to contain nothing on the topic (I can be mistaken, please point me then), except

The objective we chose to optimize was 5×QED−SAS, where QED is the Quantitative Estimation of Drug-likeness (QED),37 and SAS is the
Synthetic Accessibility score.36Thisobjective represents a rough
estimate of finding the most drug-like molecule that is also easy to
synthesize.

Other examples of the papers also leading to the same papers are 2, 3, etc.
What does penalized logP mean for drug-likeliness and why it should be maximized?

A positive logP is good (hydrophobic), but high results in ADMET problems (cytochromes break overly hydrophobic compounds) and aggregation (gunk, tar, crud are words happily thrown around). However, in most cases the logP is not greater than +3 because good binding is result of hydrogen bonds etc —some bad algorithms will suggest a polyaromatic hydrocarbon that makes lots of pi-pi interactions but these are few.
Confusingly, a high synthetic accessibility score is bad (hard). 3 is the average and 10 is basically something impossible. This score is how hard it is to make (i.e. pricier), which is one of the major drivers.
The number of cycles is number of rings. Even if you get pi-pi, cation-pi and sulfur-pi bonds you also get aggregation due to stacking (a crud). So is somewhat bad.
So basically that score is
high=good minus high=bad minus high=bad

If difference in relative Gibbs energy of binding is thrown it, negative ∆∆G is good, so that would be subtracted. As its atom count dependent ligand efficiency is better, but that ranges from 0 to -1.2 kcal/mol. So would need scaling.

You were trapped by arXiv versioning convention that shows you the latest version of the preprint as default. It seems that this experiment was removed from the final publication.
If you have a look at the v1 (https://arxiv.org/abs/1610.02415v1) instead of the latest version (v3), you find the explanation:

As a simple example, we first attempt to maximize the water-octanol
partition coefficient $logP$, as estimated by RDkit. [43] $logP$ is an
important element in characterizing the druglikeness of a molecule,
and is of interest in drug design.  To ensure that the resulting
molecules to be easy to synthesize in practice, we also incorporate
the synthetic accessibility [44] (SA) score into our objective. Our
initial experiments, optimizing only the $logP$ and SA scores, produced
novel molecules, but ones having unrealistically large rings of carbon
atoms. To avoid this problem, we added a penalty for having carbon
rings of size larger than 6 to our objective. Thus our preferred
objective is, for a given molecule m, given by:
$J(m) = logP(m) − SA(m) − ring-penalty(m)$,          (1)
where the scores $logP(m)$, $SA(m)$, and $ring-penalty(m)$ are normalized to
have zero mean.


text:
How to share a collection of sbol designs in SynBioHub
We currently work on a project that will generate a collection of recombinant plasmids targeting over 3500 genes of the cyanobacterium Synechocystis sp. PCC 6803. We would like to share these designs in SynBioHub but we are not sure how best to go about it. The plasmids have been designed and now they are starting to be produced in batches, checked and finally transformed into Synechocystis to induce a single gene deletion by homologous recombination.
During the project, we will have:

Stage 1 – the initial plasmids designs, waiting for production. They
are worth sharing as the community could influence the production
queue.
Stage 2 – the assembled plasmids. In some cases they will have
different structure than initial designs from stage 1. If nothing
else their description should reflect the status change to
'produced'.
Stage 3 – plasmids are verified by sequencing (also to
check the barcodes are correct) and sequences attached.
Stage 4 – plasmids will be transformed into Synechocystis and checked by PCR to
see if i) the target gene has been knocked out and ii) if the strain
is segregated (i.e. all copies of the genome in the strain contain
the gene deletion).  This information should also be recorded with
the description.

Our initial plan was:

Create a collection with initial designs (component definitions) Stage 1
As stage 2 progresses add new versions of udpated designs
Stage 3, 4 update plasmids descriptions (synbio desc or notes) and attach additional files (e.g. from sequencing).

However, our current understanding is the public collections are frozen so stage 2,3,4 cannot be accommodated.
We are seeking advice how to plan for sharing this resource and some clarifications.

Is there a way of making collection readable by others without freezing its content?
What operations are possible on public component definitions, can the free text descriptions/notes be updated or file attached?
Can new component definitions or their versions be added to a published collection
For 2) 3) who can perform these operations, the collection’s creator or SynBio curators?

Looking forward to input from SynBioHub community

Welcome. Typically we ask that questions only contain a single question. Can you please consider condensing your post down to a single question and/or asking the remaining question in a separate post. Thanks.

Knowing something about SynBioHub, I think this is well-answerable as one question, because it's all focused on its publication workflow. I'm not sure of the answers myself, but I'm going to try to get one of the folks who does know to answer.

There is a balancing act between allowing all edits, which defeats the purpose of having URIs point to static objects, and allowing no edits, which defeats the purpose of providing some updated information about an object.  SynBioHub resolves this issue by allowing a limited amount of edits after publication.  To answer your questions specifically:

Yes, you can share objects by generating a share link OR adding another owner to an object.  Sharing can be done by the original owner of the object.
There are some mutable fields that any owner of the object can change, these include the description, notes, source, and references in the Details section.  You can also add attachments.
Yes, when you publish (Make Public) a new object or collection, you have the option of making it public into an existing collection OR a new collection.  A person with curator privileges on the SynBioHub instance you are working on can perform this operation.
Answered above.


text:
Do proteins move inside the cell?
Do proteins move inside the cell and how fast do they move or move? Please support your answer with articles or books.

I’m voting to close this question because while it may be vaguely related to biology, it isn't related to bioinformatics.

Proteins move in cells or get transported in the cell. You can go to NCBI (https://www.ncbi.nlm.nih.gov/), and start searching for protein dynamics within cells. Also there is a technique called FRAP (https://en.wikipedia.org/wiki/Fluorescence_recovery_after_photobleaching) which is used to measure kinetics and dynamics.
Also, this differs probably a lot between the cellular system you are asking for. You can imagine that some cells have a more active transport system because stretched out (neurons) vs others with maybe less active transport systems, because of a more static nature.


text:
How to execute several processes sequentially, without channels?
Using DSL2.
I have processes A, B, C which I want to run in sequence - B waits for A to finish, C waits for B to finish.
They are not communicating via channels so I can't do A.out.collect().
I have tried making dummy channels but this is problematic and a hack. The onComplete subscription doesn't seem to let me trigger another process.
What is the correct way to run a sequence of processes which are not communicating or triggered by channels?
P.S. I am aware that this is not really the intended use case of Nextflow (which is oriented around channel communication/triggers... but... I just need to get this done.)

Could you please elaborate on how you tried to do dummy channels? Maybe by giving an example. I would suggest to use watchPath but using this type wouldn't allow your pipeline to finish.

I would guess it is not possible without channels unless you are using beforeScript and afterScript. I came up with this solution:
#!/usr/bin/env nextflow

nextflow.enable.dsl=2

workflow {
  // C(B(A("")))

  next = ""
  next = A(next)
  next = B(next)
  next = C(next)
}

process A{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

process B{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

process C{
  echo true
  input: val _
  output: val ""

  """
  echo $task.process
  """
}

By using workflows you can change the execution order. But it is still a hack and I would prefer something like a next: keyword to specify the next process.

And for those not interested in DSL2:
process A {
  echo true

  output:
  val "" into B

  """
  echo "$task.process"
  """
}

process B {
  echo true

  input:
  val _ from B

  output:
  val "" into C

  """
  echo "$task.process"
  """
}

process C {
  echo true

  input:
  val _ from C

  """
  echo "$task.process"
  """
}

Also, I am quite curious to see other solutions.


text:
How to trigger next process on completion of previous process?
Simple question but I am stumped. (I am still learning Nextflow)
I have a process which takes as input a channel of files and processes each file separately.
There is a second process which needs to wait until all files are complete, so it can merge the outputs into a single collated file.
How can I do this? I have tried subscribing to "onComplete" event to no avail as I cannot specify a process as the function to trigger.

files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    ????

    output:
    ????

    script:
    """
    echo "How to trigger this process once, after all individual files are processed?"
    """

}

??? Want to do something like this
result.subscribe onComplete: merge_results

```

I've run into the same problem but I'm interested in trying this method. Did it work for you?

So from this page...
https://groups.google.com/g/nextflow/c/jt77_-uApMs
It is the "scatter-gather" pattern. I think the answer is to collect() the channel as follows - interested to hear if there is a better way
files = Channel.from params.the_files_list

process process_files {

    input:
    val file_path from files

    output:
    stdout result

    """
    python ./process_one_file.py $file_path
    """
}

process merge_results {

    input:
    val go from result.collect()

    script:
    """
    echo "Triggered once after all files complete!"
    """

}

.collect() is the correct way to do this


text:
Nextflow: Can I parameterize similar processes inside a for loop?
I have set of ~5 processes which perform very similar tasks such that their differences can be parameterized. Rather than putting 5 items through some channel (which seems complicated), or defining 5 separate processes, it appears to be possible in Nextflow to declare the processes in a for loop, and output into an array of channels.
Below is a prototypical example and it seems to work as intended. (I understand this could be done differently with some channel cross product operator but I am just giving this as an example of a process in a for loop.)
Is there anything wrong with doing this? Is there a better way to do it?
Bonus question - the name of the process below is "foo". I would like to have the index reflected in the name i.e. "foo_$i" but this doesn't seem possible by default - any way to do this??
some_inputs = ['prot','dna', 'rna']

some_parameters = [5,3,2,7,1]

// Set up an array of channels to output to
receivers = new ArrayList<Channel>()
for (i = 0; i < 5; i++) {
    receivers.add(Channel.create())
}

// Iterate 5 times to create 5 processes each with different parameters
for (i = 0; i < 5; i++) {

    // WOW - process can be defined inside a for loop!!
    process "foo" {

      // Extract parameters here so correct $i value is used
      p = some_parameters[i]

      input:
      val x from Channel.from(some_inputs)

      output:
      val "$x $p" into receivers[i] // Note that correct index $i seems to be used here

      """
      """
    }
}

// Confirm that each of the processes outputted independently with the correct parameter
receivers[0].view { "Received[0]: $it" }
receivers[1].view { "Received[1]: $it" }
receivers[2].view { "Received[2]: $it" }
receivers[3].view { "Received[3]: $it" }
receivers[4].view { "Received[4]: $it" }

Output - note that the parameter is consistent within each input channel i.e. working correctly:
N E X T F L O W  ~  version 20.10.0
Launching `play.nf` [fabulous_keller] - revision: c89e2b2281
executor >  local (15)
[10/b477ec] process > foo (3) [100%] 3 of 3 ✔
[16/8282ba] process > foo (2) [100%] 3 of 3 ✔
[db/f7327e] process > foo (1) [100%] 3 of 3 ✔
[fb/886f6f] process > foo (1) [100%] 3 of 3 ✔
[46/c34ed3] process > foo (1) [100%] 3 of 3 ✔
Received[1]: rna 3
Received[1]: prot 3
Received[4]: rna 1
Received[2]: rna 2
Received[0]: dna 5
Received[2]: dna 2
Received[2]: prot 2
Received[0]: prot 5
Received[3]: rna 7
Received[0]: rna 5
Received[1]: dna 3
Received[3]: dna 7
Received[4]: dna 1
Received[3]: prot 7
Received[4]: prot 1

Use channels, this is what nextflow's entire design is based on. I'm honestly surprised that spawning 5 processes with the same name in a loop even works.
Turn your parameters into regular channels, use the combine operator (https://www.nextflow.io/docs/latest/operator.html#combine) to transform the 2 parameter channels into a channel that has each combination of params/inputs. Feed this channel into the  foo  process. Afterward use branch (https://www.nextflow.io/docs/latest/operator.html#branch) to split your foo output channel into 5 different channels.
typechan = Channel.from(5,3,2,7,1)
numchan = Channel.from('prot','dna', 'rna')

cartesian = numchan.combine(typechan)

process foo {
  input:
  val numtype from cartesian
  output:
  val numtype into bar
  script:
  """
  echo $numtype
  """
}

bar.branch {
  one: it[1] == 1
  two: it[1] == 2
  three: it[1] == 3
  five: it[1] == 5
  seven : it[1] == 7
}.set {result}
result.one.view()
```

You don't need to create processes in a loop when you have channels. The each qualifier allows you to repeat the execution of a process for each item in a collection. For example, the contents of 'test.nf':
some_inputs = Channel.of( 'prot', 'dna', 'rna' )
some_parameters = Channel.of( 5, 3, 2, 7, 1 )

process foo {

    tag { "task: ${task.index} - ${x},${y}" }

    input:
    val x from some_inputs
    each y from some_parameters

    output:
    tuple val(x), val(y) into receivers

    """
    echo "x: ${x}; y: ${y}"
    """
}

results = receivers.branch { x, y ->

    five: y == 5
        return tuple( x, y )

    three: y == 3
        return tuple( x, y )

    two: y == 2
        return tuple( x, y )

    seven: y == 7
        return tuple( x, y )

    one: y == 1
        return tuple( x, y )
}

results.five.view { "My result for param 5 is: $it" }
results.three.view { "My result for param 3 is: $it" }
results.two.view { "My result for param 2 is: $it" }
results.seven.view { "My result for param 7 is: $it" }
results.one.view { "My result for param 1 is: $it" }

Then running like:
nextflow run -ansi-log false test.nf

Results:
N E X T F L O W  ~  version 20.10.0
Launching `test.nf` [fabulous_franklin] - revision: 14c3290d2b
[6e/c12ce7] Submitted process > foo (task: 7 - dna,3)
[93/616f93] Submitted process > foo (task: 4 - prot,7)
[27/0bc3c4] Submitted process > foo (task: 3 - prot,2)
[0d/0d0b51] Submitted process > foo (task: 9 - dna,7)
[d4/f56d18] Submitted process > foo (task: 6 - dna,5)
[b0/d1535f] Submitted process > foo (task: 8 - dna,2)
[13/0a660a] Submitted process > foo (task: 5 - prot,1)
[f0/1c2f57] Submitted process > foo (task: 2 - prot,3)
[3f/3e893d] Submitted process > foo (task: 1 - prot,5)
[19/2921de] Submitted process > foo (task: 10 - dna,1)
My result for param 3 is: [dna, 3]
My result for param 7 is: [prot, 7]
My result for param 7 is: [dna, 7]
[a7/953a2b] Submitted process > foo (task: 11 - rna,3)
[1c/0fd4ef] Submitted process > foo (task: 12 - rna,5)
My result for param 1 is: [prot, 1]
[39/d21adf] Submitted process > foo (task: 14 - rna,7)
My result for param 5 is: [dna, 5]
My result for param 2 is: [prot, 2]
My result for param 2 is: [dna, 2]
[f6/8509a7] Submitted process > foo (task: 15 - rna,1)
My result for param 5 is: [prot, 5]
[6f/cee0ae] Submitted process > foo (task: 13 - rna,2)
My result for param 3 is: [rna, 3]
My result for param 1 is: [dna, 1]
My result for param 2 is: [rna, 2]
My result for param 3 is: [prot, 3]
My result for param 1 is: [rna, 1]
My result for param 5 is: [rna, 5]
My result for param 7 is: [rna, 7]

The answer to your bonus question is to use the task.index variable in your script. I used this to tag your job, but it can be used anywhere in the process declaration. You likely won't need this variable though.
Usually, using the branch operator is enough for most use cases. However,it's possible to implement your own branching function. The following uses the new DSL, DSL2, to automatically fork the 'receivers' channel:
nextflow.enable.dsl=2

process foo {

    tag { "task: ${task.index} - ${x},${y}" }

    input:
    val x
    each y

    output:
    tuple val(x), val(y)

    """
    echo "x: ${x}; y: ${y}"
    """
}

def branchByElement( input_ch, elements, by=0 ) {

    def targets = [:]
    for( element in elements ) {
        def key = element
        targets[ key ] = { input_ch.filter { it[by] == key } }
    }
    return targets
}

workflow {

    some_inputs = [ 'prot', 'dna', 'rna' ]
    some_parameters = [ 5, 3, 2, 7, 1 ]

    inputs_ch = Channel.fromList( some_inputs )
    parameters_ch = Channel.fromList( some_parameters )

    foo( inputs_ch, parameters_ch )

    parameter_results = branchByElement( foo.out, some_parameters, by=1 )
    input_results = branchByElement( foo.out, some_inputs, by=0 )

    parameter_results[7].call().view { "Param 7: $it" }
    input_results.dna().view { "Input dna: $it" }
}

Results:
N E X T F L O W  ~  version 20.10.0
Launching `test.nf` [distraught_meninsky] - revision: 04df251752
[bb/32325f] Submitted process > foo (task: 4 - prot,7)
[a0/19fdf5] Submitted process > foo (task: 8 - dna,2)
[e9/4c1d09] Submitted process > foo (task: 6 - dna,5)
[a3/dcb002] Submitted process > foo (task: 2 - prot,3)
[db/7e612e] Submitted process > foo (task: 5 - prot,1)
[cd/dfa84e] Submitted process > foo (task: 9 - dna,7)
[d4/6a85f2] Submitted process > foo (task: 7 - dna,3)
[f0/8e8bf4] Submitted process > foo (task: 1 - prot,5)
[89/f20a28] Submitted process > foo (task: 3 - prot,2)
Input dna: [dna, 2]
Param 7: [prot, 7]
Input dna: [dna, 5]
[53/056fe4] Submitted process > foo (task: 10 - dna,1)
[f3/a61abf] Submitted process > foo (task: 12 - rna,3)
Input dna: [dna, 7]
[f9/a568a8] Submitted process > foo (task: 13 - rna,2)
Param 7: [dna, 7]
[06/72bb6c] Submitted process > foo (task: 11 - rna,5)
Input dna: [dna, 3]
[a3/6d80d8] Submitted process > foo (task: 14 - rna,7)
[ac/acdb5a] Submitted process > foo (task: 15 - rna,1)
Param 7: [rna, 7]
Input dna: [dna, 1]

I think each is even more elegant than combine, do you also have a more elegant way to split the resulting channel than branch ?

@Pallie I think there is using a function and the new DSL2 - please see my update above.


text:
Holistic enzyme activity determination with computation
I'm currently working on a program which will determine the activity difference between an original enzyme and a variant with just 1 or 2 variants. Of course, I'm not talking about the "kinetic activity", but some indications of a probable activity difference should be given. The system has to be sensitive enough to determine the probable activity difference between unthought-of variants.
From the information that I've got from my research, there is two option. I can use one of 3 quantitive enzyme reaction analysis tools (such as EVB, FEB and LRA) or their hybrids.
EVB => Empirical Valence Bond
FEB => Free Energy Perturbation
LRA => Linear Response Approximation
Source:
https://www.researchgate.net/post/Can-the-activity-of-proteins-be-studied-using-molecular-dynamic-simulations
Now, my questions are: - Does these application be used for such purpose? - For an accurate prediction of  enzyme activity difference, are these programs good enough or require developement? - Which one of them can give the best result? - And lastly, are QM/MM methods a must for an accurate answer?
I'm a high school student so sorry if my questions sound dumb.
If you can answer my inquiry's I would appreciate it:)
Thanks

Your title says holistic. This is a tad problematic as there's layers upon layers. Say, post-translation regulation, inhibiting metabolites, interacting protein etc. That is why when talking of an enzyme inhibitor, at the biochemical level one speaks in terms of k_i (inhibition constant), while at the cellular level ("holistic") one talks of IC50.
This could be a biocatalysis/synthetic biology problem (make an enzyme better) or translational medicine problem (why is this mutation pathogenic).
In terms of in silico enzyme design, one can have all sorts of issues that need to be addressed:

catalysis
side reactions
stability
aggregation
solubility

Also there's the expression problems, which may be due to subpar promoter, insufficient cofactors and poor codon usage.
For clinical variant assessment, a variant may:

destabilise leaving to degradation
remove a post-translational modification
disrupt the active site, affecting catalysis
lose a binding partner
result in incorrect localisation
result in solubility issues

So force-field scores are not the sole solution. Furthermore, a forcefield method gives one a ∆∆G (difference in relative Gibbs free energy), not a K_M (Michaelis constant) and a k_cat (catalytic rate), although it is somewhat possible to calculate these, or a relative variant of these.

a variant with just 1 or 2 variants

A variant protein can have one or two mutations (e.g. G10R) in its sequence. A mutant and a variant are near synonyms, but variant is generally intended to be a mutant in a population.

Does these application be used for such purpose?

For enzyme design MD simulations are not generally the first step. Rosetta pmut_scan for example is a nice starting step as the mutational landscape of a protein is N^19, while a MD simulation takes weeks and 6+ core per single variant.
This calculates only an increase in stability (relative difference in Gibbs folding energy), not reactivity. For that you can use a transition state, whose results are inaccurate, but majorly affected by the set-up, i.e. junk-in junk-out plays a big part. There's a saying 'a week in the lab saves you an hour at the library': this is backwards for sarcasm as it happens too often, namely that is background research is always a lifesaver. An QM-MD simulation on the wild type reaction definitely can help in elucidating things if it is not clear.
Online tools
Generally, however, there are resources out there to do many stability calculations already.
SynBio If you are interested in picking a residue or two to make an enzyme more stable, without resorting to a thermophilic enzyme, PROSS is a very nice tool.
Clinical (Caveat: If your variant is a variant that affects you, a friend or a relative, you/they should speak to a genetic councillor and not do research on it yourself) There are several tools that can describe a variant at the protein level. Do remember that ∆∆G is a potential so more negative is more stable —except in SDM where the score is inverted.

Missense3D
Venus
Miscast
SDM


text:
Normalize RNA seq data from multiple runs for expression analysis
I have RNA samples sequenced with TruSeq Stranded Total RNA kit protocol in Illumina HiSeq (2x125bp) and NovaSeq platforms (2x150bp) - almost 100 samples altogether. I have to use the samples data for expression analysis. Question is how do we normalize samples from different runs and of read length? What is best way (with available packages) we can normalize the count (to overcome the bias) and do a DE analysis? ANy suggestion would be helpful.

There are a number of methods. If you're doing DE, you have ComBatSeq, SVAseq, RUVseq, BUSseq. You could also try Z-score normalization ($(x - \overline{x}) / \sigma$) or even quantile transformation. For the latter two, make sure you work on each batch individually, not the whole dataset at once. See more on that here.
To visually compare if transformation works, plotting PCA/UMAP/t-SNE on raw and transformed data can perhaps be of some insight.

Thanks Roman for the response. Very useful. Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or do you suggest to include read length as a factor during DE analysis? What do you suggest.

This sounds like another question. :) Feel free to post it as such.

Okay, let me do that. Don't know if this kind of experiment has been done previously (happens data is used from multiple sources for integrated analysis)

I don't think RNASEq is that sensitive to running different things on different instruments.  Differing read lengths will change things.  You could always trim the reads down so they all match.  Or include read length as a factor in the design when looking for DE genes.

Normalization and batch-effect correction are two different things.
You always have to do normalization of RNAseq data in order to be able to compare the runs and do DE analysis, for instance. Plenty of methods and you have to search for them as normalization methods indeed.
What you are asking and concerned about is "batch-effects" in your data as a function of having used two different Illumina machines. This is a valid concern, but before you try to correct for it or mitigate it, you first have to assess whether it is actually there. Maybe is not at all. To figure out if you have batch-effects you have to first do data exploration, this is why the use of PCA/UMAP/t-SNE has been already suggested. In this plots you will be trying to see whether your data cluster together by machine type instead of by the biology it represents. Read about how to assess the presence of batch effects.
If you do observe it or anyway want to correct for its presence, then indeed you can try one of the methods also already mentioned: ComBatSeq, SVAseq, RUVseq, BUSseq.
Rather than the two different machines, the batch-effect might actually be due to different people doing the lib preparation, experimental settings if experiments were done at different moments, etc.

I agree with swbarnes that you don't need to explicitly model for different read length, but you need to normalize for different production of sequencing runs.
Among the several available packages, DEseq2 is one of the most widely used.
The normalization method computes the geometric mean of expression of each gene across samples; then, expression of the gene in each sample is divided by the geometric mean. Once this has been don for all genes, in each sample it computes the median of values, and this is the "size factor". Ok, it looks convoluted, but it is indeed easy, and it is implemented in DESeq2 and in other software packages.
You can find methodological details here.
Finally, DESeq2 also offer the possibility to account for additional experimental details (e.g. read length) and this should reduce the risk of spurious results.

Thanks for your response, Fabio. If I use ether median-to-ratios (geometric mean in DEseq2) or upper-quartile normalization, I think it should work good despite different sequencing machines? Also do you have any articles where the read length is used as a factor in the design for DE analysis? How robust is this method compared to trimming of reads as trimming on large samples is a huge additional task.


text:
RNAseq DE comparison with samples of different read length
I have RNA samples but with different read lengths (Eg, HiSeq 2x125 and NovaSeq 2x150bp data). I would like to do DE analysis on these samples. What do you recommend? Do you recommend to trim the reads to 125bp for the NovaSeq to compare with the other batches from HiSeq? Or any other suggestions.

Does this answer your question? Normalize RNA seq data from multiple runs for expression analysis

You don't need to do anything, this is taken care of by the linear model - if you don't have any confounding effects.
If you do have confounding effects, e.g. all control samples on HiSeq and all treated on NovaSeq there's not a lot you can do anyway, other than be honest when you present your results.
Taken care of by the model means adding an extra 2 level factor for sequencer by sample, you need to do this anyway to take care of the batch effect between sequencer runs.


text:
polidy VCF file for Canvas somatic-wgs cnv calling
I am trying to run Illumina CANVAS cnv caller for Somatic-WGS. There is an option "--ploidy-vcf" which is mandatory to supply, but don't know what exactly that mean. I had supplied the CNV.vcf file from the germline WGS (matched normal run on Canvas smallPedigree), but it has given an error that says "Reference copy number > 2 is not supported" at the last stage of run.
I would like to know your experience when you used CANVAS on somatic WGS samples. What ploidy-file is appropriate here? Should it be only chrX and Y in it? Your suggestions please.

It seems you need to supply vcf files specifying the ploidy of the sex chromosomes.
See https://github.com/Illumina/canvas/issues/89 for explanation and example vcf files.


text:
Clustering with Seurat
My data consists of two groups control and stim. I'm concerned that my clustering in the stim group is dominated by the unregulated genes. I can remove these genes and recluster but is there a way to then project the removed genes onto the new umap space so see if they are indeed all clustered together. They are no longer in the seurat object so I can't think of how to do this. I hope this makes sense.

Yes, this can be done, particularly if your concern is with the clustering / UMAP rather than the transformation process; I've done this recently with some Rhapsody single cell data that we generated. The main trick is adding 'return.only.var.genes = FALSE', as well as calculating a variable set that is a bit larger than the desired set. Depending on the number of cells, this will balloon the size of the resulting Seurat object pretty significantly.
Note that UMAP is a technique to place cells, rather than genes. If you have the information on where the cells should go, then plotting gene expression on top of that via FeaturePlot should work without any reprojection needed.
Here's what I did (with experiment-specific variables excluded):
sct.data <- SCTransform(sct.data, vars.to.regress="percent.mt",
                        variable.features.n = 3000 + length(genesToExclude),
                        return.only.var.genes = FALSE);

sct.data <- RunPCA(sct.data, verbose = TRUE,
                   features=setdiff(VariableFeatures(sct.data), genesToExclude));
sct.data <- RunUMAP(sct.data);
sct.data <- FindNeighbors(sct.data);
sct.data <- FindClusters(sct.data);

All transformed genes should be available for querying in the results.


text:
Programmatic way of accessing the FASTA sequence similarity tool (or similar) in Python
I am looking for a tool that performs a sequence similarity algorithm for proteins. More specifically, I am looking for something that would be usable in Python (or anything else usable in the command line actually) which has a similar output to the FASTA tool, i.e:

For my use case, the 'Source' column is not super relevant, so it does not matter if it's not included in the output. What's most important is the DB:ID and the match statistics.
Have you used anything similar? Do you have any suggestions regarding such a tool?
Thanks in advance!

Hmmm ... two suggestions, you can do both via Python, Blast or a tree matrix. Its an interesting question but difficult to know what you're looking for

Will you be comparing specific, known sequences, or is part of the question finding similar sequences to your query in a large database of non-related sequences? Basicaly, do you want a sequence aligner (mafft, t_coffee, clustal etc) or tools like the various flavors of blast (psi-blast, p-blast)? Or are you looking for proteins with shared domains which means tools like HMMER or psi-blast again? Pretty much everything in bioinformatics can be done from the command line, so you need to explain exactly what question you're trying to answer for us to suggest the better options.

Apparently there is a REST API for FASTA. Looks like there are few dependencies but other than that should be straightforward to use.


text:
How to convert Foldseek .m8 alignment files to FASTA format?
I want to use sequence alignments found by Foldseek in PROSS. I know there are other ways to get alignments (also tried with results of HHPred), but I had in mind that Foldseek may find remote homologs that would still be structurally relevant.
Currently Foldseek only gives alignment results in m8 format, however PROSS only accepts FASTA alignment files.
I tried looking into BioPython but it seems it doesn't open m8 files, only m9. I also tried importing the alignment in Clustal X but it can't read it.
I could probably use Foldseek CURL API to get the alignment as json and make a script to write the results in FASTA format, but this seems like a heavy handed solution.
What would be other tools I could try to convert m8 to FASTA?

tl;dr
Use an awk script (below) or "one-liner"
awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta

Caveat emptor: I've tested this on my own Foldseek outputs and it works for my pedestrian purposes of multisequence alignment with COBALT. I haven't tested it with PROSS but expect it should work for you too.
Long answer:
Using this definition of m8 (which is a modification of the BLAST m8 format):
https://github.com/steineggerlab/foldseek#output-search
And this definition of FASTA:
https://en.wikipedia.org/wiki/FASTA_format
We see we want to take lines like this (I've added [...] here to show what I've cut off for brevity):
job.pdb_A   Some matched structure name 51.100  311 138 8   2   305 25  328 1.000   2.790E-35   1272    306 334 VDVSVSIFINKIYGVNTLEQTYKVDGYIVAQWTGKPR--KTPGDKPLIVENTQIERWINNGLWVPALEFINVVGSPDTGNKRLMLFPDGRVIYNARFLGSFSNDMDFRLFPFDRQQFVLELEPFSYNNQQLRF--SDIQVYTENIDN[...] -27.488,-8.287,-35.993,-26.838,-4.500,-36.421,-25.278,-3.958,-32.938,-22.398,-6.436,-33.555,-21.375,-4.861,-36.892,-21.373,-1.326,-35.385,-19.237,-2.585,--18.044,-1.289,8.073,-14.345,-2.072,7.367,-15.099,-5.282,5.328,-16.872,-3.358,2.522,-16.876,-5.018,-0.77 [...] MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]   1918945 Vibrio spartinae

and map them to entries like:
>some matched structure name
MIKKEILSAIASLVWFCSPALYAQYVVDTSVSINKIYGVNTVDQTYKVDGYLVTSWDEPDPTYRPASGERVYENHHADQLMNDGLWVPAFEFINIIGQRQTANRRVVISADGRVTYNERFQGLFTTTMDFRQFPFDYQTFVLMMEPFSYEQKQL[...]

We can do this with an awk script (The awk manual is your friend here https://www.gnu.org/software/gawk/manual/html_node/Long.html). Put the following in a file somewhere handy and call it using awk -f m8_to_fasta.awk my_m8_file.m8 > my_fasta_file.fasta
you may want to combine this with head, tail, and/or grep to reduce the number of sequences that end up in your output.
#! /bin/awk
{
        name = "> "
        # start on the second field, foldseek will put the job name in $1
        field = 2
        # Check each field until we see numbers
        while ( $field + 0 != $field ) {
           name = name " " $field
           field = field+1
        }
        # name now contains the name string which is the first line of a fasta entry

        field = NF
        # Now find the matched sequence
        # Sometimes an entry ENDS in a number, so we strip that
        if ( $field + 0 == $field ) { field = field -1 }
        while ( $field + 0 != $field ) {
            field = field - 1
        }
        # field is now pointing at some undocumented number. Decrement once more to get the sequence
        field = field - 1
        # field is now pointing at the sequence
        seq = $field

        #Now I/we want to sanitize the sequence to remove any gaps
        gsub(/-/, "", seq)

        print name
        print seq
    }

Or we can coerce this into a monstrous one-liner:
 awk '{name = "> "; field = 2; while ( $field + 0 != $field ) { name = name " " $field; field = field+1 } ; field = NF; if ( $field + 0 == $field ) { field = field -1 }; while ( $field + 0 != $field ) { field = field - 1 }; field = field - 1; seq = $field; gsub(/-/, "", seq); print name; print seq}' the_m8_file.m8 > the_fasta_file.fasta


text:
What are the best practices for comparing .seg copy number outputs from iChorCNA to VisCapCancer?
I have shallow whole genome sequence data that I used the tool iChorCNA to find the mean copy number segmentation per region. The tool produces a seg file and can be seen here: https://github.com/broadinstitute/ichorCNA
I have also run these same samples to higher coverage on a targeted panel and used the tool VisCap Cancer to find the mean copy number segmentation. The tool can be found here and also produces a seg file: https://github.com/pughlab/VisCapCancer
I am looking to compare the outputs of each tool to see how close they overlap with each other. The problem is that iChorCNA produces outputs in bins that span the genome per million bases, and VisCap Cancer is only looking at the sites for targetted probes. These probes are anywhere from 1.3 million basepairs to 260 million basepairs in length.
What I am looking to do is segment the VisCap output into bins that match the bins used by iChorCNA. I want to then take the sum of the mean value given by iChorCNA at each bin that overlaps with the VisCapCancer results, and divide it by the total number of bins that fall in the range of the results to find the average iChorCNA reading per VisCap reading. Then, I can do a t-test, ANOVA or other measure to see how closely these two numbers are to each other.
I was wondering if anyone has directly compared copy number results from two data sources such as this before, such as shallow whole genome sequencing against a targeted panel, and if they have any advice as to how to go about this most effectively - including if it would be better to do the comparison using R, python or another language, and best practices for setting up the comparison most effectively.
Thank you for your time


text:
Pros and cons between genomes and transcriptomes?
What are the major differences between the use cases of a shot-gun genome and an mRNA capture transcriptome? Especially when it comes to downstream analyses such as looking for orthology, and selection analysis.
From what I can visualize in my head, I can come up with the following points:

The genome has the structural data, information on chromosomes, and haplotypes while the transcriptome lacks this information.

The genome records the SNPs more directly than transcriptomes. I might be thinking of this wrong but the genomes would have the actual bases and we can then compare the bases between the genomes but for transcriptomes differences in expression would show up as presence/absence variations that have no way of being traced back?

The transcriptomes should be easier to work with when it comes to computational requirements since they tend to be 25 to 50 times smaller than genomes.

I could be completely wrong about the points I listed and could be missing something major. What do you guys think?

It really depends on your use case. If this is human data, where we already have a fairly polished reference genome, then genome sequencing adds value for the genetic variation, e.g. rare disease genetics, GWAS or cancer genomics. Transcriptome sequencing allow exploration of tissue (or cell type) specific patterns and RNA biomarkers etc. Calling variants from RNA can be interesting in the context of splicing and RNA editing.
For a non-model organism where there isn't currently a reference genome, then the genome creates a reference for future studies to map genomic/transcriptomic reads to, so it's usual to do the genome first.
Since you mention ortholgy and selection studies, well traditional comparative genomics studies have done these analyses from the genome only, but the transcriptome can add value depending on your questions.
Finally, I wouldn't consider computational requirements as you stated: it depends on sequencing depth/coverage and other factors, and is likely not to drive the choice between the two assays

Thanks for responding. So, when I mention comparative genomics what I am trying to say is that getting a good genome for a non-model organism is not easy but it seems RNA data is more abundant. I am trying to see if I can stockpile information on plant evolution from RNA data only since genomes are so rare. As for the computational part, when you feed multiple whole genomes to off-the-shelf software they usually just end up with segfaults (since whole genomes can easily be significantly larger than transcriptomes). I am sure that's not an issue in targeted research based on individual samples.


text:
Finding 2000 bp sequences that occur more than once in the genome
I am looking for a software that can identify long sequences (1500-2000 bp) that occur more than once in my genome of interest. A k-mer counter like KMC could have worked, but KMC has a max limit of 256 bp. If I can club contiguous 256-mers from KMC, that would work too, but I don't know how to do that. Any ideas?

MUMmer includes a repeat-match tool which can do exactly that:
$ repeat-match
USAGE:  repeat-match  [options]  <genome-file>

Find all maximal exact matches in <genome-file>

Options:
 -E    Use exhaustive (slow) search to find matches
 -f    Forward strand only, don't use reverse complement
 -n #  Set minimum exact match length to #
 -t    Only output tandem repeats
 -V #  Set level of verbose (debugging) printing to #

Note that a "maximal exact match" is just a maximal repeat.

This looks ideal. Will try this out and update. Thank you!

Okay I tried it out and it worked, so thank you! Ended up using mummer instead of repeat-match so as to compare between (and not within) sequences


text:
Problems about value log2(IP/Input) less than zero in ChIP-seq?
I have a question about the normalization for ChIP-seq. I used CPM to normalize my bam files of each IP and Input. Then I calculate the coverage of gene bodies for all genes on the genome. I have WT and different mutants. Take H3K4me1 as an example, I first calculate the matrix over gene bodies for H3K4me1 IP and Input. Then I used the matrix of H3K4me1 divided by the matrix of input. then log2(H3K4me1/Input) and plot this value over entire gene bodies. I found all the regions are lower than 0 and it seems that my IP does not work. In fact, this is not the case since I have checked the antibody's efficiency via ChIP-qPCR and it works well. Moreover, the distribution of H3K4me1 is consistent with the published data. But the problem is the profile of log2(H3K4me1/Input) is lower than 0.  I don't know how to deal with such problems? I would appreciate it if someone can give me some suggestions on that.  Thanks!

This is a common issue caused by the normalization method not being particularly robust (most likely it's using the mean signal for normalization, which will be a bit upward skewed when you have regions of high enrichment in the ChIP). The basic ways to fix this are:

Use a more robust normalization (e.g., the median signal, or even just manually select un-enriched regions to calculate the normalization factors)
Manually tweak the normalization factors so the background (i.e., log2 ratio in regions where you don't expect any enrichment) is 1. That will generally make the regions you're plotting more biologically realistic.

Note that there can be a negative enrichment in some areas, due to things like exclusion of histones, so don't be overly zealous if you manually tweak the normalization factors.

Dear Devon,  Thanks for your kind reply and valuable suggestions! I have tried to use the median value of the non-enriched region to calculate the scaled factor of each IP vs their input. I found it works for all modifications but failed for H3 IP. I guess this is due to the peaks for H3 are widely spread on the genome and when I use MACS2 to do the peak calling for H3, it won't work efficiently. Do you have any suggestions for H3 IP normalization? I am looking forward to your reply.  Thanks!

You may have to manually tweak normalization values to get reasonable results.

What I did: first did peak annotation, then extracted the genes(name it as gene list 1)  with enriched peaks. Then I calculated the matrix of the gene body of genes (not included in gene list 1 ), further calculate the median value for each sample (for both IP and input).  In the peak annotation step, I tried Homer and Rgmatch, unfortunately, the results of these two software are quite different.

Therefore, I recalculated the median value given the coverage of the non-peak region (just escaped the peak annotation step). It seems the results of using Homer were much closer to the results directly given the scaled factor for the non-peaks region.  I checked the default parameter of Homer, they defined TSS as from -1kb to +100bp and TES as -100bp to +1kb. Since S.cerevisiae. genome is not big like human or mouse, I am not sure how much effect using the default parameter.   I am pretty interested in the 500bp up-and downstream of gene bodies.

The Homer source code is based on perl and I haven't used Perl before.  Do you have any other suggestions? Thanks！


text:
How to choose the proper way to analyze the differential binding peaks among different samples for ChIP-seq? DiffBind or CSAW?
I have used the TMM method in CSAW package to normalize the composition bias for my different ChIP-seq samples and also got the normalized bigwig files. Next step I want to analyze the differential binding peaks between several mutant with WT for specific IP (for example, H3K4me1). But I was confused about how to choose the proper way to do the differential binding peaks analysis. What I know now, there exist two basic strategies I could use:
option 1. I use the DiffBind package to do that. The initial input file is bam files and a series of bed files generated by MACS2 call peaks (such as narrowPeak format). In this package, I will generate the normalization factors for the raw bam files before starting the differential binding analysis.
option 2. I could also use the strategy in the CSAW package to filter out the low abundance region and only keep the higher abundance region for multiple tests between different samples. (as shown in the second example mentioned in this paper(https://f1000research.com/articles/4-1080/v1). In this example, I do not need to do peak calling beforehand and I just by calculating the counts by setting the bin width as a narrow or wide range.
Based on these, I have several questions:

since I want to use TMM composition bias normalization to normalize all bam files for specific IP, if I choose the DiffBind package to get the differential binding peaks based on raw bam files, will it cause some argument since in both TMM composition bias normalization in CSAW and DiffBind will generate their own normalization factors.
is the peak calling by MACS2 is an essential step before the differential binding peak analysis? can I omit the peak calling step, just use the generate matrix generated by TMM composition bias normalization to do the differential binding peak analysis?
if choose the DiffBind to analyze differential binding peaks, for the bed files of two biological replicates, should I use the common regions or the pooling regions of two biological replicate to do the analysis?
If I pool the coordinate of peaks of all samples (macs2 callpeak) and calculate the count matrix based on the TMM composition bias-normalized bigwig files, is it possible to import the matrix directly into DiffBind to do the differential binding peak analysis?

I would appreciate it if someone would like to give me some suggestions on that! Many thanks.

Did you get the answers somewhere in the end?
What approaches did you use? Thanks and regards,
Kylie


text:
Synteny test and tools to use for mammalian organisms?
I am trying to do a synteny test for MHC genes in human and another mammalian organism
For this reason, I need to download all paralogs of rabbits from humans and then construct position the coordinates of paralog to test the synteny. However, I am unsure how to download the paralogs from Ensembl and how to perform this test?
My main aim is to research synteny between MHC genes of humans, rabbits and mice.
Any suggestions on how to do this or alternative bioinformatic methods on how to perform a synteny test would be highly appreciated!

Are you sure you mean paralogues (genes separated by a duplication event) and not orthologues (genes separated by a speciation event)? Please confirm this. I will separately answer your question assuming you meant orthologues.

This is not an easy ask. Firstly, very few species (I think only human and mouse) have alternative haplotypes annotated for MHC. Most species only have the primary assembly, only one possible haplotype of this locus. This means that the haplotype in the primary assembly is not necessarily the same haplotype as you'll find in another species. Since the MHC haplotypes can contain completely different genes, you will necessarily be able to find orthologues and synteny.
Secondly, we in Ensembl do not calculate homologues of the genes on the alternative loci, so we only have homologues of the genes that just happen to be on the primary assembly.
There is already pre-calculated synteny in Ensembl. Here's the human MHC locus compared to Rabbit. This includes all the homologous genes across this locus. However, this is just one arbitrary MHC vs another arbitrary MHC.


text:
Remove gaps from alignment?
I have an MSA (protein sequence) and I have used various programmes (Clustal, Aliview, MEGA11 etc) to align. However, in all programmes I get many gaps which is not ideal as I am trying to construct a tree from this alignment.
My question is are there any tools to delete gaps in alignments?
Your help would be greatly appreciated.

You may also find useful the BMGE tool (stands for Block Mapping and Gathering with Entropy). BMGE selects regions in a multiple sequence alignment that are suitable for phylogenetic inference by computing their entropy-like score. This allows to distinguish among biologically expected and unexpected variability for each aligned character. BMGE also implements different types of trimming (for instance, deletes gaps) to minimize phylogeny reconstruction artifacts.
Available from here and also installable from Bioconda.

If you want an ungapped alignment, you could do an ungapped alignment (the older read aligners were ungapped, e.g. bowtie). However, I suggest that you don't want an ungapped alignment, you just want a stricter one, here's a few options...

Remove/trim poorly aligning sequences. You could try GUIDANCE, which has an explicit model or GBLOCKS for a more simple filtering criteria
Tune your alignment parameters for ClustalOmega/MAFFT etc. to make it stricter
Finally, consider that your sequences must be homologous otherwise you'll trim away sequences until you have a random alignment that appears good

This new tool, CIAlign, talks about some automated removal of gaps in multiple sequence alignments
https://twitter.com/thekatybrown1/status/1306215132670889984
https://github.com/KatyBrown/CIAlign

Gblocks is a tool for removing gaps/poorly conserved positions from alignments. Available here. Here is a webserver that runs it for you.
It's a bit old, possibly there are others at this point.


text:
Is there a way to calculate the Ka/Ks value in R/unix?
I am investigating the rate of synonymous non synonymous substitutions and would like to find the Ka/Ks value also known as dn and ds between multiple alleles? I want to investigate multiple alleles because I am examining the MHC region which is highly polymorphic. Any suggestions would be very helpful!

I have also used the seqinr package for Ka/Ks

A while back, I used PAML for this purpose for a phylogenetics project. It has Windows and macOS executables but for *nix, you will need to compile it from source code.


text:
Differentiating molecules based on peptide sequence? How to annotate?
I want to differentiate between classical class I and non classical class I MHC molecules in a model organism using well conserved structural features within classical MHC I molecules (eg intradomain disulfide bridges within class I a1 domain). To do this I need to align the peptide sequences of transcripts of that model organism with transcripts of classical MHC I human genes and see which genes from those transcripts have those structural features. So I am looking for a tool to align and annotate those structural features in a replicable manner
Additionally if you have any other bioinformatic ideas /computational analyses on how to find classical MHC Class I molecules on an organism which is poorly annotated that would highly appreciated!

MAFFT is able to incorporate structural information in alignments: https://mafft.cbrc.jp/alignment/software/dash.html

Please consider the answer below and upvoting/accepting where appropriate.

Knowing almost nothing about the biology, I think that the most straightforward way to do this is:

train profile HMMs on your different classes of protein sequences (in human maybe, but ideally including MHCs from a large number of organisms)
translate your transcripts into proteins and run them against your trained profile HMMs
confirm that the profiles adequately discriminate the two classes and adjust if needed.

For an example of how to do this (on MHC), you can see here.
These approaches are then easily extensible to further organisms, as you have the profile HMMs.
It's possible that a simpler approach would consist of just finding hits against human MHC and then just checking for the existence of sequence patterns corresponding to those structural motifs, but I doubt that would work well in any significantly divergent organism.


text:
How can I speed up FASTA ssearch output?
I'm trying to predict targets for miRNA sequences, and so far after reading some papers for other tools that do this, I've found that I should be using FASTA's ssearch36 for this.
I'm starting with a small file of 500 miRNA sequences, searching a reference DB of 50k sequences. When using ssearch, I get the following output at the very bottom after ~85 minutes of running:
Total Scan time: 294.000 Total Display time: 4681.000
I am piping the output to a file, and the result file is ~13 Gb. It is taking this long even when I pipe the output to >nul so I don't think it's my disk speed.
I'm wondering if there's any way to speed this up, since it seems very unreasonable to me that the actual search itself only takes ~5 minutes but yet taking over an hour to simply generate output. Is it a bug or am I interpreting it incorrectly?

are you using the parallel or MPI version?  If so, you can add more threads with the -T option.  However, this probably won't help since the default is to use as many processors as are available


text:
How can I associate centrality score of genes(vertices) in a GRN with the probabilistic model in terms of cause and effect?
I'm learning to build Bayesian Network Based Modeling of Gene Regulatory Network of Cancerous Cellular Cells. So far, I have learned to form a Bayesian Network based on nature of genes(i.e Functional Genes and TFs). However, I cannot figure out how to give probability to each node/vertex. Is it possibile to somehow iterpret the centrality score of each vertex an associate the meaning with the probabilistic model in terms of cause and effect?
or can you atleast refer me a relevant research paper?
I believe I would be worth mentioning that I'm fairly new to this.
Here's my data of TFs and Sources.

TF
Target
Model of Regulation
References

ABL1
BCL2
Repression
11753601

AHR
FOS
Unknown
10478842

APC
BIRC5
Unknown
11751382

AR
ERBB2
Activation
21613411

AR
ERBB2
Unknown
21741601

Following is the sample of my data for network analysis, mostly containing centrality scores some graph properties.

SUID
name
BetweennessCentrality
ClosenessCentrality
ClusteringCoefficient
Eccentricity
EdgeCount
Indegree
IsSingleNode
NeighborhoodConnectivity
Outdegree
PartnerOfMultiEdgedNodePairs
selected
SelfLoops
shared name
Stress

7228
TP53
0.06611272
0.46456693
0.04141865
4
75
50
FALSE
14.40625
25
11
FALSE
0
TP53
39572

6788
AR
0.03605041
0.40410959
0.03164983
4
60
46
FALSE
13.38181818
14
5
FALSE
0
AR
20498

7029
MYC
0.03342125
0.35014837
0.03280053
5
81
74
FALSE
12.93589744
7
3
FALSE
0
MYC
19140


text:
Imputation of dog genotype
I am an undergrad working on a set of genotype matrix of dogs, it's in csv file format with information on only the each sample as rows, chromosome and its position as columns, and the rest of the matrix is filled with 0,1,2 and missing values.
So I am wondering what would be the best method to impute the genotype with all the data I get? I looked up Impute2, but it seems to require the reference data(recombination map) which I don't have...
Really appreciate whoever answers my question!!

I have never inputted genotype data, but I am pretty sure a recombination map for dogs exists. Check for example here: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5100850/

I've done imputation on dog data before - so I can hopefully help. However, I think we need some more information from you. How many individuals and SNPs do you have to begin with? Are the dogs the same breed or mixed?

Impute2 is a great, if dated choice. Another great choice is beagle. The advantage of beagle in this case is that it is not as strict about requiring a recombination map. Without a map, results will be poorer, but perhaps good enough.
Do make sure that your sample size including a reference panel is at least 200 in any case.


text:
Program for SOV calculation
Are there any programs (web server or installable) out there to calculate segment overlap score (SOV) for protein secondary structures?
Thanks

https://github.com/nanjiangshu/calSOV

Thanks a lot, that works!

Just googled it not sure it calculates the SOV in use today, seems SOV99 to me


text:
RNAseq small sample size in control group
I am doing an RNAseq experiment with affected (n=6) and control groups (n=6) in cow. However it turned out that 4 of my control samples have very low quality so I have to discard them.
My problem is that there is no way to repeat the sampling and sequencing.
Are there any solutions to perform DE,GSEA with these conditions (somehow generating bootstrapped samples from control groups or something similar)?
Thanks in advance!


text:
How to compare Log2 Fold Change values
I have 5 sets of Log2 Fold Change values pulled off DualSeqDB. Each of the 5 groups corresponds to a different microbial taxon. Within the 5 groups, there are many Log2 Fold Change values, each corresponding to a gene that was significantly affected by the presence of the microbe.
My question is: is there any way to get the average Log2 Fold Change for each group? How could I calculate that?
The end goal is to compare each microbe's effect on gene expression, hopefully to prove that the phenotype of the microbe is in some way related to how much it can change the gene expression in its host.
I have considered using p-values to compare each taxa's effect on gene expression, but because all the datasets are significant, the p-values I got using fisher's method are all 0, so I don't know how I could go ahead and compare that to categorical variables (like presence of a certain phenotype).
Thanks so much to anyone who can help!


text:
Nucleotide Frequency - HOMER?
I'm trying to find the frequency of UA dinucleotide in a bed file. I want to convert the bed file to FASTA, count all the UA dinucleotides up, then divide by the total length of the FASTA file.
Also, I probably need to compare the value I get above to the background, but I'm not sure how to go about this?
I understand homer should do the counting for me and also create a background distribution for me. It looks like I need to build a basic motif file and give homer a bedfile of locations and motif file to find. However, I'm unsure of how to set this up.
Other suggestions are welcome.


text:
How to manually curate a genome assembly for sequence variation or error?
I have a PacBio HiFi assembly of 1.1 Gb from a heterozygous species. I have aligned this assembly against a reference genome which is around 0.9 Gb. I can see that there are quite a few INDELs, Duplications, Inverted regions in my Query sequence on the Y-axis (PacBio Hifi) compared to the reference genome.
I want to check if these structural variations are real or only resulting from sequencing errors or not. So, I want to manually curate the PacBio assembly. Because I can see at which contig the structural variations are.
Can you please suggest any tools for manual curation of a genome assembly to investigate whether these structural variations are sequencing errors or real sequence variation?
I have uploaded the dot-plot below. Thank you.

based on your other question, I would suggest mapping the hi-c reads to this assembly and manually polishing in juicebox

So, you mean take the HiFi fasta and map the Hi-C fasta on it? Or do I need raw reads from the Hi-C reads?

You would need to map the raw Hi-C reads. the FASTA is the outcome of a lot of procedures that remove information. It is the Hi-C reads have the actual proximity information about which regions are close to each other on the chromosome. You can use that to correct structural errors (if any). If those are truly errors then they will look weird on the Hi-C heatmap and you can correct them manually in juicebox: https://aidenlab.gitbook.io/juicebox/assembly-tools

Ah ok. Thanks for clarifying the confusion. It is clear now.


text:
How to retrieve fasta sequence after local blast?
I have created a Blast database using a reference genome. Then, I have performed a local blast search in command line using a gene of interest. I have obtained some hits with the usual Blasting information. Now, I want to extract the exact matching sequence from the blast search with their corresponding start and stop position. I have a .gff file for the genes as well for the reference genome. Can you give me any suggestions on how to do that? Thanks a lot.

Please [edit] your question and add an example of the relevant files so we can understand what you need. Do you have multiple HSPs you want to combine? Do you need to model splicing? If you already have the gff of the target sequence, why were you blasting?

FWIW, In the Sequenceserver BLAST software, you can just click on "download alignment" to extract the aligning region.

My usual approach for extracting sequence from a reference genome is to use samtools faidx:
# create genome index; slow; only do once per reference genome
samtools faidx reference_genome.fa
# extract sequence for a specific match
samtools faidx reference_genome.fa chromosome_name1:posStart1-posEnd1 > match1.fa
samtools faidx reference_genome.fa chromosome_name2:posStart2-posEnd2 > match2.fa

You haven't mentioned how many BLAST results you have, or what format they're in, so I'm not sure if it's practical to do this manually, or if it needs to be automated in some way.


text:
How to run MD5 check for multiple fastq files in different subdirectories?
I have received Illumina sequencing reads for 100 samples. I have 8 R1.fastq.gz and 8 R2.fastq.gz files for each sample in each subfolder. I want to run an MD5 check for all the fastq files in each subfolder.
My folder structure looks like:
/mydata/sequencing/clean
/mydata/sequencing/clean/Sample1
/mydata/sequencing/clean/Sample2
......
/mydata/sequencing/clean/Sample100

I am using the following code being on the /mydata/sequencing/clean directory.
find . -type f -exec md5sum-lite {} \;

This code is printing the md5 result on Terminal.
How do I get these results printed to a .txt file? Besides, is there a more elegant way to do this MD5 check?
Thank you.

You posted the same thing in Biostars.org :) https://www.biostars.org/p/9509386/

Haha yes. Because sometimes I do not get a reply to my questions (may be they are  too naive some times) :(

Glad you got the answer.

thats because its on the wrong stack

Anyways it is always great to receive feedback on bioinformatics problems. Nice to see that many roads lead to Rome.

Small modification,
find . -type f -exec md5sum-lite {} \; >> output.txt

If you use > the output will almost certainly contain a single record.

Thank you for the prompt reply. Is there any other code that you know of for looping the md5 command over different subfolders?

Yes of course ... there are loads of ways to do this. A bash script is good (separate question IMO). However, the > and >> will still be central. In terms of elegance, your one liner looks cool. BTW >> is very likely the solution because of the find approach. One liners are cool BTW.

Thank you very much.

Your welcome. Your  second question (alternative code) is worth asking because of the variety of ways it can answered, e.g. a bash for loop, Python subprocessing, Perl one liner are the  approaches I can think of.

Yeah. totally agreed.


text:
How to get KEGG pathway names from a list of 'ko' IDs?
I have a list of ko IDs from which I want to get the pathway names. How can I get that? I have tried this link, but it does not give me any result. What to do?
  ko00010   ko00030 ko00061 ko00120 ko00121 ko00190 ko00195 ko00230 ko00250 ko00253 ko00270 ko00300 ko00330 ko00400 ko00401 ko00410 ko00430 ko00450 ko00460 ko00471 ko00472 ko00500 ko00511 ko00520 ko00521 ko00523 ko00531 ko00550 ko00562 ko00604 ko00620 ko00625 ko00640

Thanks

You can use the KEGG API:
ko_list="ko00010 ko00030 ko00061 ko00120 ko00121
         ko00190 ko00195 ko00230 ko00250 ko00253
         ko00270 ko00300 ko00330 ko00400 ko00401
         ko00410 ko00430 ko00450 ko00460 ko00471
         ko00472 ko00500 ko00511 ko00520 ko00521
         ko00523 ko00531 ko00550 ko00562 ko00604
         ko00620 ko00625 ko00640"

for i in ${ko_list}
do
    hit=$(curl -# http://rest.kegg.jp/list/pathway | grep -w "${i/ko/map}")
    echo -e "${i}\t${hit}" | cut -f1,3 >> ko_pathways.txt
done

Or if your list is too long you can change the input and use a file with a list of ids.
while read line
do
    hit=$(curl -# http://rest.kegg.jp/list/pathway | grep -w "${line/ko/map}")
    echo -e "${line}\t${hit}" | cut -f1,3 >> ko_pathways.txt
done < ko_list.txt

The result will be something like:
 ko00010    Glycolysis / Gluconeogenesis
 ko00030    Pentose phosphate pathway
 ko00061    Fatty acid biosynthesis
 ko00120    Primary bile acid biosynthesis
 ko00121    Secondary bile acid biosynthesis

Thanks a lot Zorbax, for your answer. It works. Besides that, is there any web server or something that can be used? Just wanted to know.

Hi Zorbax, I have a similar query. Is there a similar way I can use to get KEGG Module names from a list of Module IDs (e.g. M00002)? Many thanks

I don't know anything about a web server. And yes, to your second question, I think you can use the API request for the module instead pathway list.

This doesn't work anymore

Yes, it does. Use https instead http


text:
How to solve Nextflow error: "Trace file already exists"?
When trying to run epi2me-labs/wf-artic, I get the following error:
❯ nextflow run epi2me-labs/wf-artic \
 -r v0.3.18 --fastq  ~/Downloads/barcode95.fastq.gz \
 --scheme_version Midnight-ONT/V3

N E X T F L O W  ~  version 22.10.1
Launching `https://github.com/epi2me-labs/wf-artic` [ecstatic_marconi] DSL2 - revision: a60a1e1e73 [v0.3.18]
Trace file already exists: /Users/corneliusromer/Downloads/BJGMC_ONT_Run71/output/execution/trace.txt

It seems that nextflow doesn't like to rerun after it failed. Is there a convenient way to force nextflow to rerun, in a similar way to snakemake --force-rerun? Or do I need to manually delete the trace file every time? That works, but it's tedious.

Rotation of report files (trace, execution, timeline etc) was removed recently (since version 22.10.0 in a762ed5). To overwrite the existing trace file you will need to add the following to your Nextflow config:
trace.overwrite = true

See: https://github.com/nextflow-io/nextflow/issues/3317


text:
Where does Snakemake send output that's generated by shell commands?
I'm wondering how Snakemake deals with output to stdout and stderr that is generated when running a rule with a shell: command.
I ran the following rule and only got informed that there was a non-zero exit code. However, I did not get a stack trace. To get access to the stack trace I had to run the shell command manually:
# Shortened for clarity
rule parse:
    input:
        sequences = rules.join_downloads.output
    output:
        sequences = "data/sequences.fasta"
    shell:
        """
        augur parse
        """

Why is that? And can I somehow make Snakemake output more detailed information when something goes wrong to help with debugging?

You usually define a log file like below, and redirect the stdout and/or stderr to it. Stdout is >, stderr is 2> and both of them combined is &>.
rule parse:
    input:
        sequences = rules.join_downloads.output
    output:
        sequences = "data/sequences.fasta"
    log:
        my_output/my_log.log
    shell:
        """
        augur parse 2> {log}
        """