text:
Subsetting with harmony
Currently merged two Seurat objects together and then ran Harmony for batch correction. Now I want to subset out a cell type of interest do I re run harmony or just do the standard PCA?
scData <- scData %>%
    RunUMAP(reduction = "harmony", dims = 1:14) %>%
    FindNeighbors(reduction = "harmony", dims = 1:14) %>%
    FindClusters(resolution = 0.2) %>%
    identity()

numPC = 14
scData <- FindNeighbors(scData, dims = 1:numPC)
scData <- FindClusters(scData, resolution = 0.2, dims = 1:numPC)
scData <- RunUMAP(scData, dims = 1:numPC)

Both options give me different results. Not sure what the correct way to it is.

From the Broad Institute website:
https://portals.broadinstitute.org/harmony/
"Harmony is a general-purpose R package with an efficient algorithm for integrating multiple data sets. It is especially useful for large single-cell datasets such as single-cell RNA-seq."


text:
How to extract unmatched reads using bwa and samtools?
I have a single read (NOT paired) that I need to pass through the workflow described in Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.
Here in Materials and methods the procedure is described as the following:
The first step of the workflow consists of an alignment of the reads against the host genome (Fig. 2i).
This step aims at discarding reads that map to the host genome and may partially map to the viral genome after segmentation, and at reducing the working file size.
For example, MV and rMV-ΔV data sets were generated from total RNA samples of infected cells and contained mostly reads mapping the human genome (≈99%).
This step uses a combination of bwa mem and samtools view with the parameters –bS –f4.
An additional step consists of an alignment of the reads against the viral genome of interest, in order to exclude reads perfectly mapping the viral genome.
Therefore, only unmapped reads are further analyzed.
Of note, clipped reads (i.e., CIGAR motif contain S or H) are also conserved. Some of these reads may map to viral genome recombination junctions that are present in DI genomes.
I have already been suggested to use bowtie2 instead of bwa but, firstly, the output is not clear for me and secondly, I would like to test the official protocol.
Since the article suggests to use bwa and samtools on this very first step, that is what I have done so far:

(Optional) not sure if this is important but as suggested by someone I have transformed .fna --> .fa
cp GCF_000001405.39_GRCh38.p13_genomic.fna GCF_000001405.39_GRCh38.p13_genomic.fa

Indexed human genome with bwa
bwa index GCF_000001405.39_GRCh38.p13_genomic.fa

aligned my only read with indexed genome from step 1
bwa mem GCF_000001405.39_GRCh38.p13_genomic.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz > whole.sam

Transform .sam --> .bam
samtools view -S -b whole.sam -o whole.bam

Separated unmapped reads (as it is recommended in Materials and Methods using -f4)
samtools view -f4 whole.bam > sample.unmapped.sam

Converted unmapped reads into .fastq format (since this is the format used by the software later)
samtools fastq sample.unmapped.sam > unmatched.fastq

I wonder if these steps are correct?

source: Stack Exchange

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

sort: votes


text:
Why SeqIO.parse method isn't working?
I am trying to follow the tutorial from from the Biopython website here  and I am right at the beginning (2.4.1), where I am trying a Simple FASTA parsing example which is
from Bio import SeqIO
for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

For an unknown reason, I get the following error:
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-19-c64128d7dcfd> in <module>
      1 from Bio import SeqIO
----> 2 for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
      3     print(seq_record.id)
      4     print(repr(seq_record.seq))
      5     print(len(seq_record))

~\Anaconda3\lib\site-packages\Bio\SeqIO\__init__.py in parse(handle, format, alphabet)
    605     iterator_generator = _FormatToIterator.get(format)
    606     if iterator_generator:
--> 607         return iterator_generator(handle)
    608     if format in AlignIO._FormatToIterator:
    609         # Use Bio.AlignIO to read in the alignments

~\Anaconda3\lib\site-packages\Bio\SeqIO\FastaIO.py in __init__(self, source, alphabet, title2ids)
    181             raise ValueError("The alphabet argument is no longer supported")
    182         self.title2ids = title2ids
--> 183         super().__init__(source, mode="t", fmt="Fasta")
    184
    185     def parse(self, handle):

~\Anaconda3\lib\site-packages\Bio\SeqIO\Interfaces.py in __init__(self, source, alphabet, mode, fmt)
     45             raise ValueError("The alphabet argument is no longer supported")
     46         try:
---> 47             self.stream = open(source, "r" + mode)
     48             self.should_close_stream = True
     49         except TypeError:  # not a path, assume we received a stream

FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'

I am using a Jupyter Notebook.

There's something wrong with your copy of https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/ls_orchid.fasta

Is there any way not to download a .fasta file every time but say scrape it directly from pubmed genome? This string https://www.ncbi.nlm.nih.gov/nuccore/FJ211590.1?report=fasta isn't accepted as a valid argument. Thank you.

Just download the file

Thanks, this means there's no other way?

If you wanted to program in a cache you could, but that’sd be excessive IMO.

Personally I think its late and we're all tired.

"FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'"

The traceback is pointing to line 2, the SeqIO as you point out.
In my opinion, the script should be in the same directory as the fasta file. You know how to identify the path in bash (pwd or echo $PWD), although you must be using Windows (Anaconda has \ rather than /). Alternatively the file names has a typo. You know the
for seq_record in SeqIO.parse("~\path\ls_orchid.fasta", "fasta"): # where path is the dir(s) leading to ls_orchid.fasta, but obviously use / if its Linux.

I would assume you can alternatively dump the .ipynb file in the location where the fasta file is.
As personal choice I prefer a formal IDE such as 'Visual Studio Code' or "PyCharm". In my personal opinion, I am aware some will disagree, Jupyter notebook good is for running code on a remote commercial server for data science work. However, the traceback - which you clearly have activated - looks good. I can see that you're Anaconda3 installation looks good, you must be using Windows (which I've never worked with) and you're running Jupiter from Anaconda3 ... so cool.

Thanks, Michael, these are my first steps with biopython - thought it would work ok with win10 jupyter nb, and yeah checked the path - it was wrong so now everything's fine. Agree 100% that Linux is better.

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


text:
Why does the FASTA sequence for coronavirus look like DNA, not RNA?
I'm looking at a genome sequence for 2019-nCoV on NCBI. The FASTA sequence looks like this:
>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
TAATTACTGTCGTTGACAGGACACGAGTAACTCGTCTATCTTCTGCAGGCTGCTTACGGTTTCGTCCGTG
...
...
TTAATCAGTGTGTAACATTAGGGAGGACTTGAAAGAGCCACCACATTTTCACCGAGGCCACGCGGAGTAC
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

Coronavirus is an RNA virus, so I was expecting the sequence to consist of AUGC characters. But the letters here are ATGC, which looks like DNA!
I found a possible answer, that this is the sequence of a "complementary DNA". I read that

The term cDNA is also used, typically in a bioinformatics context, to refer to an mRNA transcript's sequence, expressed as DNA bases (GCAT) rather than RNA bases (GCAU).

However, I don't believe this theory that I'm looking at a cDNA. If this were true, the end of the true mRNA sequence would be ...UCUUACUGUUUUUUUUUUUU, or a "poly(U)" tail. But I believe the coronavirus has a poly(A) tail.
I also found that the start of all highlighted genes begin with the sequence ATG. This is the DNA equivalent of the RNA start codon AUG.
So, I believe what I'm looking at is the true mRNA, in 5'→3' direction, but with all U converted to T.
So, is this really what I'm looking at? Is this some formatting/representation issue? Or does 2019-nCoV really contain DNA, rather than RNA?

Excuse my ignorance, but I expected a lot more AAA (3' poly(A) tail) at the end of that sequence. Am I confusing multiple formats here or is your sequence missing a part? Is this the result of transcribing RNA as DNA?

should be 11 triplets, 33 times A

That is the correct sequence for 2019-nCov. Coronavirus is of course an RNA virus and in fact, to my knowledge, every RNA virus in Genbank is present as cDNA (AGCT, i.e. thydmine) and not RNA (AGCU, i.e. uracil).
The reason is simple, we never sequence directly from RNA because RNA is too unstable and easily degraded by RNase. Instead the genome is reverse transcribed, either by targeted reverse transcription or random amplification and thus converted to cDNA. cDNA is stable and is essentially reverse transcribed RNA.
The cDNA is either sequenced directly or further amplified by PCR and then sequenced. Hence the sequence we observe is the cDNA rather than RNA, thus we observe thymine rather than uracil and that is how it is reported.

It is really hard to sequence directly from a patient It has been done (Direct RNA Sequencing of the Complete Influenza A Virus Genome) but of course that doesn't change the correct answer here - Virtually all sequencing of RNA viruses is done through cDNA.

Hi @iayork, thanks I accept its doable, but e.g. one risk is one sloppy wet-lab mistake and the vRNA is degraded, e.g. vRNA is stored at -70oC and needs RNase inhibitors, viral cDNA is stored at -20oC and no-one cares. Sometimes there are legal issues (these are pathogens).

If this were [cDNA], the end of the true mRNA sequence would be ...UCUUACUGUUUUUUUUUUUU, or a "poly(U)" tail.

A cDNA sequence, maybe confusingly, refers to the coding strand of the cDNA (despite being called “complementary”). So while cDNA is the result of reverse transcribing RNA into DNA, by convention it has the same strandedness as the original RNA. That’s why what you’re seeing is read in 5′→3′ direction and contains a visible poly(A) tail. Having a single conventional reading direction for all archived sequences vastly simplifies data handling, and reduces errors.
In fact, since cDNA is double-stranded, there is no a priori reason why a computer-stored cDNA sequence should refer to the template strand (i.e. the opposite strand, which is synthesised from the RNA during reverse transcription).
The whole (simplified) synthesis process of cDNA is as follows:

A primer hybridises to the template RNA molecule.
The RNA template is reverse transcribed into DNA using reverse transcriptase.
The RNA template is removed.
A complementary strand is transcribed along the (currently) single-stranded cDNA, resulting in a double-stranded cDNA product.

It's not common to sequence directly from RNA because most sequencing platforms don't have that as an option. Nanopore sequencers do allow this, but I'm not aware yet of any 2019-nCov preprints involving nanopore RNA sequencing. I expect that will change in the next month or so.
Commercial kits exist; there are no insurmountable technical issues with it. Direct RNA sequencing can be done locally, on-site near the point of discovery without sample transfer or culture on a USB-powered device that fits in a pocket (RNA preparation takes about 2 hours). Flow cells that have potentially-infectious RNA inside them can be disposed as biohazard waste. However, the ease at which RNA can be quickly converted to more stable cDNA then amplified to create a much higher concentration DNA sample (which is quicker / more efficient to get results from) means that cDNA is generally preferred for sequencing unless the native RNA is needed (e.g. for looking at RNA base modifications that are destroyed when converting to cDNA).
There is a paper on coronavirus direct RNA sequencing with nanopore here; I would expect that 2019-nCoV would have a similar difficulty. The zika virus has an extremely low viral load in human blood, but has also been sequenced via direct RNA sequencing of [carefully] cultured cells (see here).
Regardless of whether or not RNA sequencing has actually been carried out, most genetic data analysis programs will only work with A/C/G/T sequences, so it's conventional to replace any U parts of an RNA sequence with T for data storage. There's no loss of information by doing this, as T replaces all Us in the RNA sequence.


text:
A new paper suggests the Corona Virus has "Uncanny similarity of unique inserts in the 2019-nCoV spike protein to HIV-1" - What does this mean?
Quote:

We found 4 insertions in the spike glycoprotein (S) which are unique to the 2019-nCoV and are not present in other coronaviruses. Importantly, amino acid residues in all the 4 inserts have identity or similarity to those in the HIV-1 gp120 or HIV-1 Gag. Interestingly, despite the inserts being discontinuous on the primary amino acid sequence, 3D-modelling of the 2019-nCoV suggests that they converge to constitute the receptor binding site. The finding of 4 unique inserts in the 2019-nCoV, all of which have identity/similarity to amino acid residues in key structural proteins of HIV-1 is unlikely to be fortuitous in nature.

Source: https://www.biorxiv.org/content/10.1101/2020.01.30.927871v1

I know this is a preprint paper and it's not yet peer reviewed, but can someone tell me what the implications are if this is true? Does this mean that the virus is artificially constructed?

UPDATE: The article has now been withdrawn with the following note:

This paper has been withdrawn by its authors. They intend to revise it
  in response to comments received from the research community on their
  technical approach and their interpretation of the results. If you
  have any questions, please contact the corresponding author.

This is very odd, and will require a rigorous investigation, but my initial reaction is one of scepticism. Considering just the 1st insert, the insert sequence is GTNGTKR, short at just 7 amino acids. A simple BLASTP vs NR did not find perfect matches to HIV sequences, but did reveal 100% identity across the full sequence length to >50 other short protein sequences, which could be spurious 'chance' hits of course. Many of the top 100 hits are against eukaryotic sequences, for example one is against Pristionchus pacificus, a type of nematode worm (see alignment below). Someone needs to do a proper peer review of this preprint before any conclusions are drawn from it.
>tank-1 [Pristionchus pacificus]
Sequence ID: PDM74036.1 Length: 2481
Range 1: 1474 to 1480

Score:24.0 bits(49), Expect:2477,
Method:,
Identities:7/7(100%), Positives:7/7(100%), Gaps:0/7(0%)

Query  1     GTNGTKR  7
             GTNGTKR
Sbjct  1474  GTNGTKR  1480

UPDATE 1: Things are moving fast. There are now 10 comments under the preprint, all of them critical of the idea that these inserts are meaningfully similar to HIV-1 sequences for the reason I outlined above (i.e. short sequences with many hits against organisms from across the tree of life) and additional points, such as small insertions being quite normal evolution for RNA viruses. There are also similar critiques on twitter and biorxiv has added the following disclaimer header to their website:

bioRxiv is receiving many new papers on coronavirus 2019-nCoV.   A
  reminder: these are preliminary reports that have not been
  peer-reviewed. They should not be regarded as conclusive, guide
  clinical practice/health-related behavior, or be reported in news
  media as established information.

I have never seen such rapid (and unanimous) post-publication peer review before. It is therefore already (just 1 day after publication) fairly clear that at least some claims in this study are completely false and there is no good evidence to support the claim that the 2019-nCoV have acquired sequences from HIV or evidence to suggest the virus is engineered.

Thank you for this. It really makes me wonder why a respected scientist would publish such an obviously flawed paper though?

Blast isn't the issue here, two proteins can been unrelated by blast but retain the same function. Functionally these proteins have similarities. The universal role of indels or what the authors call "inserts" particularly in viral surface antigens is widespread and has never been resolved. It is wrong to say the paper is "flawed", I think the observation is spurious, but has value.

@SurpriseDog I don't know the motivations of the authors but the article has now been withdrawn

Normally "inserts" used in the manuscript are "indels" in protein alignments, short for insertions and deletions.
What I think has happened is a group investigating indels in HIV env noticed indels in 2019-nCov. Essentially I think the correlation is spurious - but I haven't tested it, though the area of research in understanding indels is certainly valid and important.
What is certain is that indels induce a large structural change to a protein structure and any Gibbs free-energy style calculation will identify this.
Vaccine
The spike protein will be the primary candidate to make a 2019-nCov vaccine and that is a very important reason why the sequence was rapidly released. So it is an important protein and the structural changes indels induce mean that a SARS vaccine will probably not provide much protection against 2019-Cov, even apart from the amino acid divergence (below).
Differences 2019-nCov vs HIV
In summary, a lot. HIV env and particularly HIV gag are very different from coronaviruses, both in the mechanism of genome replication, coronavirus never leaves the cytoplasm, clinical outcomes, tissue tropism and duration of infection.
Similarities
HIV env and the glycoprotein spike of coronaviruses are the receptor binding protein to gain entry into a cell. They are called structural proteins. Entry to a cell can be blocked by antibodies and these antibodies are called "neutralizing antibodies". Neutralizing antibodies are catastrophic for a virus. Other antibody responses can be effective, such as IgM, but to clear an infection just using antibodies, you need neutralising antibodies. Both HIV env and the coronavirus spike protein are subject to neutralising antibodies. HIV gag has nothing to do with HIV env, in terms of function or antibody exposure. This is why the spike protein will be the primary vaccine candidate for a subunit vaccine.
Coincidence, law of chance
There is a large variation of indels in HIV env within HIV and what the authors are inferring is that there is a resemblence to that between SARS and 2019-nCov. In my opinion this is a coincidence, because they are comparing a large repertoire of HIV varients against a single indel pattern in the coronaviruses.
Why coronavirus indels?
That is a very good question. Generically indels in viral surface antigen genes are common, much more common in other proteins - such as those involved in virus replication (non-structural proteins). The amino acid identity between SARS and 2019-nCov is 80%, and in any virus, such as flaviviruses 80% identity means indels will be present in surface antigens between the viruses. The answer is it is not unusual in any RNA virus to see indels at a comparatively large amino acid divergence.
What function could they serve?
I've briefly looked at indel bioinformatics between flaviviruses (Zika virus, yellow fever virus etc..) notably using envelope (E) protein sequences, and they also occur between African Zika viruses in the E-protein. E-protein being the equivalent of coronavirus spike protein, the receptor-binding protein. No-one has ascribed a function to them and that is the problem with this manuscript.
Hypotheses

One theory is that a structural change in the protein will occur to stop antibody binding.
Another theory is they have functional differences, such as cell tropism

Bioinformatically separating the two theories is extremely hard without wetlab experimentation.

Really thorough explanation, clear and precise !! I learned a lot thanks ;)

"Neutralizing antibodies are catastrophic for a virus. " why because target the T cells to the virus? and if the virus can enter in T cells destroy them and reproduce too what happen???

Good common. No  HIV targets T-cells. I will update the post when I have time "neutralising" means preventing the virus binding to the cellular receptor, thus it cannot enter a cell and therefore cannot replicate

There's a new paper by a Chinese Virologist who worked at a WHO reference lab and has now fled to America claiming whistle blower status. Here is the paper: https://zenodo.org/record/4028830#.X19xByXZglR - Any thoughts? I'm hesitant to post another question about stuff like this, being so far out of my depth.

I'll have a look when I've a moment

SurpriseDog - I suspect the wording "unlikely to be fortuitous in nature" is what led you to ask "Does this mean the virus is artificially constructed?" - No!
I think the writing here is imprecise and confusing. There is no implication that these are human-generated (non-natural/engineered) insertions but rather I think they are implying that these mutations are likely to have functional implications. Whether there is any relation to HIV is questionable.
There has been some political finger pointing and threats of full-scale investigation into whether the virus escaped from a lab. Hard to say, but current thinking is that SARS-Cov-2 is a natural isolate, not engineered (see https://www.nature.com/articles/s41591-020-0820-9). It was probably transmitted from bats.


text:
a package with randomization/permutation tests for genomic analyses
please would you let me know, is there a package in R/python/C that has implemented randomization/permutation tests for joint genomic analyses
(an example of joint genomic analyses -- when jointly considering both GENE EXPRESSION and PROTEIN BINDING along the DNA).
The context of my question is the following :
let's consider 1000 UP_REGULATED genes with increased PROTEIN X and with a HISTONE MARK Y :
in order to show that PROTEIN X is related to HISTONE MARK Y for 1000 UP-REGULATED genes, what "controls" would you use for the comparison :
-- 1000 RANDOM GENES (and multiple randomization tests)
-- 1000 UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y
-- 1000 UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y
-- 1000 UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y
-- 1000 NOT-UP-REG GENES with NO PROTEIN X, NO HISTONE MARK Y
-- 1000 NOT-UP-REG GENES with PROTEIN X, and NO HISTONE MARK Y
-- 1000 NOT-UP-REG GENES with NO PROTEIN X, and with HISTONE MARK Y
-- anything else ?
any suggestions are very welcome. thanks a lot,
bogdan
ps : 've posted a similar question on BioC website : https://support.bioconductor.org/p/128828/

I would simply look at co-variance across the entire dataset and examine values >0.8 in the first instance rather than select 1000 genes.
With your approach sampling 1000 genes can be criticised for selective sampling and the way to get around this is to perform replicates, with (bootstraps) or without replacement. You could produce e.g. 1000 replicates of 1000 genes and perform co-variance on each replicate to examine Y = X.
Python has a specific permutation/combinations module call itertools, which is different to random resampling and doesn't need replicates
from itertools import combinations
X = [1,2,3]
Y = [1,2,3]
print (list(combinations(iterable, 1000) + X + Y))

The standard way is to sample is,
from random import sample
mygenes = [1, 2, 3, 4, 5]

print(sample(mygenes,3))

Output
[2, 4, 5]

There is nuance and further funcationality ... I cannot remember is "sample" is with or without replacement (its quite a big deal)

thanks a lot, Mike ! very helpful ! wondering if possible to recommend any package in R has the functions for bootstrapping, or similar functions as in itertools. thanks a lot !

Hi @bogdan, I don't use R so can't recommend anything.


text:
Can we do NGS library preparation using UMI with large amount of DNA input?
We know that in next-generation sequencing (NGS), the unique molecular identifier (UMI) can reduce or eliminate sequencing or PCR errors and result in very high accurate data. Therefore, UMI is widely used in rare mutation detection. A typical application is cancer early detection or minimal residual disease (MRD) detection that try to find the very rare DNA from cancer cells in the blood.
However, we also know that the cell-free DNA in the blood is very low, typically 5 ng in 1 mL blood. And it usually takes about 10 mL blood to get about 50 ng DNA to get sequenced with UMI. Limited DNA input leads to many false negatives.
So I came up with an idea: What if I got a sample with plenty of DNA? Can I detect rare mutations using UMI in this sample?
The UMI library preparation protocols I have met all use a small amount of DNA input, say 50ng. If I have 1μg DNA, do I have to cut it to 50ng to prepare the UMI library?

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

Do you mean that the variant allele frequencies (VAFs) of the variants will be lower than 1% for the very rare DNA from cancer cells in the blood?

There are around 3000 copies of haploid genome present in 10 ng of DNA, approximatively (50-)60 ng of cell free DNA will be required for a sensitivity of 0.01% (one rare event in 10,000 molecules),  moreover,  in blood, average cfDNA concentration in blood of cancer patients, varying between 0 and 1000 ng/mL, with an average of 180 ng/mL

Whole genome sequencing uses between 50 ng and 1 μg of high-quality DNA (A260/A280: 1.7–2.0) for PCR-amplified or non-amplified libraries if your focus is precision oncology as described in the following open access review in Seminars in Cancer Biology paper by Meggendorfer et al 2022

The question is

"Do I add 20-fold more DNA than the kit states?"

Its not my area - this is wet-lab - but there will be a tolerance however 20-fold is excessive (>1 log), i.e. no it doesn't help and could cause (expensive) downstream problems. The people to contact are the support for the given kit, there are unlikely to support your idea.
The answer isn't wet-lab but almost certainly is bioinformatics and detecting low frequency somatic driver mutations. The Broad Institute worked on this problem for years, its nothing new, i.e. identifying driver mutations in highly mixed sample comprising healthy somatic tissue and cancer cells.
I've tracked their work and have a reasonable understanding of what they do and why: oncology isn't my field. Your in-house bioinformatician should know this stuff.
One suggest is using MuTect from the Broad Institute, I assume there have been many custom solutions since then. Their famous product is GATK for VCF, I don't think thats the solution here.
Again cancer isn't my area BTW this is just common knowledge within bioinformatics, i.e. knowing what someone else is doing and why is part of the knowledge base.
Please do try to retain a bioinformatics focus for future questions.

Do you mean the following type of technology?

As it’s currently written, your answer is unclear. Please [edit] to add additional details that will help others understand how this addresses the question asked. You can find more information on how to write good answers in the help center.


text:
which NCBI tool is optimized to identify a species from a DNA fragment?
“The Iceman” was a man who lived 5300 years ago and whose body was recovered from an Alpine glacier in 1991.  Some fungal material was recovered from his clothing and sequenced.
Ice man : found as a frozen corpse in the Val Senales Glacier in Italy.
Thought to be carrying Fomitopsis betulina (formerly known as Piptoporus betulinus) commonly known as birch polypore, birch bracket, or razor strop) is a common Basidiomycota brown rot macrofungus growing on decaying birch wood.
Other species known : Fomes fomentarius (commonly known as the tinder fungus, false tinder fungus, hoof fungus,] tinder conk, tinder polypore or ice man fungus)
Given two DNA fragments from NCBI genbank, try to find what species these DNA fragments belongs to or what species these fragments are most closely related to.
>Z54169.1 Iceman fungal clone 'Sim C-a1' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAAACTTGGGCCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTTCCTTCTGGAGAACCTCATGCCCTTCACTGGGTGTGTTGGGGAA
CCAGGACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAAT
AATAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCG
GGGGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCAT
TTGCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAG
TCTTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGA
GAAATCAAAGTTTTTGG

>Z54155.1 Iceman fungal clone 'Sim C-a2' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAACCTTGGGTCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTCCTTCTGGAGAACCTCATGTCTTTTATTAGGCGTGTTGGGGACC
AGAACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAATAA
TAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCGGG
GGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCATTT
GCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAGTC
TTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGAGA
AATCAAAGTTTTTGG

Now I know I can just google what species of mushrooms (which is where the above background info is from). But say that I was unable to and had to just go off the DNA fragments in order to find the species.
I've tried BLASTn, but I'm not sure if it is the right tool or if I need to properly calibrate the parameters to use it.
so far, I've tried:
Under “Choose Search Set:
Database : Non-redundant protein sequences (nr)
Organisms excluded: Animalia (taxid:33208), Plantae, Prokaryota  (taxid:2), Protozoa
Organisms included: fungi  (taxid:4751)
Included: Models (XM/XP),   Uncultured/environmental sample sequences
Under Program Selection:  Optimized for Highly similar sequences (megablast)

Algorithm Parameters (General Parameters, Scoring Parameters, Filters and Masking,
PS/PHI?DELTA BLAST)

are all kept default.
```

So, what were the results? What makes you think this isn't the right approach? Did you not find any good hits?

For fungal species I think normally the ITS (internal transcribed spacer) is used for identification instead of 16s, 18s in your case.

My conclusion from the tree is the fungi appears to be genuine and not an extant contaminant (which is very easily done). It is not Fomitopsis betulina because this is division Basidiomycota. The fungi here is division Ascomycota and these are very distinct from Basidiomycota.
There have been many examples of contamination in ancient DNA work, but this doesn't appear to be one of them. The rationale is given below.

The first fragment was not interesting and described 100% identity to ascomycetes. The tree above is a simple distance tree from the second fragment of 18S which shows it is an outgroup to the ascomycetes (green terminal taxa). The ancient fungi contains 3 (I think) unique nucleotide residues that have not been sequenced in any extant ascomycetes to date and sequences are unique to any known species. The location of the 3 unique nucleotides is tightly clustered within the locus with 100% identities either side. This is significant because it demonstrates it is not a PCR mosaic artefact (yep its been done in ancient DNA work and subsequently proven artifactual). In the latter PCR mosaic artefact I'm referring to (which I wouldn't cite), was a fungi/bacterial mosaic from memory.
Given the homology of 18S in the extant ascomycetes I would say this appears to be a new genus, however you would need to be sure this was not sequencing error and this is very old DNA, which is subject to degradation even when frozen. The implication of the tree is quite unusual, it implies widespread extinction of this group of fungi in recent history (5000 years is slight in evolutionary terms) BUT it could be a sampling bias, i.e. no one has sequenced those fungi and NCBI is not representative.
I assume that you have performed phenotypic/morphological analysis on the spore to believe it is Basidiomycota.
The fragment in yellow is the query and the sequence next to it was the sequence already deposited on Genbank of this ancient isolate. You would need to look up the biology and ecology of ascomycete to understand its significance, but it is in a lineage of its own, so maybe difficult to interpret by comparison to extant species. You could concatenate the two alignments based on each fragment, but the problem with 18S is that is a hard to model and you'd need to perform secondary RNA structure modelling to assess whether these mutations would be viable.
In conclusion, it definitely looks genuine but it is not the species you think it is, in fact not even close.


text:
Get nucleic acids' chain names from PDB API
There is a place in my pipeline where I call PDB's API with a ribosomal id to get all the subchains of that molecule programmatically
I'd like to separate subchains that are proteins from nucleic acids in the response, or at least know how to get the chain-ids of the nucleic acids so that i can exclude them from all manually.

Wondering if there is some query parameter like is_nucleic_acid in the PDB's ADL, but that's just a hunch. I'd be happy to learn if there is a more idiomatic way to do it( from Uniprot mappings maybe? ).
My query as of now looks like when i want to get, say uniprot acccessions for subchains(in fl):
const PDBSearchSelect = "https://ww.ebi.ac.uk/pdbe/search/select?"";
const params = {
// this is of the form "3j9m", "5afi", etc...
q: `${input.toLowerCase()}`,
//wondering if there is an "fl" parameter to get all nucleic acids
fl: "molecule_name, uniprot_id, uniprot_accession, struct_asym_id,...other",
//
rows: "200"
};

axios.get(PDBSearchSelect, {params}).... => ..

I would actually be super grateful if somebody could elucidate what's the difference between

https://www.ebi.ac.uk/pdbe/api/... and

`https://www.ebi.ac.uk/pdbe/search/pdb/select?

I've been looking at sides of the PDB API for a few weeks now and it still seems like complete babylon in terms of naming. I realize that it grew organically over years.
Thanks a ton.

Not sure I follow the first question you ask in your image about 'protein vs. NA' https://www.ebi.ac.uk/pdbe/api/pdb/entry/molecules/3j9m gives "molecule_type":"polyribonucleotide" and synonym":"16S ribosomal RNA" for "in_chains":["A"].

Also, the European PDB and the San Diego PDB have their own APIs each with different parameters. For some things one is better, for other the other. (or SIFTS). So sometimes its helpful to check what the other does. The San Diego PDB API is not as good and is documented in XML and in this case it is less good, but sometimes it is.

Namely whereas https://www.rcsb.org/pdb/rest/describeMol?structureId=3j9m has the RNA chain elements in this XML clearly have the attribute type="rna", the json (changing /rest/ to /json/ gives JSON, undocumented AFAIK) does not: https://www.rcsb.org/pdb/json/describeMol?structureId=3j9m.

@MatteoFerla, this is it, I think. You are the MVP.
Remember you from my other question, actually, and grateful to learn. The opaque documentation on many of these is quite frustrating at the same time.

@MatteoFerla, are there reasons that you know of other than historical as to why it all seems fragmented? Is that even true or am i just complaining?
I frequently dream of committing all of this to a huge intelligible GraphQL schema instead of dumping hours on trying to figure the right endpoint and query string every time i need to employ a new facet of a molecule.

No idea sorry. But I totally agree with your frustration. I needed to find the offset of the PDB residue numbers (not pose id) relative to Uniprot canonical transcript. SIFTS partially did it. I ended doing it manually, which was way quicker anyway...


text:
What exactly does each of InterPro, PANTHER, Pfam bring to the table individually in classifying a protein?
I would be very grateful if somebody could sketch out the methods Pfam and PANTHER use to assign a family to a given protein and how they are different. My (cursory) understanding is that InterpProp pools databases like Pfam and PANTHER and provides a meta-classification based on some combination of their data.
In particular, I wish to understand where the discrepancy comes from described in the picture below?

I'm fairly new to the ecosystem, so feel free to throw any bone this way.

That is a low homology sequence where different thresholds are applied. Interpro is more narrow than PFam. Also, at a glance looks like a 3 residue repeating pattern (helix) featuring tryptophan and leucine, but it has two prolines in it, so it's probably a linker that sticks to the side of the protein...

Hey, @MatteoFerla! I hope you've been well. I'm following up directly on your suggestion from once upon time. I ended up using PFAM to convert to my "yet another" nomenclature  in the prototype, but as you said, some repeat proteins get assigned multiple families, which muddies the water a bit.
On top of that i've been fiding discrepancies between the families pfam's api provides and their static databases and then there also some proteins which are already classified on InterPro(with pfam as source) but have no entries on PFAM

That is all to say that I'm super interested in how to make my nomenclature converter more robust and in how these databases/classifiers work( i know pfam is based on HMMs). Say, how can it be that InterPro is narrower than the sum of the databases which it gathers under its umbrella for example.
If you know of an overview or can throw me a few clue-words to follow up on, i'd be supper happy.

Where a domain actually truly ends from a structural point of view is not trivial. The SGC (Structural Genomics Consortium) in their semi-automated pipeline for protein crystallisation actually makes multiple constructs generally cutting the ends by half a dozen residues before or after —it sounds silly way of doing it, but it works...

About the narrowness... It is just for personal observation as I do not understand how InterPro works with its "signatures" and its blending of databases. In this case the Interpro annotation is a blend of the PFam and the Panther and the latter is longer. If you get something that is more distant, you may get a PFam but not a Panther, so the InterPro won't cover it.


text:
Trouble optimizing Five Parameter Logistic (5PL) Standard Curve for ELISA data using Python
I am writing a Python script to automate a 5PL Standard Curve for ELISA data, using an XLSX file template. I am using pandas dataframes to hold the 96 wells as different rows, then in a separate dataframe, holding Standard Curve data (not all wells need a column for known concentration). I am also using scipy and numpy for analysis, then using matplotlib.pyplot for plotting the curve (to eventually be output to the XLSX), though I am open to suggestions.
So far, the program correctly takes the input, generates the aforementioned dataframes, and makes initial guesses for the 5PL's parameters. I was attempting to modify the 4PL code given here to work with my standard curve, using all actual standard replicates to optimize the curve, and obviously using 5PL instead of 4PL. However, when I use scipy's least_squares function for parameter optimization, I receive the following errors/traceback:
C:\Users\Jake\Anaconda3\lib\site-packages\ipykernel_launcher.py:6: RuntimeWarning: overflow encountered in power

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-6-70b754f01d56> in <module>
      4
      5 # Fit equation using least squares optimization
----> 6 plsq = least_squares(fun = residuals, jac = 'cs', x0 = p0, bounds = (0.0001, 100000), args=(y_meas, x))

~\Anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    793
    794     if not in_bounds(x0, lb, ub):
--> 795         raise ValueError("`x0` is infeasible.")
    796
    797     x_scale = check_x_scale(x_scale, x0)

ValueError: `x0` is infeasible.

I was initially receiving these errors with the aforementioned article's suggested scipy.least_sq function, but I assumed it was minimizing exponents to 0 and encountering an overflow, so I swapped to scipy.least_squares instead due to its bounds argument. Code in question is below. Any help would be greatly appreciated! Thank you so much! And if this question is better suited to stats.SE, just let me know!
(I did find lmfit-py, though I was hoping to complete this project on my own if possible to get a better understanding of the regression. Still open to using it instead if that seems to be a better option.)
import matplotlib.pyplot as plt
from scipy.optimize import least_squares

def logistic5(x, A, B, C, D, E):
    """5PL logistic equation"""
    return D + ((A-D)/(np.power((1 + np.power((x / C), B)), E)))

def residuals(p, y, x):
    """Deviations of data from fitted 5PL Curve"""
    A, B, C, D, E = p
    err = y - logistic5(x, A, B, C, D, E)
    return err

def peval(x, p):
    """Evaluated value at x with current parameters"""
    A, B, C, D, E = p
    return logistic5(x, A, B, C, D, E)

x = []
y_meas = []
for ind, row in std_df.iterrows(): # Iterate through rows in std dataframe.
    for rep in std_df_rep_cols: # Iterate through replicate columns.
        if not pd.isnull(row[rep]): # If there is a value in the cell...
            x.append(row['conc']) # Add an iteration of the respective conc. to the x list.
            y_meas.append(row[rep]) # Add the cell's value to the y_meas list.

x = np.asarray(x)
y_meas = np.asarray(y_meas)

print(len(x))
print(len(y_meas))

print(x)
print(y_meas)

# Initial approximations of parameters
A = np.amax([np.amin(minBlk_ser)] + [0]) # Min asymptote
D = np.amax(y_meas) # Max asymptote
B = (D - A) / np.amax(x) # Steepness
C = np.amax(x) / 2 # Inflection Point, conc at which y = (D - A) / 2
E = 1 # Asymmetry factor
print(A, B, C, D, E)

p0 = [A, B, C, D, E] # List containing initial guesses as arg for least_squares
y_true = logistic5(x, A, B, C, D, E) # Assumed "true" curve based on initial params

# Fit equation using least squares optimization
plsq = least_squares(fun = residuals, x0 = p0, bounds = (0.0001, 100000), args=(y_meas, x))

Very interesting code. I can't see where you load import pandas or the data frame, but from memory iterrows and in particular pd is obviously a pandas command. The isnull is a bit hairy because most blanks will display an OD. Anyway with regards your bug .. there is a standard package used to input data for regression analysis, it is very famous (I would need to look it up) ... this would obviate using a list, which is what you are doing to define X0.

Hi, @MichaelG.; thank you for the answer! I didn't include the pandas import/dataframe initialization in this question as they are earlier in the code. I can show them if it could help? The "isnull" is used not so much for blanks (measured and subtracted earlier in the code), but for missing replicates in the standard dataframe. It was a proactive choice in case, for instance, someone loaded wells with extra replicates for some (not all) standard concentrations and thereby left empty cells in the respective replicate column in std_df. Very nice catch! Is lmfit the package you're thinking of?


text:
How can I get the estimate of lmFit function form limma package
This is the first time that I used lmfit function from the limma package

I'm a little bit confused on how to interpret the result.
is there a way to get the estimate or the beta parameter of lmfit

Normally you run the results through eBayes() and get the coefficients from that.


text:
What is a proper way for random subsampling of metagenomic data?
Let's say we have a metagenomic sample that is paired-end FASTQ files including 10,000,000 DNA reads collected using shotgun sequencing.
How would one make a random subsample of the mentioned metagenomic sample with for example 1,000,000 reads? I know about the seqtk function that could be used in such cases, but my question is more about whether it is meaningful to represent a metagenomic sample with a small subsample generated using this way or it has to satisfy certain criteria, e.g. in terms of the number of reads and diversity to be a valid representative of the metagenomic sample?
Knowing that some species might be very rare in the sample, would it ever be proper to do subsampling on metagenomic samples?

To make more samples to be used for machine learning

Ah ... data augmentation, yep thats a common approach. I personally can't quite see what teh problem is, if you can perform ML surely you can use a randomisation module, e.g. Python's random or sample?

I am not NGS expert, I known ML. Essentially you are performing data augmentation which is essential for ML. What I question is whether a replicate of e6 from a population of e7 would really be sufficient for ML because you need sample sizes (replicates in this case) of >100000.
The answer is really simple you bootstrap the data (sample with replacement), to obtain your replicates. Numpy has a random module that will do this,
p.random.seed(345)
pop = np.random.randint(0,10000000 , size=1000000)

Use pop to select from fastaq. The other way in is to define each fastq sequence within pop then simply,
sample = np.random.choice(pop, size=1000000)

Or a combination of the two. If you can do ML you can do Python coding.

Dear Michael, thank you for the answer, but what I am asking is not how to do random subsampling but the idea whether it makes sense to do it for metagenomic data where it may remove key information of rare taxons. I don't think that subsampling should be done with replacement anyway since the subsamples would not be considered independent for training and testing.

Ok, in that case I don't understand the rationale. Sampling is used heavily in machine learning, it is also used in statistical analysis such as Monte Carlo and  null distributions. Generally it is used to ensure the observation is an outlier (except for ML)


text:
Are there any "canonical" sources of C. elegans connectome?
I found that there are a couple of projects that allow the user to simulate the neural system of C elegans, for example OpenWorm one update being here.
I don't know how to navigate sources in bioinformatics and thus struggle to evaluate their quality.
Question Are there any fundamental / most widely accepted data sources for raw data of C Elegans connectome or electron microscope scans of this organism, similar to Human Genome Project? If not, what are the ways to derive the data?

Answer from @terdon, converted from comment (please improve if possible):
Have you checked on WormBase? That's the main repository for c. elegans data. I don't know that they have connectome data, but if they exist, that's where they will be stored.
[see here for more information on the purposes of comments in the StackExchange network]


text:
Converting Pairwise single- linkage clustering distance data to "newick" format
I have pairwise distance data for single linkage cluster and I would like to convert it to newick format.
I was unable to find a conversion algorithm.
I require newick format for visualization and annotation of a phylogeny.
Here is the output.
INFO  : ======================================
INFO  : Pairwise single linkage clustering
INFO  : ======================================
INFO  : Hierarchical Tree
INFO  : ======================================
INFO  : Node     Item 1   Item 2      Distance
INFO  :     0 :       53       52        0.001  5zxf  5z02
INFO  :    -1 :       64       60        0.002  7bj6  7bir
INFO  :    -2 :       62       -1        0.002  7biv
INFO  :    -3 :       34       26        0.002  4wt2  4ode
INFO  :    -4 :       65       -2        0.003  7bmg
INFO  :    -5 :       25       13        0.003  4occ  4erf
INFO  :    -6 :       29       -3        0.009  4ogt
INFO  :    -7 :       39       38        0.010  4zyi  4zyf
INFO  :    -8 :       50       48        0.011  5oc8  5ln2
INFO  :    -9 :       54       -8        0.012  6ggn
INFO  :   -10 :       42       41        0.015  5hmi  5hmh
INFO  :   -11 :       -4       61        0.020            7bit
INFO  :   -12 :       70       40        0.028  7qdq  5c5a
INFO  :   -13 :       47       68        0.032  5laz  7na3
INFO  :   -14 :       -9       45        0.035            5law
INFO  :   -15 :       49       59        0.038  5oai  6q9o
INFO  :   -16 :       32       -5        0.039  4qo4
INFO  :   -17 :      -15        8        0.044            3uu1
INFO  :   -18 :       23      -16        0.047  4oas
INFO  :   -19 :      -14       -7        0.053
INFO  :   -20 :       24      -18        0.057  4oba
INFO  :   -21 :      -19        6        0.057            3lbl
INFO  :   -22 :       43      -10        0.058  5hmk
INFO  :   -23 :       44       19        0.059  5lav  4jvr
INFO  :   -24 :      -23      -20        0.061
INFO  :   -25 :      -11       15        0.061            4hg7
INFO  :   -26 :       57      -24        0.062  6q9h
INFO  :   -27 :       27       -6        0.065  4odf
INFO  :   -28 :      -25      -17        0.067
INFO  :   -29 :       14      -28        0.067  4hbm
INFO  :   -30 :      -13      -21        0.067
INFO  :   -31 :      -27      -29        0.067
INFO  :   -32 :      -26      -30        0.068
INFO  :   -33 :       58       21        0.068  6q9l  4mdn
INFO  :   -34 :      -31      -32        0.070
INFO  :   -35 :       17        5        0.071  4jv9  3lbk
INFO  :   -36 :       30      -34        0.074  4ogv
INFO  :   -37 :      -33        2        0.077            1t4e
INFO  :   -38 :       35      -36        0.079  4zfi
INFO  :   -39 :       20      -38        0.086  4jwr
INFO  :   -40 :      -39      -37        0.094
INFO  :   -41 :        7      -40        0.106  3tj2
INFO  :   -42 :       51      -41        0.111  5trf
INFO  :   -43 :       69      -42        0.112  7na4
INFO  :   -44 :       56      -43        0.113  6q96
INFO  :   -45 :        4      -44        0.127  3jzk
INFO  :   -46 :       33      -45        0.129  4qoc
INFO  :   -47 :       22      -46        0.137  4mdq
INFO  :   -48 :      -12      -47        0.151
INFO  :   -49 :       37      -48        0.155  4zyc
INFO  :   -50 :       31      -49        0.156  4oq3
INFO  :   -51 :       12      -50        0.158  4ere
INFO  :   -52 :      -22       46        0.161            5lay
INFO  :   -53 :      -35      -51        0.182
INFO  :   -54 :       55      -53        0.186  6i29
INFO  :   -55 :       10      -54        0.224  3w69
INFO  :   -56 :       18      -55        0.238  4jve
INFO  :   -57 :        9      -56        0.265  3vzv
INFO  :   -58 :       28      -57        0.271  4ogn
INFO  :   -59 :      -58        1        1.000            1rv1
INFO  :   -60 :       66      -59        1.000  7na1
INFO  :   -61 :        0      -60        1.000
INFO  :   -62 :      -52      -61        1.000
INFO  :   -63 :       11      -62        1.000  4dij
INFO  :   -64 :       67      -63        1.000  7na2
INFO  :   -65 :       16      -64        1.000  4jv7
INFO  :   -66 :       36      -65        1.000  4zgk
INFO  :   -67 :        3      -66        1.000  2lzg

The problem with the tree format presented is that it does not conform to a known standard. Whilst its possible to write a parser I don't think its a good idea.
This output is hierarchical clustering, UPGMA, that sort of thing. Its not phylogenetics because this is a linkage tree.
Can you access the distance matrix prior to the clustering? If so the information below the line will help.
The output should provide the raw distance matrix prior to hierarchical clustering and will certainly have generated it in the calculation. BTW that particular clustering method isn't cool compared with nj (below) for complex reasons.
The calculation will be:

Perform pairwise linkage analysis
Perform cluster analysis on the resulting matrix

You need neighbor-joining to convert the distance matrix into a newick-format. There is no possibility of a direct conversion.
If you are using R its in the ape library via
treenk <- nj(M) with M is the matrix. The treenk will be a format that can be read straight into ggtree
Briefly, a distance matrix requires a clustering algorithm and neighbor-joining is by far the best clustering method developed in phylogenetics.
The following is a phylogenetics example ... this is different from the example you have 'cause thats about linkage.
data(mynucleotides)
tr <- nj(dist.dna(mynucleotides))
plot(tr)

I your case the final line will be something like ggtree(tr) ... with various options.

Thank for the idea. It helps so much.


text:
How to calculate mutation rate and mutation sites in a genome using FASTA file?
I have 6 viral genome sequences of the same virus and 1 reference sequence in FASTA format.

How I can identify mutations and mutation sites in those genomes using FASTA sequences  but how I can do this for FASTA file?
And how I can identify the mutation rate for one genome?

Background If I've FASTQ file then I'll simply align the reads to the reference and by using variant calling tool I will get the mutate sites

The mutation rate is impossible to calculate without knowing when the samples diverged. If you have fasta sequences, you can align them and build a tree (I guess beast is the most common tool for viral genomes including time calibration).

Typically you can align FASTA or FASTQ sequences.  If your variant caller of choice doesn't accept a FASTA and you trust the bases in your FASTA you could mock up the quality scores and make a FASTQ

You can still align a FASTA file with a tool like bwa mem ( if they are short reads ) or minimap2 for long reads and run it through a variant caller like freebayes. Alternatively, if you have sufficient reads to create an assembly of all of the genomes, or if they already are you can create a MSA of them using a tool like Cactus (as it is built for genomes) or MAFFT if the sequences are short enough and then use your favorite viewer to find variable sites

I've assembled genome and the genome size is about 30-32k, so I align them and after that, I try to check variable region, It will take a lot of time.

And how I can calculate mutation rates?

I don't expect it to take too long to run as MAFFT i-INS-n ( mush slower then Cactus but probably more accurate)  can align my dateset of 208 seq at 200kb  in around 20000 clock hours  (20 hrs with 10 cores ) so 6 should be fine unless they are really gappy.  I am unclear what you mean about mutation rates in this context unless these genomes have established time points


text:
How to do BLASTP for short sequences (<20 aa) effectively?
I am trying to find sequence homology between viral sequences and my protein of interest. I have the sequences of their epitopes which varies from 5 to 500 amino acids long. For shorter sequences, it is recommended to use blastp-short but I am confused until which length it is okay to run blastp?
For shorter sequences, what are the parameters needs to be added for effective results?
blastp -task blastp-short -query POI.fasta -db virus.fasta -out output.txt -outfmt 6 -evalue 20000

You should increase the -evalue 10 default option to a higher value because this value it's adjusted to the sequence database size. Check this table, maybe this can help you to decide which parameters are the best option for you.
    OPTION      TASK          DEFAUL    DESCRIPTION
    word_size   blastp        3         Word size of initial match.
    word_size   blastp-short  2         Word size of initial match.
    gapopen     blastp        1         Cost to open a gap.
    gapopen     blastp-short  9         Cost to open a gap.
    gapextend   blastp        1         Cost to extend a gap.
    gapextend   blastp-short  1         Cost to extend a gap.
    matrix      blastp        BLOSUM62  Scoring matrix name.
    matrix      blastp-short  PAM30     Scoring matrix name.
    threshold   blastp        11        Minimum score to add a word to the BLAST lookup table.
    threshold   blastp-short  16        Minimum score to add a word to the BLAST lookup table.
    window_size blastp        40        Multiple hits window size, use 0 to specify 1-hit algorithm.
    window_size blastp-short  15        Multiple hits window size, use 0 to specify 1-hit algorithm.

Thanks! Can you please explain what is the meaning of word_size and window_size in the parameters?

The word_size parameter is used in the initial search when BLAST makes a table of all possible "words" present in the query sequence. This first search is not only for identical "words" but for all "words" that would yield an alignment with a score above a threshold, with this step, BLASTp increases speed at the expense of only a small loss in sensitivity by only trying to extend a "hit" if there is a second "hit" within a certain window_size.

How does word_size effects my BLASTp result? If I want to find 5 or 6 amino acid match at least then increasing word_size to 5 or 6 will help? It means it will try to find 5 or 6 amino acids initial match?


text:
Loop Pymol command get_area
I want to select multiple residues from different PDB files loaded on the same pymol session.
select resi LIG
get_area sele

How do I loop the above command to enable the command below:
select ligands,resn LIG

from https://www.blopig.com/blog/2019/08/how-to-iterate-in-pymol/
Example 1: Print the names of all the selected objects
You’ve selected several objects in PyMOL, and you want to see their names. By default, a new selection object called “(sele)” is created. We can create a Python variable called seleobjs, say, and then loop over them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: print obj

Example 2: Disable all the selected objects
If you still have your selection, “(sele)”, you can disable them (effectively turn them off on the right-hand sidebar):
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.disable(obj)

Example 3: Delete all the selected objects
You’re done with the objects in your selection, “(sele)”: you might be able to delete them using “delete ”, but if you can’t write a wildcard that covers all the names of the objects, you can iterate over a selection and delete them:
seleobjs = cmd.get_object_list('(sele)')
    for obj in seleobjs: cmd.delete(obj)

from https://pymolwiki.org/index.php/Get_Area
get_area calculates the surface area in square Angstroms of the selection given. Note that the accessibility is assessed in the context of the object(s) that the selection is part of. So, to get the surface areas of e.g. a component of a complex, you should make a new object containing a copy of just that component and calculate its area.
please let us know if it worked

Thanks, this code is working, but I am still not able to calculate the area of selected residues. My main hurdle is how to use the get_area pymol command and incorporate in the loop.

@shome edited my answer

I typed :
seleobjs = cmd.get_object_list('(sele)')
for obj in seleobjs: get_area obj
It is showing me invalid syntax..

cmd.get_area(str selection="(all)", int state=1, int load_b=0)     try with Parentheses (Round Berackets)

showing invalid syntax

both get_area(obj) and cmd.get_area(obj) ?

Ok went back to this question having pymol within a hand's reach this time. Question is not really clear to me so I am posting a reasonable (to me) version of the question into a new code snippet;
I am using two PDBs containing ATP as ligand '3dgl.pdb', '2p09.pdb'
code:
import pymol

list_pdbs = ['3dgl.pdb', '2p09.pdb']

pymol.finish_launching()

for _i in list_pdbs :

    pymol.cmd.load( _i, _i.rstrip('.pdb'))

objs = pymol.cmd.get_object_list(selection='(all)')

print('structure loaded : ' , objs)

pymol.cmd.remove('resn HOH')

for _i in objs:

    pymol.cmd.create(_i+'_ATP_copy' ,  _i + ' and resn ATP')

    # print('after create ')

    area = pymol.cmd.get_area(_i+'_ATP_copy', quiet = 0)

    pymol.cmd.show('dots' , _i+'_ATP_copy')

    print(_i+'ATP_copy' + ' area --->  ' + str(area))

output:
structure loaded :  ['3dgl', '2p09']
 cmd.get_area: 1013.704 Angstroms^2.
3dglATP_copy area --->  1013.7040405273438
 cmd.get_area: 1013.704 Angstroms^2.
2p09ATP_copy area --->  1013.7040405273438

picture of the ATP surfaces in Pymol itself,

you can play with  cmd.get_area(...) parameters see https://pymolwiki.org/index.php/Get_area


text:
How to analyze expression of certain group of genes in certain types of cell?
I need to analyze expression of ~ 400 genes with a certain function in the embryonic stem cells. All these genes are combined into one database with RefSeq IDs. What is a recommended workflow to quantify differential expression of these genes in RNA-Seq data of cells and visualize it?
All genes for analysis concern Human E3 Ubiquitin Ligases originally located here, hpcwebapps.cit.nih.gov/ESBL/Database/E3-ligases and the updated database is here:
https://esbl.nhlbi.nih.gov/Databases/KSBP2/Targets/Lists/E3-ligases/

By database do you mean a table? Your question is a bit broad.. have you tried anything so far?

I really apologize for it being so broad because I am a very new to bioinformatics and a little bit lost. I mean all genes I need to analyze are in a database https://hpcwebapps.cit.nih.gov/ESBL/Database/E3-ligases/

Hey no problem. Ok if you have 400+ genes, and expression values, you need to have them onto a count table, where rows are genes, and columns are your samples. If you can do some basic R, then try to run DESeq (https://www.nature.com/articles/nprot.2013.099) or edgeR (https://bioinformatics-core-shared-training.github.io/cruk-bioinf-sschool/Day3/rnaSeq_DE.pdf)

You might have a look at papers such as this one: https://www.embopress.org/doi/10.15252/msb.20188746


text:
Scooping and GeneMatcher
A colleague has uploaded about a dozen candidate disease genes to GeneMatcher, described here and here and located here,  with excellent success. However last year and again recently he has encountered the case wherein a case was uploaded and a cohort formed, but a few months down the line one of the dropped out collaborators published an n=1 paper in a low tier journal, basically using GeneMatcher to check whether the disease link of the gene was legitimate, to assess the phenotypic range and see what test are worth doing.

Is there documented precedence of this behaviour (preferably citable)? If not,
does anyone have personal experience of this and what is the ratio of incidences over total submissions (I know that is hearsay, but it would be very helpful)?
And also has anyone successfully taken action (say getting user barred)?

when you say "this behavior", do you mean an ex-collaborator trying to rush a part of the collaborative project into print to get credit? that is very common. Are you referring more narrowly to this behavior on the GeneMatcher platform?

I was not aware that is common in general for a collaboration to be formed only for the collaborator to drop out and rush to press with the intel. Turns out there are lots of collaboration agreements to protect against shenanigans in joint works which are similar to MTAs. Sickening as it sounds this is the solution.

I think that often collaborator starts with good intentions, but then the team falls apart and stop liking each other, and they still want to get something out of the work. There are a few ways it can happen. Never happened to me personally (at least not like this specifically). But anecdotally, I'd guess this is actually one of the most common paths to scooping.

Ah, no, this is more malicious. The site matches geneticists and they exchange patient details in chain emails, discuss a path to publication (agreed by all) and make the clinicians request further analyses say MRIs or blood work (within 2-3 months). Then when follow-up lab tests are planned, the collaborator suddenly goes silent and half a year later there is a paper in a low tier journal with their sole patient but with strong emphasis on the key traits discussed on in the collaborations.

Ah, I see. that is a more specific and very bad behavior.

I am having trouble finding anything that directly addresses the issue of failed collaborations leading to scooping, but here are some examples of similar behavior:

it is sort of implicit in the discussion here.

This question on academia SE is related but not exactly right.

There are some related stories here.

Here is an article about plagiarism among collaborators, which is quite similar.

I feel like a lot of the discourse around preprinting etc. is similar but again not quite there.

As for punishment, that sounds hard. Scooping is sort of built into science; fighting over precedence is a big deal and has been for a long time. My solution would be to get rid of current academic publication and have everyone's open digital lab notebooks streaming results to publication feeds in real time; that way there's no ambiguity. Preprints and micropublications are sort of a poor man's version of this.
This is not a very complete answer but there may be enough there to help.

Thanks, excellent points. Just to add some extra in terms of GeneMatcher. Most collaborations flowered and only a few sour. Someone else has had a similar experience to my colleague, so it is not uncommon. The solution is to use collaboration agreements to safeguard against this —even if these are laughably enforceable.


text:
A question on Homer normalization
When using annotatePeaks.pl script from Homer software to create histograms, the output is normalized per bp per peak (on top of normalization to 10 million tags).
What does it mean to normalize per bp per peak?

Hypergeometric Optimization of Motif EnRichment (HOMER v4.11, 10-24-2019), is used for de novo motif discovery and NGS analysis and written in Perl and C++.

Crossposted from here, as there is no reply on biostars.


text:
Find covariation of different DNA-binding protein binding between two conditions
I have ChIP-seq data of RNAPII and two other factors which we think follow RNAPII during transcription in two different conditions, and I'd like to show that the genes that lose RNAPII also lose the other two factors, and same for gains.
Do you know any methods that can calculate some sort of covariation coefficients?
I know there are methods that calculate covariation of genes between two conditions, but none for different measures on the same gene.

do you have an example dataset?

@StupidWolf Not yet, I'm still in preprocessing. I had a few issues, but I should have it soonish

Answer from @stupidwolf converted from comment:
I would try something like this, for every gene you have a logFC of RNAPII in control versus treatment. Do this for other 2 proteins, and you correlate the logFC. In very rough terms, correlation is the scaled form of covariance, so this will work.
If you have 1000 genes, do correlation between 1000 logFC from RNAPII against 1000 logFC from other proteins.
Answer from @m converted from comment:
You use a covariation stat and then perform Pearson's correlation to standardise the response between -1 to 1. An answer of >0.7 is considered robust and >0.8 highly robust. ITs an interesting question and you are right covariation is the way to answerr it


text:
Compare substitution matrices
I am looking for methods to compare substitution matrices to determine which one is better for a given set of data / type of protein.
Do you know such methods / tools?
e.g.  BLOSUM62 vs PHAT for 50k - 100k+ sequences

The technique is phylogenetics is called likelihood ratio and for a regular amino acid matrix is 2Δln which is assessed on Chi2 tables where Δln is difference in likelihood between model 1 and model 2. The degrees of freedom hmmmm... is the number of parameters that are constrained, this involves the number of taxa minus the number of parameter differences between the models. This makes performing 50 000 taxa ... tricky because the degrees of freedom will be outside the Chi2 table.
iqtree will perform a tree and provide a likelihood for 50+K taxa. There are two limitations,

No PHAT model
degrees of freedom (as mentioned)

Overcoming df is simply to say whatever is the lowest likelihood is the best model and simply state Δln without bothering with significance testing. If the difference between models is large you'll get away with that (needs to be big) and for 50+K taxa it will be need to be huge.
Getting PHAT (the model below the ---- might be PHAT) into iqtree is certainly doable because of QMaker. However, how this is done I don't know, because I never done this calculation with iqtree, you'll need to read the manual. For traditional packages such as PAML it's easy ... you make a tree then specific the custom amino acid model and it spits out the likelihood. For the number of taxa involved, PAML is not possible, it's impossibly big for PAML, even 1000 is a lot of PAML.
iqtree does have LG, WAG and JTT amino acid models and will almost certainly have all the BLOSUM matrices and you just specify those at the command line. Again a comparison with PHAT needs to specify a custom matrix and I don't know how to do that in iqtree. The only thing I know is it must be doable.
iqtree is the most efficient of the maximum likelihood phylogenetic methods by a very long way. Hence for dealing with 50+K is why I'm recommending it here.
Note The theory about whats going on and why the whole thing needs phylogenetics is very complicated. Basically its an entire lecture on phylogenetic theory.

This I believe is PHAT (needs checking)
https://github.com/biopython/biopython/pull/283
Note the code is copied and pasted from the above link, therefore it should be strictly regarded as a quote for copyright purposes.
#!/usr/bin/env python
matrix = {
('W', 'F') :  0.0, ('L', 'R') : -6.0, ('I', 'I') :  5.0, ('Q', 'Q') :  9.0,
('W', 'N') : -5.0, ('V', 'I') :  3.0, ('H', 'T') : -4.0, ('H', 'P') : -6.0,
('W', 'V') : -4.0, ('Q', 'E') :  1.0, ('W', 'R') : -7.0, ('Q', 'A') : -3.0,
('H', 'H') : 11.0, ('H', 'D') : -1.0, ('L', 'N') : -3.0, ('Y', 'M') : -2.0,
('Y', 'I') : -3.0, ('Y', 'E') : -2.0, ('E', 'S') : -3.0, ('Y', 'A') : -3.0,
('Y', 'Y') : 11.0, ('T', 'C') : -1.0, ('E', 'C') : -7.0, ('Y', 'Q') :  0.0,
('E', 'G') : -3.0, ('V', 'A') :  1.0, ('C', 'C') :  7.0, ('M', 'R') : -6.0,
('P', 'T') : -4.0, ('V', 'E') : -5.0, ('P', 'P') : 13.0, ('I', 'T') : -1.0,
('K', 'S') : -5.0, ('R', 'G') : -5.0, ('I', 'P') : -4.0, ('R', 'C') : -8.0,
('A', 'T') :  0.0, ('K', 'K') :  5.0, ('A', 'P') :  3.0, ('V', 'M') :  1.0,
('I', 'D') : -5.0, ('K', 'C') :-10.0, ('K', 'G') : -5.0, ('R', 'S') : -6.0,
('F', 'Q') : -2.0, ('F', 'A') : -1.0, ('V', 'V') :  4.0, ('M', 'N') : -2.0,
('F', 'E') : -5.0, ('D', 'N') :  2.0, ('F', 'I') :  0.0, ('F', 'M') :  0.0,
('M', 'S') : -2.0, ('S', 'S') :  6.0, ('L', 'Q') : -3.0, ('W', 'E') : -7.0,
('W', 'A') : -4.0, ('W', 'M') : -4.0, ('H', 'S') : -2.0, ('W', 'I') : -4.0,
('S', 'C') :  1.0, ('L', 'A') : -1.0, ('L', 'E') : -5.0, ('W', 'Q') :  1.0,
('H', 'G') : -4.0, ('Q', 'N') :  2.0, ('H', 'C') : -7.0, ('L', 'M') :  2.0,
('W', 'Y') :  1.0, ('Y', 'N') :  2.0, ('E', 'P') : -5.0, ('Y', 'F') :  4.0,
('E', 'T') : -5.0, ('A', 'A') :  5.0, ('I', 'N') : -3.0, ('G', 'A') :  1.0,
('Y', 'V') : -3.0, ('E', 'D') :  6.0, ('W', 'H') : -3.0, ('Y', 'R') : -6.0,
('M', 'Q') : -1.0, ('P', 'S') : -3.0, ('R', 'H') : -4.0, ('A', 'C') :  1.0,
('R', 'D') : -7.0, ('K', 'P') : -4.0, ('L', 'D') : -5.0, ('K', 'T') : -6.0,
('V', 'N') : -3.0, ('M', 'A') : -1.0, ('K', 'H') : -5.0, ('V', 'R') : -7.0,
('P', 'C') : -8.0, ('M', 'E') : -5.0, ('A', 'S') :  2.0, ('T', 'T') :  3.0,
('R', 'T') : -6.0, ('I', 'G') : -2.0, ('R', 'P') : -7.0, ('K', 'D') : -5.0,
('I', 'C') : -3.0, ('F', 'R') : -7.0, ('F', 'V') : -1.0, ('L', 'C') : -2.0,
('F', 'F') :  6.0, ('D', 'A') : -5.0, ('F', 'N') : -1.0, ('W', 'D') : -7.0,
('L', 'P') : -5.0, ('Q', 'S') : -1.0, ('N', 'C') : -2.0, ('N', 'G') : -1.0,
('H', 'N') :  4.0, ('W', 'T') : -7.0, ('Q', 'G') : -2.0, ('W', 'P') : -6.0,
('Q', 'C') : -5.0, ('N', 'S') :  1.0, ('L', 'H') : -4.0, ('L', 'L') :  4.0,
('G', 'T') : -1.0, ('M', 'M') :  6.0, ('G', 'P') : -3.0, ('Y', 'K') : -4.0,
('Y', 'G') : -3.0, ('Y', 'C') : -1.0, ('E', 'A') : -5.0, ('E', 'E') : 12.0,
('Y', 'S') : -2.0, ('M', 'P') : -5.0, ('V', 'C') : -2.0, ('M', 'T') :  0.0,
('V', 'G') : -2.0, ('R', 'E') : -6.0, ('V', 'K') : -8.0, ('K', 'Q') : -1.0,
('R', 'A') : -6.0, ('I', 'R') : -6.0, ('N', 'A') : -2.0, ('V', 'S') : -2.0,
('M', 'D') : -5.0, ('M', 'H') : -4.0, ('K', 'A') : -7.0, ('R', 'Q') : -2.0,
('K', 'E') : -4.0, ('F', 'S') : -2.0, ('I', 'K') : -7.0, ('D', 'P') : -5.0,
('D', 'T') : -5.0, ('I', 'M') :  3.0, ('F', 'C') :  0.0, ('W', 'L') : -3.0,
('F', 'G') : -2.0, ('F', 'K') : -7.0, ('F', 'T') : -2.0, ('D', 'D') : 12.0,
('Q', 'T') : -3.0, ('W', 'G') : -5.0, ('Q', 'P') : -3.0, ('W', 'C') : -4.0,
('W', 'K') : -8.0, ('H', 'Q') :  2.0, ('Q', 'D') :  0.0, ('W', 'W') : 11.0,
('V', 'L') :  1.0, ('L', 'G') : -2.0, ('W', 'S') : -5.0, ('L', 'K') : -7.0,
('N', 'P') : -4.0, ('H', 'E') : -1.0, ('N', 'T') : -1.0, ('H', 'A') : -3.0,
('Y', 'L') : -2.0, ('Y', 'H') :  3.0, ('G', 'S') :  1.0, ('Y', 'D') : -4.0,
('V', 'Q') : -3.0, ('L', 'T') : -1.0, ('G', 'G') :  9.0, ('G', 'C') : -2.0,
('E', 'N') :  0.0, ('Y', 'T') : -3.0, ('Y', 'P') : -5.0, ('R', 'N') : -3.0,
('V', 'D') : -5.0, ('K', 'R') : -1.0, ('V', 'H') : -5.0, ('I', 'Q') : -3.0,
('V', 'P') : -4.0, ('M', 'C') : -2.0, ('K', 'N') : -2.0, ('V', 'T') :  0.0,
('M', 'G') : -1.0, ('T', 'S') :  1.0, ('I', 'E') : -5.0, ('M', 'K') : -6.0,
('I', 'A') :  0.0, ('N', 'N') : 11.0, ('R', 'R') :  9.0, ('F', 'P') : -5.0,
('L', 'I') :  2.0, ('I', 'S') : -2.0, ('D', 'S') : -4.0, ('L', 'S') : -2.0,
('I', 'H') : -5.0, ('F', 'D') : -5.0, ('D', 'C') : -7.0, ('F', 'H') : -2.0,
('D', 'G') : -2.0, ('F', 'L') :  1.0
}

chars = [aa[0] for aa in matrix] + [aa[1] for aa in matrix]
chars = sorted(list(set(chars)))

output = " ,%s\n" % ",".join(chars)

for i, aa1 in enumerate(chars):
    output += aa1
    for j, aa2 in enumerate(chars):
        if (aa1, aa2) in matrix:
            output += "," + str(matrix[(aa1, aa2)])
        else:
            output += "," + str(matrix[(aa2, aa1)])
    output += "\n"
print(output)


text:
panic: division of zero by zero or infinity by infinity error when running GODON positive selection analysis
I am using godon to find genes under positive selection with the branch-site test. However, after running for a little while I get the error "panic: division of zero by zero or infinity by infinity".
Below is the full output with the error report. I can supply the phylogenetic tree and fasta file I am using if that would be helpful. Thank you.

iter=33 lnL=-73104.072
iter=66 lnL=-72956.855
iter=99 lnL=-72912.956
iter=131 lnL=-72902.900
iter=162 lnL=-72900.147
iter=194 lnL=-72898.562
iter=226 lnL=-72897.970
iter=258 lnL=-72897.651
iter=289 lnL=-72897.435
iter=321 lnL=-72897.284
iter=353 lnL=-72897.199
iter=385 lnL=-72897.153
iter=416 lnL=-72897.108
iter=448 lnL=-72897.076
iter=17 lnL=-72890.551
panic: division of zero by zero or infinity by infinity
goroutine 1 [running]:
math/big.(*Float).Quo(0xc4273a8f40, 0xc4273a8f68, 0xc4273a8f18, 0xc42e88c118)
/usr/local/go/src/math/big/float.go:1639 +0x15c
bitbucket.org/Davydov/godon/cmodel.(*BranchSite).BEBPosterior(0xc422164120, 0xc4259bb900, 0x91, 0x91)
/go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:416 +0x6e9
bitbucket.org/Davydov/godon/cmodel.(*BranchSite).Final(0xc422164120, 0x101)
/go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:455 +0x112
main.computeFinal(0x1f263a0, 0xc422164120, 0xc4201522a0, 0x101, 0x6a80, 0x0)
/go/src/bitbucket.org/Davydov/godon/godon/optimize.go:150 +0x93
main.performSingleTest(0xc420162320, 0xc428266000, 0x9399, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, ...)
/go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:335 +0x1b52
main.hypTest(0x1ee03a7, 0x1ee03a7, 0x4, 0xc4201aa501, 0x1, 0x1)
/go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:88 +0x71a
main.main()
/go/src/bitbucket.org/Davydov/godon/godon/godon.go:316 +0xefa

Hi, sorry about delay. Apparently I didn't subscribe properly to the tag.  Anyhow this looks like a bug to me, would be great if you could open an issue on bitbucket and attach an example file together with options you used to run software. https://bitbucket.org/Davydov/godon/issues

Ah wait you mentioned a branch sites test, right that makes sense. You've got two or more identical sequences in the analysis. Ensure there are no identical sequences in your alignments and ... PRESTO! It will work :-)
Below are my original thoughts before I re-read your question. I've kept the notes because it helps explain the theory a bit. There is no need to supply the tree, it should be resolved.

The tree likelihood looks to be decreasing nicely, so it can see both the tree and the codons in the alignment.
At a guess the error relates to two or more indentical nucleotide sequences in your alignment. Thus zero non-synonymous sites (amino acid change) and zero synonymous (silent) sites. This error would only occur if you were calculating selection for each 'node' in the tree (internal branch).  Alternatively, possibly if the number of synonymous sites was zero, regardless of the non-synonymous mutations this error could appear, but this doesn't seem to fit the error output.
I don't know "Godon", but recommend PAML for this analysis. Admittedly it is complicated to deploy, but it is a very robust algorithm and hugely respected.


text:
Comparison of phylogeny of SARS2 whole genome vs spike genome
I’ve built two phylogeny at minimum evolution in MEGA-X and extracted these as newick files.
My next step is to compare the two phylogeny to look at if the mutations within SARS2 are more prevalent within the spike, or as a global event across the whole genome. In order to do this I have chosen to use VisualTreeCMP. However, when it comes to choosing how to compare these I am very confused by all the possible options.
Could someone please give me some guidance on the most suitable parameters to give a good image, and comparison between the two phylogeny.
please see the newick format below:
((((((((((((((hCoV-19/Scotland/CVR140/2020|EPI_ISL_425684|2020-03-16:0.00026192,hCoV-19/England/CAMB-73211/2020|EPI_ISL_440287|2020-03-21:0.00026198):0.00000084,hCoV-19/Wuhan/HBCDC-HB-06/2020|EPI_ISL_412982|2020-02-07:0.00000000):0.00000034,hCoV-19/Gambia/NPHL-287/2020|EPI_ISL_561300|2020:0.00000000):0.00000372,hCoV-19/England/NORT-286167/2020|EPI_ISL_573313|2020-07-21:0.00025782):0.00002652,(hCoV-19/South_Africa/KRISP-0543/2020|EPI_ISL_498084|2020-07-12:0.00020599,(hCoV-19/England/ALDP-49EC3C/2020|EPI_ISL_559586|2020-05-20:0.00027741,hCoV-19/Belgium/HCM-0324151/2020|EPI_ISL_420408|2020-03-24:0.00024635):0.00005585):0.00021443):0.00002118,(hCoV-19/Colombia/MAG-INS-87968/2020|EPI_ISL_498152|2020-03-25:0.00015372,((hCoV-19/Australia/VIC8930/2020|EPI_ISL_565286|2020-07-29:0.00000000,hCoV-19/Australia/VIC8929/2020|EPI_ISL_565285|2020-07-28:0.00000000):0.00047159,(hCoV-19/Australia/VIC5819/2020|EPI_ISL_518613|2020-07-20:0.00024688,(hCoV-19/Australia/VIC4715/2020|EPI_ISL_519103|2020:0.00000000,hCoV-19/Australia/VIC2165/2020|EPI_ISL_480772|2020-03-19:0.00000000):0.00001495):0.00004719):0.00036022):0.00009570):0.00004876,((((hCoV-19/Wales/PHWC-16BD74/2020|EPI_ISL_573970|2020-09-15:0.00058675,hCoV-19/Spain/PV-IBV-001568/2020|EPI_ISL_467121|2020-03-13:0.00000000):0.00015901,hCoV-19/England/CAMB-7AAD6/2020|EPI_ISL_433837|2020-04-07:0.00009990):0.00003774,hCoV-19/Greece/HPI-6465/2020|EPI_ISL_451311|2020-03-29:0.00009055):0.00010675,((hCoV-19/USA/CA-SR0097/2020|EPI_ISL_445096|2020-03-13:0.00026193,hCoV-19/Iceland/13/2020|EPI_ISL_417765|2020-02-27:0.00026193):0.00020549,(hCoV-19/England/LIVE-A4CDD/2020|EPI_ISL_500083|2020-03-12:0.00007450,(hCoV-19/Wales/PHWC-16C441/2020|EPI_ISL_574056|2020-09-19:0.00022387,hCoV-19/England/SHEF-10B2F49/2020|EPI_ISL_537266|2020-05-21:0.00003797):0.00018736):0.00008143):0.00010981):0.00001973):0.00006712,(hCoV-19/Italy/LOM-UniMI-L134/2020|EPI_ISL_542130|2020-03-17:0.00064152,hCoV-19/Hong_Kong/HKU-200823-001/2020|EPI_ISL_516798|2020-03-26:0.00015299):0.00016120):0.00001351,(hCoV-19/USA/TX-HMH-1807/2020|EPI_ISL_542855|2020-04-20:0.00010829,hCoV-19/USA/CA-ALSR-1462/2020|EPI_ISL_494618|2020-03-20:0.00015349):0.00017136):0.00008773,(hCoV-19/USA/MA-MGH-00617/2020|EPI_ISL_460393|2020-03-28:0.00022036,(hCoV-19/USA/WA-S28/2020|EPI_ISL_417081|2020-03-02:0.00000726,hCoV-19/Singapore/746/2020|EPI_ISL_500560|2020-07-08:0.00025453):0.00022704):0.00003603):0.00001521,(hCoV-19/Spain/MD-IBV-006433/2020|EPI_ISL_510145|2020-04-03:0.00030013,hCoV-19/Austria/CeMM0009/2020|EPI_ISL_419662|2020-03-14:0.00000000):0.00025838):0.00001169,hCoV-19/Wuhan/IME-WH01/2020|EPI_ISL_529213|2019-12-30:0.00000000):0.00000016,hCoV-19/Wuhan/IPBCAMS-WH-01/2019|EPI_ISL_402123|2019-12-24:0.00000000):0.00000045,hCoV-19/USA/WA-S106/2020|EPI_ISL_417159|2020-02-29:0.00026265,NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome:0.00000000);

(((((((((((hCoV-19/England/ALDP-49EC3C/2020|EPI_ISL_559586|2020-05-20:0.00013635,hCoV-19/Belgium/HCM-0324151/2020|EPI_ISL_420408|2020-03-24:0.00006529):0.00000773,hCoV-19/South_Africa/KRISP-0543/2020|EPI_ISL_498084|2020-07-12:0.00026099):0.00013025,(hCoV-19/England/LIVE-A4CDD/2020|EPI_ISL_500083|2020-03-12:0.00016259,(hCoV-19/Gambia/NPHL-287/2020|EPI_ISL_561300|2020:0.00016807,hCoV-19/England/CAMB-73211/2020|EPI_ISL_440287|2020-03-21:0.00013706):0.00000534):0.00000008):0.00000278,(hCoV-19/England/NORT-286167/2020|EPI_ISL_573313|2020-07-21:0.00036517,((hCoV-19/Australia/VIC8930/2020|EPI_ISL_565286|2020-07-29:0.00000000,hCoV-19/Australia/VIC8929/2020|EPI_ISL_565285|2020-07-28:0.00000012):0.00009405,(hCoV-19/Australia/VIC5819/2020|EPI_ISL_518613|2020-07-20:0.00006516,(hCoV-19/Australia/VIC4715/2020|EPI_ISL_519103|2020:0.00003344,hCoV-19/Australia/VIC2165/2020|EPI_ISL_480772|2020-03-19:0.00003367):0.00000193):0.00000601):0.00017105):0.00002318):0.00000716,hCoV-19/Greece/HPI-6465/2020|EPI_ISL_451311|2020-03-29:0.00004599):0.00008629,(hCoV-19/England/SHEF-10B2F49/2020|EPI_ISL_537266|2020-05-21:0.00015654,(hCoV-19/Wales/PHWC-16C441/2020|EPI_ISL_574056|2020-09-19:0.00019758,hCoV-19/Wales/PHWC-16BD74/2020|EPI_ISL_573970|2020-09-15:0.00023814):0.00021545):0.00003087):0.00000571,((hCoV-19/USA/TX-HMH-1807/2020|EPI_ISL_542855|2020-04-20:0.00023928,hCoV-19/Colombia/MAG-INS-87968/2020|EPI_ISL_498152|2020-03-25:0.00013017):0.00003266,(hCoV-19/Scotland/CVR140/2020|EPI_ISL_425684|2020-03-16:0.00003644,hCoV-19/Italy/LOM-UniMI-L134/2020|EPI_ISL_542130|2020-03-17:0.00023391):0.00000153):0.00000121):0.00012657,(hCoV-19/Singapore/746/2020|EPI_ISL_500560|2020-07-08:0.00034258,(hCoV-19/Iceland/13/2020|EPI_ISL_417765|2020-02-27:0.00012787,(hCoV-19/England/CAMB-7AAD6/2020|EPI_ISL_433837|2020-04-07:0.00012752,(hCoV-19/USA/CA-ALSR-1462/2020|EPI_ISL_494618|2020-03-20:0.00010172,hCoV-19/Hong_Kong/HKU-200823-001/2020|EPI_ISL_516798|2020-03-26:0.00020057):0.00000680):0.00004213):0.00002475):0.00002767):0.00000390,(hCoV-19/Wuhan/HBCDC-HB-06/2020|EPI_ISL_412982|2020-02-07:0.00007999,(hCoV-19/Spain/PV-IBV-001568/2020|EPI_ISL_467121|2020-03-13:0.00019054,((hCoV-19/Wuhan/IME-WH01/2020|EPI_ISL_529213|2019-12-30:0.00003592,hCoV-19/Spain/MD-IBV-006433/2020|EPI_ISL_510145|2020-04-03:0.00020087):0.00000317,(hCoV-19/USA/CA-SR0097/2020|EPI_ISL_445096|2020-03-13:0.00011733,(hCoV-19/USA/MA-MGH-00617/2020|EPI_ISL_460393|2020-03-28:0.00013271,(hCoV-19/USA/WA-S28/2020|EPI_ISL_417081|2020-03-02:0.00003240,hCoV-19/USA/WA-S106/2020|EPI_ISL_417159|2020-02-29:0.00006806):0.00001270):0.00000336):0.00008627):0.00001038):0.00001150):0.00004587):0.00000335,hCoV-19/Austria/CeMM0009/2020|EPI_ISL_419662|2020-03-14:0.00013033):0.00000452,hCoV-19/Wuhan/IPBCAMS-WH-01/2019|EPI_ISL_402123|2019-12-24:0.00010072,NC_045512.2_Severe_acute_respiratory_syndrome_coronavirus_2_isolate_Wuhan-Hu-1_complete_genome:0.00000000);

And a link to the tree viewer:
https://eti.pg.edu.pl/treecmp/WEB

I am just requesting the basic parameters, on how one would get the most out of the programme... or failing that another programme with less in depth parameterization.

These are also rooted on the Wuhan reference genome

Interesting, will be at least until the weekend before I can respond

You need a variety of metrics in context, but like ... hmmmm, you're expecting a lot from a single tree comparison tool, albeit its parameter rich.

This is for an undergraduate dissertation, so I mean it doesn’t need to be crazy in depth analysis. Due to the pandemic, we’ve been forced away from wet lab work. So my main issue is having to self teach so much new analytic methods. So just the basic parameters, on how one would get the most out of the programme... or failing that another programme with less in depth parameters

Okay you've just made my "when I get around to doing it" that bit more "when I get around to doing it". Still its an interesting approach and it is a topical paper for sure.

I don’t know whether I should take that as positive or negative? :D

Hi @TheoJones I'll do this over next weekend, apologies for the delay.


text:
>My counter is counting genotypic combination occurences more than once, how do I ensure it counts one combination and doesnt count it again?
mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

        combinations = ['_'.join(sorted(x)) for x in combinations]
        combinations = [x.split('_') for x in list(set(combinations))]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

combination_labels = []
combination_counts = []

for lineage in lineages:

    df = pd.read_csv('Results/' + lineage).dropna()

    for combination in combinations:

            combination_df = df[list(df)]

            for mutation in combination:
                combination_df = combination_df[combination_df[mutation] == 1]

        #print(combination_df.shape[0])

            combination_labels.append('_'.join(combination))
            combination_counts.append(combination_df.shape[0])

            out_df = pd.DataFrame({'combination':combination_labels,
                               'count':combination_counts})
            out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
            out_df = out_df.sort_values('percentage', ascending = False)

            out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
                     header = True,
                     index = False)

The input CSV
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0

Output CSV
combination,count,percentage
D614G,87355,90.7084929856806

My above code is used to count all occurences of combinations of spike protein mutations

to determine the genotypes.
However, looking at the values processed into the .csv file, i can see ' D614G' is being counted even when
it is in combination with other combinations.

My question is how do i ensure once something has been counted (a row in the imported .csv file)
it is not counted for a second time?
OR possibly how can I edit this code to prevent counting of singular mutations even when presented within a combination?
Thanks in advance

Nice to see a traditional python question based on pandas. The top of your code is chopped off. Could you inset it please?

Could you also add what a row in your CSV looks like?

Input / output csv?

The input CSV in your Result folder is only two lines per file?

The input file has about 97,000 lines of data - each for a strain of SARSCOV2. The output file has around 54 lines, as genotypes but this should be quite a few less

Alright, so there are a number of problematic patterns in your code - as far as I understand what you are trying to do. Next time, try to post a reproducible example that people can use and more people will be willing to help.
L18-21:
combination_labels = []
combination_counts = []

for lineage in lineages:

Declaring these two lists before the loop, then appending all your combination labels inside the loop means that these two lists will contain all labels and counts for all runs (all your lineages). If I am understanding your intentions correctly this is probably not something you want. The simple fix here is to just move them into the for-loop.
L29-30:
for mutation in combination:
    combination_df = combination_df[combination_df[mutation] == 1]

Here you are iteratively checking if the value for each mutation is set to '1', but you are also overwriting your combination_df variable at each iteration of the loop. After this loop, combination_df will be whatever mutation last was set to '1' in your current combination. I will get back to solving this in my solution below.
L37-44:
combination_labels.append('_'.join(combination))
combination_counts.append(combination_df.shape[0])

out_df = pd.DataFrame({'combination':combination_labels,
                   'count':combination_counts})
out_df['percentage'] = (out_df['count'] / df.shape[0]) * 100
out_df = out_df.sort_values('percentage', ascending = False)

out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv',
         header = True,
         index = False)

Since this part is again scoped inside the for-loop, each iteration (each combination) will overwrite the same 'Results_2/' + lineage.replace(".csv", "") + '_3.csv' file. This you also want to move out of the for combination in ... loop.

A big concern is also performance. You have 2000-something combinations and 100,000 lines in your file. As you loop over combinations, you also process the DataFrame each time, when this is a proble you can solve in one iteration. It's not much data so it's still fine, but these habits are still good to develop early.
Now here's an idea for a rewrite that should work better. I don't have a csv from you to test with, so there might (proabably are) still things to fix, but hopefully it'll give you a starting point:
mutations = ['A222V', 'D614G', 'E484Q', 'E780Q', 'G476S', 'L18F', 'N439K',
             'S477', 'S477N', 'T478I', 'V483A']

combinations = []
for M in range(1, len(mutations)+1):
    for subset in itertools.combinations(mutations, M):
        combinations.append(subset)

combinations = ['_'.join(sorted(x)) for x in combinations]

root = "C:"

os.chdir(root)

lineages = os.listdir('Results')

for lineage in lineages:

    out_df = pd.DataFrame({"Count": 0}, index=combinations)

    df = pd.read_csv('Results/' + lineage, index_col="lineage").dropna()

    for index, row in df.iterrows():
        combination = '_'.join(sorted(df.columns[row==1]))
        out_df.loc[combination, "Count"] += 1

    out_df['Percentage'] = (out_df['Count'] / df.shape[0]) * 100
    out_df = out_df.sort_values('Percentage', ascending = False)

    out_df.to_csv('Results_2/' + lineage.replace(".csv", "") + '_3.csv', header = True,
        index = True)

Absolutely amazingly detailed reply, thank you so much! It's worked exactly how I needed it! Not sure how to mark this as answered, as I'm very new to this forum - perhaps a mod could help here?
Thank you again, sincerely!

This might help you with that


text:
Adding multiple dataframes to compare against genotypes - DATE/LOCATION - geotemporal incidences of SARS2 genotypes
I'm trying to reuse some code I have previously used to add extra parameters.
In the initial code, it was used to give a 1/0 counter if a specific mutation was present (first code below)
import re
import os
import pandas as pd
import logging
import sys
import csv

root = "C:"

os.chdir(root)

lineages = os.listdir('Lineages')
mutations = os.listdir('Mutations')

for lineage in lineages:

    df = pd.read_csv('Lineages/' + lineage).dropna()
    df.columns = ["lineage"]
    df['lineage'] = [x.split('|')[0].split('_')[-1]
                     for x in df['lineage']]

    for mutation in mutations:

        mutation_df = pd.read_csv('Mutations/' + mutation).dropna()
        mutation_header = mutation.replace(".txt", "")
        mutation_df.columns = [mutation_header]
        mutation_df[mutation_header] = [x.split('_')[-1].split('|')[0]
                          for x in mutation_df[mutation_header]]

        mutation_set = set(mutation_df[mutation_header])

        df[mutation_header] = [1 if x in mutation_set else 0
                               for x in df['lineage']]

    df.to_csv('Results9/' + lineage.replace(".txt", "") + '.csv',
                 header = True,
                 index = False)

However, I now wish to add extra parameters such as date of sample and location of sample.
The first code produces this type of counter system:
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0
417942,0,1,0,0,0,0,0,0,0,0,0
417944,0,1,0,0,0,0,0,0,0,0,0
417946,0,1,0,0,0,0,0,0,0,0,0
417947,0,1,0,0,0,0,0,0,0,0,0
417948,0,1,0,0,0,0,0,0,0,0,0
417950,0,0,0,0,0,0,0,0,0,0,0
417955,0,1,0,0,0,0,0,0,0,0,0
431011,0,1,0,0,0,0,0,0,0,0,0

This code is fine and working, but if I try and add the second parameter (date/location) with the following code:
import re
import os
import pandas as pd
import logging
import sys
import csv

root = "C:/Users/jones/Downloads/SARSCOV2/"

os.chdir(root)

lineages = os.listdir('Lineages')
mutations = os.listdir('Mutations')
dates = os.listdir('Dates')
locations = os.listdir('Locations')

for date in dates:

    date_df = pd.read_csv('Dates/' + date).dropna()
    date_df.columns = ["Date"]
    date_df['Date'] = [x.split('|')[-1].split('/')[-1]
                       for x in date_df['Date']]

    for location in locations:

        location_df = pd.read_csv('Locations/' + location).dropna()
        location_df.columns = ["Location"]
        location_df['Location'] = [x.split('/')[1].split('_')[0]
                                   for x in location_df['Location']]

        for lineage in lineages:

            df = pd.read_csv('Lineages/' + lineage).dropna()
            df.columns = ["lineage"]
            df['lineage'] = [x.split('_')[-1].split('|')[0]
                             for x in df['lineage']]

            for mutation in mutations:

                mutation_df = pd.read_csv('Mutations/' + mutation).dropna()
                mutation_header = mutation.replace(".txt", "")
                mutation_df.columns = [mutation_header]
                mutation_df[mutation_header] = [x.split('_')[-1].split('|')[0]
                                                for x in mutation_df[mutation_header]]

                mutation_set = set(mutation_df[mutation_header])

                df[mutation_header] = [1 if x in mutation_set else 0
                                       for x in df['lineage']]

            df.to_csv('Results_3/' + lineage.replace(".txt", "") + '.csv',
                 header = True,
                 index = False)

I retrieve this data, where the counting is maintaibned, but the new parameters (DATE/LOCATION) are not added
lineage,A222V,D614G,E484Q,E780Q,G476S,L18F,N439K,S477,S477N,T478I,V483A
417941,0,1,0,0,0,0,0,0,0,0,0
417942,0,1,0,0,0,0,0,0,0,0,0
417944,0,1,0,0,0,0,0,0,0,0,0
417946,0,1,0,0,0,0,0,0,0,0,0
417947,0,1,0,0,0,0,0,0,0,0,0
417948,0,1,0,0,0,0,0,0,0,0,0
417950,0,0,0,0,0,0,0,0,0,0,0
417955,0,1,0,0,0,0,0,0,0,0,0
431011,0,1,0,0,0,0,0,0,0,0,0

In short, my end goal is to:

Add Lineage ID|Location of Sample|DD/MM/YYY|Mutations| -- as the headers
still maintain the mutation counter
Use the mutation counter as a guide to determine where and when the genotypes (mutation combinations) arose.

The files where I am retrieving the parameters looks like this also:

hCoV-19/Singapore/4/2020|EPI_ISL_410535|2020-02-03

hCoV-19/USA/WA13-UW9/2020|EPI_ISL_413601|2020-03-02

hCoV-19/USA/WA-UW142/2020|EPI_ISL_416680|2020-03-11

hCoV-19/USA/WA-UW143/2020|EPI_ISL_416681|2020-03-11

hCoV-19/USA/WA-UW144/2020|EPI_ISL_416682|2020-03-09

hCoV-19/USA/WA-UW145/2020|EPI_ISL_416683|2020-03-15

hCoV-19/USA/WA-UW146/2020|EPI_ISL_416684|2020-03-14

hCoV-19/USA/WA-UW147/2020|EPI_ISL_416685|2020-03-15

hCoV-19/USA/WA-UW148/2020|EPI_ISL_416686|2020-03-14

A very laborious, and poorly explained question; I apologise in advance!

I have realised an error in my coding - however, this only allows the initiation of the mutation counter. I am still unable to add the extra dataframes. Please see the editted code.

Would you be able to share two or three of your input files?

https://docs.google.com/document/d/1fV6wiYKkTngCFOeOwcfl1JGbhpY_VpWrSpGsdTfy-bw/edit?usp=sharing
Lineage A
https://docs.google.com/document/d/1qCNtzRI6EV2GNn3BPuP6EtuB5VWkFqgmGuLT4z_XqoY/edit?usp=sharing
A222V Mutation

Sorry @BastianSchiffthaler I realise those two files may seem lacking. I was in a rush for work while sending these. So the ‘Dates’ and ‘Location’ input files are identical to the ‘Lineage’ files. I merely copied these, as the code didn’t like trying to make multiple separations from the same folder for some reason

Just to clarify, when I asked you to share some of your files, I meant to share them in a way that I can reproduce your issue from start to finish. You read from at least four different files in a quadruple for-loop and I'm pretty sure your goal could be achieved much more straightforwardly than your current implementation.
Without a reproducible example, I can offer a best-guess on what's happening: you create date_df and location_df, but after you instantiate them, you never do anything with the objects. Your final write happens from df, so you'll need to add those columns, e.g.:
df["Date"] = date_df["Date"][0]
df["Location"] = location_df["Location"][0]

The above is just an example, you'll need to make sure the correct value gets assigned to your df object instead of simply the first one.

Again, apologies I was late running for work and rather skimmed the post. So would this be implemented before the mutation_set loop or between it? Or before the df[mutation_header]? I assume the latter

I assumed the latter and it’s worked perfectly? Thanks again!