text:
How can I classify the 3 clades(S, G, V) of the coronavirus without using protein data?
On GISAID they classified the coronavirus using 4 clades(S, G, V, Other).
I downloaded around 1,000 complete genomes of the coronavirus from GISAID and I would like to classify each one as belonging to one of the 4 clades(S, G, V, Other).
On the top left of the image below is the classification of the 3 clades(S, G, V).
More specifically, the three strains were classified based on three specific variants(This answer was taken from my previous question):

Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the  gene's protein product   to a valine (V).

My problem is that using this method of classification for the coronavirus genome I am seeking information on the protein products of the genes in the genome.
So to be more specific here are my concise questions:

How can I classify a genome as "Strain S" using only the complete genome assembly, without any information about its proteins?
How can I classify a genome as "Strain G" using only the complete genome assembly, without any information about its proteins?
How can I classify a genome as "Strain V" using only the complete genome assembly, without any information about its proteins?

This picture was taken from this PDF on GISAID

Can you please link to where you find this information? I can't see it on the GISAID site and it's hard to answer if you don't tell us where you're finding this from.

The interpretation is easy at @terdon pointed out ... GISAID from my understanding is GIS mapping collaboration. A link to the phylogeny side of their work would be useful

@Michael based on the site, it looks like they get it from https://nextrstrain.org (which is an excellent resource, by the way). I'd still like to know where the OP found that image since I haven't been able to reproduce it on their site and there may be more details given there.

I didn't find any article about this topic. Only the pdf from which the picture was taken. Here it is: https://sars-cov-2.netlify.com/images/GISAID_hCoV-19_Analysis_Update_2020-04-04_1500UTC.pdf

Is there a way to classify the genome as one of the strains without using protein data? I only have the complete genomes without any protein data available.

It seems to be explained right there in the image you posted:

So, the three strains were classified based on three specific variants:

Strain S, variant ORF8-L84S: a variant in the gene "ORF8" which changes the leucine (L) residue at position 84 of the gene's protein product to a serine (S).
Strain G, variant S-D614G: a variant in the gene "S" which changes the aspartic acid (D) residue at position 614 of the gene's protein product to a glycine (G).
Strain V, variant NS3-G251V: a variant in the gene "NS3" which changes the glycine (G) residue at position 251 of the  gene's protein product   to a valine (V).

Are these the locations on the whole aligned genome or are they the location on the gene? For example is location 84 the 84th location on the whole genome or the 84th location on the "ORF8" gene in the whole genome? If the latter is the answer then where can I find the genes("ORF8", "S", "NS3") on the whole aligned genome? Please revisit my question. I reedited it.

@yuval as I say in the answer, it is the position on the protein product of the gene. So not the genome and not the gene, but the protein that is created when that gene is translated. If you now want to ask about finding the genomic positions of genes, please open a new question since that's a completely different subject.

agreed with @terdon its a standard convention, they describing the gene first, ancestral amino acid, position on the gene and the new mutation

@Michael no, they are not describing the gene or the position on the gene. That's p-notation, so everything is with respect to the protein, nothing about the gene. ORF8-L84S means the 84th amino acid of the protein, not the 84th residue of the gene that codes for it.

@terdon .. ok its typo, I erred

@terdon Is it possible that by this definition that a certain genome could be classified by multiple strains? For example Strain S and Strain G

Is there a way to classify the genome as one of the strains without using protein data? I only have the complete genomes without any protein data available.

@yuval you can look at the mutations at the DNA level and see whether they would cause one of the protein mutations mentioned.

The most important thing to say is we are grateful this information has been released. It is important to note these are very large genomes of around 9000 amino acids in total.
The clades are assigned to single amino acid mutations. The tree is built using nucleotide data therefore synonymous (silent) mutations are resulting in phylogenetic resolution and the authors have superimposed amino acid mutations over this tree structure. They have represented the tree as a circular phylogram.
So the reason why there is loads of branches (of a phylogenetic tree) of different sizes is because there are different numbers of silent (mostly 3rd codon wobble) mutations between the isolates. RNA viruses mutate rapidly. There are only 4+ amino acid changes across the entire epidemic according to this analysis.
Nomenclature
The clades are named after the amino acid mutation that arises within that group. Its a bit weird, but its okay. Hence S clade means a serine mutation, G clade a glycine mutation etc.... The thing is that S also equals Spike, which is a very important protein (below), so it is confusing because the S clade refers to ORF8. Which gene the amino acid mutation is occurring in is more important than the specific mutation. Secondly the "other" clade is dying out - which is not that informative, but others are holding their course.
Mutation scheme and protein
So @terdon is right

L->S at position 84 of ORF8 - Hmmmm I am not hugely excited by that because it is a very weird gene and is is not a radical mutation (the amino acids are related) which is absent in some closely related coronaviruses altogether. In other words some closesly related viruses can remove this gene altogether and remain funcational.
D->G (aspartic acid is D) in the Spike protein at position 614. That is an interesting mutation because S protein contains the cellular receptor and that is an unusual mutation. S protein is critical for vaccination.
G->V in the NS3 ... I dunno because glycine to valine means it will be a bit more hydrophobic. NS3 is in the structural protein region of the genome (including spike, capid, membrane and another structural protein) but I'm not sure precisely what is does

Epidemiological interpretation
Epidemiological interpretation is difficult because it doesn't really explain the possible differences in mortality rates between different countries, for example USA appears to have a low mortality rate in comparison to Italy but there is lots of mixing. It is also not easy to read the tree because it is a circular phylogram and I am struggling to read each "clade" (which has a precise phylogenetic diffinition) I have to take their word for it. Secondly is the colour scheme the 'red' swamps every other colour and gives a disproportionate red colour. This is an old epidemiologists trick in heat maps, there precise colour scheme albeit that is relative, appears to convey a "danger" message to the reader, even though mathematically its simply a consequence of the colour scheme.
Fundamentally it is not really clear whether the mutations confer phenotypic changes and these are huge genomes for an RNA virus. However, the G clade needs further investigation.
Summary
The battle continues as to whether data collection skews, social skews, age-prevelence skews, climatic skews result in variation of mortality rates between countries? Versus is their a genetic difference in mortality. On this data I'd go for the non-genetic factors.

The answer to your second question is simple, you just translate the triplet codon into amino acids and identify which amino acid is present at the above 3 positions. You then simply denote the clade within the taxon ID of your nucleotide data. MEGA X has an alignment editor that will allow this translation.

Is there a way to classify the genome as one of the strains without using protein data? I only have the complete genomes without any protein data available.

@yval I've answered your question above.


text:
Relationship between the instantaneous rate of change of a continuous time Markov model and the rate of a global clock model?
If continuous-time Markov models contain a parameter q that denotes the instantaneous rate of change and a global/strict clock model is the overall rate of evolution, how are these two parameters related? That is, for a given substitution model (e.g., JC69, F81, GTR), is there a deterministic relationship between the rate parameter(s) of these substitution models and the rate of the clock model? For instance, should the individual instantaneous rates of change (e.g., A -> C, T -> C, etc.) sum up to the global clock rate? If this is not the case, how should one understand the relationship between the rate of a continuous-time Markov model and the overall evolutionary rate?

Thats quite easy, but it'll be a few days before I can answer I'm a bit pushed for time. Could specify your MCMC method, is it Beast?

I've been using RevBayes more than Beast, but I'm familiar with the latter. Many thanks for your time.

Is that a reversible mutation matrix?

For some yes, for others no. The most basic model was an Mk model with k=2 for binary characters 0 and 1. I have two parameters, one for q_01 and the other for q_10.

Essentially you are comparing a "local clock" with a global clock and whether the instance rate relates to the global it depends on the model and in MCMC clocks in most cases the answer is no. A global clock is a very specific and very limiting model (below). Everything below is a comparison between taxa not between alignment positions.
"Local clock" in the context of MCMC clocks is in inverted commas because they are using population models. I will use the nomenclature here for convenience to describe rate heterogeneity between branches, not within a branch (last example).
A "local clock" has mutation rate heterogeniety across the tree whilst a global clock is completely homogenious. You can't say how is a "local clock" related to global clock because there are numerous different models.
Global clocks never hold for biological data because there are always some taxa evolving at a faster rate than others. This is particularly true in infectious disease.
One historic model used in traditional local clocks (no inverted commas) was autocorrelation, where a local clock would relate to a global clock on a "clade by clade" basis. That is clades are correlated by their mutation rate, so some clades evolve faster and other clades slower. Statastically an autocorrelated clock always outperforms a global clock under a likelihood ratio test.
I remember autocorrelation clearly because I implemented a variant of this approach, so I needed to understand how they did it, rather than click a few buttons on a GUI, which is what happens in Beast.
Beast has at least 3 population models and off-hand I cannot remember the population scenarios. I do remember the population scenarios made sense and secondly they can give 10-fold differences in rate in the data set I looked at, You need to read on the background theory for each model and assess what is more biologically appropriate for a given context. One model (from dim memories) reminded me of a classic population growth equation, i.e. an exponential. Population growth and collapse in phylodynamics would be informative of the rate, but again you need to go to the underlying theory.
The molecular clock theory that the rate of evolution towards the tips of trees is much faster than at the deeper branches is the sort of model that could be implemented in a MCMC clock. In the latter case the relationship between a global clock and the rate heterogeneity clock (which is a better description in that model than "local clock") is quite different from global vs. "local clock".
The fundamental issue with clocks is exemplified in the 10-fold difference issue. The lesson that is important to me is very easy to get a "wrong answer". The lesson important to you is if two MCMC models can't agree the same rate how can both obtain the same answer as a global clock? I am personally a fan of maximum heterogeneity, but in this day and age they will ask you to perform a Beast model to support it.

source: Stack Exchange

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

sort: votes


text:
Subsetting SNPs with specific exclusion criteria using bcftools
I am trying to subset SNPs from 32 cultivars. I have the following exclusion criteria:

>20% missing data
SNPs present in less than 50% of genotypes
Each genotype must have at least 50% of included SNPs
>=60% heterozygous calls
Bi-allelic SNPs only
No very rare variants (i.e. MAF<0.05)
bcftools view --output-type u --min-alleles 2 --max-alleles 2 \
  --types snps --exclude MAF[0]<0.05 --exclude ** \
    --known variants_raw_sorted.bcf > SNP_filtered_sorted_kkf.bcf

Where the ** is, I am trying to implement the remainder my exclusion criteria, (criteria 1-5), however I am having a great deal of difficulty understanding how to use the expression in bcftools. I am looking to understand how to use the information in the info fields (AN,AC, AF/MAF, NS) to implement my exclusion criteria. Any help is very much appreciated.

Further information : This is the upstream code that led to the bcf file.
fastp -i ACDC_R1.fastq.gz -I ACDC_R2.fastq.gz \
  -o fastp_ACDC_R1P.fastq.gz fastp_ACDC_R2P.fastq.gz \
  -R "ACDC_fastp_QC" -h "Fastp_ACDC" -j "ACDC_fastp" \
  -e 20 -c -a "auto" -L -r \
  -m --merged_out "ACDC_fastp_Merged" --out1 "ACDC_R1_fastp_Unmerged" \
  --out2 "ACDC_R2_fastp_Unmerged" \
  --unpaired1 "ACDC_R1Pass_R2Fail" --unpaired2 "ACDC_R2Pass_R1Fail"

bowtie2 -p 20 -x cs10 -U ACDC_fastp_Merged.fastq.gz -S ACDC_aligned.sam

samtools view -bS ACDC.sam > ACDC.bam

bcftools mpileup --output-type u -f cs10.fasta -b bam_list.txt |
    bcftools call -vmO z -o all_raw_condition_1.vcf.gz --threads=20

bcftools view --output-type u --output-file variants_raw.bcf variants_raw.vcf.gz

bcftools sort --max-mem 14000 --output-type u --output-file variants_raw_sorted.bcf

Your criteria 2, 3 and 5 don't seem to make sense to me - could you be a bit more clear about the terminology, and what you mean by SNP in this case? Generally, 'SNP' is another word used for genetic marker, or position, so it doesn't make sense for it to be present in 50% of genotypes.

Here I mean SNP in the usual sense - i.e. a mutation of a single nucleotide at a set position. Essentially, the idea is to use SNPs to discriminate between different cultivars. E.g. with a perfect set of 5 SNPs, you could discriminate between 2^5 cultivars.
To clarify the criteria, I will remove criteria 5 - misleading. For the others : Individually, we would want each SNP (i.e. mutation) in our set to be present in less than 50% of genotypes (criteria 2). Across our whole set of SNPs, we would want at least 50% of the included SNPs to be present in each genotype (crit 3).

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


text:
How to identify the reference genome assembly used to create a .hic file?
I have a bunch of .hic files of different species. But, I don't know which reference genome assembly (e.g. GRCh37 or GRCh38 in the case of human genome) the genomic co-ordinates of the .hic file correspond to. So I wonder if there is any way to figure out the version of the reference genome assembly from a .hic file.
Therefore, my question is: Does .hic file contain metadata such as the version of reference genome assembly? If yes, then how could I access it?
(I could obtain the information I need by digging into the methods sections of each paper corresponding to each .hic file. But it would be much better if I could get that right-away from the .hic file.).

Turns out HiC file does store such metadata.
straw tool from Aiden Lab contains a script (read_hic_header.py) that prints out the metadata from a given HiC file. Usage: read_hic_header.py <hic file or URL> [verbose]
In order to obtain the metadata in the form of a dictionary, one may consider using a function of the same tool that is currently available in my fork.
from straw import straw_module
straw_module.read_metadata(infile,verbose=False)

Args
    infile: str, path to the HiC file
    verbose: bool, print the metadata if True else return the metadata
Returns
    metadata: dict, the metadata is returned if verbose is False.
                Keys of the metadata:
                HiC version,
                Master index,
                Genome ID (str),
                Attribute dictionary (dict),
                Chromosomes (dict),
                chromosomes2length (dict),
                Base pair-delimited resolutions (list),
                Fragment-delimited resolutions (list).


text:
Differential Expression
I have genes differentially expressed between two groups (case and control). I would like to annotate them by classifying them according to their biological functions.
I work on parasites of Plasmodium Falciparum type Dd2 which are divided into two groups (case and control). I carried out a Differential expression analysis with cuffdiff and CummeRbund then I annotated my genes on plasmodb.org.
Example of id genes: pfDd2_14005040
Now I would like to do a pathway analysis to see which biological pathways these genes are involved in.
StupidWolf,
I resumed my annotation with the 3D7 strain of plasmoduim Falciparum then I created my geneList as you indicated but when I make the dotplot nothing is displayed
setwd("C:/cygwin64/home/DIANGO/EXCELL")
library(org.Pf.plasmo.db)
library(clusterProfiler)

gene <- c ("PF3D7_0936800","PF3D7_1478900","PF3D7_1009700","PF3D7_0508500")
geneList = keys(org.Pf.plasmo.db)
enr <- enrichGO(gene,universe=geneList,OrgDb='org.Pf.plasmo.db',keyType ="SYMBOL",ont="BP")
head(enr@result,readtable=TRUE)
dotplot(enr)

The output gives a white rectangle
                   ID                           Description GeneRatio  BgRatio     pvalue
GO:0006323 GO:0006323                         DNA packaging       1/2  21/2615 0.01599974
GO:0071103 GO:0071103               DNA conformation change       1/2  30/2615 0.02281728
GO:0051276 GO:0051276               chromosome organization       1/2  58/2615 0.04387582
GO:0006996 GO:0006996                organelle organization       1/2 144/2615 0.10712138
GO:0016043 GO:0016043       cellular component organization       1/2 203/2615 0.14925925
GO:0006464 GO:0006464 cellular protein modification process       1/2 265/2615 0.19244223
            p.adjust     qvalue        geneID Count
GO:0006323 0.1026777 0.07205456 PF3D7_0508500     1
GO:0071103 0.1026777 0.07205456 PF3D7_0508500     1
GO:0051276 0.1316275 0.09237015 PF3D7_0508500     1
GO:0006996 0.2203602 0.15463872 PF3D7_0508500     1
GO:0016043 0.2203602 0.15463872 PF3D7_0508500     1
GO:0006464 0.2203602 0.15463872 PF3D7_1009700     1
> dotplot(enr)

Once again I need your help.

EnrichR, gprofiler2, clusterprofiler, to give some keywords you can look up.

I used webtools like:
DAVID, Reactome, panther.
But they all tell me that my id genes are not in the right format. I made my annotations on plasmodb.org

I cannot help you with no information provided. How do your genes look, so first of all which species. Is it Ensembl IDs, HGNC? Please add some background.

i think david,etc might not have plasmodium annotations, so they cannot map your gene id to pathways... you use this.. http://bioconductor.org/packages/release/data/annotation/html/org.Pf.plasmo.db.html

I guess plasmodium doesn't have that many online tools dedicated to it so you can use the annotation packages from bioconductor and try it. I used clusterProfiler for analysis and org.Pf.plasmo.db for annotation below.
From this paper, grabbed table 1, a table of supposed genes involved in something:
genes = c("PF3D7_1302100", "PF3D7_0501300", "PF3D7_0302500*", "PF3D7_1030200",
"PF3D7_0302200*", "PF3D7_1327300", "PF3D7_0410000", "PF3D7_0104300",
"PF3D7_0420300", "PF3D7_0810000*", "PF3D7_0522300", "PF3D7_0935700*",
"PF3D7_0424700*", "PF3D7_1477700*", "PF3D7_0202000*", "PF3D7_1371600",
"PF3D7_1036300", "PF3D7_0422900", "PF3D7_0501400", "PF3D7_0221700",
"PF3D7_0301000", "PF3D7_1454100", "PF3D7_0831800", "PF3D7_0424900",
"PF3D7_1447500", "PF3D7_1477400", "PF3D7_0424800", "PF3D7_0823700",
"PF3D7_1477300", "PF3D7_0814100", "PF3D7_1372100", "PF3D7_0201500",
"PF3D7_0825300", "PF3D7_0201700", "PF3D7_1470000", "PF3D7_1401000",
"PF3D7_0724900", "PF3D7_0601700", "PF3D7_1410000", "PF3D7_0109300",
"PF3D7_1347300", "PF3D7_0201800", "PF3D7_1242100", "PF3D7_0403300",
"PF3D7_0114500", "PF3D7_1301200", "PF3D7_0931200", "PF3D7_0832200.1",
"PF3D7_0305900", "PF3D7_0608310", "PF3D7_0221800", "PF3D7_0511100",
"PF3D7_1001300", "PF3D7_1225900", "PF3D7_1444600", "PF3D7_1349900",
"PF3D7_1008300", "PF3D7_0919600", "PF3D7_1210300", "PF3D7_0929700",
"PF3D7_0701600", "PF3D7_0201600", "PF3D7_0801600", "PF3D7_1463600",
"PF3D7_0201900", "PF3D7_0905000", "PF3D7_1310000", "PF3D7_0315200",
"PF3D7_1139200", "PF3D7_1023500", "PF3D7_1111300", "PF3D7_0930700",
"PF3D7_1350800", "PF3D7_0412300", "PF3D7_0831900", "PF3D7_0406600",
"PF3D7_1448100", "PF3D7_0114300", "PF3D7_0526100", "PF3D7_1243200",
"PF3D7_0921100", "PF3D7_0301900", "PF3D7_0622700", "PF3D7_0811800",
"PF3D7_0114200", "PF3D7_0221500", "PF3D7_1129850", "PF3D7_1142200",
"PF3D7_0114100", "PF3D7_0512900", "PF3D7_1346200", "PF3D7_1369400",
"PF3D7_0314300", "PF3D7_0324100", "PF3D7_0107700", "PF3D7_1122000",
"PF3D7_0618800", "PF3D7_0601200", "PF3D7_1360300", "PF3D7_1318500",
"PF3D7_1111700", "PF3D7_0901900", "PF3D7_1313900", "PF3D7_1101700",
"PF3D7_1313200", "PF3D7_1039700", "PF3D7_0222100", "PF3D7_1371500",
"PF3D7_1344700", "PF3D7_1408300", "PF3D7_0929800", "PF3D7_0115000",
"PF3D7_1462900", "PF3D7_1352600", "PF3D7_0626200", "PF3D7_0831750",
"PF3D7_0602300", "PF3D7_0521100", "PF3D7_1452100", "PF3D7_0718900",
"PF3D7_0101300", "PF3D7_0815700", "PF3D7_0113600", "PF3D7_0937200",
"PF3D7_1108800", "PF3D7_0221100", "PF3D7_0921500", "PF3D7_1038300",
"PF3D7_0201400", "PF3D7_0822200", "PF3D7_0937000", "PF3D7_0301100"
)

We can for example using clusterProfiler:
library(org.Pf.plasmo.db)
library(clusterProfiler)

Besides the list of interest, you need to get a "universe" of all genes tested, in your situation it will be genes that you tested for differential gene expression, below I used all the genes annotated:
geneList = keys(org.Pf.plasmo.db)

enr <- enrichGO(genes,universe=geneList, OrgDb='org.Pf.plasmo.db',
keyType ="SYMBOL",ont="BP")

Above is a GO term enrichment, biological processes, you can look at the result like this:
head(enr@result,readtable=TRUE)
                       ID
    GO:0042000 GO:0042000
    GO:0044417 GO:0044417
    GO:0044766 GO:0044766
    GO:0051808 GO:0051808
    GO:0051836 GO:0051836
    GO:1902579 GO:1902579
                                                                                Description
    GO:0042000                                             translocation of peptides or proteins into host
    GO:0044417                                                        translocation of molecules into host
    GO:0044766                                                                    multi-organism transport
    GO:0051808 translocation of peptides or proteins into other organism involved in symbiotic interaction
    GO:0051836            translocation of molecules into other organism involved in symbiotic interaction
    GO:1902579                                                                 multi-organism localization

And you can plot it:
dotplot(enr)

You can read more about the analysis you can do at the package website

genes = c("PfDd2_090042100*", "PfDd2_060005300*","PfDd2_020024700*","PfDd2_050038200*", "PfDd2_130078600*","PfDd2_050037500*",
"PfDd2_060005700*","PfDd2_120045400*" )                                                                                --> No gene can be mapped....
--> Expected input gene ID: PF3D7_1405400,PF3D7_0623400,PF3D7_1025600,PF3D7_1439000,PF3D7_1423000,PF3D7_0716100
--> return NULL...

Apparently in this database there is no emotion that matches my gene_id (pdDd2_xxxxx) instead of PF3D7_xxxxxx. ps: I am working on the Dd2 strain of the parasite and not the 3D7 strain

Anyway you can map it to the PF3D7? For pathway etc analysis, you are as good as your annotation. I am not so familiar with plasmodium, but surely somehow must have encountered this problem?

Okey thank you very much for your suggestions I will resume the annotation with 3D7


text:
Why do we scale the rate matrix of a substitution model to make the average rate of substitution equal 1?
The matrix of transition probabilities for a substitution model Q over time t is found as follows: $e^{Qt}$. Since Q and t only show up as a product, we cannot distinguish time and rate of change. So we multiple Q by a constant so that the average rate of substitution is 1. Time t is then measured by the expected number of substitutions per site.
But why does forcing the average rate of substitution to equal 1 have this effect?

Anyway because this is assumed to be per alignment position without rescaling you cannot derrive the number of substitutions per site per unit time across an alignment.

@Michael's comment is correct but possibly a little difficult to understand, so I'm going to rephrase it here with an example (please let me know if I'm misrepresenting):
This matrix (pretending for a moment that we care only about purine/pyrimidine transitions and thus only two states) is directly interpretable in terms of the number of transitions per unit time:
$$ Q = \begin{matrix}
0 & 1.5 \\
0.5 & 0 \\
\end{matrix}
$$
Whereas this matrix is harder to think about because it is talking about the number of transitions per (I think?) $\frac{1}{0.37}$ units of time:
$$ Q = \begin{matrix}
0 & 0.37 \\
0.37 & 0 \\
\end{matrix}
$$
You will have to rescale any result anyways to talk about it in terms of "per unit time", so you might as well do it in the matrix itself so that at each point we are dealing with a more or less intuitively comprehensible numbers, given that the scaling doesn't actually change anything because time is dimensionless as you point out.
Of course, if you have a time-calibrated phylogeny then time is not dimensionless, so it might make sense to use a different scaling, such that the average rate in $Q$ is instead "per year" or "per million years" or whatever works within the allowances of floating-point precision and brains.
So it's not that rescaling has any particular effect on the math, but rather that rescaling makes it easier for human brains to think about, because we like to represent time in units.


text:
How to filter the PDB for single-domain proteins with full length 3D-structures solved
I am trying to reproduce a machine learning model that has been developed here. As one of the datasets, they use single-domain proteins with full length 3D-structures solved. Since I don't have a background in biology, I am not sure how to filter the PDB for these kind of proteins. Any help would be appreciated.

If you want to reproduce the model in that paper, you should email the corresponding author asking for the exact set of PDB IDs/chain IDs/residue ranges they used. There is too much flexibility in the criteria as given to get the exact set any other way.
If you want to assemble a similar set for your own work, the RCSB PDB advanced search is very powerful and would let you filter e.g. by max sequence length or for only monomers.
However the PDB doesn't deal with domains (regions that evolve and fold independently) directly. A PDB chain may contain part of a domain, one domain, or multiple domains.
If you want a dataset of protein domains you could use CATH, SCOP or ECOD.
As a side issue, please remember to remove overlap of homologous proteins between the training and test sets. People entering the field often neglect to do that. See for example this comment piece.

That is very helpful. Thank you for your answer.


text:
Set new Idents based on gene expression in Seurat and mix n match identities to compare using FindAllMarkers
I am relatively new to Bioinformatics and scRNA-seq data analysis. I am using Seurat V3 to analyze a scRNA-seq dataset in R. Currently, I have merged three scRNA-seq samples from the same donor into one Seurat object, All_Samples. We'll call them:

Uninfected
Virus1
Virus2
All_Samples <- merge(x = Uninfected, y = c(Virus1, Virus2))

Uninfected cells did not receive virus (i.e. express no viral genes), the other two samples were infected with a virus introducing new genes into the cells. I would like to accomplish two things.
First: I want to create a new Ident for Virus1 and Virus2 samples based on expression of a viral gene  of interest: "GeneA". All cells from Virus1 and Virus2 samples that have "GeneA" expression > 0.5 would be labeled "Pos", those with "GeneA" < 0.5 would be labeled "Neg" in a  the new Column called GeneA.
I am able to subset the objects based on GeneA expression that applies to all samples in the object. For example:
All_Samples_GeneA_Pos <- subset(All_Samples, subset = GeneA > 0.5)
All_Samples_GeneA_Neg <- subset(All_Samples, subset = GeneA < 0.5)

I assume I would have to modify All_Samples@meta.data based on this post and this post, but I admit I'm not entirely sure how to implement the suggested answers into my data.
Second: Once I have separated the data from Virus1 and Virus2 into "Pos" and "Neg" cells for "GeneA", I want to look for differentially expressed genes between all Uninfected cells and "Pos" cells from Virus1 or Virus2. I have been using something like this to compare samples:
 FindMarkers(object = All_Samples, group.by = 'virus', ident.1 = "1", ident.2 = "2")

But how would I do something equivalent to this:
 FindMarkers(object = All_Samples, ident.1 = "Uninfected", ident.2 = "Virus1_GeneA_Pos")

Or maybe
 FindMarkers(object = All_Samples, ident.1 = "Uninfected_GeneA_Neg", ident.2 = "Virus1_GeneA_Pos")

Let me know if I can clarify any points. I appreciate any help!

For the first question, you can use ifelse() to create a new column in the meta.data slot:
All_Samples@meta.data$new_column <-
  ifelse(
    rownames(All_Samples@meta.data) %in% colnames(All_Samples_GeneA_Pos),
    "GeneA_Pos",
    "GeneB_Pos"
  )

colnames(seurat_object) provides a vector of cell names in a given Seurat object. Here whatever cell that is in the All_Samples_GeneA_Pos object would be GeneA_Pos and whatever is not GeneB_Pos. To better control the behavior, you can use a "nested" ifelse(); you can have another ifelse() instead of the "GeneB_Pos" bit above.
For the second question, the rationale is the same as above, you can define an extra column with the "identities of interest", then use these identities with SetIdent() and then in turn using two identities of interest in the differential expression call.

Thank you for your response! It was very helpful.


text:
Extract autodocked protein-ligand connections programatically
I have 2 crystal structures, one for my protein and one for my ligand (I have several protein-ligand pairs). I am using AutoDock Vina to simulate docking, which returns another file of the ligand with updated 3D coordinates. Everything is in pdbqt format, but I can use OpenBabel to convert them if needed.
I am looking for a programmatic way (mainly python, but other languages work as well) to load both the protein and the docked ligand, then extract only the acting atoms which connect the 2 structures together (and the bonds between them).
If a direct answer is too complex, I would be happy if someone could suggest a tool/framework I could use to manipulate this sort of data easily so I can implement a solution on my own.

See a related post at Biostars for examples of how to do this with BioStructures.jl and LigPlot+: https://www.biostars.org/p/425276

Thanks. That helps.

With pymol as a python module (conda install -c schrodinger pymol, not the GUI one), it is very easy.
import pymol2
with pymol2.PyMOL() as pymol:
     pymol.cmd.load('protein.pdbqt', 'protein') # the second is the name within pymol instance.
     pymol.cmd.load('ligand.sdf', 'ligand')
     pymol.cmd.save('neighbours.pdb','ligand expand 3') # get atoms 3 Å around

Done!
However, if you are doing stuff in python in the first place you might want the last line as a string: pymol.cmd.get_pdbstr('ligand expand 3'). Changing the selection 'ligand expand 3' to 'byres ligand expand 3' will save the full residues.
Alternatively you can get the atom coordinates themselves etc. with it. To see what attributes Atom has you have you can do this:
for atom in pymol.cmd.get_model('ligand expand 3').atom:
    print(dir(atom))
    break

Which will give a large amount of attributes, such as atom.coord, which is a 3-element list.


text:
How to import data from cell ranger to R (Seurat)?
I will have some scRNA-seq data. The goal of the experiment will be to see if there is any difference in gene expression between treatment groups using the package Seurat from R.
I have read a tutorial how to do the analyze, but this tutorial does not explain how to import data.
I am a total beginner with R, I searched the internet,but the majority of the explanations related to csv/excel/txt... common files.
Data are from Cell ranger and spread in 3 files with following file extensions : .tsv and .mtx (barcodes.tsv, genes.tsv and matrix.mtx).
Collaborators ran Cell Ranger and gave these cell ranger output files : barcodes.tsv, genes.tsv and matrix.mtx.
Can someone give me the code to import these kind of data to R ?
Thank you in advance.
update: By "directory" you means "folder" on my computer for example ? So I need to save the 3 files under one folder and then to put name of the folder into brackets after the Read10x function ? I am a total beginner with R (I have used to R commander for statistical tests but that's it...).
update 2: I successed to import data into Seurat with the explanation of StupidWolf. I do not know C language ( I do not know almost any language ..). My system crashed a little further in the analyze after importing, not enough RAM... so now I am learning Linux commands to switch on the cluster.

Answer from StupidWolf, converted from comment:
There is a function call Read10X in seurat, you can see how to use it with the example in https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/rkit
You place the three files into a directory / folder, then you specify Read10X()... hey how big is your data? you might need to run this on a machine with enough RAM.
You can try downloading the data following instructions here:
https://satijalab.org/seurat/v3.1/pbmc3k_tutorial.html

Have a look here and the entire tutorial:
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md

To convert data from one format into another, you may find some software tool sometimes. But mostly no such tool can be found to satisfy all the requirements that we want.
In such situations, I recommend that you write a C program that reads data from the input file and writes it into the output file. In this way, you have complete flexibility so that you can determine the exact output format and what part of data to ignore or include. To do that, you should be familiar with some concepts in C language such as arrays, loops, strings, and file operations.
For an example of such conversions, you can get ideas from the following page:
https://stackoverflow.com/questions/55347322/c-program-to-convert-a-text-file-into-a-csv-file
Of course, such a converter program can be written in other languages too such as C# and Python.
The program can be developed in R language too. You need to be familiar with import functions, loops, strings, and arrays in this language.
We cannot expect that custom tools exist for all our file convertions. Every bioinformatician must be skilled in a standard programming language such as C, C#, Java, Python, or R. Then, he/she will be able to develop any custom program that reads an input file and generates a custom output file.

It's a rather onerous workaround to write C programs as interfaces between R packages and data (R generally has excellent data import functionality).

This is arguably the worst possible solution as it requires some expertise with the formats to capture edge cases and avoid messing up the data.


text:
Snakemake on cluster produces empty output files
Problem:
When I run my Snakefile on our cluster, all the processes complete successfully. I can see from the logs that my shell scripts are being run correctly and are claiming to complete successfully. About half-way through when analysis steps begins, the workflow crashes due to empty input files (which is normal for PLINK-1.9, what's not normal is that the input files are empty). I went back and checked, most of the data-containing files (*.map and *.ped, in some cases/processes Plinks binary file set of *.bed, *.bim and *.fam) are completely empty since the beginning of the workflow? (Logs still claim processes all completed and created output files successfully?)
Nuances:

When I was playing around, I accidentally messed up the dynamic naming variables (escaped the variable interpret syntax, causing literal usage) and caused creation of single file with actual variable names, not dataset names. When I did this however, the file was not empty?
Some 'side' files produced by the software (E.g. PLINK-1.9 produces a *.nosex file indicating samples with missing sex information) are not empty but the actuall data-containing output is empty?
Had a problem earlier using cut -f 4 {params.prefix}.unlifted | sed '/^#/d' >> {params.prefix3} (Intention was to collate list of entries ID's to exclude in next workflow step), where the use of >> with Snakemake on cluster environment simply would not be allowed (Permission denied error despite it being created line before by the same script?). Had to fix by adding in chmod a+rwx {params.prefix3} just before. Might this be a symptom of the same problem?
This is being executed on a PBS/Torque cluster environment where the scripts are run as user who submitted them so permissions are not supposed to be an issue. Cluster profile was built for Snakemake using this PBS/Torque Profile. (I just stripped out dependency flags because they were breaking things)

Code:
Here is the first three processes. So far, it crashes at the third CLEAN rule as the software PLINK-1.9 does not play nice with empty inputs.
rule TRIM:
    input:
        expand("rawData/{{sample}}.{extension}", extension=['bed', 'bim', 'fam'])

    output:
        expand("Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}", fileExtensions=["map", "ped"])

    params:
        fromBP = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["from"],
        toBP = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["to"],
        chr = lambda wildcards: config["locations"][wildcards.location]["GRCh37"]["chromosome"]

    shell:
        """
        module load plink-1.9
        plink --bfile rawData/{wildcards.sample} --chr {params.chr} --set-missing-var-ids @_# --recode --keep-allele-order --from-bp {params.fromBP} --to-bp {params.toBP} --out Results/{wildcards.sample}_{wildcards.location}_TRIMMED
        """

#/********* LiftOverPlink (LIFTOVER) *********/

rule LIFTOVER:
    input:
        expand("Results/{{sample}}_{{location}}_TRIMMED.{fileExtensions}", fileExtensions=["map", "ped"])

    output:
        outMap="Results/{sample}_{location}_LIFTED.map",
        outPed="Results/{sample}_{location}_LIFTED.ped",
        exclusion="Results/{sample}_{location}_EXCLUDE.dat"

    params:
        prefix=lambda wildcards: "Results/{}_{}_PRE_FILTER".format(wildcards.sample, wildcards.location),
        prefix2=lambda wildcards: "Results/{}_{}_POST_FILTER".format(wildcards.sample, wildcards.location),
        prefix3=lambda wildcards: "Results/{}_{}_EXCLUDE.dat".format(wildcards.sample, wildcards.location),
        chainFile="Binaries/hg19ToHg38.over.chain",
        LiftOver="Binaries/liftOverPlink.py",
        rmBadLifts="Binaries/rmBadLifts.py",
        sexes="rawData/1000g.sexes"

    run:
        print("Determining Liftover requirements now...")
        if config['samples'][wildcards.sample]['refGenome'] != "GRCh38":
            print("Liftover required. Dataset {} is mapped to {}".format(wildcards.sample, config['samples'][wildcards.sample]['refGenome'])),
            shell("module load liftover"),
            if config['samples'][wildcards.sample]['refGenome'] == "GRCh37" or config['samples'][wildcards.sample]['refGenome'] == "Hg19":
                print("Lifting from GRCh37 to GRCh38."),
                shell("python {params.LiftOver} -m Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map -p Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped -o {params.prefix} -c {params.chainFile}"),
                shell("ls Results"),
                shell("python {params.rmBadLifts} --map {params.prefix}.map --out {output.outMap} --log {params.prefix2}.log"),
                shell("cut -f 2 {params.prefix2}.log > {params.prefix3}"),
                shell("chmod a+rwx {params.prefix3}"),
                shell("cut -f 4 {params.prefix}.unlifted | sed '/^#/d' >> {params.prefix3}"),
                shell("less {params.prefix3} | {output.exclusion}"),
                shell("awk 'NR==FNR{{a[$1]=$4;next}}{{ if (a[$1]!= 'NULL') {{ if (a[$1] == 'female') {{ $5=2; print }} else if (a[$1] == 'male') {{ $5=1; print }} else {{$5=0; print}}}}}}' {params.sexes} {params.prefix}.ped > Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped")
                # ToDo: Add conditionals for other human reference genome builds
        else:
            print("No liftover required. Dataset {} is already mapped to GRCh38.".format(wildcards.sample)),
            shell("cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.map Results/{wildcards.sample}_{wildcards.location}_LIFTED.map"),
            shell("cp Results/{wildcards.sample}_{wildcards.location}_TRIMMED.ped Results/{wildcards.sample}_{wildcards.location}_LIFTED.ped"),
            shell("touch Resuts/{wildcards.sample}_{wildcards.location}_EXCLUDE.dat")

rule CLEAN:
    input:
        mapFile="Results/{sample}_{location}_LIFTED.map",
        pedFile="Results/{sample}_{location}_LIFTED.ped",
        exclude="Results/{sample}_{location}_EXCLUDE.dat"

    output:
        expand("Results/{{sample}}_{{location}}_CLEANED.{extension}", extension=['bed', 'bim', 'fam'])

    shell:
        """
        module load plink-1.9
        plink --map {input.mapFile} --ped {input.pedFile} --out Results/{wildcards.sample}_{wildcards.location}_CLEANED --make-bed --keep-allele-order --exclude {input.exclude}
        """
```

#Update:
I just re-enabled the dependency flags and now all outputs are looking normal again? (When I don't remove it, my cluster management software (PBS/Torque) kicks up an error, saying Invalid Job Dependancy) Why would dependancy flags on qsub matter or impact output file write operations?

https://stackoverflow.com/questions/54964785/the-touch-of-death
This is extremely enlightening. Turns out it may be a Snakemake problem?

Ok, so I have done some more debugging. It seems if I interrupt my Workflow between the rule TRIM and LIFTOVER, I can see that it is in fact creating valid outputs. As soon as I re-run the pipeline and it carries on, it wipes the input files and then as a result, generates empty output files? Why does Snakemake do this?

What command are you using to submit your jobs to cluster?

It’s a PBS/Torque cluster so uses the qsub command

snakeMake isn't generating empty output files unless you have it touch() them. It will delete them if it thinks an input file may have changes (e.g., due to being temp() and therefore having been previously removed).

Ok but my problem is the inputs. It’s wiping the inputs? Also I didn’t use touch or temp (mostly bc I didn’t know touch existed and bc I want the intermediate files as well XD)

Ok, so I feel like a right PawPaw. This is NOT a snakemake problem. The issue was multi-part but it was me, python and good ol' PBS/Torque Cluster environments.
Empty input files:
I was using a wrapper python script by the name of liftOverPlink.py to convert my files from PLINK-1.9's binary file format into the format needed by the software liftOver by UCSC. In my infinite wisdom, I did make some custom alterations. Turns out in python, there is a big difference between "w+" and "a+" for the open() function ("w+" = write mode and create if needed, "a+" = append mode and create if needed). "w+" opens or creates a file but places the cursor at the beginning of the file, effectively wiping the file if it already had content inside, wheras "a+" puts the cursor at the end of the files content (which is beginning of file if it had to create file). I used "w+" which was effectively erasing my input files content...
Why alter the script in the first place?
As it turns out, our PBS/Torque cluster was not properly serving up the software at runtime (software in question being liftOver). From what I have been able to gather, normally liftOverPlink would create intermediate files, liftOver would run and create its files where necessary, and liftOverPlink would then step in and perform format conversion back into PLINK-1.9 format. Since liftOver was not firing, the files did not exist and liftOverPlink kicked up errors since it was never designed to handle that (functionality to detect missing files would have been redundant). I attempted to ugly-fix it by using
filename = open(filename, "w+")
filename.close()

to open and close the file, effectively forcefully creating it if needed. Unfortunately that let to the above problem. Providing an absolute path to the liftOver executable fixed runtime issues and my custom code was no-longer necessary.


text:
Is there a command line tool to split a SAM/BAM file by CB (cell barcode) tag?
I have a BAM file from a single cell sequencing experiment. Each read has had the cell barcode annotated in the CB tag. Some reads do not have a CB tag. For each of an arbitrary set of barcodes, I would like to split out the reads for each of those barcodes into a separate file. Is there a tool that does this?

I wrote a package called sinto (https://github.com/timoast/sinto) that contains a function to split a BAM file into different files based on the cell barcodes.
The installation can be done via anaconda:
conda create -p conda/envs/sinto python=3.8 sinto

And a usage example from the documentation is:
sinto filterbarcodes -b my_bam.bam -c my_barcodes.txt --barcodetag "CB"

Where the my_barcodes.txt file contains cluster names in a tab separated list format:
TGGCAATGTTGAAGCG-1  A
GACCAATCACCATTCC-1  A
CAGGATTCAGAACTTC-1  B
GAACCTAAGAGAGGTA-1  B,A
..

I wrote a python script based on pysam that is free for anyone to use:
"""
MIT License
Copyright (c) 2020 Warren W. Kretzschmar
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
"""

from pathlib import Path

import pysam

class BamWriter:
    def __init__(self, alignment, barcodes, prefix):
        self.alignment = alignment
        self.prefix = prefix
        self.barcodes = set(barcodes)
        self._out_files = {}

    def write_record_to_barcode(self, rec, barcode):
        if barcode not in self.barcodes:
            return
        if barcode not in self._out_files:
            self._open_file_for_barcode(barcode)
        self._out_files[barcode].write(rec)

    def _open_file_for_barcode(self, barcode):
        self._out_files[barcode] = pysam.AlignmentFile(
            f"{self.prefix}_{barcode}.bam", "wb", template=self.alignment
        )

def main(input_bam, barcodes_file, output_prefix, contigs):
    """Split a 10x barcoded sequencing file into barcode-specific BAMs

    input:
    barcodes_file: a file containing barcodes or a single barcode
    contigs: '.' for all contigs, 'chr1' for the contig 'chr1',
    or '1-5' for chromosomes 1, 2, 3, 4, and 5
    """
    alignment = pysam.AlignmentFile(input_bam)
    if Path(barcodes_file).is_file():
        with open(barcodes_file, "r") as fh:
            barcodes = [l.rstrip() for l in fh.readlines()]
    else:
        barcodes = [barcodes_file]
        print(f"Extracting single barcode: {barcodes}")
    writer = BamWriter(alignment=alignment, barcodes=barcodes, prefix=output_prefix)
    if contigs == ".":
        print("Extracting reads from all contigs")
        recs = alignment.fetch()
    else:
        if "-" in contigs:
            start, end = contigs.split("-")
            print(f"Extracting reads from contigs {start} to {end}")
            recs = (alignment.fetch(str(contig)) for contig in range(start, end + 1))
        elif "," in contigs:
            contigs = contigs.split(",")
            print(f"Extracting reads from contigs {contigs}")
            recs = (alignment.fetch(str(contig)) for contig in contigs)
        else:
            print("Extracting reads for one contig: {contigs}")
            recs = (alignment.fetch(c) for c in [contigs])

    for rec in recs:
        try:
            barcode = rec.get_tag("CB")
            writer.write_record_to_barcode(rec=rec, barcode=barcode)
        except KeyError:
            pass

if __name__ == "__main__":
    import sys

    main(
        input_bam=sys.argv[1],
        barcodes_file=sys.argv[2],
        output_prefix=sys.argv[3],
        contigs=sys.argv[4],
    )
```

You can do it with samtools:
# Save the header lines
samtools view -H $BAM_FILE > SAM_header

# Filter alignments using filter.txt. Use LC_ALL=C to set C locale instead of UTF-8
samtools view $BAM_FILE | LC_ALL=C grep -F -f filter.txt > filtered_SAM_body

# Combine header and body
cat SAM_header filtered_SAM_body > filtered.sam

# Convert filtered.sam to BAM format
samtools view -b filtered.sam > filtered.bam

And 10x Genomics has developed subset-bam a tool that can split a BAM file by barcode.

If you have access to a lot of cores:
10x has subset-bam a pretty fast way to pull out all the cells with a specific CB/BX if you only need a small number.


text:
How are the scores of GeneSplicer, MaxEntScan and SpliceRegion interpreted from VEP annotations?
I am using VEP from Ensembl to annotate my VCF files with the extra plugins of GeneSplicer, MaxEntScan and SpliceRegion. However, I don't fully understand the output of these scores.
I know that they are used for splicing predictions on trained sets.

What are the range of these values ?
What does a high or low value mean ?
What is the threshold ?

Can someone please explain them please ?


text:
Platanus-allee phasing fail: Error(13): Error, SolveDBG exception!
I am using Platanus-allee 2.2.2 for heterozygous genome (~500mb) assembly with Illumina short reads and PacBio reads input data.
I have the file contigs.fa from short reads but phasing step with PacBio reads gives the following error after running for couple of hours

Error(13): Error, SolveDBG exception!! platanus_allee solve_DBG
command failed.

Here is the script used for contig assembly and phasing
num_threads=20
tmp_dir=/path-to-/tem/platanus_allee
memory_GB=100
out_prefix=out
minimap2_path=minimap2

# Input reads

PE=trimmed_L1_*.fq
PacBio=pacbio???.fastq.gz

# Contig assembly commands
platanus_allee assemble -t ${num_threads} -tmp ${tmp_dir} -m ${memory_GB} -o ${out_prefix} \
-f $PE \
2>${out_prefix}.assembleLog

# phasing commands
platanus_allee phase -t ${num_threads} -tmp ${tmp_dir} -o ${out_prefix} -mapper ${minimap2_path} \
-c out_contig.fa \
-IP1 $PE -p $PacBio \
2>${out_prefix}.phaseLog

cat out.phaseLog
platanus_allee version: 2.2.2
platanus_allee phase -t 20 -tmp /path-to-/tem/platanus_allee -o out -mapper minimap2 -c out_contig.fa -IP1 trimmed_L1_1.fq trimmed_L1_2.fq -p /path-to-/pacbio*.fastq.gz

Error(13): Error, SolveDBG exception!!
platanus_allee solve_DBG command failed.

tail -n 15 path-to-/out_intermediateResults/round1/out.solveDBGLog
dividing nodes based on bubbles ...
dividing scaffolds to adjust bubble-breakpoints...
dividing indonsistent bubble-ends...
pairing bubble nodes...
dividing scaffolds to adjust bubble-breakpoints distiguishing strands ...
NUM_DIVISION = 1832
dividing scaffolds to adjust bubble-breakpoints...
dividing indonsistent bubble-ends...
pairing bubble nodes...
dividing scaffolds to adjust bubble-breakpoints distiguishing strands ...
NUM_DIVISION = 0
TOTAL_NUM_DIVISIONS =1832

dividing erroneous scaffolds based on base-level coverages ...
Segmentation fault (core dumped)

A bug causing the segmentation fault seems issue on Platanus-allee v2.0.0 and earlier

I would suggest talking to the platanus authors. A segfault without any accompanying error message or informative logging suggests an unforeseen bug, and this is really not much information for us to work from. It does not look like it is related to the previously fixed bug that causes a segfault. As usual of course it is a good idea to check all your input files but it does not appear that that is an issue.

Thanks @MaximilianPress. I wrote to correspondent author of Platanus-allee and didn't get response!

That isn't good! I hope they get back to you. In the meantime, you might look into alternatives for assemblers/scaffolders. I hesitate to recommend because I don't know what your goals are, though.


text:
Getting sequences of one transcript per gene out of annotation and genome
I got a genome (data/genome/genome.fasta) and braker-based genome annotation (data/genome/annotation.gff3), now I would like to get sequences of one transcript per gene (let's say the longest one).
Cufflinks has a gffread program to get the transcriptome out of the genome and the annotation. Something like
gffread data/genome/annotation.gff3 -g data/genome/genome.fasta -w data/genome/transcripts.fasta -y data/genome/proteins.faa

Will do the job. However, it is pulling multiple transcripts per gene. I read the help page twice, found parameters that could potentially do what I would like to (collapsing alternative transcripts into single one), but even when I specify my candidate parameters -M -K I get still .t1 and .t2 sequences for some genes.
The task is surprisingly simple, and I got already a custom script solution, but it just feels wrong. There must be a better way!
What is a reasonable/elegant way to get nucleotide sequences of genes out of a genome and annotation? Am I missing something with cufflinks?
Related, not exactly identical:

How to get the longest fasta sequence including all possible switching isoforms of a gene out of isoforms (dealing more with terminology of Trinity)
When I used Kallisto/Salmon for differential expression analysis, I collapsed transcripts into genes while importing to R via tximport.

By gene sequences, you mean the entire sequence including the introns? I have been in similar situations before and my approach was to create a BED file from the GFF3 using some set of awk commands followed by bedtools getfasta for sequence. I am curious what your custom script solution was.

I made an edit, I mean to have one transcript per gene. So I suppose the longest transcript would be the most appropriate term. I just want to avoid multiple isoforms from the same gene.

I posted the script here too...

Here is the custom python script solution:
#!/usr/bin/env python3

from collections import defaultdict
from Bio import SeqIO

# read file into dictionary of transcrpt->gene
# jg25690   jg25690.t1
# jg25690   jg25690.t2
# jg25691   jg25691.t1
transcripts2genes = dict()
gene2transcript_mapfilename = 'data/genome/transcripts2genes.map'
with open(gene2transcript_mapfilename, 'r') as g2t_file:
    for line in g2t_file:
        gene, trainscript = line.rstrip('\n').split('\t')
        transcripts2genes[trainscript] = gene

# construct a dictionary gene->list of transcripts
# this fasta file is output of default gffread -w
genes2transcript_seuqneces = defaultdict(list)
transcripts_fasta_filename = 'data/genome/transcripts.fasta'
transcripts_fasta_file = SeqIO.parse(transcripts_fasta_filename, 'fasta')
for seq_record in transcripts_fasta_file:
    gene = transcripts2genes[seq_record.name]
    genes2transcript_seuqneces[gene].append(seq_record)

genes_fasta_filename = 'data/genome/genes.fasta'
with open(genes_fasta_filename, 'w') as out_file:
    # iterate through genes
    for gene in genes2transcript_seuqneces.keys():
        # get all the associated transcripts
        transcripts = genes2transcript_seuqneces[gene]
        # find the longest, by default the one with index 0
        lengest = 0
        for index, tr in enumerate(transcripts):
            if len(transcripts[lengest].seq) > len(tr.seq):
                lengest = index
        # this would be nicer with Bio interface (with wrapping of sequence lines)
        # but I was lazy and blast does not care, only thing I do is converting sequence to uppercase (not even sure if necessary)
        fasta_string = '>' + gene + '\n' + str(transcripts[lengest].upper().seq) + '\n'
        out_file.write(fasta_string)

Do you know which one of the 3 steps takes the longest in your case? UCSC Kent Utils have some tools that you can use as substitution. One change I would suggest is that when you are constructing the genes2transcript_seuqneces dictionary, just keep the longest transcript and nothing else. That way, you are not storing the entire transcripts.fasta in memory.

Not sure which of the steps takes the longest, but all of this takes only a few seconds, so I don't think it's really worth any optimisation effort. I will take a look at UCSC Kent Utils.


text:
Phylogenetics analysis between diverse species
Hope everyone is well.
So i was appointed with a task, and i need some guidance. Attention, the one that gave me the task may not have explained this very clearly. I will write what he exactly asked of me, and what he gave me as primary data.
So he gave me 4 bed files (human genome) with chr and coordinates(start/stop), and asked me to do: "sequencing/coordinates conservation between diverse species".
Is that doable with coordinates only or do i need sequences in fasta format (phast uses fasta format)? I just dont understand how i will find "conservation between diverse species" only with coordinates. It seems weird to me.
I have found these programs: TargetOrtho2, Phast and BUSCO(?). Are they any good? Should i continue with one of them that you may have experience with?
Thanks

It is indeed a very unclear task. Are all the 4 bed files human bed files? What distant species you are interested in?

Hi @JimHat That is cryptic, I think you need clarification! Literally that means non-humans.

I would suggest looking at those BED regions in the UCSC 100-way conservation tracks. Just doing that might just suffice.

hi all, and thanks for the answers! @KamilSJaron, yes they are bed files. He didnt clarify which distant species but i am guessing species that represent some taxonomy (or sub catergory of one, also maybe a model organism).

@DevonRyan i wll try that. Thanks

I second commenter Devon Ryan's suggestion, I'm just going to give a little more detail about how to do that. You can plug the BED files into the UCSC browser following these instructions: https://genome.ucsc.edu/goldenPath/help/customTrack.html
Now you can look at the overlap between your bed tracks (black boxes I believe) and conservation (red box).

Assuming that is all you want to do, that should do it. If you want to do further data analysis or processing, you can download the conservation track as a .wib file from UCSC. You can convert this file to BED as described here, and use BEDTools to analyze it together with your BED of regions of interest.

Hi and thanks for the answer. By doing this, i will find whatever correlation there is between these bed files (which by doing what you suggested i will find the regions that are conserved(?) ) and -some- diverse species with bedtools?

Hi @JimHat, it is a little hard to tell exactly what it is that you are supposed to do based on the question, the simple answer is "yes", you will look for overlaps between conserved regions and your regions of interest using commands like bedtools intersect. But I don't have a clear idea of what your goals are beyond that. Also, I looked a little deeper and the browser doesn't have BED files of the conserved regions but it does have .wib files, which can yield bed files. See updated answer.


text:
What RNA-Seq expression value would be closest to Microarray equivalent?
I know this question may seem strange.
I'm using Spearman correlation between gene expression profiles for various reasons (I won't go into details here). As a result, I often compare RNA-Seq and Microarray samples.
For preliminary analysis, I usually grab what version of the data is easily accessible (RPKM, FPKM,...), but I'd like to dig a bit more.
Intuitively, I'd think a value similar to RPKM or such would make more sense than raw counts, which is why I usually convert raw counts to RPKM (I know they're considered obsolete for statistical analysis of the data, what I'm doing relies on the correlation)
So I ask what RNA-Seq gene expression value would be closest to a microarray equivalent ?
i.e. What would (theoretically) maximize correlation between gene expression of the  same sample profiled with RNA-seq and microarray. RPKM ? FPKM ? TPM ? etc.

“I know they're considered obsolete for statistical analysis of the data” — No, they’re considered obsolete for every purpose. Just don’t use them. Read https://bioinformatics.stackexchange.com/a/69/29.

I think it is very hard to say which are the closest because they are not really comparable. But since you are using Spearman correlation, I guess RPKM, FPKM, and TPM do not change the order of gene expression levels. You might also want to normalize RNA-seq and microarray data so that they are more comparable.

I did a comparison of cDNA count data against microarray data that was published a few years ago:

For comparisons to published data (Fig. S2; Miller et al., 2012), a generalized linear model was fitted to the relationship between log-transformed microarray and VSTPk expression levels obtained from the ImmGen Project database, and was used to transform the microarray data into values comparable to our VSTPks.

I found that the Variance-Stabilizing Transformation that was carried out by DESeq2 was close to what I wanted, but there still seemed to be a length-based bias to the reads. I corrected this by dividing by the length of the longest gene isoform in kilobases (creating something that I called VSTPk).
After doing this, there were range differences between the microarray and cDNA data, so I did an additional linear transformation to get the data fitting as close as possible.


text:
How to visualise the proportion of clonal genes to sub clonal genes
I have a data frame
Response_TO_CHEMOTHERAPY    Gene    CCF
Responders                  POLQ    1
Responders                 CDKN2A   1
Responders                  TP53    1
Responders              CCDC102B    0.906729075
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders           CCDC102B   1
Non-responders             TSHZ3    1
Non-responders             PIK3CA   1
Responders                 SLIT2    1
Non-responders              SLIT2   1

I have two groups of patients; Responders to chemotherapy and non-responders. For each group I have calculated cancer cell faction (CCF) for a fixed set of genes so that an individual gene has several reads in each group
If CCF = 1 the gene clonal and if CCF < 1 gene is sub clonal, for a given gene for each group, I therefore would have the number of clonal and sub clonal for this gene, for example 2 clonal 7 sub clonal which would be 2/7 something like the below figure

This is another potential example which reveals what is in my mind

Which show the clonal frequency for each group
I don't know how to use dput() but this is my complete data
https://www.dropbox.com/s/ol7m4nlkuyuy85a/c1.txt?dl=0

Sorry why you are unvoting my question? Is not it clear enough? When I am asking a question that means I need to solve that not that I am posting a question for fun. I invest time to prepare input and to know what exactly I am looking for

I've used some dummy data for simplicity. It's fairly easy to do this with ggplot, but you have to format your matrix-like data for ggplot input.
eg :
# make data
set.seed(1)
ngenes <- 20
nsamps <- 10
my_dat <- matrix(rnorm(ngenes*nsamps), ncol = nsamps, nrow = ngenes)
colnames(my_dat) <- paste("samp", seq_len(nsamps))
rownames(my_dat) <- paste("gene", seq_len(ngenes))

library(ggplot2)
library(reshape2) # for melt

# format data for plotting
dat_toplot <- melt(my_dat)

#plot
ggplot(dat_toplot, aes(x = Var1, y = Var2, fill = value)) + geom_tile()

plots

If you need specifics on plot tweaking, google's your friend. There are plenty of tutorials on how to use ggplot for heatmaps. (eg. here, or here)
EDIT
I loaded your data into the dat object.
library(tidyverse)
dat <- as_tibble(dat)

dat_counts <- dat %>% group_by(Sample, Gene) %>% summarise(count_subclonal = sum(CCF<1), count_clonal = sum(CCF>=1))
dat_counts$ratio <- dat_counts$count_subclonal/dat_counts$count_clonal

# set infinite values (/0) to NA
dat_counts$ratio[is.infinite(dat_counts$ratio)] <- NA

ggplot(dat_counts, aes(x = Sample, y = Gene, fill = ratio)) + geom_tile() +
  geom_text(aes(label = paste(count_subclonal, "/", count_clonal)), col = 'white') # add ratio text

Outputs :

Thank you but this is not what I need, I need the proportion of clonal to sub clonal genes in each box not a heat map of cancer cell fraction

@Exhausted Ah, then your question is not on generating a heatmap. It is "how to compute the proportion of clonal to subclonal genes"

@Exhausted See my updated post


text:
How can I find the position of every mutation where the Allele Frequency is greater than X, in regards to a reference such as Hg19?
I have a Human genome sequenced in a BAM file (along with other files with the indels, snps, cnvs).
I want to find every mutation with regards to the reference Human genome. However the majority of "variants" are noise/read errors. I am therefore required to filter the variants via a minimum threshold based on allele freq (eg. 0,5).
Using a reference (eg Hg19), how can I find the position of every mutation where the allele frequency is greater than the threshold?
PS. I am just starting, so I don't know which are the best (command line?) tools for the job nor how to build a pipeline for this.
Help would be greatly appreciated

Hello and welcome to the world of bioinformatics,
BAM files are compressed formats of SAM and you can read more about it here.
To answer your question, I presume that your BAM is generated from a Whole Genome library preparation. Calling mutations on a Whole Genome takes a longer time than a Whole Exome since you are looking at every position instead of the exons of your genome.
What you are looking for is a special format of files called Variant Caller Format (VCF) generated by Variant Callers tools listed below in comparison studies. There isn't a tool better than another, they all have their advantages and disadvantages. I suggest that you read about them as they also describe their pipeline approaches.

Comparison of three variant callers for human whole genome sequencing
Comparing the performance of selected variant callers using synthetic data and genome segmentation
Performance assessment of variant calling pipelines using human whole exome sequencing and simulated data
Systematic comparison of germline variant calling pipelines cross multiple next-generation sequencers
Variant Callers for Next-Generation Sequencing Data: A Comparison Study

Once the VCF is generated, you will also have a column named AF which corresponds to the allele frequency also computed when using these tools. There are tools that allow you to filter VCF files to extract the needed information with a wished criteria like SnpSift,  vcftools, VariantFiltration. So, I suggest reading the above articles as a start.
As for building pipelines, it depends on what language you feel at ease at. However, if I could suggest, Nextflow is a good start. You can also use Snakemake if you like Python. Or maybe simply bash ? In any case, there are many pipelines already made to help you through your needs.
(If this helps, don't forget to either accept the answer, vote the answer or both).


text:
Concerning the current outbreak, using bioinformatics
Dear Researches community
A question concerning the current outbreak, recently we attempting to contribute to the study of the Covid19, with a tight submission time ( 4 weeks).
Inquiring help because I'm not a virologist with an insufficient biological background. However, I assume able to do in silico analysis or ( dry lab).
Moreover, new sample sequences from local strain will be provided. Assembly and annotation could be done at the first stage. Then perhaps a comparative genome investigation between the old SARS and the Cuvid19.
Furthermore, The virus estimation length that causes all this is only around 29,829 base pairs of RNA. The genome contains at least 10 predicted open reading frames (ORFs): ORF1a,  ORF1b, S,  3,  4a, 4b,  5,  E,  M  and  N,  with sixteen predicted nonstructural proteins being encoded by ORF1a/b.
I am wondering if someone can assist by considering a particular section to focus on or a study suggests. Perhaps we can contribute at least 0.001 % to the communities.
I have been reading articles around, and seems are thousand of them every day.
so I have come up with starting list:

Genome assemblies, annotation
epitome studies
Here you predict the protein sequences based on the genome that will sequence, then there are a number of tools that analyze the immunogenicity of different proteins
That can help you identify parts of the COVID-19 proteins that are likely to be detected by the immune system and they, of course, can be targeted by vaccine developers.
comparison between old saras Vs newCoVID
The M protein is a good target. It is extremely abundant in the nucleocapsid and is the protein responsible for giving the virus it's the shape, Almost everyone is looking at the S protein, which is justified. But there is not as much focus on the M protein.

Your response is highly appreciated if you could attach the any reference related to the list above.

Is there a precise scientific question?

I would advise editting your question, firstly grammar (it is difficult understand) and then to a technical question. One vote for closure has already been initiated (not by me). At present you are outside the remit of the site for asking a question.

It is bioinformatics community where we all exchange ideas, students, doctors, researchers if you don't have an answer, please escape the page and keep the comment to your self @ATpoint

@Michael  An apology dear for my writing, excuse my language it's not my mother tongue still learning :(, anyway thank you I have edited again I hope looks better

@BioInfo Calm down. My question intended to know what you want to analyse in order to give a productive suggestion on how to achieve that. The answer would've essentially been the four points you added to your question by now. I will leave this thread then as you request. Good luck with your research.

My personal approach to something like this is to first try to reproduce existing work in order to learn to use relevant tools and understand the concepts involved. Then as you get more comfortable with the subject matter, you may spot an opportunity to add value.
If you'd like to get started with the phylogenomics side of things, you can use this public workspace developed by the viral genomics group at the Broad Institute, which includes viral genome sequences from NCBI SRA as well as some workflows for doing genome assembly and quality control: https://app.terra.bio/#workspaces/pathogen-genomic-surveillance/COVID-19

FWIW we recently updated the viral genome assembly workflow and added a workflow for running Augur (from NextStrain) to do the phylogenomic analysis of multiple strains.

Please note that Cuvid19 is COVID-19, the WHO nomenclature. SARS2 is specifically called SARS-CoV-2 and is the virologists name for COVID-19 because that is the nomenclature of the International Committee on the Taxonomy of Viruses (ICTV) for the same virus following formal taxonomic investigation. It is important not to get the names confused.
Technically you are asking to compare COVID-19 with COVID-19. Whilst that is a standard investigation this would be technically rewritten "to investigate the molecular epidemiology of COVID-19".  If you are working within the virological community then the virus is termed SARS2-CoV-2. Phylogeny is a key approach, rather than "comparative genomics".
I wish you luck with your study.

Dear Thank you, so you suggest after sequencing and assembly, I construct a phylogenetic tree

SARS2 and  SARS-CoV-2 would be the names of the virus (as synonyms). COVID-19 is the name of the patology, the name of illness


text:
Use a specific python.exe instance with PyMOL/Windows?
I want to use RDKit with a python script in PyMOL. RDKit requires a conda environment, so I think I need to be able to specify that the instance of python that I want to use would be in that environment, instead of the default for PyMOL, so that I can run it through PyMOL. The script I made is working for me through PyCharm when I test it, by using the RDkit environment "rdkit-env" as my python interpreter, with all of the additional modules I need installed there as well.
I would like to run it through PyMOL so I can add some PyMOL specific functionality as well. I've tried using shebang to specify the RDkit environment that I made, but that does not seem to be working for me with Windows.
Could someone please help me figure out what other approach might work for this? I'm new to this area, so I would appreciate any insights. Thank you!

So it is really hard to make a given python instance run a library from a different instance. It is not just one dist-packages library folder, but multiple that you would need to link, and all that assuming they use the same version of python.
So the best bet is to install a new instance of pymol the python module.
The module can be installed in a python that uses conda with
conda install -c schrodinger pymol

This is the incentive version of PyMOL 2 and does not require a licence.
You installed rdkit, so I am assuming you used conda as conda install -c rdkit rdkit. But just in case you went through the pain of compiling rdkit for a given python instance without conda and want the same for PyMOL. There are guides of instructions required to compile open source. Such as https://pymolwiki.org/index.php/Windows_Install. On a barebone OS like a Pi's Raspsian, I had to do a lot more than claimed say sudo apt-get install python3-dev libglm-dev freeglut3-dev libglew-dev libpng12-dev libfreetype6-dev build-essential libxml++3.6-dev and link a few things manually. So just like for rdkit, conda make life a breeze.
PS. RDKit and PyMOL? Add PyRosetta and you just my favourite modules!


text:
Unmapping Alternative reads in BAM file
I want to do some HLA typing, most of the tools require the bam file is aligned to primary genome without alternative read handling. From the International Genome Sample Resource project, I can get the high coverage file, but they are all aligned to primary genome with hla, decoy and alternative reads. I want to make the alternative reads into unmapped reads.
I think about unmapping the whole file to fastq, and then realigning according to different assembly, but it's time-consuming and also costs a lot of storage, so I don't want to follow that approach.
I also tried this script, but so far it's not working (I have all requirement packages installed).
Any help, tool or insight for this problem is warmly welcome.
Thank you

hmmm are you familiar with samtools? you can convert unmapped reads on the fly (on the console) and then remap it and output a bam file

the only thing i am not sure about where do your HLA reads go in the genome, if they are multimapping or you know that they map to a specific location, you can also do something similar to above

thank you for your comment. I'm not so familiar with samtools, actually, any answer for me is valuable. For the HLA allele, I use package like optitype and kourami, they can handle it if I can give them the correct input.


text:
BERT Language Model and Gene Sequences - How Do I Relate Clusters of Sequences?
I hope you'll indulge a question from a computer scientist with limited bioinformatics knowledge.   I've been working with the Google tool for language modeling called BERT.   It's generally regarded as state of the art when encoding language.   It allows you to say that the sentences "I went to the bank to cash my check." and "My wife is going to the Credit Union to pay the mortgage" are more similar than if compared to "he was fishing and fell off the bank."   In short, BERT takes sentences and encodes them so that related concepts are near each other and unrelated concepts are farther apart.
BERT is an amazing tool and I wondered if it could successfully model the language of the genome if it were encoded into "sentences" (I chose 7-mer words and sentences of length 64, arbitrarily).  Below is the result of encoding the full genomes of human, bonobo, gorilla, chimp, and orangutan.  There are about 270K "sentences" here.   You can see there are four major blobs, with indications of small tight clusters.   BERT has decided that the points indicated in red are somehow semantically related.  I think is is significant, or at the very least interesting.

My first thought was that clustering was trivial degenerate cases.  Maybe all the "AAAAAAA..." were put near each other, and all the "TTTTTTTTT...." etc.  I printed out samples from the clusters and visually that is not the case.  Also I've made histograms of the average Levenshtein distance in groups.  Within a tight group, the LD is statistically much lower (but not trivially low) as compared to unclustered regions (the big blob is all unclustered).

If there is a biological basis for the clustering, then new sequences could be similarly transformed, and classified according to cluster membership (and related biological function).  This exhausts my biological knowledge.  I've taken sequences from clusters and used BLAST to try and get some insight but I'm way out of my comfort zone with that.
My primary question is, given a set of sequences (in my case len=448), how do you figure out what they all have in common (given my hypothesis is that BERT says they are related based on proximity)?    And my backup question is, do you think this is interesting or a waste of time?

I love the question because I used to work on this topic. However, even an incomplete answer is a thesis (namely Titouan Bougouin-Laessle's thesis). The association can be biological, we used a much more naive clustering approach and got quite a decent association. I will answer at length later!

blasting doesn't tell you anything?

StupidWolf - I can BLAST the sequences in each cluster, but I don't have the background/knowledge to interpret the results.   Is there something you would look for specifically?

Can you share the code/data ?

The UMAP results look very pretty :)
If I were to try to make sense of this blob, I'd try to do the following.

Pull out Gene Ontology (GO) terms for a list of candidate genes (you can try to subset Uniprot / Swissprot for some examples). Then I feed through some genes through your fragmented BERT kmer approach to see where they land. I'd be particularly interested in seeing where transmembrane proteins and enzymes land here.
I'd try to color this by introns vs exons (you should be able to extract this sort of information from your genome files). There is a good chance that the difference between introns / exons could explain some of these differences.
It maybe worth looking into genes that only encode RNAs vs genes that encode proteins.

In truth deep neural networks quite mainstream now, albeit they are a pain because of their reliance on very expensive graphics cards to process the data (GPU) rather than traditional CPU.
There are two issues,

The algorhithms need an application to a biological phenomena and this is dependent on the model. DNA applications (usually) use a different model to language applications (below)
The model is only as good what it is trained against. In languages this would be people saying stuff, like in your example. In biology it is a phenotype. Thus you need to feed it the right answer and all variants of the right answer. It will then spot the pattern and recognise this pattern in de novo inputs.

What we are talking is deep learning using a neural network  yes we use them, but again the models are (usually)  different. Convolution neural networks (CNN) are common in DNA/genomics because it assess the physical linkage of DNA. CNN is also used in image recognition, so like when you type 'dog' into google and see dog images, chances are the search is a consequence of CNN classification.
In languages recurrent neural networks particularly long short term memory (LSTM) models are common because they need to learn the patterns and idioms of a language. Google of course will only use their development called TensorFlow to build their models, which is very geeky, for example they have a programming construct called a "placeholder". In my view a placeholder is something found in Adobe PhotoShop rather than a programming language etc... TensorFlow exponents outside Google tend to use the Keras interface. In Python there is Natural Language Toolkit — NLTK 3.5 to do this.
Summary So, to answer your question, by dumping DNA sequence into a language RNN its fun, but it wouldn't mean very much because the evolution of language, or construction of a given language (i.e. the training set) is a very different data to mechanism(s) involved in the evolution of DNA.
Fundamentally, we don't speak DNA.
I think it would be fun dumping whale songs and dolphin squeaks into these programs to see whether there is any convergence with our own speech.

My advice is to use an unsupervised training approach to identify your clusters. This would be model-free and training-free and would help further analysis in "supervised-learning" OR just stick with the clusters it generates.


text:
how do i change a variable with iterator under a for loop in python
This question comes up in my bioinfo class.
I have a csv file containing gene names and the cluster they belong to. I want to run GESA on each individual cluster. So I would like to extract the list of gene names corresponding to each cluster. So I wrote code like this
import pandas as pd
gene_list = pd.read_csv('marker_genes.csv')
for i in range(12):
    "cluster_"+str(i) = gene_list.gene[gene_list.cluster==i]

I know this raises syntax error, so what is the correct way of doing it.
In the end, I would want cluster_0 to cluster_12 each being a list of gene names.

You're looking for the groupby() function https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.groupby.html

Groupby is probably right, not always easy to get the object to spit out the values. Difficult without knowing the table format

Reshaping via melt is another possible way in or pivoting might work. There's probably lots of solutions. The pd.read_csv is a neat part of pandas, however you rarely use loops, only when you are stuck and its a different syntax.

cluster_0 to cluster_12 are just strings. they cannot be the names of arrays or matrices. Strings can only be used in print statements and string processing statements.
The correct way to define an array or matrix is as follows: call the arrays as: cluster[0] to cluster[12]
then, assign:
  for i in range(12):
        cluster[i] = gene_list.gene[gene_list.cluster==i]

then whenever you need to print the names of these arrays, you use:  "cluster_"+str(i)

Okay .... yeah, this is better. There is an issue with the "cluster" dictionary (??), so I still think it will throw an error, but it is a definite improvement. The pandas command that is the formal solution is "iterrows" for this style of loop. – its a fun question.


text:
Visualizing protein-ligand interactions with LIGPLOT or any other suitable software?
I have used DOCK software to dock 10 ligands into the protein structure.
I want to generate a pdb file that contains all 10 binding modes of the ligands together with the protein. So all together in one pdb file.
Now, what I want to know is when I give this file to LIGPLOT, is it possible to generate something like:
%polar interactions XX
%non-polar interactions XX
%hydrogen interactions XX
Is there a way to do it with LIGPLOT, if now is there some other software that can be recommended for this type of assessment?
I have been having hard time, trying to find a way or software that can perform similar type of calculation and give me similar type of output that I want.
Best

"binding modes" --> "poses".
You will need to do it programmatically. This post for example use Python pymol for querying atoms. Can you use Python, Perl or Julia?

Yeah, I assumed so, cus it tends to be a bit more complicated. Unfortunately, I cant say for my self that I know any of the three languages :(

You can visualize with Discovery Studio Visualizer and Schrodinger.
Arpeggio and PLIP can generate output files of the different interactions.

Hey, thanks! already found that PLIP can do it and it works pretty well so I went for it.


text:
How to download RNAseq gene expression data from GTEx
I would like to download RNAseq gene expression data from GTEx from the liver only. What is the best way to do this?
I have tried a few things. I downloaded several files from the GTEx website (https://www.gtexportal.org/home/datasets, under "RNA-Seq Data"), but the problems I run into are:

Not sure what is the best file: "Gene read counts" or "Gene TPMs"? What is the difference?
I would like 1 value per sample per gene. There are many duplicate genes in these 2 files. How can I collapse them to 1 value per gene?
How can I download the data for liver tissue only, instead of ALL data for ALL tissues?

Thank you for suggestions.

What do you eventually want to do with the data? That helps deciding what to use.

What do you mean by duplicates? I assume you mean multiple entries with the same gene name. Gene names are not unique identifiers. Also, whether you want raw counts or TPMs depends on what you want to do with the data. If you want to calculate differential expression between 2 samples, you will need counts. Honestly, comparing expression values of different genes is apples to oranges. Plus, you need replicates, ideally from the same batch. RNA-Seq is very noisy data. Finally, if you want a specific tissue, you have to find samples of that tissue.

@ATpoint and hepcat72  Thank you. I have found that subsetting the data worked to get data from one tissue only. The gene names start with ENSGxxxx. I found out that they are different transcripts, sometimes of the same gene, so I suppose there are no real duplicates.

Gene read counts are the number of reads or fragments that are mapped to each gene. So it is an integer and the values are not directly comparable between samples due to differences in sequencing depth. TPM stands transcripts per million. It is the normalized gene expression level. Basically, you first normalize read counts by gene length and then normalize the library size to 1e6 in all libraries. As for which one should you use, it depends on what analysis you want to do and what tool you are going to use.
I have not noticed that there are duplicated genes in those files, can you show an example?
I do not think you can download expression from one tissue, but you can always subset. The whole dataset is not very large.

Thank you very much. Subsetting worked for me.


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:
Loading all PPI information from Uniprot into a graph in NetworkX
I am looking for guidance to load all SwissProt PPI data for Homo sapiens into a graph (from NetworkX). In particular, the part that I'm struggling with is gathering the data from SwissProt for each human protein.

The data that I'm looking to extract is protein-protein interaction. This is not included in the files found here, which @rtviii (answer below) suggested to use. Rather, when say you click on an entry in SwissProt, for each entry, there is a section devoted to "Interaction" as shown in the image attached below. I would like to extract that information for each human protein, and load it into a graph data structure using NetworkX

EDIT 2:
How I ended up doing it:

    query = "reviewed:yes+AND+organism:9606"
    args = {'query': query, 'format': 'xls', 'columns':'id,entry name,interactor'}

    uniprot_url = "https://www.uniprot.org/uniprot/?{}".format(urllib.parse.urlencode(args))

    with open(<DESIRED FILE PATH>, "wb") as target:
          target.write(requests.get(url).content)

What are you trying to do exactly? Are you starting with a list of proteins and want to find their interactions? If so, why not go straight to Mint or Intact (but be careful with Intact, there's a lot of dirty data in there)? Why would you start from uniprot? If you want to build a high quality human PPI network, this really isn't the simplest or best approach. You might want to ask a new question, taking a step back to basics.

The query link for all interactors would be thus : https://www.uniprot.org/uniprot/?query=organism%3A9606+AND+interactor%3A*&sort=score

There is always an option of downloading the static version of the service if you find it incovnient to deal with their API for your use case from here and have it all to yourself, pipe data into your dataframes locally. I find it useful to just broswe their repositories sometimes instead of reading the specifications.
In general, you can get a variety of datashapes  from uniprot that you cannot get from their static database-image by specifying the appropriate parameters in the request. You can browse and explore parameters by expanding their GUI's search fields as described in this little guide. Then you can request your data like so, in python for example.
import requests
URL = "https://www.uniprot.org/uniprot/"
PARAMS = {'organism': 9606, 'format': 'tab', ....#Additional parameters#...}
response = requests.get(url=URL, params=PARAMS)

Thanks for the reply. So I was already able to download the database data, I'll explain in more detail what it is exactly that I need in the post.

@AnIgnorantWanderer Added query for interactors. Tag on &format=tab to make your life easier once you get to python.