text:
Missing ',' in line when Biopython reads a nexus tree
I want to edit a tree that I got from BEAST2 treeannotator in nexus-format.
Usually I use the module Phylo from Biopython for such work but Phylo.read(r"filename.tree", "nexus") gave me the next exception:
---------------------------------------------------------------------------
NexusError                                Traceback (most recent call last)
Input In [29], in <cell line: 1>()
----> 1 Phylo.read(r"filename.tree", "nexus")

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:60, in read(file, format, **kwargs)
     58 try:
     59     tree_gen = parse(file, format, **kwargs)
---> 60     tree = next(tree_gen)
     61 except StopIteration:
     62     raise ValueError("There are no trees in this file.") from None

File ~\miniconda3\lib\site-packages\Bio\Phylo\_io.py:49, in parse(file, format, **kwargs)
     34 """Parse a file iteratively, and yield each of the trees it contains.
     35
     36 If a file only contains one tree, this still returns an iterable object that
   (...)
     46
     47 """
     48 with File.as_handle(file) as fp:
---> 49     yield from getattr(supported_formats[format], "parse")(fp, **kwargs)

File ~\miniconda3\lib\site-packages\Bio\Phylo\NexusIO.py:40, in parse(handle)
     32 def parse(handle):
     33     """Parse the trees in a Nexus file.
     34
     35     Uses the old Nexus.Trees parser to extract the trees, converts them back to
   (...)
     38     eventually change Nexus to use the new NewickIO parser directly.)
     39     """
---> 40     nex = Nexus.Nexus(handle)
     42     # NB: Once Nexus.Trees is modified to use Tree.Newick objects, do this:
     43     # return iter(nex.trees)
     44     # Until then, convert the Nexus.Trees.Tree object hierarchy:
     45     def node2clade(nxtree, node):

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:668, in Nexus.__init__(self, input)
    665 self.options["gapmode"] = "missing"
    667 if input:
--> 668     self.read(input)
    669 else:
    670     self.read(DEFAULTNEXUS)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:718, in Nexus.read(self, input)
    716     break
    717 if title in KNOWN_NEXUS_BLOCKS:
--> 718     self._parse_nexus_block(title, contents)
    719 else:
    720     self._unknown_nexus_block(title, contents)

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:759, in Nexus._parse_nexus_block(self, title, contents)
    757 for line in block.commandlines:
    758     try:
--> 759         getattr(self, "_" + line.command)(line.options)
    760     except AttributeError:
    761         raise NexusError("Unknown command: %s " % line.command) from None

File ~\miniconda3\lib\site-packages\Bio\Nexus\Nexus.py:1144, in Nexus._translate(self, options)
   1142         break
   1143     elif c != ",":
-> 1144         raise NexusError("Missing ',' in line %s." % options)
   1145 except NexusError:
   1146     raise

NexusError: Missing ',' in line 1 AB298157.1_2015_-7.9133750332192605_114.8086828279248, 2 AB298158.1_2007_-8.41698974207…

Using Nexus.read(Nexus(), input=r"filename.tree") gave the same result. Please could anyone help with this? I cannot understand the reason of this error because nexus file looks correct.

Nexus format
A NEXUS tree file is a rigid but very useful format which can be used to store a phylogenetic tree, along with associated metadata, for example, the alignment.
It will always comprise a:

TAXA block

The TAXA block contains information about taxa.
... can comprise a,

DATA block

The DATA block contains the data matrix (e.g. sequence alignment).
... in this context importantly will have a,

TREES block

The TREES block contains trees described via the format, e.g. ((A,B),C);:
Each block will open with a BEGIN statement and terminate with an END statement. Again its versatility to store specific information within a block is valuable property of the format.

NB! I encountered some problems with my following work. So I would like to warn readers about them.
It is allowed in the newick format to surround the label with quotes, & some programmes or scripts add them to those names that have ambiguous characters. But it can lead to exceptions during the following phylogenetic analysis, for instance, in BEAST. I wish you would be careful with this.

That's a complicated output.
Two solutions, the easiest is simply to use Andrew Rambaut's figtree http://tree.bio.ed.ac.uk/software/figtree/ . Figtree is a very, very good tree drawing program and it will work; it works on Mac, Windows and Linux.
You can then export the file from Figure - from memory in numerous formats, e.g. phylip format, and import it into Biopython if that is where it is needed. If you are parsing Beast output, Beast has A LOT of internal tools and you might be attempting to replace those via Biopython. I would avoid it. Also there is the package for processing Beast output (which I forget), I think its called Tracer.
With Biopython try parse, it might help if there are multiple trees:
from Bio import Phylo
trees = Phylo.parse("examples.nex", "nexus")
for tree in trees:
    print(tree.name)

I would avoid the 'r'. I agree it makes no difference, but really you'd need to be using unusual file names for 'r' to be needed. Its possible you are using it to get around \s, i.e. spaces in the file name. You could just think about simply not including \ and space type characters in file names. In regex it is used for certain but that's because its 'inescapable' (no pun intended, but I thought it was funny).

Just to wrap up and clarify, I recommended:

To open the tree in Figtree;
Find the export menu and output as either 'newick' or preferably 'nexus' or both;
Re-run Biopython importing as 'nexus' or 'newick' (if nexus fails).

The OP reported ...

I have just compared them (nexus->FigTree->nexus).

Good so it opened up in FigTree and exported, this was expected.

Fig Tree compose all this actually in newick & stash it in nexus block decor. While in the previous case newick tree contains IDs, & there is a seclude block with ID-name pairs

I missed this bit, "newick stashed in a nexus block". That is potentially really cool because I suspect that is within the TREES block. If that's correct then that is nexus format. This should now open Biopython.
Why the bug?
One of the issues with Nexus is there are two variants regarding how the TAXA block relates to the TREES block. In one format they recode the taxa - I definitely understand why. Recoding is possibly (very likely) be being used in Beast to increase algorithm efficiency - it makes no sense slushing a load of taxa names around a convergence algorithm. Writing all that stuff out is not cool either.
However, when you are looking at a tree that's not what you want at all - efficiency is an irrelevance and readability is important. Thus recoding the taxa labels is definitely not wanted. This is a possible explanation for the read error.
Whatever the reason for the bug I am certain Biopython has a work around because in phylip alignment format this is certainly true. There is a custom module that handles the specific issues the format generates.
It was fun but time to move on.

Could I use Figtree algorithms through command line? If it has to be converted to another format, I would like to do this automatically.

My advice is to get it working, even if point and click, - it will work then compare the input formats to understand what the issue is. You can then script the parsing.

So do you suggest to convert nexus to nexus & compare the difference?

I have just compared them (nexus->FigTree->nexus). Fig Tree compose all this actually in newick & stash it in nexus block decor. While in the previous case newick tree contains IDs, & there is a seclude block with ID-name pairs.

Good, should work. I've rewritten the answer to accommodate this.


text:
How do I use bash for loops to loop through a WDL Array?
Within the example below a WDL array of strings, test is defined as "a","b","c","d" the value total is a variable that signifies the length of the WDL array. Using the bash for loop, which will iterate from 1 up to length of the array, I seek to print out the values at each indicy in the array and output it to another array. For example when the loop is on 1 the value at indicy 1 within the array will be printed out in this case "a" and stored in the array out, if the loop is on 2 then the value "b" will be printed out and stored in the array out, etc. I tried the below example, however, I get an error when running. Can anyone help me with this problem?
How I tried to do this in WDL:
task A {
    input {
        Array[String] test = ["a","b","c","d"]
        Int total = length(test)
    }

    command {
       for (( c = 1; c <= ${total}; c++ ))
       do
           echo ${test[$c]}
       done
    }
    output {
        Array[String] out = read_lines(stdout())
    }
}

I suspect the issue here is that you are trying to use a bash variable inside of a WDL placeholder. The ${} are interpreted by WDL as places where some value should should be inserted into the command. So when you say echo ${test[$c]}, test[$c] is seen a WDL expression. This isn't a valid WDL expression, causing a parsing error.
The choice to use ${} as placeholder syntax in WDL was a rather inconvenient decision by the WDL team, because it clashes with the bash syntax. There is an alternative syntax available for the command section in WDL 1.0, which is preferred for this very reason. This alternative uses ~{} as placeholder syntax, instead of ${}, and wraps the command in <<< >>> (referred to as the "heredoc-like syntax"), rather than the usual { }. In this heredoc-like syntax the ${} isn't seen as a WDL placeholder and can, thus, be used in bash commands like normal. You can read about it in the WDL SPEC.
Using this syntax your WDL would look something like this:
version 1.0

workflow test {
    call A
}

task A {
    input {
        Array[String] test = ["a","b","c","d"]
        Int total = length(test)
    }

    command <<<
       ARRAY=(~{sep=" " test}) # Load array into bash variable
       for (( c = 0; c < ~{total}; c++ )) # bash array are 0-indexed ;)
       do
           echo ${ARRAY[$c]}
       done
    >>>

    output {
        Array[String] out = read_lines(stdout())
    }
}

Note that if the values in your array are more complex and might include characters with special meaning in bash (like spaces), it would probably be better to wrap the values in single quotes like so:
ARRAY=('~{sep="' '" test}')
And just in case you weren't aware, bash for-loops can also loop over arrays directly:
command <<<
       for c in ~{sep=" " test}
       do
           echo $c
       done
>>>

This is a comprehensive answer, giving a clear insight into WDL and its integration with bash.

Thank you for this. In the first answer you gave, is it possible to load the array into a bash variable i.e. ARRAY=(~{sep=" " test}) within the workflow definition so it would like the bellow example? Can I then take that as an input into my task definition?
workflow test {
Array[String] test = ["a","b","c","d"]
Int total = length(test)
ARRAY=(~{sep=" " test})
call A  {
      inputs:  testarray =  ARRAY
 }

}
task A {
inputs {
Array[String] = testarray
}
}

Currently, WDL doesn't have a way of defining environment variables in WDL itself. You have to do so in the command section. There has been some discussion on introducing such a feature though: https://github.com/openwdl/wdl/discussions/458

source: Stack Exchange

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

sort: votes


text:
Target genes for piRNA
Where I can find a database or tool to give me the target genes of PIWI (piRNAs) in human?
I found one but works for worm like piRTarBase
There are a lot for miRNAs around but not for piRNAs

Update
It does indeed appear that while piRBase etc. have human piRNAs, they do not have target information.
Looking at related questions, recommendations are made to map piRNA sequences to the genome to find the targets. One tool recommended is PiPipes, which might assume read input but otherwise might give some ideas.
Original answer
piRBase has this. See "target mRNA" field.

Thank you. Can I use that for human? because all are for non human likely

@Exhausted You are right- how odd. THey have human piRNAs in the db but none of them seem to have targets. I missed that, sorry. I will see if I can figure it out and if

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


text:
What can be the reason of getting negative branches lengths after BEAST analysis?
BEAST2 is currently being used for tree reconstruction prior phylogeographic analysis. The sample size and loci are described below.
I thought that BEAST/BEAST2 does not allow negative lengths of branches on prior distribution level so I was surprised when I have looked at my consensus tree using FigTree.

Question What could be the reason of this and how to wrestle with this phenomenon?
I am analysing about 1500 taxa of HAV (subtype: IA, region: VP1) under the next parameters:

GTR+Г,
Strict Clock,
Coalescent constant population.

The tree was generated by TreeAnnotator through
~/beast/bin/treeannotator -heights median -burnin 20 -limit 0.5 VP1_test_bt_ExpPop.trees VP1_test_bt_ExpPop.tree

Is this an issue of the tree that was generated or an issue of the algorithm that makes the layout for plotting the tree?

Not either/or @pallie. The tree drawing program should be excluded.

What I mean is, does the topology of the tree as saved on disk match the topology of the tree drawn in the graph? @M_

@Pallie It is the problem of the tree. For instance, =0.008838086259721895]:0.6701304093024589)[&posterior=0.005774216356351638]:-0.3146123701838377,(((26[&height=7.0,height_95%_HPD=

@M__ 20 is not a number. This is percentage.

I found the answer on ResearchGate from Santiago Sanchez-Ramirez & let me quote it here:

Try using the option -heights ca when you run TreeAnnotator. This stands for "common ancestor trees" and aims at summarizing clade ages across all posterior trees and not only the values for subset of trees that have that clade. If that subset of trees is low in frequency the average node height might end up being older than the direct ancestor.
Here is the explanation for negative branches:
"MCC trees produced by TreeAnnotator can have a descendent node that is older than its direct ancestor (a negative branch length). This may seem like an error but is actually the correct behaviour. The MCC tree is, by default, generated with average node heights across all trees in the sample which contain that clade. The negative branch lengths result when a clade is at low frequency and tends not to occur in those trees that have the MCC tree's ancestral clade (or vice versa). This means the average heights are for the adjacent nodes are derived from different sets of trees and may not have any direct ancestor-descendent relationship."
In this paper you will find more information about the ca option: Looking for trees in the forest: Summary tree from posterior samples

I am avoiding directly answering the question because in summary, I feel, this given issue is akin to project supervision and is simply too much time.
One issue raised however is something called 'burnin'.
Burnin is foundational to MCMC.
In the comments I stated the burnin is too low. The burnin used ...
~/beast/bin/treeannotator -heights median -burnin 20 ...

The OP believes the burnin represents a percentage, i.e. 20%.
Citing the documentation on TreeAnnotator at https://beast.community/treeannotator

BurnIn (as trees)
In this case you just specify the actual number of
trees you wish to discard rather than the number of states.

In this given case, considering this is 1500 taxa, I would estimate the actual burnin  used in comparison to an initial assessment of completion time (technically called 'convergence') to be 0.4% maximum, i.e. well over 0.5 log smaller than the burnin believed to be used.
It is easy to assess empirically e.g. comparing wc infile and wc outfile, where infile and outfield are MCMC treefiles.
A standard burnin is 10% of the total MCMC. The threshold would be independently assessed and if inadequate a higher burnin would be used, of which 20% is reasonable.
To explain the documentation, in an MCMC only a sample of the trees are written to disk (often 1 in 1000, but it can be anything) and the burnin is 10% of the total (unsampled).
Thus factoring the sampling is important in calculating the precise burnin in any MCMC output analysis. Otherwise the burnin that was intended, isn't the burnin used in practice.
Perhaps you can now see my point?  Not the specific technical matter above, but the issue concerning complexity and time commitment.

To address the matter concerning Beast1 vs Beast2 raised as a comment below, i.e. the proposal the original MCMC was performed in Beast2 and Beast2 is exempt. The issue under discussion is TreeAnnotator.
The TreeAnnotator section of the Beast2 documentation is here https://www.beast2.org/treeannotator/

TREEANNOTATOR ....

BurnIn
This option allows you to select the amount of burn-in, i.e., the
number of samples that will be discarded at the start of the run, so
that you are only analysing the part of the trace that is in
equilibrium.

Again the Beast2 documentation is clear, 'the number of samples', not the percentage or proportion. Furthermore note it states 'samples', not 'states', i.e. the number of trees resulting from the sample frequency that was assigned.
A burnin can be performed prior to TreeAnnotator analysis and that is good practice, albeit not essential. This might be where the report output the OP cites was obtained. Furthermore, if this was done here upstream, then no burnin is required for TreeAnnotator. Moreover if a correct upstream burnin was performed, burnin should either be; avoided for TreeAnnotator, or else the total combined burnin calculated and correctly reported. The final output ain't just a tree, there's a lot of important analysis/ 'meta data'.
Note  I agree with the OP's second comment below this answer using information from the comment below the question: a value of -0.31 can't be ignored, i.e. its a notable value for any branch length and very large in the given context.
Its definitely time to move on
OP upvoted, good luck with the project and all the best.

This is correct for BEAST 1, if I memorise this well. In BEAST 2 the parameter -burnin or -b is percentage. Because it is specified in the -help message & you cannot specify any number above 100 without an exception. Moreover, there are numbers of treated & total trees in stdout of treeannotator. For example, from my another analysis with 10% burning: Processing 18669 trees from file after ignoring first 10% = 2074 trees..

So the question is about possible common reasons of such phenomenon.

Dear @Vovin I have addressed the first comment you raised above as an append to previous answer below the separator.

I have addressed the second point you raised in the section called Note

This is not productive for seeking the solution of the question but I think that I should point this out so others do not get confused by your answer. This is certainly percentage & the fact is based on the next evidence: the help message of the programme where it is named as percentage, the log message where the number is equal to my burnin times the total trees number, & the source code where this variable is called as burninPercentage & used in the appropriate way.

Moreover, for instance, treeannotator -heights median -burnin 101 -limit 1 example.trees example.tree gives beast.app.util.Arguments$ArgumentException: Argument, -burnin has a bad integer value: 101. But it works with 99. The example from exeption message: Example: treeannotator -burnin 10 -heights mean test.trees out.txt. It is obvious that this is not a number because it would be too small. I was myself confused in the past by the blog article that you cited but it works as percentage.

Time to move on


text:
What is coarse-grained modeling of proteins?
I need beginner-level study materials for coarse-grained modeling of proteins.
Can anyone supply me with entry level resources besides the Wikipedia article?

Background
Coarse grain protein modelling could be very broadly defined as an approximation of all-atom protein models and find current application in folding studies but are early approaches extend back to 1975.
Citing a recent review where comments in brackets are for elucidation and are not part of the citation.

Native structures (proteins) are not completely fixed, but they change when
proteins perform their biological function, interact with other
biomacromolecules, or undergo unfolding–folding transitions (for example conformational changes).
Computational modeling of these processes is crucial for ... (examples truncated). Classical atom-level
molecular modeling can address many of the these tasks, but its
practical applications are still limited by its algorithmic efficiency
and the available computing power.

Kmiecik et al "Coarse-Grained Protein Models and Their Applications"  Chem. Rev. 2016, 116, 14, 7898–7936
https://doi.org/10.1021/acs.chemrev.6b00163

http://cgmartini.nl/index.php/about


text:
Running a list of apis via Python's request module
I am running the following Python code with python requests in my Jupiter notebook to get cancer oncology id's from EBI-EMBL.
import requests
cancer_all = requests.get('https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000376/children')
obi = cancer_all.json()

The above code is for one link. However I want to get results for 30 such links together.
Other examples of links:
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000621/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000637/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000653/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002402/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002516/children
https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002813/children

Questions

Is there any way I can run all these api in one go?
Could I fetch the results (in json format) of these links in a single download. Thus there would be one json file for all the links?

import requests
import json #??

apiString = "https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000621/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000637/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000653/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002402/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002516/children\nhttps://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0002813/children"

def requestIt():
    cancer_all = requests.get('https://www.ebi.ac.uk/ols/api/ontologies/efo/terms/http%253A%252F%252Fpurl.obolibrary.org%252Fobo%252FMONDO_0000376/children')
    return cancer_all.json()

if __name__ == '__main__':
callList = split('\n')
obi = [requestIt(x) for x in callList]

# e.g.
print (obi[0]

# an alternative approach to list comprehension
# computationally faster but will not make much difference here
obi = list(map(requestIt, callList)

Note If you want to parallelise it, which might be what you meant by 'all at once', you instantiate Pool() from multiprocessor via map and use the same comprehension but this time the function requestIt will spread this over your available CPUs. Thats a little bit more complicated - namely getting the result back.
If the apis are in a file they will load into a list nicely:
with open ('apis.txt', 'r') as fin:
    callList = fin.readlines()

However you want to combine the json into a single json ... you import json then use json.loads to convert the json into a dictionary, combine the dictionaries and then use json.dumps to convert it back to json. The code for that is a separate question altogether and the comprehension might be better expressed as a dictionary comprehension (not sure).


text:
Can I use Aho-Corasick for epidemic spread simulation?
I found that with this string-finding algorithm, this algo can determine if a string(person) holds something(if he has sick, already got vaccinated, and et cetera.).
I am trying to simulate an epidemic spread with simulation parameters and i can't find a simple mathematical method or algorithm that can help me to understand how i will do it with the simulation process.
ps: its my first time to ask a question here so dont judge if i cant follow the format..

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

The problem with the question is there isn't sufficient information to understand for example how you simulated the epidemic, spatial models, temporal, incorporating immunity? Is it the traditional R0 style deterministic models etc ... Thus it is not clear where Aho-Corasick would fit here. At a guess I can't see how it could drive a simulation.


text:
BAM files with no RNAME and POS, how to map contents to SNPs?
I have a set of 4 .bam files containing the exome of an individual, around 400 MB each. I used samtools to generate a 2.4 GB .sam file out of one of the .bam files, and I found it contains lines with 101 base sequences but seemingly no RNAME or POS information on them. The first columns of the first lines in the read like this (I won't include here the SEQ, QUAL and optional fields):
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:18818    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:18818    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:32033    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:32033    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:12712    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:12712    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:21261    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:21261    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:34882    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:34882    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:07764    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:07764    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:11569    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:11569    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:26537    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:26537    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:21010    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:21010    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:27085    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:27085    141 *   0   0   *   *   0   0

Most fields are nullified. As far as I can understand, subsequent lines with the 77 and 141 flags form corresponding pairs of bases and the QNAME field first element is an individual identifier, as it won't change through the head and tail of the file. Same goes for the second element (1). The last three elements change and might hold the information I need, but I'm not sure how to read them.
What I want to do is to map bases to SNPs (or a nucleotide identifier like NC_000006.11).
I have no real experience with bioinformatics besides enriching a file from MyHeritage (which contains around 609,000 base pairs with a corresponding rSNP number and chromossome) and crossed it with a reference human genome database I downloaded from NBCI. As it's a personal project, I somewhat validated the final result by correctly spotting a few known mutations/polymorphisms for that individual tested through usual lab tests, like the one output by the code I've wrote:
 RS1799945 (Chromosome 6 ): CG (NC_000006.11:g.26091179C>G, snv), https://www.ncbi.nlm.nih.gov/snp/rs1799945
          Reference Allele: C
                 Frequency:  9.6% (0.09631744626441) (~ 1:10)
              Loci (Genes): HFE  (homeostatic iron regulator)
                          : LOC108783645  (HFE antisense RNA)
        Detected Mutations: c.187C>G
             Variation URL: https://www.ncbi.nlm.nih.gov/clinvar/variation/10/
                  Diseases: not provided (conflicting-interpretations-of-pathogenicity), None (other), Hereditary hemochromatosis (pathogenic), Cystic fibrosis (risk-factor), Alzheimer disease (pathogenic), Familial porphyria cutanea tarda (None), Hemochromatosis type 1 (None), Microvascular complications of diabetes 7 (None), Transferrin serum leve, Microvascular complications of diabetes 7 (risk-factor), Hemochromatosis type 1 (pathogenic), Hemochromatosis type 2 (pathogenic)
                    PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8696333, https://www.ncbi.nlm.nih.gov/pubmed/8896550, https://www.ncbi.nlm.nih.gov/pubmed/8943161, https://www.ncbi.nlm.nih.gov/pubmed/9106528, https://www.ncbi.nlm.nih.gov/pubmed/9326341, https://www.ncbi.nlm.nih.gov/pubmed/9462220, https://www.ncbi.nlm.nih.gov/pubmed/10673304, https://www.ncbi.nlm.nih.gov/pubmed/11358905, https://www.ncbi.nlm.nih.gov/pubmed/11423500, https://www.ncbi.nlm.nih.gov/pubmed/11812557, https://www.ncbi.nlm.nih.gov/pubmed/11874997, https://www.ncbi.nlm.nih.gov/pubmed/11904676, https://www.ncbi.nlm.nih.gov/pubmed/12429850, https://www.ncbi.nlm.nih.gov/pubmed/14673107, https://www.ncbi.nlm.nih.gov/pubmed/15347835, https://www.ncbi.nlm.nih.gov/pubmed/15858186, https://www.ncbi.nlm.nih.gov/pubmed/18194558, https://www.ncbi.nlm.nih.gov/pubmed/18199861, https://www.ncbi.nlm.nih.gov/pubmed/18566337, https://www.ncbi.nlm.nih.gov/pubmed/18603647, https://www.ncbi.nlm.nih.gov/pubmed/18795173, https://www.ncbi.nlm.nih.gov/pubmed/19084217, https://www.ncbi.nlm.nih.gov/pubmed/19159930, https://www.ncbi.nlm.nih.gov/pubmed/19165391, https://www.ncbi.nlm.nih.gov/pubmed/19237423, https://www.ncbi.nlm.nih.gov/pubmed/19401444, https://www.ncbi.nlm.nih.gov/pubmed/19554541, https://www.ncbi.nlm.nih.gov/pubmed/19862010, https://www.ncbi.nlm.nih.gov/pubmed/19879291, https://www.ncbi.nlm.nih.gov/pubmed/19884647, https://www.ncbi.nlm.nih.gov/pubmed/20110814, https://www.ncbi.nlm.nih.gov/pubmed/20301613, https://www.ncbi.nlm.nih.gov/pubmed/20435227, https://www.ncbi.nlm.nih.gov/pubmed/20556870, https://www.ncbi.nlm.nih.gov/pubmed/20659343, https://www.ncbi.nlm.nih.gov/pubmed/20846924, https://www.ncbi.nlm.nih.gov/pubmed/20876667, https://www.ncbi.nlm.nih.gov/pubmed/21138591, https://www.ncbi.nlm.nih.gov/pubmed/21149283, https://www.ncbi.nlm.nih.gov/pubmed/21208937, https://www.ncbi.nlm.nih.gov/pubmed/21540342, https://www.ncbi.nlm.nih.gov/pubmed/21553046, https://www.ncbi.nlm.nih.gov/pubmed/21679129, https://www.ncbi.nlm.nih.gov/pubmed/21909115, https://www.ncbi.nlm.nih.gov/pubmed/21978626, https://www.ncbi.nlm.nih.gov/pubmed/22074419, https://www.ncbi.nlm.nih.gov/pubmed/22100073, https://www.ncbi.nlm.nih.gov/pubmed/22125427, https://www.ncbi.nlm.nih.gov/pubmed/22232660, https://www.ncbi.nlm.nih.gov/pubmed/22408404, https://www.ncbi.nlm.nih.gov/pubmed/22611049, https://www.ncbi.nlm.nih.gov/pubmed/22735619, https://www.ncbi.nlm.nih.gov/pubmed/22815867, https://www.ncbi.nlm.nih.gov/pubmed/23389292, https://www.ncbi.nlm.nih.gov/pubmed/23446634, https://www.ncbi.nlm.nih.gov/pubmed/23792061, https://www.ncbi.nlm.nih.gov/pubmed/23935582, https://www.ncbi.nlm.nih.gov/pubmed/23996192, https://www.ncbi.nlm.nih.gov/pubmed/24033266, https://www.ncbi.nlm.nih.gov/pubmed/24121126, https://www.ncbi.nlm.nih.gov/pubmed/24398642, https://www.ncbi.nlm.nih.gov/pubmed/24755831, https://www.ncbi.nlm.nih.gov/pubmed/24798971, https://www.ncbi.nlm.nih.gov/pubmed/24825461, https://www.ncbi.nlm.nih.gov/pubmed/25071582, https://www.ncbi.nlm.nih.gov/pubmed/25112650, https://www.ncbi.nlm.nih.gov/pubmed/25287020, https://www.ncbi.nlm.nih.gov/pubmed/25741868, https://www.ncbi.nlm.nih.gov/pubmed/26024779, https://www.ncbi.nlm.nih.gov/pubmed/26025379, https://www.ncbi.nlm.nih.gov/pubmed/26159428, https://www.ncbi.nlm.nih.gov/pubmed/26365338, https://www.ncbi.nlm.nih.gov/pubmed/26416403, https://www.ncbi.nlm.nih.gov/pubmed/26597663, https://www.ncbi.nlm.nih.gov/pubmed/26695521, https://www.ncbi.nlm.nih.gov/pubmed/26716680, https://www.ncbi.nlm.nih.gov/pubmed/26833210, https://www.ncbi.nlm.nih.gov/pubmed/27221532, https://www.ncbi.nlm.nih.gov/pubmed/27255824, https://www.ncbi.nlm.nih.gov/pubmed/27280446, https://www.ncbi.nlm.nih.gov/pubmed/27317329, https://www.ncbi.nlm.nih.gov/pubmed/27437086, https://www.ncbi.nlm.nih.gov/pubmed/27547017, https://www.ncbi.nlm.nih.gov/pubmed/27584680, https://www.ncbi.nlm.nih.gov/pubmed/27657935, https://www.ncbi.nlm.nih.gov/pubmed/27661980, https://www.ncbi.nlm.nih.gov/pubmed/28151915, https://www.ncbi.nlm.nih.gov/pubmed/28684612, https://www.ncbi.nlm.nih.gov/pubmed/28820869, https://www.ncbi.nlm.nih.gov/pubmed/29167213, https://www.ncbi.nlm.nih.gov/pubmed/29201641, https://www.ncbi.nlm.nih.gov/pubmed/30409984, https://www.ncbi.nlm.nih.gov/pubmed/30536387, https://www.ncbi.nlm.nih.gov/pubmed/30798813, https://www.ncbi.nlm.nih.gov/pubmed/31220083, https://www.ncbi.nlm.nih.gov/pubmed/31226389, https://www.ncbi.nlm.nih.gov/pubmed/31469255, https://www.ncbi.nlm.nih.gov/pubmed/31970519
             Polymorphisms: c.118=, c.187=, c.118C>G, c.187C>G, g.8671=, g.8671C>G, g.26090951=, g.26091179=, g.26090951C>G, g.26091179C>G, p.His40=, p.His63=, p.His40Asp, p.His63Asp, n.84=, n.281=, n.309=, n.84G>C, n.281C>G, n.309C>G
     Clinical Significance: RCV000000026.18, RCV000000027.5, RCV000175607.6, RCV000394716.6, RCV000763144.1, RCV000844708.1, RCV000991133.1
                      HGVS:

This question was also asked on Reddit

"What I want to do is to map bases to SNPs (or a nucleotide identifier like NC_000006.11)." - can you please ask this as a separate question; it's inconsistent with the title of this question.

These lines (with no REF or POS) are unmapped reads. If those lines are the entirety of the BAM file, then it is simply another file format to store the raw sequenced reads (as an alternative to FASTQ files). It's somewhat common for BAM files to include all reads, with mapped reads sorted by chromosome and position at the start, and unmapped reads appearing at the end of the file. This means that only the BAM file needs to be transferred and stored, rather than a BAM file and a separate 'unmapped' file (which can change size/content depending on the reference genome and mapper used).
FWIW, the lines of a BAM file can be viewed using samtools view; there's no need to convert to SAM first:
samtools view -h input.bam | less

If you only want to look at the unmapped reads, then you can filter for inclusion using the read flag 4:
samtools view -f 0x04 -h input.bam | less

If you want to exclude unmapped reads, you can filter that flag for exclusion (capital F instead of lower-case f):
samtools view -F 0x04 -h input.bam | less

Its a simple and very clear explanation.

This is an unmapped .bam file.  Someone just converted the fastq to bam format, but no one has aligned these to anything.  So you need to find an aligner program and a genome and align them.


text:
Is there an R package to calculate global FST, with confidence intervals, for tetraploid populations
I'm interested in estimating global FST among several populations of a tetraploid plant species. I'm looking for an R package that can do this, and can preferably also produce an upper and lower bound for a 95% confidence interval for FST. Does anyone know of a package that can be used for this purpose? I'll post a brief overview of my understanding of some of the options, in case I am not understanding something 100% correctly:
The StAMPP package seems really useful for this sort of thing with tetraploid or mixed populations, but unfortunatley it calculates pairwise FST and not a global measure. I'm aware that adegenet can calculate global FST for tetraploid populations, but it doesn't produce the confidence interval, and hierfstat can only be used with haploid or diploid populations. I've also came across the assigner package https://thierrygosselin.github.io/assigner/index.html, but am not sure if this is valid for tetraploid populations.

Its an interesting question. I'm skeptical about the value of population specific Fst (global), because it incorporates strong temporal perspective, rather than representing the extant population. Thus its quite different from pairwise Fst. Upvoted.

Personally, I would use adegenet and incorporate resampling with replacement or a simulated replication using your extant data set. This would introduce a confidence interval. It would need some insight into how tetraploid plants interbreed however.


text:
What does it mean by "topology" in the case of secondary and tertiary structures of proteins?
What does it mean by "topology" in the case of secondary and tertiary structures of proteins?
N.B. I am not talking about DNA/RNA/genomics. I am talking about protein folding, as it relates to protein modeling and simulation.

The word topology gets used in two different contexts.
Common case: the geometric configuration of the protein (=fold)
In common non-computational biochemistry parlance it is just the standard English meaning of topology and just geometry, in the case of protein it means protein fold:
https://en.wikipedia.org/wiki/Protein_fold_class
Wiktionary (https://en.wiktionary.org/wiki/topology) say:

The branch of mathematics dealing with those properties of a geometrical object (of arbitrary dimensionality) that are unchanged by continuous deformations (such as stretching, bending, etc., without tearing or gluing)

In the of protein, unchanged by small changes in length or amino acid sequence identity —no proteases or meat-glue ligase required!
Computational (Bio)chemistry
Topology has a technical jargon definition in computational biochemistry and computational chemistry.
In computational chemistry how the connections of a molecule are configured, expressed as a graph network https://en.wikipedia.org/wiki/Graph_theory, is referred to as topology. For example isopropanol and n-propanol have same number of heavy atom (nodes) and these will have the same hybridisation, but the network will have different connectivities (vertices).
The topology definition generally runs off atom types more so than elements, these combine the element and hybridisation and other properties specific to a certain element in a molecule —e.g. amine vs. amide nitrogen vs. tryptophan ring nitrogens are different.
The 3D embedding of said molecule is not part of its topology. Things get a bit problematic as ideal bond lengths, angles and dihedrals in dihedral space (i.e. not in cartesian space) go into the 'configuration' file (topology file) for several forcefield systems. And in others still there may be constraints or restraints specified.
In computational biochemistry, this holds true for ligands. But with proteins there is less. The connectivity within each amino acids is preset in stone: where one to use the allo isomer of threonine as opposed to the threo one would need to define a non-canonical amino acid. So by topology of a protein it is intended the connectivity of the amino acids in chains, e.g. primary sequence devoid of cartesian coordinates.
The word topology appears often in the context of ligands: these require custom definitions (=topologies) required to use a ligand not parameterised in the force field used, a common source of trouble.
See also

https://manual.gromacs.org/documentation/2019/reference-manual/topologies.html
https://mdtraj.org/1.9.4/api/generated/mdtraj.Topology.html

I am actually talking about protein modeling and simulation. :)

I have amended my answer, but I should warn that by virtue of being used just to mean configuration/definition in compbiochem, there are a hundred and one cases where my answer is 100% wrong due to the layers of complexity.


text:
Questions about CreateGeneActivityMatrix function in Seurat
Seurat V3 has a function CreateGeneActivityMatrix, which can convert ATAC peaks to gene-level's matrix.
Here is the prototype:
CreateGeneActivityMatrix(
  peak.matrix,
  annotation.file,
  seq.levels = c(1:22, "X", "Y"),
  include.body = TRUE,
  upstream = 2000,
  downstream = 0,
  verbose = TRUE
)

Here are my questions:

I want to know what is the format of the first parameter peak.matrix ? For example, what are the row names of this matrix? (Is there a strict rule of the row names?)
The parameter seq.levels. Besides 1...22,X,Y , Can I add other chromosome ? like hs37d5


text:
How do I intrepret the generated motifs of a dataset of 5-methylcytosine (5mC) sequences?
I'm a computer science person currently working on a project which uses deep learning to identify 5mC sites of Nipponbare Rice. Within the project, there is an option to view what the model 'discovers' through the generated motifs. Here are the motifs generated by most of the models

with the most common motif of being ACGGC. I've also read this article which states that

Alignment of the 5mC methylation sites indicated the conserved dinucleotide CG in the center with less conserved flanking bases.

and

The 5mC motif was similar to the methylation motif 5′-CmCGG-3′ of NgoAXIV in Neisseria gonorrhoeae.

Am I correct in assuming that the model has 'discovered' something that has been in the literature, specifically the CG motif? Also, how do I make sense of the second statement in terms of the connection of the identified motif in rice and in a bacteria?


text:
Perform protein structure-based sequence alignment in Python
I am looking for a Python package that performs pairwise structural alignment of protein structures (i.e., PDB files) and returns a sequence alignment. PyMOL is able to do this through the GUI, for example:

For two protein PDBs, one can be aligned to the other with the "align" menu or align command:

With it, the sequences appear aligned in the sequence viewer according to the structural alignment, and an "alignment" object is created, which can be saved with the command:
save A_aligned_with_B.aln, aln_structure_1_to_structure_0
which produces the A_aligned_with_B.aln alignment file in Clustal format (nevermind the ?s, which correspond to special HIS resnames):
CLUSTAL

structure_0  YIYQDSISLPWKVLLVMLLALITLATTLSNAFVIATVYRTRKL?TPANYLIASLAVTDLLVSI
structure_1  EEQGNKL?WAALLILMVIIPTIG-----GNTLVILAVSLEKKLQYATNYFLMSLAVADLLVGL
                                         .*..**..*....**....**...****.****..

structure_0  LVMPISTMYTVT-GRWTLGQVVCDFWLSSDITCCTASIL?LCVIALDRYWAITDAVEYSAKRT
structure_1  FVMPIALLTIMFEAMWPLPLVLCPAWLFLDVLFSTASIM?LCAISVDRYIAIKKPIQANQYNS
             .****......  ..*.*..*.*..**..*....****.***.*..***.**........ ..

structure_0  PKRAAVMIALVWVFSISISLPP-FFWRQAKAEEEVSE-CVVNTD?IL-YTVYSTVGAFYFPTL
structure_1  RATAFIKITVVWLISIGIAIPVPIKGIETDVDNPNNITCVLTKERFGDFMLFGSLAAFFTPLA
             ...*...*..**..**.*..*  .              **....    ........**..*..

structure_0  LLIALYGRIYVEARSRILKQTPNRTGKLLEKKKLMAAR--ERKATKTLGIILGAFIVCWLPFF
structure_1  IMIVTYFLTI?ALQKKAYLVKNKPPQRLTIGKKSVQTISNEQRASKVLGIVFFLFLLMWCPFF
             ..*..*....                              * .*.*.***....*...*.***

structure_0  IISLVMPICKDACW--F?LAIFDFFTWLGYLNSLINPIIYTMSNEDFKQAF?KLIRFK--
structure_1  ITNITLVLCDS-CNQTTLQMLLEIFVWIGYVSSGVNPLVYTLFNKTFRDAFGRYITCNYR
             *.......*   *   .  .....*.*.**..*..**..**..*..*..**

I need a Python package that is easily installable/distributable that accepts two PDBs or protein structures and returns the structural alignment of the sequences somehow. All within a Python script, and not creating intermediate files, e.g. for the alignment.
There are plenty of Python tools that allow structural superimposition and RMSD calculation (e.g., Biopython), but I haven't found any that are able to return a sequence alignment from the structural alignment.
Edit:
Possible answer, thanks to the answer and comments from Matteo Ferla.
import pymol2
with pymol2.PyMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar', object='aln')
    aln = pymol2.cmd2.pymol.exporting.get_alnstr('aln', _self=pymol.cmd) # I don't know if this line could be cleaner
    # aln is itself a string

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Please provide more details to help us in answering your question. For example, showing the alignment that you would like to see from Chimera or PyMOL. Showing the first few lines of any input files is also helpful (or, at least, providing two demonstrative PDB IDs that you are interested in).

This was asked before here —I know because I answered it. You want the primary sequence alignment based on PyMOL's superimposition of tertiary structure.

Ignore my comment. I saw your Q was just updated...

"align" performs a sequence alignment followed by a structural
superposition, and then carries out zero or more cycles of
refinement in order to reject structural outliers found during
the fit. "align" does a good job on proteins with decent sequence
similarity (identity >30%). For comparing proteins with lower
sequence identity, the "super" and "cealign" commands perform
better.

cealign aligns two proteins using the CE algorithm. It is very robust for proteins with little to no sequence similarity (twilight zone). For proteins with decent structural similarity, the super command is preferred and with decent sequence similarity, the align command is preferred, because these commands are much faster than cealign.

yoiu can save a clustal object as a result of cealign, cealign target, mobile [, target_state [, mobile_state
[, quiet [, guide [, d0 [, d1 [, window [, gap_max
[, transform [, object ]]]]]]]]]]   --> object = string: name of alignment object to create {default: (no alignment object)}  https://pymolwiki.org/index.php/Cealign

PyMOL in itself can be used as python package.
The GUI app has its own Python, but in a conda environment you can install it as a package.
Any PyMOL app command can be used as a cmd command. In PyMOL Wiki you see both. They follow a similar pattern and actually have more functionality.
import pymol2
with pymol2.PYMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar')
    pymol.cmd.save('foobar.aln')

For more discussion of how the alignment output object is used see Assessing PyMol sequence alignment object

Hi! Thanks so much for your answer. Sorry it took me so long to say anything, I was testing all the options. Also, I'm sorry I don't have enough reputation to upvote your answer. Is there any way I could save the alignment in a Python text stream, as in an io.StringIO() stream? So that I don't save any file in the process. I made some attempts but I couldn't hack it and I didn't see anything related in your previous answer that you linked. Thanks in advance!

Nah, the important thing is that marked as solved, thus helping the healthiness of the site, which struggled to leave beta because of high unsolved cases —so thanks!
I don't know off the top of my head, but https://docs.python.org/3/library/tempfile.html may be of use. There is likely a function that is not imported by the top of cmd module, but is present in one of the submodules. Maybe one of the searches in the open-source PyMOL repo might help https://github.com/schrodinger/pymol-open-source/search?q=.aln (sorry I'm tight for time atm)

Thanks for the pointers! I finally found out how to retrieve the alignment string, I'm going to edit my question to put the final answer in it.

@MatteoFerla any way to use it in python without anaconda enviroment ? https://stackoverflow.com/questions/63958085/how-to-embed-the-pymol-module-into-python

Hi @pippo1980 maybe a separate question? Generally not advised, you could use venv and pip, but historically there were compatibility issues between conda and pip (probably resolved)

Anyway thanks @FranchoNerínFonz for finalising the post. Very interesting (and useful) question and answer indeed.

@pippo1980 The classic solution was installing the OS pymol —which is not installed via conda but with wizardry and  the deprecated python setup.py install (not pip install .) https://github.com/schrodinger/pymol-open-source/blob/master/INSTALL —alternatively, there's always the nightmarish über-hack of copypasting of site-packages folders from a conda installation to a vanilla Python... But that's a big no no due to irreproducibility.


text:
How deep variant and others preprocessed insertion and deletion for extracting candidate variant?
I am trying to build a model for variant calling,but I am a bit confused how to pre-process the reads that has insertion and deletion and align them with reference?
Should I add something to the location of insertion or deletion to reference or should I just delete from read?

Can you provide some reference for why you think preprocessing is necessary?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Are you building a model specifically for deep variant, or trying to create your own model via a different algorithm?


text:
How to explain disparities in Flavonoid Docking results between Vina and Schrodinger?
I am planning to dock a protein to Selleck Chem's Flavonoid Library. I have designed two workflows using Autodock Vina and Schrodinger as follows:
Schrodinger

Use their proprietary tools for protein preparation followed by grid preparation under default settings.
Run LigPrep on the SelleckChem library to generate one ligand (at most) per compound under default settings.

I executed another run (let's call it 'special run'; see below) generating eight ligands (at most) per compound under default settings.

Run Glide in SP mode under default settings + add a reward for H-bonds.

Autodock Vina

Run every module of UCSF Chimera's Dock Prep on the Protein under default settings.

Except use Reduce in place of AddH - suggested by Chimera as a better alternative.

Calculate grid-center to be the centroid of the default ligand and use the same grid size as used in Schrodinger workflow.
Run prepare_receptor4.py but in non-default settings to only remove lone pairs and non-standard residues.

We do not remove waters or add charges, as done under default settings; they have been already executed by Chimera.

Download 3D files of each compound in the Selleckchem library from Pubchem in sdf format.

Use Obabel to convert them to pdb format.
Run prepare_ligand4.py under default settings to get output in pdbqt format.
Use Obabel to replace Gasteigger charges with eem2015ba.
Use Obabel to minimize energy using MMFF94 forcefield under default settings.

Run Vina 1.2 to dock at exhaustiveness of 32.

The results from the two workflows exhibit vast disparity with the correlation coefficient of the rank being about 0.2. For example, the best result of Schrodinger is ranked 150 by Vina (the library has about 200 molecules in total) and so on. If I take the best 20 results from Vina and Schrodinger each, there is a single common compound.
In my 'special run', I saw that Schrodinger assigns vastly different overall docking scores to various conformations of the same ligand. This might answer the above disparity unless I am doing something acutely wrong. But that raises the question about how to choose the initial conformation of the ligand if I am constrained to use open-source software and starting conformations influence the final result so severely?

Please show the first few lines of results, or whatever other appropriate way can demonstrate this desparity.

@gringer I have added a couple of lines. Does this suffice?

Yes, thanks. I think that's enough.


text:
syndip dataset for benchmark variant
I have a question about syndip dataset : https://github.com/lh3/CHM-eval .
I'm struggling to find the syndip vcf.
In the release ( https://github.com/lh3/CHM-eval/releases ), we have a file named : rep2.37.broad.hc.raw.vcf.gz, that i don't know what it is. And we have a file named CHM-evalkit-20180222.tar which contains full.37m.vcf and other files (bed, eval ...).
So i did my search and according to this file: they mentioned that full.37m.vcf is the truth dataset. ( https://www.biorxiv.org/content/biorxiv/suppl/2018/10/29/456103.DC1/456103-1.pdf Page 16).
The problem is that the file rep2.37.broad.hc.raw.vcf.gz contain variants with MQ, DP, GQ ... that I need to extract. But the full.37m.vcf doesn't contain this information.. ( just Chrom pos ref alt and QUAL.)
So I tried to intersect rep2.37.broad.hc.raw.vcf.gz with full.37m.vcf and take the variant that present in two files, with the DP MQ GQ in rep2.37.broad.hc.raw.vcf.gz.
Is that okay? Since I don't know what is rep2.37.broad.hc.raw.vcf.gz.
And I also noticed that the QUAL in the full.37m.vcf is always 30. Is it normal?

If you want the truth VCF, it is full.37m.vcf. rep2.37.broad.hc.raw.vcf.gz is an example call set by GATK.

thank you for the response  !


text:
recover files after groupTuple
I have a process that emits BAM files with a key containing the sample name:
process source{
  output:
    tuple val(sample), path("sample.bam") into ch

  [...]

}

1st question: I have several files per sample, let's assume 2 for now. I want to merge them together with samtools merge. I can regroup that channel with groupTuple(), but how can I recover them in the next process? Considering they all have the same name:
grouped_ch = ch.groupTuple(by: 0, size: 2)

process merge{
  input:
    grouped_ch

  shell:
  '''
  samtool merge sample.bam sample.bam > !{sample}.bam
  '''
}

Is there perhaps a way in the input: block to say that I expect a list, and name the list components?
2nd question: In practice, I have a variable number of "sample.bam" files for each sample value, does that change the way to define the input:? It also appears I will need to use the built-in groupKey(), is there documentation available? I have read this Github issue, but it would help if a general description was available.

My preference is to avoid working with files with the same name as much as possible. However, sometimes this might be unavoidable. Fortunately, file name collisions can be avoided when working with multiple input files by using the * and ? wildcards:

When a target file name is defined in the input parameter and a collection of files is received by the process, the file name will be appended by a numerical suffix representing its ordinal position in the list.
The target input file name can contain the * and ? wildcards, that can be used to control the name of staged files.

Note that there is also the stageAs path option (see input of type 'path'), which can also take a name pattern (using the above wildcards) to let you use a variable in your script. For example:
process samtools_merge {

    conda 'samtools=1.15.1'
    debug true
    stageInMode 'rellink'

    input:
    tuple val(sample), path(bam_files, stageAs: 'mysample*.bam')

    output:
    tuple val(sample), path("${sample}.bam")

    """
    samtools merge "${sample}.bam" ${bam_files}
    ls -nl
    """
}

workflow {

    Channel.fromPath( ['./path/to/files/*.bam', './path/to/more_files/*.bam'] ) \
        | map { bam -> tuple( bam.baseName, bam ) } \
        | groupTuple(by: 0, size: 2) \
        | samtools_merge
}

Results:
N E X T F L O W  ~  version 22.04.0
Launching `script.nf` [elated_montalcini] DSL2 - revision: 1427d64a1d
executor >  local (3)
[83/25654c] process > samtools_merge (3) [100%] 3 of 3 ✔
total 184
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample1.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample1.bam
-rw-r--r-- 1 1000 985 187911 May 25 12:18 sample1.bam

total 128
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample2.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample2.bam
-rw-r--r-- 1 1000 985 129947 May 25 12:18 sample2.bam

total 196
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample3.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample3.bam
-rw-r--r-- 1 1000 985 197225 May 25 12:18 sample3.bam

If each group contains a variable number of files, you can keep the same input declaration (as above) but you will want use a special groupKey object so that the collected values can be streamed as soon as possible. It's not well documented, but here's a simple example that sorts the input BAM files prior to merging. Once all of the 'samtools_sort' jobs are finished for a particular sample, that sample can proceed to merging without us having to wait for all the other 'samtools_sort' jobs to be finished (which would otherwise be the case if the 'size' parameter had not been specified):
process samtools_sort {

    tag { bam.baseName }

    conda 'samtools=1.15.1'

    input:
    tuple val(sample), path(bam)

    output:
    tuple val(sample), path("${bam.baseName}.sorted.bam")

    """
    samtools sort -o "${bam.baseName}.sorted.bam" "${bam}"
    """
}

process samtools_merge {

    tag { sample }

    conda 'samtools=1.15.1'

    input:
    tuple val(sample), path(bam_files, stageAs: 'mysample*.bam')

    output:
    tuple val(sample), path("${sample}.bam")

    """
    samtools merge "${sample}.bam" ${bam_files}
    """
}

workflow {

    Channel.fromPath( params.bam_files ) \
        | map { bam -> tuple( bam.baseName, bam ) } \
        | groupTuple() \
        | map { sample, files -> tuple( groupKey(sample, files.size()), files ) } \
        | transpose() \
        | set { bam_files }

    samtools_sort( bam_files ) \
        | groupTuple() \
        | samtools_merge \
        | view()
}

Results:
$ nextflow run script.nf -ansi-log false --bam_files './path/to/*/*.bam'
N E X T F L O W  ~  version 22.04.0
Launching `script.nf` [cheesy_swanson] DSL2 - revision: da1a82a059
[04/643129] Submitted process > samtools_sort (sample3)
[1a/93a939] Submitted process > samtools_sort (sample1)
[93/fde0a0] Submitted process > samtools_sort (sample3)
[db/208bac] Submitted process > samtools_sort (sample3)
[8f/f2a231] Submitted process > samtools_sort (sample2)
[f3/8f0f9f] Submitted process > samtools_sort (sample2)
[bb/cca23e] Submitted process > samtools_merge (sample1)
[sample1, /home/steve/testing/work/bb/cca23e016a641b201b7e42f9db74ef/sample1.bam]
[ff/bd3a2a] Submitted process > samtools_merge (sample2)
[sample2, /home/steve/testing/work/ff/bd3a2ac83a0691d2b80b04ef97b755/sample2.bam]
[94/a4f8eb] Submitted process > samtools_merge (sample3)
[sample3, /home/steve/testing/work/94/a4f8eb483602e9335b28f5d40058c8/sample3.bam]

Note that the val specified as input to 'samtools_sort' is an object of class nextflow.extension.GroupKey. In the example above, the input to 'samtools_merge' is also a groupKey - but it doesn't need to be. If you need a string value for whatever reason, you can just call '.toString()' to get it back:
    ...

    samtools_sort( bam_files ) \
        | groupTuple() \
        | map { sample_key, files -> tuple( sample_key.toString(), files ) } \
        | samtools_merge \
        | view()

It's awesome how much effort you put into these nf answers!

Thank you so much for this very comprehensive answer! I had totally missed that section of the input documentation. And as ATpoint notes, this is the second time you help me on a Nextflow question, this has been extremely helpful!


text:
Generate consensus protein sequence from relatively gappy alignment?
I would like to know what would be the best way to generate a full-length consensus sequence?
I want to obtain a single representative sequence from a relatively gappy multiple sequence alignment of peptides.
See example below.
So far I have attempted Consensus.py from https://github.com/jrjhealey/bioinfo-tools but it just gives me gaps, even with a threshold close to zero.
>HX5_33766:7:1222:1022:72139
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------SEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFKDKVYFNVT--------------------------------------------
>HX5_33766:7:1211:6319:29841
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIIINCDGY------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1222:1022:72139.2
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------SEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFKDKVYFNVT--------------------------------------------
>HX5_33766:7:1219:12723:20928.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-ISTSYKDIQVNGGELILIHQNPGGFCAY--------------
>HX5_33766:7:1119:11069:37612
----------------------------------------------------------------------------------------------------------------------------------------VLLGAAIIQGLGLLLC--F----------------------TYI----------------------CLHFYASQ-V----------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1219:12723:20928
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-ISTSYKDIQVNGGELILIHQNPGGFCAY--------------
>HX5_33766:7:1211:6319:29841.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIIINCDGY------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1219:12723:20928.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-ISTSYKDIQVNGGELILIHQNPGGFCAY--------------
>HX5_33766:7:1107:2564:39757
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------VFSLKGYFSEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFK----------------------------------------------------
>HX5_33766:7:1211:6319:29841.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIIINCDGY------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1212:30371:5686
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIXLGCN--------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1119:11069:37612.2
----------------------------------------------------------------------------------------------------------------------------------------VLLGAAIIQGLGLLLC--F----------------------TYI----------------------CLHFYASQ-V----------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1212:30371:5686.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIXLGCN--------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:10013:25165
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMK-------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:10013:25165.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMK-------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1106:31832:24286
----------------------------------------------------------------------SVYLRL--IVFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1117:3691:48881
-------------------------------------------------------------------------------VFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:9364:63314
-------------------------------GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:10013:25165.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMK-------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:9364:63314.2
-------------------------------GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1212:30371:5686.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIXLGCN--------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:9364:63314.1
-------------------------------GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1119:11069:37612.1
----------------------------------------------------------------------------------------------------------------------------------------VLLGAAIIQGLGLLLC--F----------------------TYI----------------------CLHFYASQ-V----------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1107:2564:39757.1
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------VFSLKGYFSEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFK----------------------------------------------------
>HX5_33766:7:1106:31832:24286.1
----------------------------------------------------------------------SVYLRL--IVFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1117:3691:48881.1
-------------------------------------------------------------------------------VFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1222:1022:72139.1
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------SEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFKDKVYFNVT--------------------------------------------
>HX5_33766:7:1107:2564:39757.2
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------VFSLKGYFSEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFK----------------------------------------------------
>HX5_33766:7:1106:31832:24286.2
----------------------------------------------------------------------SVYLRL--IVFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1117:3691:48881.2
-------------------------------------------------------------------------------VFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1121:4320:48213
MTVISQSLKSELECHKNNSLLAQ-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1121:4320:48213.1
MTVISQSLKSELECHKNNSLLAQ-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1121:4320:48213.2
MTVISQSLKSELECHKNNSLLAQ-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1214:16701:50287
----------------------------------------------------------------------------------------------------------------------------------------------R--KQR--EXPKIEQLSAFTRSRSPESAFSNNHSVTYFHTLWGD------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1220:24545:3085
------------------------------------------------------------------------------------------------------------------------------------------------------------------RSRSPESAFSNNHSVTYFHTLWGD-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------R

Given the alignment you've presented, there are obvious groups of sequence similarity, but I'm seeing a lot of different peptide sequences in this (e.g. CENEKGCIITSPSKDETMKVQDNSIXLGCN and GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR as two obviously different peptides). Have you done a pre-filter for similarity to make sure you're looking at the same proteins for alignment?
Multiple alignment will not work properly with peptides that are very different from each other. You need to group proteins based on similarity first (e.g. with Linclust), then carry out multiple alignments within the groups.

from Bio import AlignIO
from pathlib import Path
from Bio.Align import AlignInfo

path = '/my/path/'
test1 = 'testset.fa'
mypath = Path(path,test1)

myalign = AlignIO.read(mypath, "fasta")
summary = AlignInfo.SummaryInfo(myalign)
consensus = summary.dumb_consensus()
print (consensus)

This will give a single consensus and 'X' low prevalence ambiguity. Easy to check whether its what is required because the code is ready to go.


text:
How can I improve or otherwise investigate an unreliable genome tree?
Summary My genome tree doesn't agree with my gene trees and I get the feeling that my genome tree might be wrong, possibly due to long branch attraction, but I don't know how to check/fix it.
Background
I have a set of genes of interest recovered from genomes and metagenomes. I also have some bins built from the metagenomes of interest.
Previous work I've built multiple gene trees, which tend to agree with each other, as well as with other phylogenetic methods (16S rRNA, initial genome trees where available, conserved orientation of the genes of interest). So at this point I though I had a pretty good idea of how my genes / species of interest are linked to the species tree.
Concern
I've assembled the metagenomes and recovered bins. I've refined bins of interest with anvio, based on sequence composition / differential coverage, and kept only the bins which appeared quite high quality (>70% completion*, <10 - often <5% - redundancy). Thus contamination in the metagenome assemblies (MAGs) is unlikely. The contamination would have to happen multiple times in the same way to different samples from different environments, at the same coverage, to get the same effect.

with the exception of a bin of particular interest which I accepted with 62%

I have de-replicated these bins, added some known genomes, and tried to build a species tree with Orthofinder (with the proteomes). There is a problem occurs with a clade of 3 bins I will call 'X'.

Trees Based on the gene-of-interest phylogenies, X looks like probably a sister-group to the rest of known genomes from order 'A'. Which is very interesting.
With the 16S rRNA, it's not quite as clear (there are no 16S rRNA sequences in the bins); however, some 16S rRNA amplified from one of the original samples fit the bill as 'sister-clade to basically everything else in group A'). The first (100 or so) BLAST hits for this 16S sequence are all from group 'A', still.
Tree conflict However, in the proteome tree, X as well as Y (Y being a sequenced genome that I suspected was related to X based on 16S data) end up in order 'B' instead of 'A'. GTDBtk also labels X as belonging to clade 'B', not 'A'.
Now,I know that genes and genomes don't always have the same evolutionary trajectories. But multiple sets of data lead me to a different conclusion than the proteome tree, so I have to at least question whether the proteome tree is correct; what's likely to be the problem (I suspect long-branch attraction), and, if it is long-branch attraction (or anything else) how can I uncover the correct signal?
Thank you for your time!

Thanks @Laura ... I ain't sure LBR is occurring but would need to think. What is the tree building algorithm and parameters?

@M__ Thank you for having a look! I used Orthofinder (https://github.com/davidemms/OrthoFinder) which appears to use an algorithm called STAG (https://www.biorxiv.org/content/10.1101/267914v1) for the species tree building. I'm not sure on the details. The taxonomy tool GTDB-Tk (which appears to use HMMER /pplacer) is consistent with this proteome tree (but not with everything else).

Where is taxa 'Y' in the far left tree? I can see 'X' ... but no 'Y'

The difficulty in answering the question is:

The absence of "Y" in the "Gene of interest" tree
The sampling bias in the number of taxa mostly between the "genome tree" and the "gene of interest" and "16S" tree
The "genome" tree is unrooted so its difficult to compare it with the "16S" tree, thus its possible that if "16S" was unrooted it falls within monophyly "B"
The trees don't have any branch or statistical support
We don't really know the tree building method (I haven't looked it up)
The bacteria might have significant recombination (we don't know the bacteria).

The answer is there is a lot of room for different tree artefacts here, however it could be a biologically real result. I agree "the tree" is simply a tool which uses a series of assumptions and obtaining confidence in a tree requires a lot of analysis. Vast amounts of computational time are dedicated to obtaining a robust tree.
The principle step is to compare like for like:

same number of taxa and
same rooting position (or unrooted).
second step is to obtain statistical robustness for the observation.

There is a formal process for identifying conflicting trees and it needs code - I have written stuff that does this and maybe I should release it. Essentially, you would data chunk gene by gene to produce a genome of gene trees and use e.g. ETE3 or similar to plot the tree conflicts. It ain't rocket science. Thus the amount of user time is trivial. The idea is that once the "problem trees" are identified, e.g. rapid evolution.
The one observation I would say is that the "Y" branch looks a little bit weird and this could indicate part of the data is has resulted from different evolutionary rate and that is the basis for long-branch attraction.

Scenario 1 It could be that the 16S and "genome" tree are part of the same tree (monophyly B), they are just rooted differently. Thus the conflict is just a small number of genes and this could lead to ...
Scenario 2 Keep in mind bacteria do undergo 'horizontal transmission' and if 'gene of interest' is under strong selective pressure recombination is entirely possible.

Hello. My mistake - I get these results with the rooted trees. HGT would be possible, but hasn't been shown for that particular set of genes before, so I would have to provide fairly strong evidence. I also think it would be highly unusual (assuming that organism X really belonged to B) that when blasting the 16S of an organism, the results I'd get would mostly be from another clade. Thank you very much for clarifying this and telling me about ETE3. The program I used for the genome tree has the gene trees stored, so I can compare them : )

Great @Laura glad its useful and its worked out. ETE3 toolkit BTW, bit different from ETE3 python module is stuck in Python 3.6. The all-singing all-dancing ETE4 was supposed to be released already - it got full funding, but maybe never.

It seems like you've got the right approach, in trying out a few different genes and looking for consistency. As you say, different genes can have different selective pressures, which influences how much they change in relation to the whole genome background.
Certainly don't put too much trust in a single gene, and especially not 16S: it has a low mutation rate, can jump around and be repeated in the genome, and is an occasional target for horizontal gene transfer.
The answer may be that there isn't any single correct phylogeny. A primary assumption of phylogenetics is that evolutionary progression has a tree-like structure (with roots and tips), but that model is broken by horizontal gene transfer and plasmid sharing. It may be the case that you need to break the genome up into fragments, each with its own inherited history.


text:
Is this RNA seq data good based on the MDS plot
I am analyzing some RNA sequencing data from a collaborator where the effect of some ligands on the transcriptome is being looked at. I am using DESeq2 for my analysis. I am looking to compare each of the ligand treatments vs the DMSO. I have attached the MDS plot with this post. There seems to be a lot of variation between the 3 replicates across all the treatments. Based on this plot, is it worth it to use this data to run a differential expression?
It is also at this link below.
https://i.stack.imgur.com/ZXhk5.jpg

You definitely see separation by treatment in PC1 and PC2, so that is an indication that the treatments are doing something. DMSO and TCDD cluster close, so these are probably similar. Still, you see a variation along PC1 and it seems that there is within-group separation in all groups that could indicate a batch effect. Check whether the samples per group that are more leftish have been processed as one batch, if so one could include that information into the design, see DESeq2 vignette on batch effects. Alternatively, check tools such as RUVseq or SVA to detect unknown batch effects to correct for this within-grup separation.


text:
To create a new column with Pandas on the basis of two types of column values
I have a  data frame df. The below table only represent top 8 rows of the data frame. There are multiple combination of LIBRARY_NAME and ANCHOR_NAME
BARCODE ANCHOR_NAME     LIBRARY_NAME    ANCHOR_CONC Synergy HSA     Bliss
14482   5-Fluorouracil  Cisplatin       2.5         FALSE   FALSE   FALSE
                                        10          FALSE   FALSE   FALSE
                        Docetaxel       2.5         FALSE   TRUE    TRUE
                                        10          FALSE   FALSE   FALSE
        AZD7762         Cisplatin       0.0625      TRUE    TRUE    TRUE
                                        0.25        FALSE   TRUE    TRUE
                        Docetaxel       0.0625      FALSE   FALSE   FALSE
                                        0.25        FALSE   TRUE    FALSE

I want to create a new column 'new_synergy 'with pandas which is based on column ANCHOR_CONC and Synergy. There are two ANCHOR_CONC for each anchor_name and library_name combination. If for any concentration of ANCHOR_CONC the synergy is TRUE, for other concentration also synergy becomes TRUE
Eg: new_synergy: Condition for TRUE (for AZD7762 and Cisplatin), IF at 0.0625 Synergy is TRUE, it will be TRUE for 0.25 also.
Expected data_frame
new_df
BARCODE ANCHOR_NAME     LIBRARY_NAME    ANCHOR_CONC Synergy HSA     Bliss  new_synergy
14482   5-Fluorouracil  Cisplatin       2.5         FALSE   FALSE   FALSE    FALSE
                                        10          FALSE   FALSE   FALSE    FALSE
                        Docetaxel       2.5         FALSE   TRUE    TRUE     FALSE
                                        10          FALSE   FALSE   FALSE    FALSE
        AZD7762         Cisplatin       0.0625      TRUE    TRUE    TRUE     TRUE
                                        0.25        FALSE   TRUE    TRUE     TRUE
                        Docetaxel       0.0625      FALSE   FALSE   FALSE    FALSE
                                        0.25        FALSE   TRUE    FALSE    FALSE

I tried the following code
import pandas as pd
import numpy as np
dd = data[["BARCODE", "LIBRARY_NAME", "ANCHOR_NAME", "ANCHOR_CONC", "Synergy", "Bliss", "HSA"]]
df = pd.pivot_table(data=dd,index=["BARCODE", "ANCHOR_NAME", "LIBRARY_NAME", "ANCHOR_CONC"])
## setting conditions as pivot table step converted true to 1 and false to 0.
conditions = [(dd['Synergy'] == 1.0), (dd['Synergy'] == 0.0)]
values = ['True', 'False']
dd['new_synergy'] = np.select(conditions, values)
## However it is giving me the same results in both new_synergy and synergy column. I have no idea how to fit `ANCHOR_CONC` condition in this step.

Could you supply your code please?

You can use groupby() as following :
dd["new_synergy"] = dd[["ANCHOR_NAME", "LIBRARY_NAME", "Synergy"]].groupby(by=["ANCHOR_NAME", "LIBRARY_NAME"]).any()
Explenation : you want to group your dataframe by both ANCHOR_NAME and LIBRARY_NAME, and return True if any of the synergyvalue is true, which is what any() does.
Note : you can check for all methods available for GroupBy objects here.

To elaborate on this answer, here is a full working example where I assume that "BARCODE", "ANCHOR_NAME" and "LIBRARY_NAME" are used as MultiIndex:
Creating the test data:
from io import StringIO
import pandas as pd
data_str = """BARCODE\tANCHOR_NAME\tLIBRARY_NAME\tANCHOR_CONC\tSynergy\tHSA\tBliss
14482\t5-Fluorouracil\tCisplatin\t2.5\tFALSE\tFALSE\tFALSE
14482\t5-Fluorouracil\tCisplatin\t10\tFALSE\tFALSE\tFALSE
14482\t5-Fluorouracil\tDocetaxel\t2.5\tFALSE\tTRUE\tTRUE
14482\t5-Fluorouracil\tDocetaxel\t10\tFALSE\tFALSE\tFALSE
14482\tAZD7762\tCisplatin\t0.0625\tTRUE\tTRUE\tTRUE
14482\tAZD7762\tCisplatin\t0.25\tFALSE\tTRUE\tTRUE
14482\tAZD7762\tDocetaxel\t0.0625\tFALSE\tFALSE\tFALSE
14482\tAZD7762\tDocetaxel\t0.25\tFALSE\tTRUE\tFALSE
"""
# StringIO is to simulate an import from a file
data = pd.read_table(StringIO(data_str), sep="\t", index_col=[0, 1, 2])
# reset_index is here just for displaying purposes
# (MultiIndex does not seem well supported in (pandas?) markdown export)
print(data.reset_index().to_markdown(index=False))

BARCODE
ANCHOR_NAME
LIBRARY_NAME
ANCHOR_CONC
Synergy
HSA
Bliss

14482
5-Fluorouracil
Cisplatin
2.5
False
False
False

14482
5-Fluorouracil
Cisplatin
10
False
False
False

14482
5-Fluorouracil
Docetaxel
2.5
False
True
True

14482
5-Fluorouracil
Docetaxel
10
False
False
False

14482
AZD7762
Cisplatin
0.0625
True
True
True

14482
AZD7762
Cisplatin
0.25
False
True
True

14482
AZD7762
Docetaxel
0.0625
False
False
False

14482
AZD7762
Docetaxel
0.25
False
True
False

Creating a DataFrame with the new column:
new_synergy = data.reset_index().groupby(
    ["BARCODE", "ANCHOR_NAME", "LIBRARY_NAME"])[["Synergy"]].any()
# Rename the resulting "Synergy" column
new_synergy.columns = ["new_synergy"]
print(new_synergy.reset_index().to_markdown(index=False))

BARCODE
ANCHOR_NAME
LIBRARY_NAME
new_synergy

14482
5-Fluorouracil
Cisplatin
False

14482
5-Fluorouracil
Docetaxel
False

14482
AZD7762
Cisplatin
True

14482
AZD7762
Docetaxel
False

Merging the two DataFrames based on their indices:
new_data = data.merge(new_synergy, left_index=True, right_index=True)
print(new_data.reset_index().to_markdown(index=False))

BARCODE
ANCHOR_NAME
LIBRARY_NAME
ANCHOR_CONC
Synergy
HSA
Bliss
new_synergy

14482
5-Fluorouracil
Cisplatin
2.5
False
False
False
False

14482
5-Fluorouracil
Cisplatin
10
False
False
False
False

14482
5-Fluorouracil
Docetaxel
2.5
False
True
True
False

14482
5-Fluorouracil
Docetaxel
10
False
False
False
False

14482
AZD7762
Cisplatin
0.0625
True
True
True
True

14482
AZD7762
Cisplatin
0.25
False
True
True
True

14482
AZD7762
Docetaxel
0.0625
False
False
False
False

14482
AZD7762
Docetaxel
0.25
False
True
False
False


text:
Weird data in VCF file
I am using GATK HaploTypeCaller for SNP calling. I have 300 samples. As far as I understand that if a sample is missing a SNP it should be denoted by ./., however, I am getting a lot of values like 0/0:0,0:0:0:0,0,0.
As a result it is considered as homozygous and downstream processing detecting no missing data at all, like when I am using vcftools. This is weird as my sequencing coverage is only 1.0X, so I should have a lot of missing data per SNP. What is going wrong here? Is it okay if I convert all the 0/0:0,0:0:0:0,0,0 values to ./.:.,.:.:.:.,.,.
The commands I used:
gatk HaplotypeCaller
    -I P1_q10_sorted_rg_dedup.bam
    -R Zm-B73-REFERENCE-NAM-5.0.fa
    -ERC GVCF
    -O P1.g.vcf

gatk GenomicsDBImport \
    -V P1.g.vcf
    --genomicsdb-workspace-path /fs/project/PMIU0157/Mahmood_Hasan/Second/vcf_files/database/mop1/mop1_db_chr1
    --intervals chr1

gatk GenotypeGVCFs
    --max-alternate-alleles 300
    -R Zm-B73-REFERENCE-NAM-5.0.fa
    -V gendb://mop1_db_chr1
    -O mop1_variants_chr1.vcf

gatk SelectVariants
    -R Zm-B73-REFERENCE-NAM-5.0.fa
    -V mop1_variants_chr1.vcf
    --select-type-to-include SNP
    -O raw_mop1_SNPs_chr1.vcf

gatk VariantFiltration \
-R /fs/project/PMIU0157/Mahmood_Hasan/B73_ref_genomeV5/Zm-B73-REFERENCE-NAM-5.0.fa \
-V raw_wt_SNPs_chr6.vcf \
-O filtered_wt_SNPs_chr6.vcf \
--filter-name "QD2" \
--filter-expression "QD < 2.0" \
--filter-name "FS60" \
--filter-expression "FS > 60.0" \
--filter-name "MQ40" \
--filter-expression "MQ < 20.0" \
--filter-name "MQRankSum-12.5" \
--filter-expression "MQRankSum < -12.5" \
--filter-name "ReadPosRankSum-8" \
--filter-expression "ReadPosRankSum < -8.0"

#The output looks like the following
chr1    321     .       T       A       526.30  PASS    AC=10;AF=0.030;AN=336;DP=42;ExcessHet=0.0000;FS=0.000;InbreedingCoeff=0.3471;MLEAC=43;MLEAF=0.128;MQ=35.97;QD=25.36;SOR=0.892   GT:AD:DP:GQ:PGT:PID:PL:PS
       0/0:1,0:1:3:.:.:0,3,36  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,26  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,12  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,33  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,12  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:2,0:2:6:.:.:0,6,49  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,14  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:2,0:2:6:.:.:0,6,39  1|1:0,2:2:6:1|1:321_T_A:90,6,0:321      0/0:2,0:2:6:.:.:0,6,40  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,16  0/0:0,0:0:0:.:.:0,0,0   0/0:2,0:2:6:.:.:0,6,45  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,25  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,24  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,42  0/0:2,0:2:6:.:.:0,6,41  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,40  0/0:0,0:0:0:.:.:0,0,0
   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,15  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,11  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   1|1:0,2:2:6:1|1:321_T_A:90,6,0:321      0/0:0,0:0:0:.:.:0,0,0   0/0:2,0:2:6:.:.:0,6,49  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,22  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   1|1:0,2:2:6:1|1:321_T_A:90,6,0:321      0/0:1,0:1:3:.:.:0,3,20  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   1|1:0,2:2:6:1|1:321_T_A:90,6,0:321      0/0:0,0:0:0:.:.:0,0,0   1|1:0,1:1:3:1|1:321_T_A:45,3,0:321      0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,14  0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,26  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0
   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,39  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:2,0:2:6:.:.:0,6,49  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:1,0:1:3:.:.:0,3,16  0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0   0/0:0,0:0:0:.:.:0,0,0


text:
Checking that two multiline FASTAs are identical (allowing for different order and lower/upper case)
Is there a good CLI tool that can diff two multiline fasta files, ignoring order, and comparing sequences with identical names with each other pairwise.
Bonus points if lower/upper case is tolerated (as in, it's considered identical) and if also indels can be tolerated, maybe through a flag?
I could script this in Python myself, but I was wondering whether there was a custom tool for this purpose already out there. Maybe something more performant than a Biopython script.
If there's no custom tool, I'd also be happy with a short script that uses a number of tools, e.g. seqkit samtools diff etc.

Why are indels tolerated but not SNVs?

Try md5fa from samtools. It computes two MD5s with and without ordering, respectively. Alternatively, there is seqkit sum, which is perhaps more powerful.

There are a lot of ambiguities in the question, but I'll provide one approach by means of a fasta/fastq sorting program (see here), using process substitution to avoid intermediary files. This will only work for small-ish fasta files, which can fit into computer memory:
diff -q <(fastx-sort.pl a.fasta) <(fastx-sort.pl b.fasta)

Note that this will compare the entire file. If there are any sequences in one file but not the other, then the result will be that the files are different.
The -q option for diff will just output whether or not the two files differ; you didn't specify in your question whether you wanted any more detail. This option can be removed to show differences, or replaced with -u to add context. Context and/or line changes are not necessarily going to be helpful for multi-line fasta sequences, because the compared files may have a different order from the input, and the context lines may not include the name of the sequence. This could be fixed by changing the files to single-line fasta format (e.g. by piping through fasta_formatter from the FASTX-Toolkit):
diff -u <(fastx-sort.pl a.fasta | fasta_formatter) <(fastx-sort.pl b.fasta | fasta_formatter)

Note that this approach will output the entire sequence for any different sequences, together with adjacent context (including the sequence name). There are ways round that as well... but it would be far better if you could make your question more specific so that we don't need to explore the universe when answering.
If you want to ignore case, convert everything to upper (or lower) case:
diff -q <(fastx-sort.pl a.fasta | perl -pe 'tr/a-z/A-Z/') <(fastx-sort.pl b.fasta | perl -pe 'tr/a-z/A-Z/')

Excellent thanks! I was thinking of using seqkit, will try in a bit.
Does this output the diffing lines?